Merge pull request #2137 from jtclemm/ptemp

Adding new ptemp keyword to fix/nh

doc/src/fix_nh.rst

@@ -46,7 +46,7 @@ Syntax
 
   .. parsed-literal::
 
-     keyword = *temp* or *iso* or *aniso* or *tri* or *x* or *y* or *z* or *xy* or *yz* or *xz* or *couple* or *tchain* or *pchain* or *mtk* or *tloop* or *ploop* or *nreset* or *drag* or *dilate* or *scalexy* or *scaleyz* or *scalexz* or *flip* or *fixedpoint* or *update*
+     keyword = *temp* or *iso* or *aniso* or *tri* or *x* or *y* or *z* or *xy* or *yz* or *xz* or *couple* or *tchain* or *pchain* or *mtk* or *tloop* or *ploop* or *nreset* or *drag* or *ptemp* or *dilate* or *scalexy* or *scaleyz* or *scalexz* or *flip* or *fixedpoint* or *update*
        *temp* values = Tstart Tstop Tdamp
          Tstart,Tstop = external temperature at start/end of run
          Tdamp = temperature damping parameter (time units)
@@ -69,6 +69,8 @@ Syntax
        *nreset* value = reset reference cell every this many timesteps
        *drag* value = Df
          Df = drag factor added to barostat/thermostat (0.0 = no drag)
+       *ptemp* value = Ttarget
+         Ttarget = target temperature for barostat
        *dilate* value = dilate-group-ID
          dilate-group-ID = only dilate atoms in this group due to barostat volume changes
        *scalexy* value = *yes* or *no* = scale xy with ly
@@ -208,6 +210,28 @@ a timespan of (roughly) 10 time units (e.g. :math:`\tau` or fs or ps
    time units, and that timesteps are NOT the same as time units for most
    :doc:`units <units>` settings.
 
+The relaxation rate of the barostat is set by its inertia :math:`W`:
+
+.. math::
+
+   W = (N + 1) k T_{\rm target} P_{\rm damp}^2
+
+where :math:`N` is the number of atoms, :math:`k` is the Boltzmann constant,
+and :math:`T_{\rm target}` is the target temperature of the barostat :ref:`(Martyna) <nh-Martyna>`.
+If a thermostat is defined, :math:`T_{\rm target}` is the target temperature
+of the thermostat. If a thermostat is not defined, :math:`T_{\rm target}`
+is set to the current temperature of the system when the barostat is initialized.
+If this temperature is too low the simulation will quit with an error.
+Note: in previous versions of LAMMPS, :math:`T_{\rm target}` would default to
+a value of 1.0 for *lj* units and 300.0 otherwise if the system had a temperature
+of exactly zero.
+
+If a thermostat is not specified by this fix, :math:`T_{\rm target}` can be
+manually specified using the *Ptemp* parameter. This may be useful if the
+barostat is initialized when the current temperature does not reflect the
+steady state temperature of the system. This keyword may also be useful in
+athermal simulations where the temperature is not well defined.
+
 Regardless of what atoms are in the fix group (the only atoms which
 are time integrated), a global pressure or stress tensor is computed
 for all atoms.  Similarly, when the size of the simulation box is

doc/utils/sphinx-config/false_positives.txt

@@ -135,6 +135,7 @@ atan
 atc
 AtC
 ATC
+athermal
 athomps
 atm
 atomeye

src/fix_nh.cpp

@@ -39,6 +39,7 @@ using namespace FixConst;
 
 #define DELTAFLIP 0.1
 #define TILTMAX 1.5
+#define EPSILON 1.0e-6
 
 enum{NOBIAS,BIAS};
 enum{NONE,XYZ,XY,YZ,XZ};
@@ -83,6 +84,7 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) :
   flipflag = 1;
   dipole_flag = 0;
   dlm_flag = 0;
+  p_temp_flag = 0;
 
   tcomputeflag = 0;
   pcomputeflag = 0;
@@ -259,6 +261,12 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) :
       drag = force->numeric(FLERR,arg[iarg+1]);
       if (drag < 0.0) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       iarg += 2;
+    } else if (strcmp(arg[iarg],"ptemp") == 0) {
+      if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
+      p_temp = force->numeric(FLERR,arg[iarg+1]);
+      p_temp_flag = 1;
+      if (p_temp <= 0.0) error->all(FLERR,"Illegal fix nvt/npt/nph command");
+      iarg += 2;
     } else if (strcmp(arg[iarg],"dilate") == 0) {
       if (iarg+2 > narg) error->all(FLERR,"Illegal fix nvt/npt/nph command");
       if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
@@ -466,6 +474,10 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) :
       (p_flag[5] && p_period[5] <= 0.0))
     error->all(FLERR,"Fix nvt/npt/nph damping parameters must be > 0.0");
 
+  // check that ptemp is not defined with a thermostat
+  if (tstat_flag && p_temp_flag)
+    error->all(FLERR,"Thermostat in fix nvt/npt/nph is incompatible with ptemp command");
+
   // set pstat_flag and box change and restart_pbc variables
 
   pre_exchange_flag = 0;
@@ -760,16 +772,19 @@ void FixNH::setup(int /*vflag*/)
   } else if (pstat_flag) {
 
     // t0 = reference temperature for masses
+    // set equal to either ptemp or the current temperature
     // cannot be done in init() b/c temperature cannot be called there
     // is b/c Modify::init() inits computes after fixes due to dof dependence
-    // guesstimate a unit-dependent t0 if actual T = 0.0
+    // error if T less than 1e-6
     // if it was read in from a restart file, leave it be
 
     if (t0 == 0.0) {
-      t0 = temperature->compute_scalar();
+      if(p_temp_flag) {
-      if (t0 == 0.0) {
+        t0 = p_temp;
-        if (strcmp(update->unit_style,"lj") == 0) t0 = 1.0;
+      } else {
-        else t0 = 300.0;
+        t0 = temperature->compute_scalar();
+        if(t0 < EPSILON)
+          error->all(FLERR, "Current temperature too close to zero, consider using ptemp setting");
       }
     }
     t_target = t0;

src/fix_nh.h

@@ -71,6 +71,9 @@ class FixNH : public Fix {
   char *id_dilate;                 // group name to dilate
   class Irregular *irregular;      // for migrating atoms after box flips
 
+  double p_temp;                   // target temperature for barostat
+  int p_temp_flag;
+
   int nlevels_respa;
   double *step_respa;
 
@@ -226,6 +229,10 @@ E: Fix nvt/npt/nph damping parameters must be > 0.0
 
 Self-explanatory.
 
+E: Thermostat in fix nvt/npt/nph is incompatible with ptemp command
+
+Self-explanatory.
+
 E: Cannot use fix npt and fix deform on same component of stress tensor
 
 This would be changing the same box dimension twice.
@@ -238,6 +245,10 @@ E: Pressure ID for fix npt/nph does not exist
 
 Self-explanatory.
 
+E: Current temperature too close to zero, consider using ptemp setting
+
+The current temperature is close to zero and may cause numerical instability. The user may want to specify a different target temperature using the ptemp setting.
+
 E: Non-numeric pressure - simulation unstable
 
 UNDOCUMENTED