From 9ece9d499630385f360e5c4225b44eade51aaa44 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Thu, 16 Dec 2021 15:24:54 -0500
Subject: [PATCH] sPuReMD: add support for custom charge constraints with EEM
for QM/MM simulations.
---
sPuReMD/src/forces.c | 51 ++++++++++-
sPuReMD/src/init_md.c | 26 ++++--
sPuReMD/src/reax_types.h | 18 ++++
sPuReMD/src/spuremd.c | 181 ++++++++++++++++++++++++++++++++++-----
sPuReMD/src/spuremd.h | 6 +-
5 files changed, 252 insertions(+), 30 deletions(-)
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index 7a1258a9..d8c5b6aa 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -361,7 +361,8 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
break;
case EE_CM:
- if ( system->num_molec_charge_constraints == 0 )
+ if ( system->num_molec_charge_constraints == 0
+ && system->num_custom_charge_constraints == 0 )
{
H->start[system->N] = *Htop;
H_sp->start[system->N] = *H_sp_top;
@@ -439,6 +440,47 @@ static void Init_Charge_Matrix_Remaining_Entries( reax_system const * const syst
}
++(*H_sp_top);
}
+
+ for ( i = system->num_molec_charge_constraints;
+ i < system->num_molec_charge_constraints + system->num_custom_charge_constraints; ++i )
+ {
+ H->start[system->N + i] = *Htop;
+ H_sp->start[system->N + i] = *H_sp_top;
+
+ for ( j = system->custom_charge_constraint_start[i - system->num_molec_charge_constraints];
+ j <= system->custom_charge_constraint_start[i - system->num_molec_charge_constraints + 1]; ++j )
+ {
+ /* custom charge constraint on atoms */
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = system->custom_charge_constraint_atom_index[j] - 1;
+ H->val[*Htop] = system->custom_charge_constraint_coeff[j];
+ }
+ ++(*Htop);
+
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = system->custom_charge_constraint_atom_index[j] - 1;
+ H_sp->val[*H_sp_top] = system->custom_charge_constraint_coeff[j];
+ }
+ ++(*H_sp_top);
+ }
+
+ /* explicit zeros on diagonals */
+ if ( *Htop < H->m )
+ {
+ H->j[*Htop] = system->N + i;
+ H->val[*Htop] = 0.0;
+ }
+ ++(*Htop);
+
+ if ( *H_sp_top < H_sp->m )
+ {
+ H_sp->j[*H_sp_top] = system->N + i;
+ H_sp->val[*H_sp_top] = 0.0;
+ }
+ ++(*H_sp_top);
+ }
}
break;
@@ -1367,7 +1409,8 @@ static void Estimate_Storages_CM( reax_system const * const system,
break;
case EE_CM:
- if ( system->num_molec_charge_constraints == 0 )
+ if ( system->num_molec_charge_constraints == 0
+ && system->num_custom_charge_constraints == 0 )
{
Htop += system->N_cm;
}
@@ -1378,6 +1421,10 @@ static void Estimate_Storages_CM( reax_system const * const system,
Htop += system->molec_charge_constraint_ranges[2 * i + 1]
- system->molec_charge_constraint_ranges[2 * i] + 1;
}
+ for ( i = 0; i < system->num_custom_charge_constraints; ++i )
+ {
+ Htop += system->custom_charge_constraint_count[i];
+ }
}
break;
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 8563ff82..29db6775 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -397,11 +397,13 @@ static void Init_Workspace( reax_system * const system,
break;
case EE_CM:
system->N_cm = system->N
- + (system->num_molec_charge_constraints == 0 ? 1 : system->num_molec_charge_constraints);
+ + (system->num_molec_charge_constraints == 0 ? 1 : system->num_molec_charge_constraints)
+ + (system->num_custom_charge_constraints == 0 ? 1 : system->num_custom_charge_constraints);
if ( realloc == TRUE || system->N_cm > system->N_cm_max )
{
system->N_cm_max = system->N_max
- + (system->num_molec_charge_constraints == 0 ? 1 : system->num_molec_charge_constraints);
+ + (system->num_molec_charge_constraints == 0 ? 1 : system->num_molec_charge_constraints)
+ + (system->num_custom_charge_constraints == 0 ? 1 : system->num_custom_charge_constraints);
}
break;
case ACKS2_CM:
@@ -467,15 +469,27 @@ static void Init_Workspace( reax_system * const system,
workspace->b_s[i] = -system->reax_param.sbp[ system->atoms[i].type ].chi;
}
- if ( system->num_molec_charge_constraints == 0 )
+ if ( system->num_molec_charge_constraints == 0
+ && system->num_custom_charge_constraints == 0 )
{
workspace->b_s[system->N] = control->cm_q_net;
}
else
{
- for ( i = 0; i < system->num_molec_charge_constraints; ++i )
+ if ( system->num_molec_charge_constraints > 0 )
{
- workspace->b_s[system->N + i] = system->molec_charge_constraints[i];
+ for ( i = 0; i < system->num_molec_charge_constraints; ++i )
+ {
+ workspace->b_s[system->N + i] = system->molec_charge_constraints[i];
+ }
+ }
+ if ( system->num_custom_charge_constraints > 0 )
+ {
+ for ( i = 0; i < system->num_custom_charge_constraints; ++i )
+ {
+ workspace->b_s[system->N + system->num_molec_charge_constraints + i]
+ = system->custom_charge_constraint_rhs[i];
+ }
}
}
break;
@@ -1329,6 +1343,8 @@ static void Finalize_System( reax_system *system, control_params *control,
system->max_num_molec_charge_constraints = 0;
system->num_molec_charge_constraints = 0;
+ system->max_num_custom_charge_constraints = 0;
+ system->num_custom_charge_constraints = 0;
reax = &system->reax_param;
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index 49975296..61da0a81 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -898,6 +898,24 @@ struct reax_system
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;
+ /* num. of custom charge constraints */
+ unsigned int num_custom_charge_constraints;
+ /* max. num. of custom charge constraints */
+ unsigned int max_num_custom_charge_constraints;
+ /* num. of custom charge constraint entries (atom indices, coefficients) */
+ unsigned int num_custom_charge_constraint_entries;
+ /* max. num. of custom charge constraint entries (atom indices, coefficients) */
+ unsigned int max_num_custom_charge_constraint_entries;
+ /* entry counts for each custom charge constraint */
+ int *custom_charge_constraint_count;
+ /* indices for custom charge constraint info (atom indices, coefficients) */
+ int *custom_charge_constraint_start;
+ /* 1-based atom numbers used to identify entries for each custom charge constraint */
+ int *custom_charge_constraint_atom_index;
+ /* coefficients of entries for each custom charge constraint */
+ real *custom_charge_constraint_coeff;
+ /* right-hand side (RHS) constant values for each custom charge constraint */
+ real *custom_charge_constraint_rhs;
/* atomic interaction parameters */
reax_interaction reax_param;
/* grid specifying domain (i.e., spatial) decompisition
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index 81feaa42..c778d81d 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -227,6 +227,9 @@ void * setup( const char * const geo_file, const char * const ffield_file,
spmd_handle->system->N_max = 0;
spmd_handle->system->max_num_molec_charge_constraints = 0;
+ spmd_handle->system->num_molec_charge_constraints = 0;
+ spmd_handle->system->max_num_custom_charge_constraints = 0;
+ spmd_handle->system->num_custom_charge_constraints = 0;
Read_Input_Files( geo_file, ffield_file, control_file,
spmd_handle->system, spmd_handle->control,
@@ -1021,18 +1024,27 @@ int set_control_parameter( const void * const handle, const char * const keyword
* sim_box_info: simulation box information, where the entries are
* - box length per dimension (3 entries)
* - angles per dimension (3 entries)
- * num_charge_constraints: num. of charge constraints for charge model
- * charge_constraint_start: starting atom num. (1-based) of atom group for a charge constraint
- * charge_constraint_end: ending atom num. (1-based) of atom group for a charge constraint
- * charge_constraint_value: charge constraint value for atom group
+ * num_charge_constraint_contig: num. of contiguous charge constraints for charge model
+ * charge_constraint_contig_start: starting atom num. (1-based) of atom group for a charge constraint
+ * charge_constraint_contig_end: ending atom num. (1-based) of atom group for a charge constraint
+ * charge_constraint_contig_value: charge constraint value for atom group
+ * num_charge_constraint_custom: num. of custom charge constraints for charge model
+ * charge_constraint_custom_count: counts for each custom charge constraint
+ * charge_constraint_custom_atom_index: atom indices (1-based) for custom charge constraints
+ * charge_constraint_custom_coeff: coefficients for custom charge constraints
+ * charge_constraint_custom_rhs: right-hand side (RHS) constants for custom charge constraints
* ffield_file: file containing force field parameters
* control_file: file containing simulation parameters
*/
void * setup_qmmm( int qm_num_atoms, const char * const qm_symbols,
const double * const qm_pos, int mm_num_atoms, const char * const mm_symbols,
const double * const mm_pos_q, const double * const sim_box_info,
- int num_charge_constraints, const int * const charge_constraint_start,
- const int * const charge_constraint_end, const double * const charge_constraint_value,
+ int num_charge_constraint_contig, const int * const charge_constraint_contig_start,
+ const int * const charge_constraint_contig_end, const double * const charge_constraint_contig_value,
+ int num_charge_constraint_custom, const int * const charge_constraint_custom_count,
+ const int * const charge_constraint_custom_atom_index,
+ const double * const charge_constraint_custom_coeff,
+ const double * const charge_constraint_custom_rhs,
const char * const ffield_file, const char * const control_file )
{
int i;
@@ -1051,8 +1063,8 @@ void * setup_qmmm( int qm_num_atoms, const char * const qm_symbols,
spmd_handle->data, spmd_handle->workspace,
spmd_handle->out_control, FALSE );
- spmd_handle->system->max_num_molec_charge_constraints = num_charge_constraints;
- spmd_handle->system->num_molec_charge_constraints = num_charge_constraints;
+ spmd_handle->system->max_num_molec_charge_constraints = num_charge_constraint_contig;
+ spmd_handle->system->num_molec_charge_constraints = num_charge_constraint_contig;
if ( spmd_handle->system->num_molec_charge_constraints > 0 )
{
@@ -1065,9 +1077,54 @@ void * setup_qmmm( int qm_num_atoms, const char * const qm_symbols,
for ( i = 0; i < spmd_handle->system->num_molec_charge_constraints; ++i )
{
- spmd_handle->system->molec_charge_constraint_ranges[2 * i] = charge_constraint_start[i];
- spmd_handle->system->molec_charge_constraint_ranges[2 * i + 1] = charge_constraint_end[i];
- spmd_handle->system->molec_charge_constraints[i] = charge_constraint_value[i];
+ spmd_handle->system->molec_charge_constraint_ranges[2 * i] = charge_constraint_contig_start[i];
+ spmd_handle->system->molec_charge_constraint_ranges[2 * i + 1] = charge_constraint_contig_end[i];
+ spmd_handle->system->molec_charge_constraints[i] = charge_constraint_contig_value[i];
+ }
+ }
+
+ spmd_handle->system->max_num_custom_charge_constraints = num_charge_constraint_custom;
+ spmd_handle->system->num_custom_charge_constraints = num_charge_constraint_custom;
+ spmd_handle->system->num_custom_charge_constraint_entries = 0;
+
+ if ( spmd_handle->system->num_custom_charge_constraints > 0 )
+ {
+ spmd_handle->system->custom_charge_constraint_count = smalloc(
+ sizeof(int) * spmd_handle->system->num_custom_charge_constraints,
+ __FILE__, __LINE__ );
+ spmd_handle->system->custom_charge_constraint_start = smalloc(
+ sizeof(int) * (spmd_handle->system->num_custom_charge_constraints + 1),
+ __FILE__, __LINE__ );
+ spmd_handle->system->custom_charge_constraint_rhs = smalloc(
+ sizeof(real) * spmd_handle->system->num_custom_charge_constraints,
+ __FILE__, __LINE__ );
+
+ for ( i = 0; i < spmd_handle->system->num_custom_charge_constraints; ++i )
+ {
+ spmd_handle->system->custom_charge_constraint_count[i] = charge_constraint_custom_count[i];
+ spmd_handle->system->custom_charge_constraint_start[i] = (i == 0 ? 0 :
+ charge_constraint_custom_count[i - 1] + charge_constraint_custom_count[i - 1]);
+ spmd_handle->system->num_custom_charge_constraint_entries += charge_constraint_custom_count[i];
+ spmd_handle->system->custom_charge_constraint_rhs[i] = charge_constraint_custom_rhs[i];
+ }
+ }
+
+ spmd_handle->system->max_num_custom_charge_constraint_entries
+ = spmd_handle->system->num_custom_charge_constraint_entries;
+
+ if ( spmd_handle->system->num_custom_charge_constraint_entries > 0 )
+ {
+ spmd_handle->system->custom_charge_constraint_atom_index = smalloc(
+ sizeof(int) * spmd_handle->system->num_custom_charge_constraint_entries,
+ __FILE__, __LINE__ );
+ spmd_handle->system->custom_charge_constraint_coeff = smalloc(
+ sizeof(real) * spmd_handle->system->num_custom_charge_constraint_entries,
+ __FILE__, __LINE__ );
+
+ for ( i = 0; i < spmd_handle->system->num_custom_charge_constraint_entries; ++i )
+ {
+ spmd_handle->system->custom_charge_constraint_atom_index[i] = charge_constraint_custom_atom_index[i];
+ spmd_handle->system->custom_charge_constraint_coeff[i] = charge_constraint_custom_coeff[i];
}
}
@@ -1176,10 +1233,15 @@ void * setup_qmmm( int qm_num_atoms, const char * const qm_symbols,
* sim_box_info: simulation box information, where the entries are
* - box length per dimension (3 entries)
* - angles per dimension (3 entries)
- * num_charge_constraints: num. of charge constraints for charge model
- * charge_constraint_start: starting atom num. (1-based) of atom group for a charge constraint
- * charge_constraint_end: ending atom num. (1-based) of atom group for a charge constraint
- * charge_constraint_value: charge constraint value for atom group
+ * num_charge_constraint_contig: num. of contiguous charge constraints for charge model
+ * charge_constraint_contig_start: starting atom num. (1-based) of atom group for a charge constraint
+ * charge_constraint_contig_end: ending atom num. (1-based) of atom group for a charge constraint
+ * charge_constraint_contig_value: charge constraint value for atom group
+ * num_charge_constraint_custom: num. of custom charge constraints for charge model
+ * charge_constraint_custom_count: counts for each custom charge constraint
+ * charge_constraint_custom_atom_index: atom indices (1-based) for custom charge constraints
+ * charge_constraint_custom_coeff: coefficients for custom charge constraints
+ * charge_constraint_custom_rhs: right-hand side (RHS) constants for custom charge constraints
* ffield_file: file containing force field parameters
* control_file: file containing simulation parameters
*
@@ -1189,8 +1251,12 @@ int reset_qmmm( const void * const handle, int qm_num_atoms,
const char * const qm_symbols, const double * const qm_pos,
int mm_num_atoms, const char * const mm_symbols,
const double * const mm_pos_q, const double * const sim_box_info,
- int num_charge_constraints, const int * const charge_constraint_start,
- const int * const charge_constraint_end, const double * const charge_constraint_value,
+ int num_charge_constraint_contig, const int * const charge_constraint_contig_start,
+ const int * const charge_constraint_contig_end, const double * const charge_constraint_contig_value,
+ int num_charge_constraint_custom, const int * const charge_constraint_custom_count,
+ const int * const charge_constraint_custom_atom_index,
+ const double * const charge_constraint_custom_coeff,
+ const double * const charge_constraint_custom_rhs,
const char * const ffield_file, const char * const control_file )
{
int i, ret;
@@ -1219,7 +1285,7 @@ int reset_qmmm( const void * const handle, int qm_num_atoms,
spmd_handle->data, spmd_handle->workspace,
spmd_handle->out_control, TRUE );
- spmd_handle->system->num_molec_charge_constraints = num_charge_constraints;
+ spmd_handle->system->num_molec_charge_constraints = num_charge_constraint_contig;
if ( spmd_handle->system->num_molec_charge_constraints
> spmd_handle->system->max_num_molec_charge_constraints )
@@ -1247,9 +1313,82 @@ int reset_qmmm( const void * const handle, int qm_num_atoms,
{
for ( i = 0; i < spmd_handle->system->num_molec_charge_constraints; ++i )
{
- spmd_handle->system->molec_charge_constraint_ranges[2 * i] = charge_constraint_start[i];
- spmd_handle->system->molec_charge_constraint_ranges[2 * i + 1] = charge_constraint_end[i];
- spmd_handle->system->molec_charge_constraints[i] = charge_constraint_value[i];
+ spmd_handle->system->molec_charge_constraint_ranges[2 * i] = charge_constraint_contig_start[i];
+ spmd_handle->system->molec_charge_constraint_ranges[2 * i + 1] = charge_constraint_contig_end[i];
+ spmd_handle->system->molec_charge_constraints[i] = charge_constraint_contig_value[i];
+ }
+ }
+
+ spmd_handle->system->num_custom_charge_constraints = num_charge_constraint_custom;
+
+ if ( spmd_handle->system->num_custom_charge_constraints
+ > spmd_handle->system->max_num_custom_charge_constraints )
+ {
+ if ( spmd_handle->system->max_num_custom_charge_constraints > 0 )
+ {
+ sfree( spmd_handle->system->custom_charge_constraint_count,
+ __FILE__, __LINE__ );
+ sfree( spmd_handle->system->custom_charge_constraint_start,
+ __FILE__, __LINE__ );
+ sfree( spmd_handle->system->custom_charge_constraint_rhs,
+ __FILE__, __LINE__ );
+ }
+
+ spmd_handle->system->custom_charge_constraint_count = smalloc(
+ sizeof(int) * spmd_handle->system->num_custom_charge_constraints,
+ __FILE__, __LINE__ );
+ spmd_handle->system->custom_charge_constraint_start = smalloc(
+ sizeof(int) * (spmd_handle->system->num_custom_charge_constraints + 1),
+ __FILE__, __LINE__ );
+ spmd_handle->system->custom_charge_constraint_rhs = smalloc(
+ sizeof(real) * spmd_handle->system->num_custom_charge_constraints,
+ __FILE__, __LINE__ );
+
+ spmd_handle->system->max_num_custom_charge_constraints
+ = spmd_handle->system->num_custom_charge_constraints;
+ }
+
+ spmd_handle->system->num_custom_charge_constraint_entries = 0;
+ if ( spmd_handle->system->num_custom_charge_constraints > 0 )
+ {
+ for ( i = 0; i < spmd_handle->system->num_custom_charge_constraints; ++i )
+ {
+ spmd_handle->system->custom_charge_constraint_count[i] = charge_constraint_custom_count[i];
+ spmd_handle->system->custom_charge_constraint_start[i] = (i == 0 ? 0 :
+ charge_constraint_custom_count[i - 1] + charge_constraint_custom_count[i - 1]);
+ spmd_handle->system->num_custom_charge_constraint_entries += charge_constraint_custom_count[i];
+ spmd_handle->system->custom_charge_constraint_rhs[i] = charge_constraint_custom_rhs[i];
+ }
+ }
+
+ if ( spmd_handle->system->num_custom_charge_constraint_entries
+ > spmd_handle->system->max_num_custom_charge_constraint_entries )
+ {
+ if ( spmd_handle->system->max_num_custom_charge_constraint_entries > 0 )
+ {
+ sfree( spmd_handle->system->custom_charge_constraint_atom_index,
+ __FILE__, __LINE__ );
+ sfree( spmd_handle->system->custom_charge_constraint_coeff,
+ __FILE__, __LINE__ );
+ }
+
+ spmd_handle->system->custom_charge_constraint_atom_index = smalloc(
+ sizeof(int) * spmd_handle->system->num_custom_charge_constraint_entries,
+ __FILE__, __LINE__ );
+ spmd_handle->system->custom_charge_constraint_coeff = smalloc(
+ sizeof(real) * spmd_handle->system->num_custom_charge_constraint_entries,
+ __FILE__, __LINE__ );
+
+ spmd_handle->system->max_num_custom_charge_constraint_entries
+ = spmd_handle->system->num_custom_charge_constraint_entries;
+ }
+
+ if ( spmd_handle->system->num_custom_charge_constraint_entries > 0 )
+ {
+ for ( i = 0; i < spmd_handle->system->num_custom_charge_constraint_entries; ++i )
+ {
+ spmd_handle->system->custom_charge_constraint_atom_index[i] = charge_constraint_custom_atom_index[i];
+ spmd_handle->system->custom_charge_constraint_coeff[i] = charge_constraint_custom_coeff[i];
}
}
diff --git a/sPuReMD/src/spuremd.h b/sPuReMD/src/spuremd.h
index 5b4d44f0..a644bb21 100644
--- a/sPuReMD/src/spuremd.h
+++ b/sPuReMD/src/spuremd.h
@@ -77,13 +77,15 @@ void * setup_qmmm( int, const char * const,
const double * const, int, const char * const,
const double * const, const double * const,
int, const int * const, const int * const, const double * const,
- const char * const, const char * const );
+ int, const int * const, const int * const, const double * const,
+ const double * const, const char * const, const char * const );
int reset_qmmm( const void * const, int, const char * const,
const double * const, int, const char * const,
const double * const, const double * const,
int, const int * const, const int * const, const double * const,
- const char * const, const char * const);
+ int, const int * const, const int * const, const double * const,
+ const double * const, const char * const, const char * const);
int get_atom_positions_qmmm( const void * const, double * const,
double * const );
--
GitLab