From 9e08fb6bcf8ec5d6fa53bd3a435f568e7a27166a Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnk@msu.edu>
Date: Sun, 15 Apr 2018 09:12:20 -0400
Subject: [PATCH] sPuReMD: remove debug prints. Default to no good initial
guess for out of range values for spline extrapolation.
---
sPuReMD/src/charges.c | 16 ++++++++++++++--
sPuReMD/src/ffield.c | 3 +++
sPuReMD/src/init_md.c | 13 ++++++-------
sPuReMD/src/integrate.c | 17 +++++++++++------
sPuReMD/src/spuremd.c | 3 +++
sPuReMD/src/system_props.c | 6 +++---
6 files changed, 40 insertions(+), 18 deletions(-)
diff --git a/sPuReMD/src/charges.c b/sPuReMD/src/charges.c
index b61db4b8..17c013e6 100644
--- a/sPuReMD/src/charges.c
+++ b/sPuReMD/src/charges.c
@@ -47,7 +47,7 @@ static void Extrapolate_Charges_QEq( const reax_system * const system,
#endif
for ( i = 0; i < system->N_cm; ++i )
{
- /* no extrapolation */
+ /* no extrapolation, previous solution as initial guess */
if ( control->cm_init_guess_extrap1 == 0 )
{
s_tmp = workspace->s[0][i];
@@ -74,8 +74,12 @@ static void Extrapolate_Charges_QEq( const reax_system * const system,
s_tmp = 5 * (workspace->s[0][i] - workspace->s[3][i]) +
10 * (-workspace->s[1][i] + workspace->s[2][i] ) + workspace->s[4][i];
}
+ else
+ {
+ s_tmp = 0.0;
+ }
- /* no extrapolation */
+ /* no extrapolation, previous solution as initial guess */
if ( control->cm_init_guess_extrap2 == 0 )
{
t_tmp = workspace->t[0][i];
@@ -102,6 +106,10 @@ static void Extrapolate_Charges_QEq( const reax_system * const system,
t_tmp = 5.0 * (workspace->t[0][i] - workspace->t[3][i]) +
10.0 * (-workspace->t[1][i] + workspace->t[2][i] ) + workspace->t[4][i];
}
+ else
+ {
+ t_tmp = 0.0;
+ }
workspace->s[4][i] = workspace->s[3][i];
workspace->s[3][i] = workspace->s[2][i];
@@ -160,6 +168,10 @@ static void Extrapolate_Charges_EE( const reax_system * const system,
s_tmp = 5.0 * (workspace->s[0][i] - workspace->s[3][i]) +
10.0 * (-workspace->s[1][i] + workspace->s[2][i] ) + workspace->s[4][i];
}
+ else
+ {
+ s_tmp = 0.0;
+ }
workspace->s[4][i] = workspace->s[3][i];
workspace->s[3][i] = workspace->s[2][i];
diff --git a/sPuReMD/src/ffield.c b/sPuReMD/src/ffield.c
index a7b05662..ed0251a7 100644
--- a/sPuReMD/src/ffield.c
+++ b/sPuReMD/src/ffield.c
@@ -233,6 +233,7 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
val = atof(tmp[7]);
reax->sbp[i].acore2 = val;
+ /* inner-wall */
if ( reax->sbp[i].rcore2 > 0.01 && reax->sbp[i].acore2 > 0.01 ) // Inner-wall
{
/* shielding vdWaals */
@@ -392,6 +393,7 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
/* calculating combination rules and filling up remaining fields. */
for (i = 0; i < reax->num_atom_types; i++)
+ {
for (j = i; j < reax->num_atom_types; j++)
{
reax->tbp[i][j].r_s = 0.5 * (reax->sbp[i].r_s + reax->sbp[j].r_s);
@@ -437,6 +439,7 @@ void Read_Force_Field( FILE* fp, reax_interaction* reax )
reax->tbp[i][j].rcore = SQRT( reax->sbp[i].rcore2 * reax->sbp[j].rcore2 );
reax->tbp[j][i].rcore = SQRT( reax->sbp[j].rcore2 * reax->sbp[i].rcore2 );
}
+ }
/* next line is number of 2-body offdiagonal combinations and some comments */
/* these are two body off-diagonal terms that are different from the
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 30e9ef4c..d098007c 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -169,8 +169,8 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
*Evolve = Velocity_Verlet_Nose_Hoover_NVT_Klein;
break;
-
- case NPT: // Anisotropic NPT
+ /* anisotropic NPT */
+ case NPT:
fprintf( stderr, "THIS OPTION IS NOT YET IMPLEMENTED! TERMINATING...\n" );
exit( UNKNOWN_OPTION );
data->N_f = 3 * system->N + 9;
@@ -186,14 +186,14 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
*Evolve = Velocity_Verlet_Berendsen_Isotropic_NPT;
break;
-
- case sNPT: // Semi-Isotropic NPT
+ /* semi-isotropic NPT */
+ case sNPT:
data->N_f = 3 * system->N + 4;
*Evolve = Velocity_Verlet_Berendsen_SemiIsotropic_NPT;
break;
-
- case iNPT: // Isotropic NPT
+ /* isotropic NPT */
+ case iNPT:
data->N_f = 3 * system->N + 2;
*Evolve = Velocity_Verlet_Berendsen_Isotropic_NPT;
break;
@@ -201,7 +201,6 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
case bNVT:
data->N_f = 3 * system->N + 1;
*Evolve = Velocity_Verlet_Berendsen_NVT;
- fprintf (stderr, " Initializing Velocity_Verlet_Berendsen_NVT .... \n");
break;
default:
diff --git a/sPuReMD/src/integrate.c b/sPuReMD/src/integrate.c
index fbcc1667..4e948406 100644
--- a/sPuReMD/src/integrate.c
+++ b/sPuReMD/src/integrate.c
@@ -447,7 +447,6 @@ void Velocity_Verlet_Berendsen_SemiIsotropic_NPT( reax_system* system,
/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/
/*!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*/
/************************************************/
-
#ifdef ANISOTROPIC
void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system,
control_params* control, simulation_data *data,
@@ -657,8 +656,10 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
E_kin = 0;
for ( i = 0; i < system->N; ++i )
+ {
E_kin += (0.5 * system->reaxprm.sbp[system->atoms[i].type].mass *
rvec_Dot( system->atoms[i].v, system->atoms[i].v ) );
+ }
P_int = inv_3V * ( 2.0 * E_kin + P_int_const );
@@ -684,13 +685,12 @@ void Velocity_Verlet_Isotropic_NPT( reax_system* system,
fprintf(out_control->log, "xi: \tG- %8.3f v- %8.3f xi - %8.3f\n",
therm->G_xi, therm->v_xi, therm->xi);
}
-
#endif
-/* uses Berendsen-type coupling for both T and P.
- All box dimensions are scaled by the same amount,
- there is no change in the angles between axes. */
+/* Uses Berendsen-type coupling for both T and P.
+ * All box dimensions are scaled by the same amount,
+ * there is no change in the angles between axes. */
void Velocity_Verlet_Berendsen_NVT( reax_system* system,
control_params* control, simulation_data *data,
static_storage *workspace, reax_list **lists,
@@ -705,6 +705,7 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
#if defined(DEBUG_FOCUS)
fprintf( stderr, "step%d\n", data->step );
#endif
+
dt = control->dt;
steps = data->step - data->prev_steps;
renbr = (steps % control->reneighbor == 0);
@@ -734,7 +735,9 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
Reset( system, control, data, workspace, lists );
if ( renbr )
+ {
Generate_Neighbor_Lists( system, control, data, workspace, lists, out_control );
+ }
Compute_Forces( system, control, data, workspace,
lists, out_control, interaction_functions );
@@ -764,7 +767,9 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
}
lambda = SQRT( lambda );
- fprintf (stderr, "step:%d lambda -> %f \n", data->step, lambda);
+#if defined(DEBUG_FOCUS)
+ fprintf( stderr, "step:%d lambda -> %f \n", data->step, lambda );
+#endif
/* Scale velocities and positions at t+dt */
for ( i = 0; i < system->N; ++i )
diff --git a/sPuReMD/src/spuremd.c b/sPuReMD/src/spuremd.c
index 7517edb2..968847ce 100644
--- a/sPuReMD/src/spuremd.c
+++ b/sPuReMD/src/spuremd.c
@@ -242,6 +242,9 @@ int simulate( const void * const handle )
Output_Results( spmd_handle->system, spmd_handle->control, spmd_handle->data,
spmd_handle->workspace, &spmd_handle->lists, spmd_handle->out_control );
}
+
+ Check_Energy( spmd_handle->data );
+
if ( spmd_handle->callback != NULL )
{
spmd_handle->callback( spmd_handle->system->atoms, spmd_handle->data,
diff --git a/sPuReMD/src/system_props.c b/sPuReMD/src/system_props.c
index 5b6ee58f..a7bf0a2f 100644
--- a/sPuReMD/src/system_props.c
+++ b/sPuReMD/src/system_props.c
@@ -227,9 +227,9 @@ void Compute_Total_Energy( reax_system* system, simulation_data* data )
q = system->atoms[i].q;
type_i = system->atoms[i].type;
- e_pol += ( system->reaxprm.sbp[ type_i ].chi * q +
- (system->reaxprm.sbp[ type_i ].eta / 2.0) * SQR( q ) ) *
- KCALpMOL_to_EV;
+ e_pol += ( system->reaxprm.sbp[ type_i ].chi * q
+ + (system->reaxprm.sbp[ type_i ].eta / 2.0) * SQR( q ) )
+ * KCALpMOL_to_EV;
}
data->E_Pol = e_pol;
--
GitLab