From 72c9f2ac25ca23b65ffde9cd3ef6f63feac08952 Mon Sep 17 00:00:00 2001
From: "Kurt A. O'Hearn" <ohearnku@msu.edu>
Date: Thu, 10 Dec 2020 11:24:58 -0500
Subject: [PATCH] sPuReMD: fixes to Berendsen thermostat for NPT ensembles.
---
sPuReMD/src/integrate.c | 26 ++++++++++++++++++--------
sPuReMD/src/reax_types.h | 13 ++++++++++---
2 files changed, 28 insertions(+), 11 deletions(-)
diff --git a/sPuReMD/src/integrate.c b/sPuReMD/src/integrate.c
index 12e53a4a..766ea766 100644
--- a/sPuReMD/src/integrate.c
+++ b/sPuReMD/src/integrate.c
@@ -83,7 +83,7 @@ void Velocity_Verlet_NVE( reax_system *system, control_params *control,
}
-/* Velocity Verlet integrator for constant volume and constant temperature
+/* Velocity Verlet integrator for constant volume and temperature
* with Berendsen thermostat. */
void Velocity_Verlet_Berendsen_NVT( reax_system* system,
control_params* control, simulation_data *data,
@@ -133,9 +133,9 @@ void Velocity_Verlet_Berendsen_NVT( reax_system* system,
rvec_ScaledAdd( system->atoms[i].v, scalar1 * inv_m, system->atoms[i].f );
}
- /* temperature scaler */
Compute_Kinetic_Energy( system, data );
+ /* temperature scaler */
lambda = 1.0 + ((control->dt * 1.0e-12) / control->Tau_T)
* (control->T / data->therm.T - 1.0);
@@ -334,16 +334,20 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
mu = (mu_3[0] + mu_3[1] + mu_3[2]) / 3.0;
/* temperature scaler */
- lambda = 1.0 + ((dt * 1.0e-10) / control->Tau_T) * (control->T / data->therm.T - 1.0);
+ lambda = 1.0 + ((dt * 1.0e-12) / control->Tau_T)
+ * (control->T / data->therm.T - 1.0);
+
if ( lambda < MIN_dT )
{
lambda = MIN_dT;
}
- else if (lambda > MAX_dT )
+
+ lambda = SQRT( lambda );
+
+ if ( lambda > MAX_dT )
{
lambda = MAX_dT;
}
- lambda = SQRT( lambda );
/* Scale velocities and positions at t+dt */
for ( i = 0; i < system->N; ++i )
@@ -358,6 +362,7 @@ void Velocity_Verlet_Berendsen_Isotropic_NPT( reax_system* system,
rvec_Scale( system->atoms[i].x, mu, system->atoms[i].x );
}
+ /* update kinetic energy and temperature based on new velocities */
Compute_Kinetic_Energy( system, data );
Update_Box_Isotropic( &system->box, mu );
@@ -435,16 +440,20 @@ void Velocity_Verlet_Berendsen_Semi_Isotropic_NPT( reax_system* system,
}
/* temperature scaler */
- lambda = 1.0 + ((dt * 1.0e-10) / control->Tau_T) * (control->T / data->therm.T - 1.0);
+ lambda = 1.0 + ((dt * 1.0e-12) / control->Tau_T)
+ * (control->T / data->therm.T - 1.0);
+
if ( lambda < MIN_dT )
{
lambda = MIN_dT;
}
- else if ( lambda > MAX_dT )
+
+ lambda = SQRT( lambda );
+
+ if ( lambda > MAX_dT )
{
lambda = MAX_dT;
}
- lambda = SQRT( lambda );
/* Scale velocities and positions at t+dt */
for ( i = 0; i < system->N; ++i )
@@ -459,6 +468,7 @@ void Velocity_Verlet_Berendsen_Semi_Isotropic_NPT( reax_system* system,
rvec_Multiply( system->atoms[i].x, mu, system->atoms[i].x );
}
+ /* update kinetic energy and temperature based on new velocities */
Compute_Kinetic_Energy( system, data );
Update_Box_Semi_Isotropic( &system->box, mu );
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index df75dcfb..b030ac67 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -229,12 +229,19 @@
/* ensemble type */
enum ensemble
{
+ /* microcanonical ensemble */
NVE = 0,
+ /* Berendsen NVT ensemble */
bNVT = 1,
+ /* Nose-Hoover NVT ensemble */
nhNVT = 2,
+ /* Parrinello-Rehman-Nose-Hoover semi-isotropic NPT ensemble */
sNPT = 3,
+ /* Parrinello-Rehman-Nose-Hoover isotropic NPT ensemble */
iNPT = 4,
+ /* Parrinello-Rehman-Nose-Hoover anisotropic NPT ensemble */
aNPT = 5,
+ /* total number of ensemble types */
ens_N = 6,
};
@@ -860,12 +867,12 @@ struct control_params
* (smallest point) in the simulation box */
int reposition_atoms;
/* ensemble values:
- * 0 : NVE
+ * 0 : microcanonical ensemble (NVE)
* 1 : Berendsen NVT (bNVT)
* 2 : Nose-Hoover NVT (nhNVT)
* 3 : Parrinello-Rehman-Nose-Hoover semi-isotropic NPT (sNPT)
- * 4 : Parrinello-Rehman-Nose-Hoover isotropic (iNPT)
- * 5 : Parrinello-Rehman-Nose-Hoover anisotropic (aNPT) */
+ * 4 : Parrinello-Rehman-Nose-Hoover isotropic NPT (iNPT)
+ * 5 : Parrinello-Rehman-Nose-Hoover anisotropic NPT (aNPT) */
int ensemble;
/* number of simulation time steps */
int nsteps;
--
GitLab