diff --git a/sPuReMD/src/analyze.c b/sPuReMD/src/analyze.c
index fa83135365692fc8be8e90c8a86ca0b69ca4276b..33e5e1cdf0644add0906051fd8de36fa5002b7bd 100644
--- a/sPuReMD/src/analyze.c
+++ b/sPuReMD/src/analyze.c
@@ -337,6 +337,7 @@ static void Analyze_Molecules( reax_system *system, control_params *control,
}
+#if defined(DEBUG_FOCUS)
static void Report_Bond_Change( reax_system *system, control_params *control,
static_storage *workspace, reax_list *old_bonds,
reax_list *new_bonds, int a1, int a2, int flag, FILE *fout )
@@ -436,9 +437,11 @@ static void Report_Bond_Change( reax_system *system, control_params *control,
fprintf( fout, "\n" );
}
}
+#endif
/* ASSUMPTION: Bond lists are sorted */
+#if defined(DEBUG_FOCUS)
static void Compare_Bonding( int atom, reax_system *system, control_params *control,
static_storage *workspace, reax_list *old_bonds,
reax_list *new_bonds, FILE *fout )
@@ -563,6 +566,7 @@ static void Compare_Bonding( int atom, reax_system *system, control_params *cont
}
}
}
+#endif
static void Visit_Bonds( int atom, int *mark, int *type, reax_system *system,
@@ -642,7 +646,7 @@ static void Analyze_Fragments( reax_system *system, control_params *control,
{
/* it is a new one, add to the fragments list */
strncpy( fragments[num_fragment_types], fragment_out,
- sizeof(fragments[num_fragment_types]) - 1 );
+ sizeof(fragments[num_fragment_types]) );
fragments[num_fragment_types][sizeof(fragments[num_fragment_types]) - 1] = '\0';
fragment_count[num_fragment_types] = 1;
++num_fragment_types;
@@ -671,6 +675,7 @@ static void Analyze_Fragments( reax_system *system, control_params *control,
}
+#if defined(DEBUG_FOCUS)
static void Analyze_Silica( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list **lists, FILE *fout )
@@ -798,6 +803,7 @@ static void Analyze_Silica( reax_system *system, control_params *control,
fflush( fout );
}
+#endif
static int Get_Type_of_Molecule( molecule *m )
@@ -932,6 +938,7 @@ static void Calculate_Drift( reax_system *system, control_params *control,
}
+#if defined(DEBUG_FOCUS)
static void Calculate_Density_3DMesh( reax_system *system, simulation_data *data,
FILE *fout )
{
@@ -1009,8 +1016,10 @@ static void Calculate_Density_3DMesh( reax_system *system, simulation_data *data
fprintf( stderr, "AMU_to_GRAM : %g\n", AMU_to_GRAM );
fprintf( stderr, "ANG_to_CM : %g\n", ANG_to_CM );
}
+#endif
+#if defined(DEBUG_FOCUS)
static void Calculate_Density_Slice( reax_system *system, simulation_data *data, FILE *fout )
{
real slice_thickness = 0.5;
@@ -1049,11 +1058,12 @@ static void Calculate_Density_Slice( reax_system *system, simulation_data *data,
}
}
}
+#endif
void Analysis( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- reax_list **lists, output_controls *out_control )
+ simulation_data *data, static_storage *workspace,
+ reax_list **lists, output_controls *out_control )
{
int steps;
@@ -1061,6 +1071,12 @@ void Analysis( reax_system *system, control_params *control,
if ( steps == 1 )
{
+ if ( lists[OLD_BONDS]->allocated == FALSE )
+ {
+ Make_List( lists[BONDS]->n, lists[BONDS]->n_max, lists[BONDS]->total_intrs,
+ TYP_BOND, lists[OLD_BONDS] );
+ }
+
if ( control->molec_anal == REACTIONS )
{
Copy_Bond_List( system, control, lists );
@@ -1072,25 +1088,26 @@ void Analysis( reax_system *system, control_params *control,
}
/****** Molecular Analysis ******/
- if ( control->molec_anal && steps % control->freq_molec_anal == 0 )
+ if ( control->molec_anal
+ && steps % control->freq_molec_anal == 0 )
{
if ( control->molec_anal == FRAGMENTS )
{
/* discover molecules */
Analyze_Fragments( system, control, data, workspace, lists,
- out_control->mol, 0 );
+ out_control->mol, 0 );
/* discover fragments without the ignored atoms */
if ( control->num_ignored )
{
Analyze_Fragments( system, control, data, workspace, lists,
- out_control->ign, 1 );
+ out_control->ign, 1 );
}
}
else if ( control->molec_anal == REACTIONS )
{
/* discover molecular changes - reactions */
- Analyze_Molecules( system, control, data, workspace,
- lists, out_control->mol );
+ Analyze_Molecules( system, control, data, workspace, lists,
+ out_control->mol );
}
}
diff --git a/sPuReMD/src/bond_orders.c b/sPuReMD/src/bond_orders.c
index 294dac11b06dc768464b7ccd09666000388dde7f..bc60d19321ad3f0715bfa5817756edfe615a9390 100644
--- a/sPuReMD/src/bond_orders.c
+++ b/sPuReMD/src/bond_orders.c
@@ -1072,10 +1072,10 @@ void BO( reax_system *system, control_params *control,
workspace->vlpex[j] = workspace->Delta_e[j]
- 2.0 * (int)(workspace->Delta_e[j] / 2.0);
- explp1 = EXP(-p_lp1 * SQR(2.0 + workspace->vlpex[j]));
+ explp1 = EXP( -1.0 * p_lp1 * SQR(2.0 + workspace->vlpex[j]) );
workspace->nlp[j] = explp1 - (int)(workspace->Delta_e[j] / 2.0);
workspace->Delta_lp[j] = sbp_j->nlp_opt - workspace->nlp[j];
- workspace->Clp[j] = 2.0 * p_lp1 * explp1 * (2.0 + workspace->vlpex[j]);
+ workspace->Clp[j] = 2.0 * p_lp1 * (2.0 + workspace->vlpex[j]) * explp1;
/* Adri uses different dDelta_lp values than the ones in notes... */
workspace->dDelta_lp[j] = workspace->Clp[j];
// workspace->dDelta_lp[j] = workspace->Clp[j] + (0.5 - workspace->Clp[j])
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index f3311c1a185370e9bc8d046d2d85a0043de50c42..998357a81743362acc038fceddef5dcfb7f33b11 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -1518,6 +1518,9 @@ static void EE( reax_system * const system, control_params * const control,
workspace->mask_qmmm[i] = system->atoms[i].qmmm_mask;
}
workspace->mask_qmmm[system->N_cm - 1] = 1;
+
+ /* Mask the b vector as well */
+ Vector_Mask_qmmm( workspace->b_s, workspace->mask_qmmm, system->N_cm );
#endif
switch ( control->cm_solver_type )
@@ -1646,8 +1649,11 @@ static void ACKS2( reax_system * const system, control_params * const control,
{
workspace->mask_qmmm[i] = system->atoms[i - system->N].qmmm_mask;
}
- workspace->mask_qmmm[2 * system->N_cm] = 1;
- workspace->mask_qmmm[2 * system->N_cm + 1] = 1;
+ workspace->mask_qmmm[2 * system->N] = 1;
+ workspace->mask_qmmm[2 * system->N + 1] = 1;
+
+ /* Mask the b vector as well */
+ Vector_Mask_qmmm( workspace->b_s, workspace->mask_qmmm, system->N_cm );
#endif
switch ( control->cm_solver_type )
diff --git a/sPuReMD/src/ffield.c b/sPuReMD/src/ffield.c
index 569b6844168fa3ce2b2d4dd09ebea28018c65853..39d78b2026a56c217ca3ab32dc55e8bd93f09f89 100644
--- a/sPuReMD/src/ffield.c
+++ b/sPuReMD/src/ffield.c
@@ -257,7 +257,8 @@ void Read_Force_Field( const char * const ffield_file,
fgets( s, MAX_LINE, fp );
Tokenize( s, &tmp, MAX_TOKEN_LEN );
- strncpy( reax->sbp[i].name, tmp[0], sizeof(reax->sbp[i].name) );
+ strncpy( reax->sbp[i].name, tmp[0], sizeof(reax->sbp[i].name) - 1 );
+ reax->sbp[i].name[sizeof(reax->sbp[i].name) - 1] = '\0';
for ( j = 0; j < strlen( reax->sbp[i].name ); ++j )
{
reax->sbp[i].name[j] = toupper( reax->sbp[i].name[j] );
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 6e2fc638af30e704e420b227204fb5e6d123a471..4b42eaa0bf677acd230ee6e337445aee275ed382 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -493,49 +493,101 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system *system,
break;
case EE_CM:
- H->start[system->N_cm - 1] = *Htop;
- H_sp->start[system->N_cm - 1] = *H_sp_top;
-
- for ( i = 0; i < system->N_cm - 1; ++i )
+ if ( system->num_molec_charge_constraints == 0 )
{
-#if defined(QMMM)
- // total charge constraint on QM atoms
- if ( system->atoms[i].qmmm_mask == TRUE )
+ H->start[system->N_cm - 1] = *Htop;
+ H_sp->start[system->N_cm - 1] = *H_sp_top;
+
+ for ( i = 0; i < system->N_cm - 1; ++i )
{
+#if defined(QMMM)
+ /* total charge constraint on QM atoms */
+ if ( system->atoms[i].qmmm_mask == TRUE )
+ {
+ H->j[*Htop] = i;
+ H->val[*Htop] = 1.0;
+
+ H_sp->j[*H_sp_top] = i;
+ H_sp->val[*H_sp_top] = 1.0;
+ }
+ else
+ {
+ H->j[*Htop] = i;
+ H->val[*Htop] = 0.0;
+
+ H_sp->j[*H_sp_top] = i;
+ H_sp->val[*H_sp_top] = 0.0;
+ }
+#else
H->j[*Htop] = i;
H->val[*Htop] = 1.0;
H_sp->j[*H_sp_top] = i;
H_sp->val[*H_sp_top] = 1.0;
- }
- else
- {
- H->j[*Htop] = i;
- H->val[*Htop] = 0.0;
+#endif
- H_sp->j[*H_sp_top] = i;
- H_sp->val[*H_sp_top] = 0.0;
+ *Htop = *Htop + 1;
+ *H_sp_top = *H_sp_top + 1;
}
+
+ H->j[*Htop] = system->N_cm - 1;
+ H->val[*Htop] = 0.0;
*Htop = *Htop + 1;
+
+ H_sp->j[*H_sp_top] = system->N_cm - 1;
+ H_sp->val[*H_sp_top] = 0.0;
*H_sp_top = *H_sp_top + 1;
+ }
+ else
+ {
+ for ( i = 0; i < system->num_molec_charge_constraints; ++i )
+ {
+ H->start[system->N + i] = *Htop;
+ H_sp->start[system->N + i] = *H_sp_top;
+
+ for ( j = system->molec_charge_constraint_ranges[2 * i];
+ j <= system->molec_charge_constraint_ranges[2 * i + 1]; ++j )
+ {
+#if defined(QMMM)
+ /* molecule charge constraint on QM atoms */
+ if ( system->atoms[i].qmmm_mask == TRUE )
+ {
+ H->j[*Htop] = j - 1;
+ H->val[*Htop] = 1.0;
+
+ H_sp->j[*H_sp_top] = j - 1;
+ H_sp->val[*H_sp_top] = 1.0;
+ }
+ else
+ {
+ H->j[*Htop] = j - 1;
+ H->val[*Htop] = 0.0;
+
+ H_sp->j[*H_sp_top] = j - 1;
+ H_sp->val[*H_sp_top] = 0.0;
+ }
#else
- H->j[*Htop] = i;
- H->val[*Htop] = 1.0;
- *Htop = *Htop + 1;
+ H->j[*Htop] = j - 1;
+ H->val[*Htop] = 1.0;
- H_sp->j[*H_sp_top] = i;
- H_sp->val[*H_sp_top] = 1.0;
- *H_sp_top = *H_sp_top + 1;
+ H_sp->j[*H_sp_top] = j - 1;
+ H_sp->val[*H_sp_top] = 1.0;
#endif
- }
- H->j[*Htop] = system->N_cm - 1;
- H->val[*Htop] = 0.0;
- *Htop = *Htop + 1;
+ *Htop = *Htop + 1;
+ *H_sp_top = *H_sp_top + 1;
+ }
- H_sp->j[*H_sp_top] = system->N_cm - 1;
- H_sp->val[*H_sp_top] = 0.0;
- *H_sp_top = *H_sp_top + 1;
+ /* explicit zeros on diagonals */
+ H->j[*Htop] = system->N + i;
+ H->val[*Htop] = 0.0;
+ *Htop = *Htop + 1;
+
+ H_sp->j[*H_sp_top] = system->N + i;
+ H_sp->val[*H_sp_top] = 0.0;
+ *H_sp_top = *H_sp_top + 1;
+ }
+ }
break;
case ACKS2_CM:
@@ -793,6 +845,11 @@ static void Init_Forces( reax_system *system, control_params *control,
flag = FALSE;
flag_sp = FALSE;
+#if defined(QMMM)
+ if ( system->atoms[i].qmmm_mask == TRUE
+ || system->atoms[j].qmmm_mask == TRUE )
+ {
+#endif
/* check if reneighboring step --
* atomic distances just computed via
* Verlet list, so use current distances */
@@ -1024,11 +1081,14 @@ static void Init_Forces( reax_system *system, control_params *control,
Set_End_Index( j, btop_j + 1, bonds );
}
+#if defined(QMMM)
+ }
+#endif
}
- }
#if defined(QMMM)
}
#endif
+ }
}
/* diagonal entry */
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index 3b361e20f487a213891a2f1975c319421d617459..32e0fe573c97062b9b89904fd4b2f37f2ce6c132 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -359,7 +359,7 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
char s_x[9], s_y[9], s_z[9];
char occupancy[7], temp_factor[7];
char seg_id[5], element[3], charge[3];
- char alt_loc, chain_id, icode;
+// char alt_loc, chain_id, icode;
rvec x;
reax_atom *atom;
@@ -409,13 +409,13 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
serial[sizeof(serial) - 1] = '\0';
strncpy( atom_name, s1 + 12, sizeof(atom_name) - 1 );
atom_name[sizeof(atom_name) - 1] = '\0';
- alt_loc = s1[16];
+// alt_loc = s1[16];
strncpy( res_name, s1 + 17, sizeof(res_name) - 1 );
res_name[sizeof(res_name) - 1] = '\0';
- chain_id = s1[21];
+// chain_id = s1[21];
strncpy( res_seq, s1 + 22, sizeof(res_seq) - 1 );
res_seq[sizeof(res_seq) - 1] = '\0';
- icode = s1[26];
+// icode = s1[26];
strncpy( s_x, s1 + 30, sizeof(s_x) - 1 );
s_x[sizeof(s_x) - 1] = '\0';
strncpy( s_y, s1 + 38, sizeof(s_y) - 1 );
@@ -441,13 +441,13 @@ void Read_PDB( const char * const pdb_file, reax_system* system, control_params
serial[sizeof(serial) - 1] = '\0';
strncpy( atom_name, s1 + 12, sizeof(atom_name) - 1 );
atom_name[sizeof(atom_name) - 1] = '\0';
- alt_loc = s1[16];
+// alt_loc = s1[16];
strncpy( res_name, s1 + 17, sizeof(res_name) - 1 );
res_name[sizeof(res_name) - 1] = '\0';
- chain_id = s1[21];
+// chain_id = s1[21];
strncpy( res_seq, s1 + 22, sizeof(res_seq) - 1 );
res_seq[sizeof(res_seq) - 1] = '\0';
- icode = s1[26];
+// icode = s1[26];
strncpy( s_x, s1 + 30, sizeof(s_x) - 1 );
s_x[sizeof(s_x) - 1] = '\0';
strncpy( s_y, s1 + 38, sizeof(s_y) - 1 );
@@ -589,7 +589,12 @@ void Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
/ (system->box.box_norms[2] * system->box.box_norms[1]) );
/* write header */
- snprintf( fname, sizeof(fname), "%s-%d.pdb", control->sim_name, data->step );
+ snprintf( fname, sizeof(fname) - 1, "%.*s-%.*d.pdb",
+ (int) MIN( strnlen(control->sim_name,
+ sizeof(fname) - snprintf(NULL, 0, "%d", data->step) - 6 ),
+ sizeof(control->sim_name) ),
+ control->sim_name, snprintf(NULL, 0, "%d", data->step), data->step );
+ fname[sizeof(fname) - 1] = '\0';
pdb = sfopen( fname, "w" );
fprintf( pdb, PDB_CRYST1_FORMAT_O,
"CRYST1",
@@ -625,6 +630,13 @@ void Write_PDB( reax_system* system, reax_list* bonds, simulation_data *data,
}
+/* Parser for geometry files in BGF format
+ *
+ * system: struct containing atom-related information
+ * control: struct containing simulation parameters
+ * data: struct containing information on active simulations
+ * workspace: struct containing intermediate structures used for calculations
+ * */
void Read_BGF( const char * const bgf_file, reax_system* system, control_params *control,
simulation_data *data, static_storage *workspace )
{
@@ -636,9 +648,9 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
char s_x[11], s_y[11], s_z[11];
char occupancy[4], temp_factor[3];
char element[6], charge[9];
- char chain_id;
+// char chain_id;
char s_a[12], s_b[12], s_c[12], s_alpha[12], s_beta[12], s_gamma[12];
- int i, n, atom_cnt, token_cnt, bgf_serial, ratom, crystx_found;
+ int i, n, num_mcc, atom_cnt, token_cnt, bgf_serial, ratom, crystx_found;
rvec x;
ratom = 0;
@@ -651,6 +663,7 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
/* count number of atoms in the BGF file */
n = 0;
+ num_mcc = 0;
line[0] = 0;
while ( fgets( line, MAX_LINE, bgf ) )
@@ -663,6 +676,10 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
{
++n;
}
+ else if ( strncmp( tokens[0], "MOLCHARGE", 9 ) == 0 )
+ {
+ ++num_mcc;
+ }
line[0] = 0;
}
@@ -672,20 +689,50 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
exit( INVALID_INPUT );
}
- sfclose( bgf, "Read_BGF::bgf" );
-
if ( system->prealloc_allocated == FALSE || n > system->N_max )
{
PreAllocate_Space( system, control, workspace, (int) CEIL( SAFE_ZONE * n ) );
+
+ if ( system->prealloc_allocated == FALSE && num_mcc > 0 )
+ {
+ system->molec_charge_constraints = smalloc(
+ sizeof(real) * num_mcc, "Read_BGF::molec_charge_constraints" );
+ system->molec_charge_constraint_ranges = smalloc(
+ sizeof(int) * 2 * num_mcc, "Read_BGF::molec_charge_constraint_ranges" );
+
+ system->max_num_molec_charge_constraints = num_mcc;
+ }
+ else if ( num_mcc > system->max_num_molec_charge_constraints )
+ {
+ if ( system->max_num_molec_charge_constraints > 0 )
+ {
+ sfree( system->molec_charge_constraints, "Read_BGF::molec_charge_constraints" );
+ sfree( system->molec_charge_constraint_ranges, "Read_BGF::molec_charge_constraint_ranges" );
+ }
+
+ system->molec_charge_constraints = smalloc(
+ sizeof(real) * num_mcc, "Read_BGF::molec_charge_constraints" );
+ system->molec_charge_constraint_ranges = smalloc(
+ sizeof(int) * 2 * num_mcc, "Read_BGF::molec_charge_constraint_ranges" );
+
+ system->max_num_molec_charge_constraints = num_mcc;
+ }
}
system->N = n;
+ system->num_molec_charge_constraints = num_mcc;
+ num_mcc = 0;
+
+#if defined(DEBUG_FOCUS)
+ fprintf( stderr, "[INFO] num_atoms = %d, num_mcc = %d\n", system->N,
+ system->num_molec_charge_constraints );
+#endif
for ( i = 0; i < MAX_ATOM_ID; ++i )
{
workspace->map_serials[i] = -1;
}
- bgf = sfopen( bgf_file, "r" );
+ fseek( bgf, 0, SEEK_SET );
atom_cnt = 0;
token_cnt = 0;
@@ -748,7 +795,7 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
atom_name[sizeof(atom_name) - 1] = '\0';
strncpy( res_name, backup + 19, sizeof(res_name) - 1 );
res_name[sizeof(res_name) - 1] = '\0';
- chain_id = backup[23];
+// chain_id = backup[23];
strncpy( res_seq, backup + 25, sizeof(res_seq) - 1 );
res_seq[sizeof(res_seq) - 1] = '\0';
strncpy( s_x, backup + 30, sizeof(s_x) - 1 );
@@ -838,6 +885,24 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
}
}
}
+ else if ( strncmp( tokens[0], "MOLCHARGE", 9 ) == 0 )
+ {
+ assert( token_cnt == 4 );
+
+ system->molec_charge_constraint_ranges[2 * num_mcc] = sstrtol( tokens[1], __FILE__, __LINE__ );
+ system->molec_charge_constraint_ranges[2 * num_mcc + 1] = sstrtol( tokens[2], __FILE__, __LINE__ );
+ system->molec_charge_constraints[num_mcc] = sstrtod( tokens[3], __FILE__, __LINE__ );
+
+#if defined(DEBUG_FOCUS)
+ fprintf( stderr,
+ "[INFO] num_mcc = %d, mcc = %f, mcc_range = (%d, %d)\n", num_mcc,
+ system->molec_charge_constraints[num_mcc],
+ system->molec_charge_constraint_ranges[2 * num_mcc],
+ system->molec_charge_constraint_ranges[2 * num_mcc + 1] );
+#endif
+
+ ++num_mcc;
+ }
/* clear previous input line */
line[0] = '\0';
diff --git a/sPuReMD/src/grid.c b/sPuReMD/src/grid.c
index 25f9591c29d6625cb0bb587c8d79b15313040bc3..ad4e4fb6feff108d44825fb9bd4f86c905c06cae 100644
--- a/sPuReMD/src/grid.c
+++ b/sPuReMD/src/grid.c
@@ -26,7 +26,7 @@
#include "vector.h"
-static int Estimate_GCell_Population( reax_system* system )
+static int Estimate_GCell_Population( reax_system * const system )
{
int i, j, k, l;
int max_atoms;
@@ -69,7 +69,7 @@ static int Estimate_GCell_Population( reax_system* system )
}
-static void Allocate_Space_for_Grid( reax_system *system, int alloc )
+static void Allocate_Space_for_Grid( reax_system * const system, int alloc )
{
int i, j, k, l, max_atoms;
grid *g;
@@ -256,7 +256,7 @@ static void Allocate_Space_for_Grid( reax_system *system, int alloc )
}
-static void Deallocate_Grid_Space( grid *g )
+static void Deallocate_Grid_Space( grid * const g )
{
int i, j, k;
@@ -305,7 +305,7 @@ static void Deallocate_Grid_Space( grid *g )
}
-static inline int Shift( int p, int dp, int dim, grid *g )
+static inline int Shift( int p, int dp, int dim, grid const * const g )
{
int dim_len, newp;
@@ -339,7 +339,7 @@ static inline int Shift( int p, int dp, int dim, grid *g )
/* finds the closest point between two grid cells denoted by c1 and c2.
* periodic boundary conditions are taken into consideration as well. */
-static void Find_Closest_Point( grid *g, int c1x, int c1y, int c1z,
+static void Find_Closest_Point( grid const * const g, int c1x, int c1y, int c1z,
int c2x, int c2y, int c2z, rvec closest_point )
{
int i, d;
@@ -385,7 +385,7 @@ static void Find_Closest_Point( grid *g, int c1x, int c1y, int c1z,
}
-static void Find_Neighbor_Grid_Cells( grid *g )
+static void Find_Neighbor_Grid_Cells( grid * const g )
{
int i, j, k;
int di, dj, dk;
@@ -460,7 +460,7 @@ static void Find_Neighbor_Grid_Cells( grid *g )
}
-void Setup_Grid( reax_system* system )
+void Setup_Grid( reax_system * const system )
{
int d, alloc;
grid *g;
@@ -530,7 +530,7 @@ void Setup_Grid( reax_system* system )
}
-void Update_Grid( reax_system* system )
+void Update_Grid( reax_system * const system )
{
int d, i, j, k, x, y, z, itr;
ivec ncell;
@@ -614,7 +614,7 @@ void Update_Grid( reax_system* system )
}
-void Bin_Atoms( reax_system* system, static_storage *workspace )
+void Bin_Atoms( reax_system * const system, static_storage * const workspace )
{
int i, j, k, l;
int max_atoms;
@@ -666,7 +666,7 @@ void Bin_Atoms( reax_system* system, static_storage *workspace )
}
-void Finalize_Grid( reax_system* system )
+void Finalize_Grid( reax_system * const system )
{
if ( system->g.allocated == TRUE )
{
@@ -675,7 +675,7 @@ void Finalize_Grid( reax_system* system )
}
-static inline void reax_atom_Copy( reax_atom *dest, reax_atom *src )
+static void reax_atom_Copy( reax_atom * const dest, reax_atom * const src )
{
dest->type = src->type;
strncpy( dest->name, src->name, sizeof(dest->name) - 1 );
@@ -691,9 +691,10 @@ static inline void reax_atom_Copy( reax_atom *dest, reax_atom *src )
}
-static void Copy_Storage( reax_system *system, static_storage *workspace,
- control_params *control, int top, int old_id, int old_type, int *num_H,
- real **v, real **s, real **t, int *orig_id, rvec *f_old )
+static void Copy_Storage( reax_system const * const system,static_storage * const workspace,
+ control_params const * const control, int top, int old_id, int old_type,
+ int * const num_H, real ** const v, real ** const s, real ** const t,
+ int * const orig_id, rvec * const f_old )
{
int i;
@@ -726,7 +727,8 @@ static void Copy_Storage( reax_system *system, static_storage *workspace,
}
-static void Free_Storage( static_storage *workspace, control_params * control )
+static void Free_Storage( static_storage * const workspace,
+ control_params const * const control )
{
int i;
@@ -759,8 +761,9 @@ static void Assign_New_Storage( static_storage *workspace,
}
-void Cluster_Atoms( reax_system *system, static_storage *workspace,
- control_params *control )
+/* Reorder atom list to improve cache performance */
+void Reorder_Atoms( reax_system * const system, static_storage * const workspace,
+ control_params const * const control )
{
int i, j, k, l, top, old_id, num_H;
reax_atom *old_atom, *new_atoms;
@@ -774,23 +777,23 @@ void Cluster_Atoms( reax_system *system, static_storage *workspace,
top = 0;
g = &system->g;
- new_atoms = scalloc( system->N, sizeof(reax_atom), "Cluster_Atoms::new_atoms" );
- orig_id = scalloc( system->N, sizeof(int), "Cluster_Atoms::orig_id" );
- f_old = scalloc( system->N, sizeof(rvec), "Cluster_Atoms::f_old" );
+ new_atoms = scalloc( system->N, sizeof(reax_atom), "Reorder_Atoms::new_atoms" );
+ orig_id = scalloc( system->N, sizeof(int), "Reorder_Atoms::orig_id" );
+ f_old = scalloc( system->N, sizeof(rvec), "Reorder_Atoms::f_old" );
- s = scalloc( 5, sizeof(real *), "Cluster_Atoms::s" );
- t = scalloc( 5, sizeof(real *), "Cluster_Atoms::t" );
+ s = scalloc( 5, sizeof(real *), "Reorder_Atoms::s" );
+ t = scalloc( 5, sizeof(real *), "Reorder_Atoms::t" );
for ( i = 0; i < 5; ++i )
{
- s[i] = scalloc( system->N_cm, sizeof(real), "Cluster_Atoms::s[i]" );
- t[i] = scalloc( system->N_cm, sizeof(real), "Cluster_Atoms::t[i]" );
+ s[i] = scalloc( system->N_cm, sizeof(real), "Reorder_Atoms::s[i]" );
+ t[i] = scalloc( system->N_cm, sizeof(real), "Reorder_Atoms::t[i]" );
}
v = scalloc( control->cm_solver_restart + 1, sizeof(real *),
- "Cluster_Atoms::v" );
+ "Reorder_Atoms::v" );
for ( i = 0; i < control->cm_solver_restart + 1; ++i )
{
- v[i] = scalloc( system->N_cm, sizeof(real), "Cluster_Atoms::v[i]" );
+ v[i] = scalloc( system->N_cm, sizeof(real), "Reorder_Atoms::v[i]" );
}
for ( i = 0; i < g->ncell[0]; i++ )
@@ -818,7 +821,7 @@ void Cluster_Atoms( reax_system *system, static_storage *workspace,
}
}
- sfree( system->atoms, "Cluster_Atoms::system->atoms" );
+ sfree( system->atoms, "Reorder_Atoms::system->atoms" );
Free_Storage( workspace, control );
system->atoms = new_atoms;
diff --git a/sPuReMD/src/grid.h b/sPuReMD/src/grid.h
index d7e9a1986af74e0568010d1650c3a53b4209e8c8..1be1a4d7a852145a4a34368635f5c51ee2ba2ce3 100644
--- a/sPuReMD/src/grid.h
+++ b/sPuReMD/src/grid.h
@@ -25,15 +25,16 @@
#include "reax_types.h"
-void Setup_Grid( reax_system * );
+void Setup_Grid( reax_system * const );
-void Update_Grid( reax_system * );
+void Update_Grid( reax_system * const );
-void Bin_Atoms( reax_system *, static_storage * );
+void Bin_Atoms( reax_system * const, static_storage * const );
-void Finalize_Grid( reax_system * );
+void Finalize_Grid( reax_system * const );
-void Cluster_Atoms( reax_system *, static_storage *, control_params * );
+void Reorder_Atoms( reax_system * const, static_storage * const,
+ control_params const * const );
#endif
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index dd1b8291313824d2ef3cf944dd64af789d6eef48..7c9fc2e97302986609bd231604c77be11d58dd72 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -377,8 +377,16 @@ static void Init_Workspace( reax_system *system, control_params *control,
system->N_cm_max = system->N_max;
break;
case EE_CM:
- system->N_cm = system->N + 1;
- system->N_cm_max = system->N_max + 1;
+ if ( system->num_molec_charge_constraints == 0 )
+ {
+ system->N_cm = system->N + 1;
+ system->N_cm_max = system->N_max + 1;
+ }
+ else
+ {
+ system->N_cm = system->N + system->num_molec_charge_constraints;
+ system->N_cm_max = system->N_max + system->num_molec_charge_constraints;
+ }
break;
case ACKS2_CM:
system->N_cm = 2 * system->N + 2;
@@ -440,7 +448,17 @@ static void Init_Workspace( reax_system *system, control_params *control,
workspace->b_s[i] = -system->reax_param.sbp[ system->atoms[i].type ].chi;
}
- workspace->b_s[system->N] = control->cm_q_net;
+ if ( system->num_molec_charge_constraints == 0 )
+ {
+ workspace->b_s[system->N] = control->cm_q_net;
+ }
+ else
+ {
+ for ( i = 0; i < system->num_molec_charge_constraints; ++i )
+ {
+ workspace->b_s[system->N + i] = system->molec_charge_constraints[i];
+ }
+ }
break;
case ACKS2_CM:
@@ -796,7 +814,7 @@ static void Init_Lists( reax_system *system, control_params *control,
lists[FAR_NBRS]->n = system->N;
}
- Generate_Neighbor_Lists( system, control, data, workspace, lists, out_control );
+ Generate_Neighbor_Lists( system, control, data, workspace, lists );
Htop = 0;
hb_top = scalloc( system->N, sizeof(int), "Init_Lists::hb_top" );
@@ -1339,6 +1357,15 @@ static void Finalize_System( reax_system *system, control_params *control,
system->prealloc_allocated = FALSE;
system->ffield_params_allocated = FALSE;
+ if ( system->max_num_molec_charge_constraints > 0 )
+ {
+ sfree( system->molec_charge_constraints, "Read_BGF::molec_charge_constraints" );
+ sfree( system->molec_charge_constraint_ranges, "Read_BGF::molec_charge_constraint_ranges" );
+ }
+
+ system->max_num_molec_charge_constraints = 0;
+ system->num_molec_charge_constraints = 0;
+
reax = &system->reax_param;
Finalize_Grid( system );
diff --git a/sPuReMD/src/integrate.c b/sPuReMD/src/integrate.c
index f23591a31afc1b6a7152e0812072a4df7e6d46dd..23e725ba35e7bcae684bfdafa3bc5f24495e7cb2 100644
--- a/sPuReMD/src/integrate.c
+++ b/sPuReMD/src/integrate.c
@@ -67,8 +67,7 @@ void Velocity_Verlet_NVE( reax_system *system, control_params *control,
if ( renbr == TRUE )
{
- Generate_Neighbor_Lists( system, control, data, workspace,
- lists, out_control );
+ Generate_Neighbor_Lists( system, control, data, workspace, lists );
}
Compute_Forces( system, control, data, workspace, lists, out_control, FALSE );
@@ -119,7 +118,7 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
if ( renbr == TRUE )
{
- Generate_Neighbor_Lists( system, control, data, workspace, lists, out_control );
+ Generate_Neighbor_Lists( system, control, data, workspace, lists );
}
Compute_Forces( system, control, data, workspace,
@@ -208,8 +207,7 @@ void Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system* system, control_params*
if ( renbr == TRUE )
{
- Generate_Neighbor_Lists( system, control, data, workspace,
- lists, out_control );
+ Generate_Neighbor_Lists( system, control, data, workspace, lists );
}
/* Calculate Forces at time (t + dt) */
@@ -298,8 +296,8 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
if ( renbr == TRUE )
{
Update_Grid( system );
- Generate_Neighbor_Lists( system, control, data, workspace,
- lists, out_control );
+
+ Generate_Neighbor_Lists( system, control, data, workspace, lists );
}
Compute_Forces( system, control, data, workspace, lists, out_control, FALSE );
@@ -406,8 +404,8 @@ void Velocity_Verlet_Berendsen_Semi_Isotropic_NPT( reax_system* system,
if ( renbr == TRUE )
{
Update_Grid( system );
- Generate_Neighbor_Lists( system, control, data, workspace,
- lists, out_control );
+
+ Generate_Neighbor_Lists( system, control, data, workspace, lists );
}
Compute_Forces( system, control, data, workspace, lists, out_control, FALSE );
@@ -518,8 +516,8 @@ void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system, control_params* contr
fprintf(out_control->log, "reset-");
fflush( out_control->log );
- Generate_Neighbor_Lists( system, control, data, workspace,
- lists, out_control );
+ Generate_Neighbor_Lists( system, control, data, workspace, lists );
+
fprintf(out_control->log, "nbrs-");
fflush( out_control->log );
@@ -622,8 +620,8 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system, control_params* control
fprintf(out_control->log, "reset-");
fflush( out_control->log );
- Generate_Neighbor_Lists( system, control, data, workspace,
- lists, out_control );
+ Generate_Neighbor_Lists( system, control, data, workspace, lists );
+
fprintf( out_control->log, "nbrs-" );
fflush( out_control->log );
diff --git a/sPuReMD/src/list.c b/sPuReMD/src/list.c
index 4dd71e25922e02897ae59aba097684cfe93748f2..21dd6a69e64ad0e742f939e8d31440c684b1ccd0 100644
--- a/sPuReMD/src/list.c
+++ b/sPuReMD/src/list.c
@@ -49,18 +49,6 @@ void Make_List( int n, int n_max, int total_intrs, int type, reax_list* l )
switch ( type )
{
- case TYP_VOID:
- if ( l->total_intrs > 0 )
- {
- l->v = smalloc( l->total_intrs * sizeof(void),
- "Make_List::l->v" );
- }
- else
- {
- l->v = NULL;
- }
- break;
-
case TYP_THREE_BODY:
if ( l->total_intrs > 0 )
{
@@ -174,13 +162,6 @@ void Delete_List( int type, reax_list* l )
switch ( type )
{
- case TYP_VOID:
- if ( l->v != NULL )
- {
- sfree( l->v, "Delete_List::l->v" );
- }
- break;
-
case TYP_THREE_BODY:
if ( l->three_body_list != NULL )
{
diff --git a/sPuReMD/src/lookup.c b/sPuReMD/src/lookup.c
index a63b676372dcbb2b69f72f6ad058e4562b1d2e7f..a7b0b99c139e65034ea1591796649becb4a029bf 100644
--- a/sPuReMD/src/lookup.c
+++ b/sPuReMD/src/lookup.c
@@ -184,6 +184,7 @@ static void Complete_Cubic_Spline( const real *h, const real *f, real v0, real v
}
+#if defined(DEBUG_FOCUS)
static void LR_Lookup( LR_lookup_table *t, real r, LR_data *y )
{
int i;
@@ -212,6 +213,7 @@ static void LR_Lookup( LR_lookup_table *t, real r, LR_data *y )
y->H = y->e_ele * EV_to_KCALpMOL / C_ELE;
//y->H = ((t->H[i].d*dif + t->H[i].c)*dif + t->H[i].b)*dif + t->H[i].a;
}
+#endif
void Make_LR_Lookup_Table( reax_system *system, control_params *control,
diff --git a/sPuReMD/src/multi_body.c b/sPuReMD/src/multi_body.c
index b9b1e67a52f6403393f5c40db2dec939724fa642..dae6b9039bae12ea1e35dc82ef90b924a0cbd492 100644
--- a/sPuReMD/src/multi_body.c
+++ b/sPuReMD/src/multi_body.c
@@ -40,7 +40,7 @@ void Atom_Energy( reax_system *system, control_params *control,
real exp_ovun2n, exp_ovun6, exp_ovun8;
real inv_exp_ovun1, inv_exp_ovun2, inv_exp_ovun2n, inv_exp_ovun8;
real e_un, CEunder1, CEunder2, CEunder3, CEunder4;
- real p_lp1, p_lp2, p_lp3;
+ real p_lp2, p_lp3;
real p_ovun2, p_ovun3, p_ovun4, p_ovun5, p_ovun6, p_ovun7, p_ovun8;
single_body_parameters *sbp_i;
two_body_parameters *twbp;
@@ -49,7 +49,6 @@ void Atom_Energy( reax_system *system, control_params *control,
reax_list *bonds;
bonds = lists[BONDS];
- p_lp1 = system->reax_param.gp.l[15];
p_lp3 = system->reax_param.gp.l[5];
p_ovun3 = system->reax_param.gp.l[32];
p_ovun4 = system->reax_param.gp.l[31];
diff --git a/sPuReMD/src/neighbors.c b/sPuReMD/src/neighbors.c
index d44330b2f5c2baa300ff0b0a6da959e268d6e089..8d6f59966fa4b4052434fb434d418f647705ceca 100644
--- a/sPuReMD/src/neighbors.c
+++ b/sPuReMD/src/neighbors.c
@@ -48,7 +48,8 @@ typedef int (*find_far_neighbors_function)( rvec, rvec, int, int,
simulation_box*, real, far_neighbor_data* );
-static void Choose_Neighbor_Counter( reax_system *system, control_params *control,
+static void Choose_Neighbor_Counter( reax_system const * const system,
+ control_params const * const control,
count_far_neighbors_function *Count_Far_Neighbors )
{
if ( control->periodic_boundaries == TRUE )
@@ -75,7 +76,8 @@ static void Choose_Neighbor_Counter( reax_system *system, control_params *contro
}
-static void Choose_Neighbor_Finder( reax_system *system, control_params *control,
+static void Choose_Neighbor_Finder( reax_system const * const system,
+ control_params const * const control,
find_far_neighbors_function *Find_Far_Neighbors )
{
if ( control->periodic_boundaries == TRUE )
@@ -103,7 +105,7 @@ static void Choose_Neighbor_Finder( reax_system *system, control_params *control
#if defined(DEBUG_FOCUS)
-static int compare_far_nbrs(const void *v1, const void *v2)
+static int compare_far_nbrs( const void *v1, const void *v2 )
{
return ((*(far_neighbor_data *)v1).nbr - (*(far_neighbor_data *)v2).nbr);
}
@@ -129,8 +131,10 @@ static inline real DistSqr_to_CP( rvec cp, rvec x )
}
-int Estimate_Num_Neighbors( reax_system *system, control_params *control,
- static_storage *workspace, reax_list **lists )
+/* Estimate the storage requirements for the far neighbor list */
+int Estimate_Num_Neighbors( reax_system * const system,
+ control_params const * const control, static_storage * const workspace,
+ reax_list ** const lists )
{
int i, j, k, l, m, itr;
int x, y, z;
@@ -189,6 +193,11 @@ int Estimate_Num_Neighbors( reax_system *system, control_params *control,
{
atom2 = nbr_atoms[m];
+#if defined(QMMM)
+ if ( system->atoms[atom1].qmmm_mask == TRUE
+ || system->atoms[atom2].qmmm_mask == TRUE )
+ {
+#endif
if ( atom1 >= atom2 )
{
count = Count_Far_Neighbors( system->atoms[atom1].x,
@@ -197,6 +206,9 @@ int Estimate_Num_Neighbors( reax_system *system, control_params *control,
num_far += count;
}
+#if defined(QMMM)
+ }
+#endif
}
}
@@ -216,9 +228,10 @@ int Estimate_Num_Neighbors( reax_system *system, control_params *control,
}
-void Generate_Neighbor_Lists( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace,
- reax_list **lists, output_controls *out_control )
+/* Generate the far neighbor list */
+void Generate_Neighbor_Lists( reax_system * const system,
+ control_params const * const control, simulation_data * const data,
+ static_storage * const workspace, reax_list ** const lists )
{
int i, j, k, l, m, itr;
int x, y, z;
@@ -243,9 +256,15 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
Bin_Atoms( system, workspace );
#if defined(REORDER_ATOMS)
- //Cluster_Atoms( system, workspace, control );
+ Reorder_Atoms( system, workspace, control );
#endif
+ for ( i = 0; i < far_nbrs->n; ++i )
+ {
+ Set_Start_Index( i, 0, far_nbrs );
+ Set_End_Index( i, 0, far_nbrs );
+ }
+
/* for each cell in the grid along the 3
* Cartesian directions: (i, j, k) => (x, y, z) */
for ( i = 0; i < g->ncell[0]; i++ )
@@ -287,6 +306,11 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
{
atom2 = nbr_atoms[m];
+#if defined(QMMM)
+ if ( system->atoms[atom1].qmmm_mask == TRUE
+ || system->atoms[atom2].qmmm_mask == TRUE )
+ {
+#endif
if ( atom1 >= atom2 )
{
nbr_data = &far_nbrs->far_nbr_list[num_far];
@@ -297,6 +321,9 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
num_far += count;
}
+#if defined(QMMM)
+ }
+#endif
}
}
diff --git a/sPuReMD/src/neighbors.h b/sPuReMD/src/neighbors.h
index 06ee65dc93ab64a8cfbcfda7f66c519e572f9cce..45292a9dd34f47cf0250dd64af97d07704a3a1fa 100644
--- a/sPuReMD/src/neighbors.h
+++ b/sPuReMD/src/neighbors.h
@@ -25,11 +25,11 @@
#include "reax_types.h"
-void Generate_Neighbor_Lists( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list**, output_controls* );
+void Generate_Neighbor_Lists( reax_system * const, control_params const * const,
+ simulation_data * const, static_storage * const, reax_list ** const );
-int Estimate_Num_Neighbors( reax_system*, control_params*,
- static_storage*, reax_list** );
+int Estimate_Num_Neighbors( reax_system * const, control_params const * const,
+ static_storage * const, reax_list ** const );
#endif
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index 74d03c2519cdd09d631af7ae4a1b9910fb068df1..4f006123ac87d560ac717a355fe45917764c12ea 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -48,8 +48,8 @@
//#define USE_REF_FORTRAN_EREAXFF_CONSTANTS
/* constants defined in LAMMPS ReaxFF code (useful for comparisons) */
//#define USE_LAMMPS_REAXFF_CONSTANTS
-/* enables reordering atoms after neighbor list generation (optimization) */
-#define REORDER_ATOMS
+/* enables reordering atoms after neighbor list generation for improved cache performance */
+//#define REORDER_ATOMS
/* enables support for small simulation boxes (i.e. a simulation box with any
* dimension less than twice the Verlet list cutoff distance, vlist_cut),
* which means multiple periodic images must be searched for interactions */
@@ -273,15 +273,14 @@ enum interaction_list_offets
/* interaction type */
enum interaction_type
{
- TYP_VOID = 0,
- TYP_THREE_BODY = 1,
- TYP_BOND = 2,
- TYP_HBOND = 3,
- TYP_DBO = 4,
- TYP_DDELTA = 5,
- TYP_FAR_NEIGHBOR = 6,
- TYP_NEAR_NEIGHBOR = 7,
- TYP_N = 8,
+ TYP_THREE_BODY = 0,
+ TYP_BOND = 1,
+ TYP_HBOND = 2,
+ TYP_DBO = 3,
+ TYP_DDELTA = 4,
+ TYP_FAR_NEIGHBOR = 5,
+ TYP_NEAR_NEIGHBOR = 6,
+ TYP_N = 7,
};
/* error codes for simulation termination */
@@ -878,6 +877,16 @@ struct reax_system
int N_cm;
/* max. dimension of the N x N sparse charge method matrix H across all simulations */
int N_cm_max;
+ /* molecular charge constraints
+ * NOTE: these constraints are currently only supported in BGF files using EEM */
+ real *molec_charge_constraints;
+ /* molecular charge constraints encoded as pairs of 1-based atom numbers indicating a range of atoms
+ * NOTE: these constraints are currently only supported in BGF files using EEM */
+ int *molec_charge_constraint_ranges;
+ /* num. of charge constraints on groups of atoms (i.e., molecules) */
+ unsigned int num_molec_charge_constraints;
+ /* max. num. of charge constraints on groups of atoms (i.e., molecules) */
+ unsigned int max_num_molec_charge_constraints;
/* atom info */
reax_atom *atoms;
/* atomic interaction parameters */
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index efbb5ad1c5c752aa06322d4d49194ce220c54462..7658e1a93f860fab1279c8b8ea89e758d3b78531 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -232,9 +232,8 @@ void * setup( const char * const geo_file, const char * const ffield_file,
Allocate_Top_Level_Structs( &spmd_handle );
Initialize_Top_Level_Structs( spmd_handle );
- /* note: assign here to avoid compiler warning
- * of uninitialized usage in PreAllocate_Space */
spmd_handle->system->N_max = 0;
+ spmd_handle->system->max_num_molec_charge_constraints = 0;
Read_Input_Files( geo_file, ffield_file, control_file,
spmd_handle->system, spmd_handle->control,
@@ -376,9 +375,6 @@ int simulate( const void * const handle )
Reset( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists );
- Generate_Neighbor_Lists( spmd_handle->system, spmd_handle->control, spmd_handle->data,
- spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control );
-
Compute_Forces( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, spmd_handle->lists, spmd_handle->out_control,
spmd_handle->realloc );