From b09550ea317c7165c4f17816e72fd08bbcf4b085 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Tue, 7 Dec 2021 14:51:18 -0500
Subject: [PATCH] sPuReMD: add support for setting charge computation frequency
(charge_freq).
---
sPuReMD/src/control.c | 5 +++
sPuReMD/src/forces.c | 85 +++++++++++++++++++++++++------------------
2 files changed, 54 insertions(+), 36 deletions(-)
diff --git a/sPuReMD/src/control.c b/sPuReMD/src/control.c
index da6c1ff..03f9a79 100644
--- a/sPuReMD/src/control.c
+++ b/sPuReMD/src/control.c
@@ -137,6 +137,10 @@ int Set_Control_Parameter( const char * const keyword,
{
control->charge_method = sstrtol( values[0], __FILE__, __LINE__ );
}
+ else if ( strncmp(keyword, "charge_freq", MAX_LINE) == 0 )
+ {
+ control->charge_freq = sstrtol( values[0], __FILE__, __LINE__ );
+ }
else if ( strncmp(keyword, "cm_q_net", MAX_LINE) == 0 )
{
val = sstrtod( values[0], __FILE__, __LINE__ );
@@ -512,6 +516,7 @@ void Set_Control_Defaults( reax_system * const system,
control->restrict_type = 0;
control->charge_method = QEQ_CM;
+ control->charge_freq = 1;
control->cm_q_net = 0.0;
control->cm_solver_type = GMRES_S;
control->cm_solver_max_iters = 100;
diff --git a/sPuReMD/src/forces.c b/sPuReMD/src/forces.c
index cb4d77c..7a1258a 100644
--- a/sPuReMD/src/forces.c
+++ b/sPuReMD/src/forces.c
@@ -140,8 +140,6 @@ static void Compute_NonBonded_Forces( reax_system *system, control_params *contr
simulation_data *data, static_storage *workspace,
reax_list** lists, output_controls *out_control, int realloc )
{
- real t_start, t_elapsed;
-
#if defined(TEST_ENERGY)
fprintf( out_control->evdw, "step: %d\n%6s%6s%12s%12s%12s\n",
data->step, "atom1", "atom2", "r12", "evdw", "total" );
@@ -149,38 +147,6 @@ static void Compute_NonBonded_Forces( reax_system *system, control_params *contr
data->step, "atom1", "atom2", "r12", "q1", "q2", "ecou", "total" );
#endif
- t_start = Get_Time( );
- Compute_Charges( system, control, data, workspace, out_control, realloc );
- t_elapsed = Get_Timing_Info( t_start );
- data->timing.cm += t_elapsed;
-
- if ( control->cm_solver_pre_comp_refactor == -1 )
- {
- if ( data->step <= 4 || is_refactoring_step( control, data ) )
- {
- if ( is_refactoring_step( control, data ) )
- {
- data->timing.cm_last_pre_comp = data->timing.cm_solver_pre_comp;
- }
-
- data->timing.cm_optimum = data->timing.cm_solver_pre_app
- + data->timing.cm_solver_spmv
- + data->timing.cm_solver_vector_ops
- + data->timing.cm_solver_orthog
- + data->timing.cm_solver_tri_solve;
- data->timing.cm_total_loss = 0.0;
- }
- else
- {
- data->timing.cm_total_loss += data->timing.cm_solver_pre_app
- + data->timing.cm_solver_spmv
- + data->timing.cm_solver_vector_ops
- + data->timing.cm_solver_orthog
- + data->timing.cm_solver_tri_solve
- - data->timing.cm_optimum;
- }
- }
-
if ( control->tabulate <= 0 )
{
vdW_Coulomb_Energy( system, control, data, workspace, lists, out_control );
@@ -1274,7 +1240,9 @@ static int Init_Forces( reax_system * const system,
dist_done = TRUE;
}
- if ( cm_done == FALSE )
+ if ( (control->charge_freq > 0
+ && (data->step - data->prev_steps) % control->charge_freq == 0)
+ && cm_done == FALSE )
{
if ( control->tabulate <= 0 )
{
@@ -1597,9 +1565,19 @@ int Compute_Forces( reax_system * const system, control_params * const control,
simulation_data * const data, static_storage * const workspace,
reax_list ** const lists, output_controls * const out_control, int realloc )
{
- int ret;
+ int charge_flag, ret;
real t_start, t_elapsed;
+ if ( control->charge_freq > 0
+ && (data->step - data->prev_steps) % control->charge_freq == 0 )
+ {
+ charge_flag = TRUE;
+ }
+ else
+ {
+ charge_flag = FALSE;
+ }
+
t_start = Get_Time( );
ret = Init_Forces( system, control, data, workspace, lists, out_control );
t_elapsed = Get_Timing_Info( t_start );
@@ -1618,6 +1596,41 @@ int Compute_Forces( reax_system * const system, control_params * const control,
t_elapsed = Get_Timing_Info( t_start );
data->timing.bonded += t_elapsed;
+ if ( charge_flag == TRUE )
+ {
+ t_start = Get_Time( );
+ Compute_Charges( system, control, data, workspace, out_control, realloc );
+ t_elapsed = Get_Timing_Info( t_start );
+ data->timing.cm += t_elapsed;
+ }
+
+ if ( control->cm_solver_pre_comp_refactor == -1 )
+ {
+ if ( data->step <= 4 || is_refactoring_step( control, data ) )
+ {
+ if ( is_refactoring_step( control, data ) )
+ {
+ data->timing.cm_last_pre_comp = data->timing.cm_solver_pre_comp;
+ }
+
+ data->timing.cm_optimum = data->timing.cm_solver_pre_app
+ + data->timing.cm_solver_spmv
+ + data->timing.cm_solver_vector_ops
+ + data->timing.cm_solver_orthog
+ + data->timing.cm_solver_tri_solve;
+ data->timing.cm_total_loss = 0.0;
+ }
+ else
+ {
+ data->timing.cm_total_loss += data->timing.cm_solver_pre_app
+ + data->timing.cm_solver_spmv
+ + data->timing.cm_solver_vector_ops
+ + data->timing.cm_solver_orthog
+ + data->timing.cm_solver_tri_solve
+ - data->timing.cm_optimum;
+ }
+ }
+
t_start = Get_Time( );
Compute_NonBonded_Forces( system, control, data, workspace,
lists, out_control, realloc );
--
GitLab