From 1ec5e509bf3611be21dc44252410de4c7ba57632 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Tue, 8 Jun 2021 12:48:50 -0400
Subject: [PATCH] sPuReMD: fix issue with atoms exactly on the upper boundary
of the simulation box not being correctly binned into grid cells for the cell
method. Refine conversion constants for kinetic energy and temperature
calculations to obtain better agreement with the Fortran ReaxFF codebase.
---
sPuReMD/src/box.c | 32 ++++++++++++--------------------
sPuReMD/src/grid.c | 32 ++++++++++++++++++++++++--------
sPuReMD/src/init_md.c | 16 ++++++++--------
sPuReMD/src/integrate.c | 23 +++++++++++++----------
sPuReMD/src/io_tools.c | 12 +++++++-----
sPuReMD/src/io_tools.h | 5 +++--
sPuReMD/src/reax_types.h | 7 ++++---
sPuReMD/src/system_props.c | 9 +++++----
sPuReMD/src/tool_box.c | 9 ++++++---
sPuReMD/src/traj.c | 4 ++--
10 files changed, 84 insertions(+), 65 deletions(-)
diff --git a/sPuReMD/src/box.c b/sPuReMD/src/box.c
index 588aa340..7a20fa9f 100644
--- a/sPuReMD/src/box.c
+++ b/sPuReMD/src/box.c
@@ -230,7 +230,7 @@ void Update_Box_Semi_Isotropic( simulation_box *box, rvec mu )
* is within the simulation box
*
* Assumption: (0, 0, 0) is the minimum box coordinate,
- * the the maximum coordinate is strictly positive in each dimension
+ * the maximum coordinate is strictly positive in each dimension
*
* Inputs:
* x: current atom position
@@ -253,30 +253,22 @@ void Update_Atom_Position_Periodic( rvec x, rvec dx, ivec rel_map, simulation_bo
/* new atomic position after update along dimension i */
tmp = x[i] + dx[i];
- /* outside of simulation box boundary [0.0, d_i],
+ /* outside of simulation box boundary [0.0, d_i),
* where d_i is the simulation box length along dimesion i */
- if ( tmp < 0.0 || tmp > box->box_norms[i] )
+ if ( tmp < 0.0 || tmp >= box->box_norms[i] )
{
remapped = FALSE;
- /* re-map the position to be in the range [-d_i, d_i] */
- if ( tmp <= -box->box_norms[i] || tmp >= box->box_norms[i] )
+ /* re-map the position to be in the range (-d_i, d_i) */
+ if ( tmp <= -1.0 * box->box_norms[i] || tmp >= box->box_norms[i] )
{
- if ( tmp >= 0.0 )
- {
- rel_map[i] += ((int) CEIL( tmp / box->box_norms[i] )) - 1;
- }
- else
- {
- rel_map[i] += (int) FLOOR( tmp / box->box_norms[i] );
- }
-
+ rel_map[i] += (int) (tmp / box->box_norms[i]);
tmp = FMOD( tmp, box->box_norms[i] );
remapped = TRUE;
}
/* if the new position after re-mapping is negative,
- * translate back to the range [0, d_i] (i.e., inside the box) */
+ * translate back to the range [0, d_i) (i.e., inside the box) */
if ( tmp < 0.0 )
{
tmp += box->box_norms[i];
@@ -288,7 +280,7 @@ void Update_Atom_Position_Periodic( rvec x, rvec dx, ivec rel_map, simulation_bo
}
}
- assert( tmp >= 0.0 && tmp <= box->box_norms[i] );
+ assert( tmp >= 0.0 && tmp < box->box_norms[i] );
x[i] = tmp;
}
@@ -327,7 +319,7 @@ void Update_Atom_Position_Non_Periodic( rvec x, rvec dx, ivec rel_map, simulatio
}
else if ( tmp > box->box_norms[i] )
{
- tmp = box->box_norms[i];
+ tmp = box->box_norms[i] - 1.0E-10;
}
x[i] = tmp;
@@ -665,17 +657,17 @@ int Find_Periodic_Far_Neighbors_Small_Box( rvec x1, rvec x2, int atom, int nbr_a
kmax = (int) CEIL( vlist_cut / box->box_norms[2] );
/* assumption: orthogonal coordinates */
- for ( i = -imax; i <= imax; ++i )
+ for ( i = -1 * imax; i <= imax; ++i )
{
d_i = (x2[0] + i * box->box_norms[0]) - x1[0];
sqr_d_i = SQR( d_i );
- for ( j = -jmax; j <= jmax; ++j )
+ for ( j = -1 * jmax; j <= jmax; ++j )
{
d_j = (x2[1] + j * box->box_norms[1]) - x1[1];
sqr_d_j = SQR( d_j );
- for ( k = -kmax; k <= kmax; ++k )
+ for ( k = -1 * kmax; k <= kmax; ++k )
{
d_k = (x2[2] + k * box->box_norms[2]) - x1[2];
sqr_d_k = SQR( d_k );
diff --git a/sPuReMD/src/grid.c b/sPuReMD/src/grid.c
index ad4e4fb6..f5bfde66 100644
--- a/sPuReMD/src/grid.c
+++ b/sPuReMD/src/grid.c
@@ -26,6 +26,11 @@
#include "vector.h"
+/* Counts the number of atoms binned to each cell within the grid
+ *
+ * NOTE: this assumes the atom positions are within the
+ * simulation box boundaries [0, d_i) where d_i is the simulation
+ * box length along a particular dimension */
static int Estimate_GCell_Population( reax_system * const system )
{
int i, j, k, l;
@@ -38,8 +43,10 @@ static int Estimate_GCell_Population( reax_system * const system )
for ( l = 0; l < system->N; l++ )
{
- /* NOTE: this assumes the atom positions
- * are within the simulation box boundaries */
+ assert( system->atoms[l].x[0] >= 0.0 && system->atoms[l].x[0] < system->box.box_norms[0] );
+ assert( system->atoms[l].x[1] >= 0.0 && system->atoms[l].x[1] < system->box.box_norms[1] );
+ assert( system->atoms[l].x[2] >= 0.0 && system->atoms[l].x[2] < system->box.box_norms[2] );
+
i = (int) (system->atoms[l].x[0] * g->inv_len[0]);
j = (int) (system->atoms[l].x[1] * g->inv_len[1]);
k = (int) (system->atoms[l].x[2] * g->inv_len[2]);
@@ -406,15 +413,15 @@ static void Find_Neighbor_Grid_Cells( grid * const g )
stack_top = 0;
/* choose an unmarked neighbor cell */
- for ( di = -g->spread[0]; di <= g->spread[0]; di++ )
+ for ( di = -1 * g->spread[0]; di <= g->spread[0]; di++ )
{
x = Shift( i, di, 0, g );
- for ( dj = -g->spread[1]; dj <= g->spread[1]; dj++ )
+ for ( dj = -1 * g->spread[1]; dj <= g->spread[1]; dj++ )
{
y = Shift( j, dj, 1, g );
- for ( dk = -g->spread[2]; dk <= g->spread[2]; dk++ )
+ for ( dk = -1 * g->spread[2]; dk <= g->spread[2]; dk++ )
{
z = Shift( k, dk, 2, g );
@@ -614,6 +621,11 @@ void Update_Grid( reax_system * const system )
}
+/* Utilize a cell method approach to bin atoms to each cell within the grid
+ *
+ * NOTE: this assumes the atom positions are within the
+ * simulation box boundaries [0, d_i) where d_i is the simulation
+ * box length along a particular dimension */
void Bin_Atoms( reax_system * const system, static_storage * const workspace )
{
int i, j, k, l;
@@ -626,9 +638,13 @@ void Bin_Atoms( reax_system * const system, static_storage * const workspace )
for ( l = 0; l < system->N; l++ )
{
- i = (int)(system->atoms[l].x[0] * g->inv_len[0]);
- j = (int)(system->atoms[l].x[1] * g->inv_len[1]);
- k = (int)(system->atoms[l].x[2] * g->inv_len[2]);
+ assert( system->atoms[l].x[0] >= 0.0 && system->atoms[l].x[0] < system->box.box_norms[0] );
+ assert( system->atoms[l].x[1] >= 0.0 && system->atoms[l].x[1] < system->box.box_norms[1] );
+ assert( system->atoms[l].x[2] >= 0.0 && system->atoms[l].x[2] < system->box.box_norms[2] );
+
+ i = (int) (system->atoms[l].x[0] * g->inv_len[0]);
+ j = (int) (system->atoms[l].x[1] * g->inv_len[1]);
+ k = (int) (system->atoms[l].x[2] * g->inv_len[2]);
/* atom index in grid cell => atom number */
g->atoms[i][j][k][g->top[i][j][k]] = l;
/* count of current atoms in this grid cell */
diff --git a/sPuReMD/src/init_md.c b/sPuReMD/src/init_md.c
index 7c9fc2e9..c97368cb 100644
--- a/sPuReMD/src/init_md.c
+++ b/sPuReMD/src/init_md.c
@@ -193,24 +193,24 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
switch ( control->ensemble )
{
case NVE:
- data->N_f = 3 * system->N;
+ data->N_f = 3.0 * system->N;
*Evolve = &Velocity_Verlet_NVE;
break;
case bNVT:
- data->N_f = 3 * system->N;
+ data->N_f = 3.0 * system->N;
*Evolve = &Velocity_Verlet_Berendsen_NVT;
break;
case nhNVT:
- data->N_f = 3 * system->N + 1;
+ data->N_f = 3.0 * system->N + 1.0;
*Evolve = &Velocity_Verlet_Nose_Hoover_NVT_Klein;
if ( control->restart == FALSE
|| (control->restart == TRUE && control->random_vel == TRUE) )
{
data->therm.G_xi = control->Tau_T * (2.0 * data->E_Kin
- - data->N_f * K_B * control->T);
+ - data->N_f * K_B / F_CONV * control->T);
data->therm.v_xi = data->therm.G_xi * control->dt;
data->therm.v_xi_old = 0.0;
data->therm.xi = 0.0;
@@ -222,13 +222,13 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
fprintf( stderr, "[ERROR] THIS OPTION IS NOT YET IMPLEMENTED! TERMINATING...\n" );
exit( UNKNOWN_OPTION );
- data->N_f = 3 * system->N + 9;
+ data->N_f = 3.0 * system->N + 9.0;
*Evolve = &Velocity_Verlet_Berendsen_Isotropic_NPT;
if ( control->restart == FALSE )
{
data->therm.G_xi = control->Tau_T * (2.0 * data->E_Kin -
- data->N_f * K_B * control->T);
+ data->N_f * K_B / F_CONV * control->T);
data->therm.v_xi = data->therm.G_xi * control->dt;
data->iso_bar.eps = 1.0 / 3.0 * LOG( system->box.volume );
// data->inv_W = 1.0 / (data->N_f * K_B * control->T * SQR(control->Tau_P));
@@ -238,13 +238,13 @@ static void Init_Simulation_Data( reax_system *system, control_params *control,
/* semi-isotropic NPT */
case sNPT:
- data->N_f = 3 * system->N + 4;
+ data->N_f = 3.0 * system->N + 4.0;
*Evolve = &Velocity_Verlet_Berendsen_Semi_Isotropic_NPT;
break;
/* isotropic NPT */
case iNPT:
- data->N_f = 3 * system->N + 2;
+ data->N_f = 3.0 * system->N + 2.0;
*Evolve = &Velocity_Verlet_Berendsen_Isotropic_NPT;
break;
diff --git a/sPuReMD/src/integrate.c b/sPuReMD/src/integrate.c
index 23e725ba..fc1a95f1 100644
--- a/sPuReMD/src/integrate.c
+++ b/sPuReMD/src/integrate.c
@@ -32,7 +32,8 @@
#include "vector.h"
-/* Velocity Verlet integrator for microcanonical ensemble. */
+/* Perform simulation using the microcanonical ensemble with
+ * the Velocity Verlet integrator. */
void Velocity_Verlet_NVE( reax_system *system, control_params *control,
simulation_data *data, static_storage *workspace,
reax_list **lists, output_controls *out_control )
@@ -83,8 +84,8 @@ void Velocity_Verlet_NVE( reax_system *system, control_params *control,
}
-/* Velocity Verlet integrator for constant volume and temperature
- * with Berendsen thermostat. */
+/* Perform simulation using constant volume and temperature (NVT ensemble)
+ * with a Berendsen thermostat and the Velocity Verlet integrator. */
void Velocity_Verlet_Berendsen_NVT( reax_system* system,
control_params* control, simulation_data *data,
static_storage *workspace, reax_list **lists,
@@ -162,8 +163,8 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
}
-/* Velocity Verlet integrator for constant volume and constant temperature
- * with Nose-Hoover thermostat.
+/* Perform simulation using constant volume and temperature (NVT ensemble)
+ * with a Nose-Hoover thermostat and the Velocity Verlet integrator.
*
* Reference: Understanding Molecular Simulation, Frenkel and Smit
* Academic Press Inc. San Diego, 1996 p. 388-391 */
@@ -258,7 +259,8 @@ void Velocity_Verlet_Nose_Hoover_NVT_Klein( reax_system* system, control_params*
}
-/* Velocity Verlet integrator for constant pressure and constant temperature.
+/* Perform simulation using constant pressure and temperature (NPT ensemble)
+ * with a Berendsen thermostat and the Velocity Verlet integrator.
*
* NOTE: All box dimensions are scaled by the same amount, and
* there is no change in the angles between axes. */
@@ -369,7 +371,8 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
}
-/* Velocity Verlet integrator for constant pressure and constant temperature. */
+/* Perform simulation using constant pressure and temperature (NPT ensemble)
+ * with a Berendsen thermostat and the Velocity Verlet integrator. */
void Velocity_Verlet_Berendsen_Semi_Isotropic_NPT( reax_system* system,
control_params* control, simulation_data *data, static_storage *workspace,
reax_list **lists, output_controls *out_control )
@@ -510,7 +513,7 @@ void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system, control_params* contr
/* Compute zeta(t + dt/2), E_Kininetic(t + dt/2)
* IMPORTANT: What will be the initial value of zeta? and what is g? */
data->therm.xi += 0.5 * dt * control->Tau_T
- * (2.0 * data->E_Kin - data->N_f * K_B * control->T);
+ * (2.0 * data->E_Kin - data->N_f * K_B / F_CONV * control->T);
Reset( system, control, data, workspace );
fprintf(out_control->log, "reset-");
@@ -544,10 +547,10 @@ void Velocity_Verlet_Nose_Hoover_NVT( reax_system* system, control_params* contr
/* Compute zeta(t + dt) */
data->therm.xi += 0.5 * dt * control->Tau_T
- * (2.0 * data->E_Kin - data->N_f * K_B * control->T);
+ * (2.0 * data->E_Kin - data->N_f * K_B / F_CONV * control->T);
fprintf( out_control->log, "Xi: %8.3f %8.3f %8.3f\n",
- data->therm.xi, data->E_Kin, data->N_f * K_B * control->T );
+ data->therm.xi, data->E_Kin, data->N_f * K_B / F_CONV * control->T );
fflush( out_control->log );
}
diff --git a/sPuReMD/src/io_tools.c b/sPuReMD/src/io_tools.c
index 819d7678..24fda648 100644
--- a/sPuReMD/src/io_tools.c
+++ b/sPuReMD/src/io_tools.c
@@ -460,8 +460,10 @@ void Print_Near_Neighbors2( reax_system *system, control_params *control,
/* far nbrs contain only i-j nbrhood info, no j-i. */
-void Print_Far_Neighbors( reax_system *system, control_params *control,
- simulation_data *data, static_storage *workspace, reax_list **lists )
+void Print_Far_Neighbors( reax_system const * const system,
+ control_params const * const control,
+ simulation_data const * const data,
+ static_storage const * const workspace, reax_list ** const lists )
{
int i, pj, id_i, id_j;
char fname[MAX_STR];
@@ -482,7 +484,7 @@ void Print_Far_Neighbors( reax_system *system, control_params *control,
id_j = workspace->orig_id[far_nbrs->far_nbr_list[pj].nbr];
fprintf( fout, "%6d%6d%3d%3d%3d %12.7f\n",
- id_i, id_j,
+ id_j, id_i,
far_nbrs->far_nbr_list[pj].rel_box[0],
far_nbrs->far_nbr_list[pj].rel_box[1],
far_nbrs->far_nbr_list[pj].rel_box[2],
@@ -595,7 +597,7 @@ void Output_Results( reax_system *system, control_params *control,
#if defined(TEST_ENERGY) || defined(DEBUG) || defined(DEBUG_FOCUS)
fprintf( out_control->out,
"%-6d%24.15f%24.15f%24.15f%13.5f%13.5f%16.5f%13.5f%13.5f\n",
- data->step, data->E_Tot, data->E_Pot, E_CONV * data->E_Kin,
+ data->step, data->E_Tot, data->E_Pot, data->E_Kin,
data->therm.T, control->T, system->box.volume, data->iso_bar.P,
(control->P[0] + control->P[1] + control->P[2]) / 3.0 );
@@ -610,7 +612,7 @@ void Output_Results( reax_system *system, control_params *control,
#else
fprintf( out_control->out,
"%-6d%16.2f%16.2f%16.2f%11.2f%11.2f%13.2f%13.5f%13.5f\n",
- data->step, data->E_Tot, data->E_Pot, E_CONV * data->E_Kin,
+ data->step, data->E_Tot, data->E_Pot, data->E_Kin,
data->therm.T, control->T, system->box.volume, data->iso_bar.P,
(control->P[0] + control->P[1] + control->P[2]) / 3.0 );
diff --git a/sPuReMD/src/io_tools.h b/sPuReMD/src/io_tools.h
index 742a070c..7b37aeb2 100644
--- a/sPuReMD/src/io_tools.h
+++ b/sPuReMD/src/io_tools.h
@@ -32,8 +32,9 @@ char *Get_Atom_Name( reax_system*, int );
void Print_Near_Neighbors( reax_system*, control_params*, static_storage*,
reax_list** );
-void Print_Far_Neighbors( reax_system*, control_params*, simulation_data*,
- static_storage*, reax_list** );
+void Print_Far_Neighbors( reax_system const * const, control_params const * const,
+ simulation_data const * const, static_storage const * const,
+ reax_list ** const );
void Print_Total_Force( reax_system*, control_params*, simulation_data*,
static_storage*, reax_list**, output_controls* );
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index 4f006123..e95616c3 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -73,8 +73,8 @@
/* unit conversion from ??? to kcal / mol */
#define C_ELE (332.0638)
/* Boltzmann constant, AMU * A^2 / (ps^2 * K) */
- #define K_B (0.831687)
-// #define K_B (0.8314510)
+// #define K_B (0.831687)
+ #define K_B (0.8314510)
/* unit conversion for atomic force to AMU * A / ps^2 */
#define F_CONV (4.184e2)
/* energy conversion constant from electron volts to kilo-calories per mole */
@@ -122,7 +122,8 @@
/* unit conversion for atomic energy */
#if !defined(E_CONV)
/* AMU * Angstroms^2 / ps^2 --> kcal / mol */
- #define E_CONV (0.002391)
+// #define E_CONV (0.002391)
+ #define E_CONV (1.0 / 418.40)
#endif
/* energy conversion constant from electron volts to kilo-calories per mole */
#if !defined(EV_to_KCALpMOL)
diff --git a/sPuReMD/src/system_props.c b/sPuReMD/src/system_props.c
index 05f82000..57a060a8 100644
--- a/sPuReMD/src/system_props.c
+++ b/sPuReMD/src/system_props.c
@@ -196,10 +196,10 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data )
rvec_Scale( p, m, system->atoms[i].v );
data->E_Kin += rvec_Dot( p, system->atoms[i].v );
}
- data->E_Kin *= 0.5;
+ data->E_Kin *= 0.5 * E_CONV;
/* temperature scalar */
- data->therm.T = (2.0 * data->E_Kin) / (data->N_f * K_B);
+ data->therm.T = (2.0 * data->E_Kin) / (data->N_f * K_B / F_CONV);
/* avoid T being an absolute zero! */
if ( FABS( data->therm.T ) < ALMOST_ZERO )
@@ -209,17 +209,18 @@ void Compute_Kinetic_Energy( reax_system* system, simulation_data* data )
}
+/* Compute potential and total energies */
void Compute_Total_Energy( simulation_data* data )
{
data->E_Pot = data->E_BE + data->E_Ov + data->E_Un + data->E_Lp +
data->E_Ang + data->E_Pen + data->E_Coa + data->E_HB +
data->E_Tor + data->E_Con + data->E_vdW + data->E_Ele + data->E_Pol;
-
- data->E_Tot = data->E_Pot + data->E_Kin * E_CONV;
+ data->E_Tot = data->E_Pot + data->E_Kin;
}
+/* Check for numeric breakdowns in the energies */
void Check_Energy( simulation_data* data )
{
if ( IS_NAN_REAL(data->E_Pol) )
diff --git a/sPuReMD/src/tool_box.c b/sPuReMD/src/tool_box.c
index 45411f1b..0483e048 100644
--- a/sPuReMD/src/tool_box.c
+++ b/sPuReMD/src/tool_box.c
@@ -101,7 +101,10 @@ void Transform_to_UnitBox( rvec x1, simulation_box *box, int flag, rvec x2 )
}
-/* determine whether point p is inside the box */
+/* Check and remap (if necessary) an atom position to fall
+ * within the boundaries of a periodic simulation box where
+ * the boundaries are [0, d_i) with d_i being the length
+ * of the simulation box in a particular dimension */
void Fit_to_Periodic_Box( simulation_box *box, rvec p )
{
int i;
@@ -116,10 +119,10 @@ void Fit_to_Periodic_Box( simulation_box *box, rvec p )
p[i] += box->box_norms[i];
}
}
- else if ( p[i] > box->max[i] )
+ else if ( p[i] >= box->max[i] )
{
/* handle higher coords */
- while ( p[i] > box->max[i] )
+ while ( p[i] >= box->max[i] )
{
p[i] -= box->box_norms[i];
}
diff --git a/sPuReMD/src/traj.c b/sPuReMD/src/traj.c
index a19d4a89..896c0400 100644
--- a/sPuReMD/src/traj.c
+++ b/sPuReMD/src/traj.c
@@ -371,7 +371,7 @@ int Append_Custom_Frame( reax_system *system, control_params *control,
/* calculate total frame length*/
snprintf( buffer, SIZE, FRAME_GLOBALS,
data->step, data->time,
- data->E_Tot, data->E_Pot, E_CONV * data->E_Kin, data->therm.T,
+ data->E_Tot, data->E_Pot, data->E_Kin, data->therm.T,
P, system->box.volume,
system->box.box_norms[0],
system->box.box_norms[1],
@@ -620,7 +620,7 @@ int Append_xyz_Frame( reax_system *system, control_params *control,
out_control->write( out_control->trj, "%d\t%8.3f\t%8.3f\t%8.3f\t%8.3f\n",
data->step, data->E_Tot, data->E_Pot,
- E_CONV * data->E_Kin, data->therm.T );
+ data->E_Kin, data->therm.T );
for ( i = 0; i < system->N; ++i )
{
--
GitLab