From 4e9bbf009770eb08fec0f75060fe882e2275020c Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Tue, 1 Jun 2021 12:48:06 -0400
Subject: [PATCH] sPuReMD: add support for user-defined molecular charge
constraints with EEM (using MOLCHARGE keyword in BGF geometry files). Fix
issue with far neighbor list not being initialized for certain edge cases.
---
sPuReMD/src/forces.c | 106 +++++++++++++++++++++++++++++----------
sPuReMD/src/geo_tools.c | 68 +++++++++++++++++++++++--
sPuReMD/src/init_md.c | 33 ++++++++++--
sPuReMD/src/neighbors.c | 6 +++
sPuReMD/src/reax_types.h | 10 ++++
sPuReMD/src/spuremd.c | 3 +-
6 files changed, 190 insertions(+), 36 deletions(-)
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 6e2fc638..1b9cdc68 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:
diff --git a/sPuReMD/src/geo_tools.c b/sPuReMD/src/geo_tools.c
index 3b361e20..a31b4326 100644
--- a/sPuReMD/src/geo_tools.c
+++ b/sPuReMD/src/geo_tools.c
@@ -625,6 +625,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 )
{
@@ -638,7 +645,7 @@ void Read_BGF( const char * const bgf_file, reax_system* system, control_params
char element[6], charge[9];
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 +658,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 +671,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 +684,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;
@@ -838,6 +880,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/init_md.c b/sPuReMD/src/init_md.c
index dd1b8291..e6c20c76 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:
@@ -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/neighbors.c b/sPuReMD/src/neighbors.c
index d44330b2..c2cf87a8 100644
--- a/sPuReMD/src/neighbors.c
+++ b/sPuReMD/src/neighbors.c
@@ -246,6 +246,12 @@ void Generate_Neighbor_Lists( reax_system *system, control_params *control,
//Cluster_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++ )
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index 74d03c25..5829d839 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -878,6 +878,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 efbb5ad1..7f9408b1 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,
--
GitLab