diff --git a/sPuReMD/src/allocate.c b/sPuReMD/src/allocate.c
index c235913acf00b81a11737c009e9f31d50a1a1019..d2bf99717f4d0f4f9b5c1fda4db9026469699507 100644
--- a/sPuReMD/src/allocate.c
+++ b/sPuReMD/src/allocate.c
@@ -152,33 +152,6 @@ static void Reallocate_Matrix( sparse_matrix *H, int n, int n_max, int m )
}
-void Initialize_HBond_List( int n, int const * const h_index,
- int * const hb_top, reax_list * const hbond_list )
-{
- int i;
-
- /* find starting indexes for each H and the total number of hbonds */
- for ( i = 1; i < n; ++i )
- {
- hb_top[i] += hb_top[i - 1];
- }
-
- for ( i = 0; i < n; ++i )
- {
- if ( h_index[i] == 0 )
- {
- Set_Start_Index( 0, 0, hbond_list );
- Set_End_Index( 0, 0, hbond_list );
- }
- else if ( h_index[i] > 0 )
- {
- Set_Start_Index( h_index[i], hb_top[i - 1], hbond_list );
- Set_End_Index( h_index[i], hb_top[i - 1], hbond_list );
- }
- }
-}
-
-
static void Reallocate_List( reax_list * const list, int n, int n_max,
int max_intrs, int type )
{
@@ -238,7 +211,6 @@ void Reallocate_Part2( reax_system const * const system,
{
Reallocate_Neighbor_List( lists[FAR_NBRS],
system->N, system->N_max, realloc->total_far_nbrs );
- Init_List_Indices( lists[FAR_NBRS] );
realloc->far_nbrs = FALSE;
}
@@ -255,7 +227,6 @@ void Reallocate_Part2( reax_system const * const system,
{
Reallocate_List( lists[BONDS], system->N, system->N_max,
realloc->total_bonds, TYP_BOND );
- Init_List_Indices( lists[BONDS] );
realloc->bonds = FALSE;
realloc->thbody = TRUE;
@@ -263,9 +234,8 @@ void Reallocate_Part2( reax_system const * const system,
if ( control->hbond_cut > 0.0 && workspace->num_H > 0 && realloc->hbonds == TRUE )
{
- Reallocate_List( lists[HBONDS], workspace->num_H, workspace->num_H_max,
+ Reallocate_List( lists[HBONDS], system->N, system->N_max,
realloc->total_hbonds, TYP_HBOND );
- Init_List_Indices( lists[HBONDS] );
realloc->hbonds = FALSE;
}
@@ -275,7 +245,6 @@ void Reallocate_Part2( reax_system const * const system,
Reallocate_List( lists[THREE_BODIES], realloc->total_bonds,
(int) CEIL( realloc->total_bonds * SAFE_ZONE ),
realloc->total_thbodies, TYP_THREE_BODY );
- Init_List_Indices( lists[THREE_BODIES] );
realloc->thbody = FALSE;
}
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index dfe47bb0c150426bda671abdb901367f89d939aa..9e653002539ff61940e84222524cf9d7df1108cc 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -1021,7 +1021,7 @@ static void Init_Bond_Full( reax_system const * const system,
if ( ihb == H_ATOM )
{
- ihb_top = End_Index( workspace->hbond_index[i], hbonds );
+ ihb_top = End_Index( i, hbonds );
}
else
{
@@ -1082,11 +1082,11 @@ static void Init_Bond_Full( reax_system const * const system,
{
if ( num_hbonds < hbonds->total_intrs )
{
- jhb_top = End_Index( workspace->hbond_index[j], hbonds );
+ jhb_top = End_Index( j, hbonds );
hbonds->hbond_list[jhb_top].nbr = i;
hbonds->hbond_list[jhb_top].scl = -1;
hbonds->hbond_list[jhb_top].ptr = nbr_pj;
- Set_End_Index( workspace->hbond_index[j], jhb_top + 1, hbonds );
+ Set_End_Index( j, jhb_top + 1, hbonds );
}
++num_hbonds;
}
@@ -1224,19 +1224,31 @@ static void Init_Bond_Full( reax_system const * const system,
Set_End_Index( i, btop_i, bonds );
if ( ihb == H_ATOM )
{
- Set_End_Index( workspace->hbond_index[i], ihb_top, hbonds );
+ Set_End_Index( i, ihb_top, hbonds );
}
}
- if ( num_bonds > bonds->total_intrs )
+ for ( i = 0; i < bonds->n; ++i )
{
- workspace->realloc.bonds = TRUE;
+ if ( Num_Entries( i, bonds ) > system->bonds[i] )
+ {
+ workspace->realloc.bonds = TRUE;
+ break;
+ }
+
}
- if ( control->hbond_cut > 0.0 && workspace->num_H > 0
- && num_hbonds > hbonds->total_intrs )
+ if ( control->hbond_cut > 0.0 && workspace->num_H > 0 )
{
- workspace->realloc.hbonds = TRUE;
+ for ( i = 0; i < hbonds->n; ++i )
+ {
+ if ( Num_Entries( i, hbonds ) > system->hbonds[i] )
+ {
+ workspace->realloc.hbonds = TRUE;
+ break;
+ }
+
+ }
}
}
@@ -1299,6 +1311,9 @@ static int Init_Forces( reax_system * const system,
rvec_MakeZero( workspace->dDeltap_self[i] );
}
+ Init_List_Indices( lists[BONDS], system->bonds );
+ Init_List_Indices( lists[HBONDS], system->hbonds );
+
// if ( lists[FAR_NBRS]->format == HALF_LIST )
// {
// Init_Bond_Half( system, control, workspace, lists );
@@ -1528,13 +1543,13 @@ static void Estimate_Storages_Bonds( reax_system const * const system,
}
workspace->realloc.total_thbodies = 0;
- for ( i = 0; i < system->N; ++i )
+ for ( i = 0; i < system->N_max; ++i )
{
//TODO: the factor of 2 depends on the bond list format (and also effects thbody list), revisit later
system->bonds[i] = MAX( (int) CEIL( system->bonds[i] * 2.0 * SAFE_ZONE ), MIN_BONDS );
}
- for ( i = 0; i < system->N; ++i )
+ for ( i = 0; i < system->N_max; ++i )
{
system->hbonds[i] = MAX( (int) CEIL( system->hbonds[i] * SAFE_HBONDS ), MIN_HBONDS );
}
diff --git a/sPuReMD/src/grid.c b/sPuReMD/src/grid.c
index 3c6eb47bdfe8348e4d821364254125e8d82df2d1..49accae1128f0f1c2a4875fb9baf1ade9bc364fe 100644
--- a/sPuReMD/src/grid.c
+++ b/sPuReMD/src/grid.c
@@ -730,15 +730,6 @@ static void Copy_Storage( reax_system const * const system,static_storage * cons
workspace->b_s[top] = -system->reax_param.sbp[ old_type ].chi;
workspace->b_t[top] = -1.0;
- if ( system->reax_param.sbp[ old_type ].p_hbond == H_ATOM )
- {
- workspace->hbond_index[top] = (*num_H)++;
- }
- else
- {
- workspace->hbond_index[top] = -1;
- }
-
rvec_Copy( f_old[top], workspace->f_old[old_id] );
}
diff --git a/sPuReMD/src/hydrogen_bonds.c b/sPuReMD/src/hydrogen_bonds.c
index ccbff7f6849b32e46a2b32c33c39aec7621b7661..14d30b2489cd11dceb1fce0838eb7242914477ae 100644
--- a/sPuReMD/src/hydrogen_bonds.c
+++ b/sPuReMD/src/hydrogen_bonds.c
@@ -95,8 +95,8 @@ void Hydrogen_Bonds( reax_system *system, control_params *control,
type_j = system->atoms[j].type;
start_j = Start_Index( j, bonds );
end_j = End_Index( j, bonds );
- hb_start_j = Start_Index( workspace->hbond_index[j], hbonds );
- hb_end_j = End_Index( workspace->hbond_index[j], hbonds );
+ hb_start_j = Start_Index( j, hbonds );
+ hb_end_j = End_Index( j, hbonds );
#if defined(_OPENMP)
f_j = &workspace->f_local[tid * system->N + j];
#else
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index cd9d60759b4f0263a7aac60f93d3e8db0e9539a0..81377aea620625d77249b022694953f54d0dfb45 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -345,10 +345,6 @@ static void Init_Workspace( reax_system * const system,
if ( realloc == TRUE )
{
- /* hydrogen bond list */
- workspace->hbond_index = smalloc( system->N_max * sizeof( int ),
- __FILE__, __LINE__ );
-
/* bond order related storage */
workspace->total_bond_order = smalloc( system->N_max * sizeof( real ),
__FILE__, __LINE__ );
@@ -903,16 +899,8 @@ static void Init_Lists( reax_system * const system,
{
if ( system->reax_param.sbp[ system->atoms[i].type ].p_hbond == H_ATOM )
{
- workspace->hbond_index[i] = workspace->num_H++;
+ workspace->num_H++;
}
- else
- {
- workspace->hbond_index[i] = -1;
- }
- }
- for ( i = system->N; i < system->N_max; ++i )
- {
- workspace->hbond_index[i] = -1;
}
}
@@ -920,34 +908,25 @@ static void Init_Lists( reax_system * const system,
{
if ( lists[HBONDS]->allocated == FALSE )
{
- workspace->num_H_max = (int) CEIL( SAFE_ZONE * workspace->num_H );
-
- Make_List( workspace->num_H, workspace->num_H_max,
+ Make_List( system->N, system->N_max,
(int) CEIL( SAFE_ZONE * workspace->realloc.total_hbonds ),
TYP_HBOND, lists[HBONDS] );
}
- else if ( workspace->num_H_max < workspace->num_H
+ else if ( system->N_max < system->N
|| lists[HBONDS]->total_intrs < workspace->realloc.total_hbonds )
{
- if ( workspace->num_H_max < workspace->num_H )
- {
- workspace->num_H_max = (int) CEIL( SAFE_ZONE * workspace->num_H );
- }
-
if ( lists[HBONDS]->allocated == TRUE )
{
Delete_List( TYP_HBOND, lists[HBONDS] );
}
- Make_List( workspace->num_H, workspace->num_H_max,
+ Make_List( system->N, system->N_max,
MAX( workspace->realloc.total_hbonds, lists[HBONDS]->total_intrs ),
TYP_HBOND, lists[HBONDS] );
}
else
{
- lists[HBONDS]->n = workspace->num_H;
+ lists[HBONDS]->n = system->N;
}
-
- Init_List_Indices( lists[HBONDS] );
}
/* bonds list */
@@ -971,8 +950,6 @@ static void Init_Lists( reax_system * const system,
lists[BONDS]->n = system->N;
}
- Init_List_Indices( lists[BONDS] );
-
/* 3bodies list */
if ( lists[THREE_BODIES]->allocated == FALSE )
{
@@ -1407,7 +1384,6 @@ static void Finalize_Workspace( reax_system *system, control_params *control,
{
int i;
- sfree( workspace->hbond_index, __FILE__, __LINE__ );
sfree( workspace->total_bond_order, __FILE__, __LINE__ );
sfree( workspace->Deltap, __FILE__, __LINE__ );
sfree( workspace->Deltap_boc, __FILE__, __LINE__ );
diff --git a/sPuReMD/src/list.c b/sPuReMD/src/list.c
index cd20c0894c3017310f8859285522542a5caa32cc..9731fd657e0368404fd0d47eeb620c79872bf52e 100644
--- a/sPuReMD/src/list.c
+++ b/sPuReMD/src/list.c
@@ -224,17 +224,23 @@ void Delete_List( int type, reax_list * const l )
/* Initialize list indices
*
* l: pointer to list
+ * max_intrs: max. num. of interactions for each list element
* */
-void Init_List_Indices( reax_list * const l )
+void Init_List_Indices( reax_list * const l, int * const max_intrs )
{
int i;
assert( l != NULL );
assert( l->n > 0 );
+ assert( max_intrs > 0 );
- for ( i = 0; i < l->n; ++i )
+ /* exclusive prefix sum of max_intrs replaces start indices,
+ * set end indices to the same as start indices for safety */
+ Set_Start_Index( 0, 0, l );
+ Set_End_Index( 0, 0, l );
+ for ( i = 1; i < l->n; ++i )
{
- Set_Start_Index( i, 0, l );
- Set_End_Index( i, 0, l );
+ Set_Start_Index( i, Start_Index( i - 1, l ) + max_intrs[i - 1], l );
+ Set_End_Index( i, Start_Index( i, l ), l );
}
}
diff --git a/sPuReMD/src/list.h b/sPuReMD/src/list.h
index 43b5dca66090c84571bbb2bbe94ca784049860fa..1a8948b0aa03a2a8c534eab99c00cd6c1ecddbce 100644
--- a/sPuReMD/src/list.h
+++ b/sPuReMD/src/list.h
@@ -31,7 +31,7 @@ void Make_List( int, int, int, int, reax_list * const );
/* Deallocates any space allocated for a reax_list instance */
void Delete_List( int, reax_list * const );
-void Init_List_Indices( reax_list * const );
+void Init_List_Indices( reax_list * const, int * const );
/* Return the size for the i-th internal list
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index 7cf2fabfb5f1ca0c96145fbc07210d42d2a5071a..4bf2c52ae79da6dfab61d47bacf3a6be4b8b433b 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -1641,10 +1641,6 @@ struct static_storage
/* num. hydrogen atoms for this simulation */
int num_H;
- /* max. num. hydrogen atoms across all simulations */
- int num_H_max;
- /* for hydrogen bonds */
- int *hbond_index;
rvec *a; // used in integrators
rvec *f_old;
diff --git a/sPuReMD/src/reset_tools.c b/sPuReMD/src/reset_tools.c
index 6756bb879e2e383b23f458c8f12c9922bb7bced6..ec58b9fd1b48d35dbf15df5dfecf694c94fac8ca 100644
--- a/sPuReMD/src/reset_tools.c
+++ b/sPuReMD/src/reset_tools.c
@@ -178,60 +178,6 @@ void Reset_Workspace( reax_system *system, static_storage *workspace )
}
-void Reset_Neighbor_Lists( reax_system *system, control_params *control,
- static_storage *workspace, reax_list **lists )
-{
- int i, tmp;
- reax_list *bonds, *hbonds;
-#if defined(TEST_FORCES)
- reax_list *dDeltas, *dBOs;
-#endif
-
- bonds = lists[BONDS];
- hbonds = lists[HBONDS];
-#if defined(TEST_FORCES)
- dDeltas = lists[DDELTA];
- dBOs = lists[DBO];
-#endif
-
- for ( i = 0; i < bonds->n; ++i )
- {
- tmp = Start_Index( i, bonds );
- Set_End_Index( i, tmp, bonds );
- }
-
- if ( control->hbond_cut > 0.0 && workspace->num_H > 0 )
- {
- for ( i = 0; i < system->N; ++i )
- {
- if ( system->reax_param.sbp[system->atoms[i].type].p_hbond == H_ATOM )
- {
- tmp = Start_Index( workspace->hbond_index[i], hbonds );
- Set_End_Index( workspace->hbond_index[i], tmp, hbonds );
-
-// fprintf( stderr, "i:%d, hbond: %d-%d\n",
-// i, Start_Index( workspace->hbond_index[i], hbonds ),
-// End_Index( workspace->hbond_index[i], hbonds ) );
- }
- }
- }
-
-#if defined(TEST_FORCES)
- for ( i = 0; i < dDeltas->n; ++i )
- {
- tmp = Start_Index( i, dDeltas );
- Set_End_Index( i, tmp, dDeltas );
- }
-
- for ( i = 0; i < dBOs->n; ++i )
- {
- tmp = Start_Index( i, dBOs );
- Set_End_Index( i, tmp, dBOs );
- }
-#endif
-}
-
-
void Reset( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace, reax_list **lists )
{
@@ -246,8 +192,6 @@ void Reset( reax_system *system, control_params *control,
}
Reset_Workspace( system, workspace );
-
- Reset_Neighbor_Lists( system, control, workspace, lists );
}
diff --git a/sPuReMD/src/reset_tools.h b/sPuReMD/src/reset_tools.h
index f725d356880df739018a0fe1a7291a0062247bce..e5f7474273ee913843bb3b2d4ec3f4cdf2dcf346 100644
--- a/sPuReMD/src/reset_tools.h
+++ b/sPuReMD/src/reset_tools.h
@@ -33,9 +33,6 @@ void Reset_Energies( simulation_data* );
void Reset_Workspace( reax_system*, static_storage* );
-void Reset_Neighbor_Lists( reax_system*, control_params*,
- static_storage*, reax_list** );
-
void Reset( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list** );
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index 91ea395f8525cc4ad11d2f729f2449e268b64711..f8cb8ec6e14520ef7501638da087dc5c8430eb30 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -386,6 +386,10 @@ int simulate( const void * const handle )
spmd_handle->workspace->realloc.bonds
|| spmd_handle->workspace->realloc.hbonds );
+ Reallocate_Part2( spmd_handle->system, spmd_handle->control,
+ spmd_handle->data, spmd_handle->workspace,
+ spmd_handle->lists );
+
ret_forces = Compute_Forces( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
spmd_handle->realloc );