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,16 @@ the :doc:`units <units>` command).
    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_\mbox{target} P_\mbox{damp}^2
+  
+where :math:`N` is the number of atoms, :math:`k` is the Boltzmann constant, and :math:`T_\mbox{target}` is the target temperature of the barostat :ref:`(Martyna) <nh-Martyna>`. If a thermostat is defined, :math:`T_\mbox{target}` is the target temperature of the thermostat. If a thermostat is not defined, :math:`T_\mbox{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_\mbox{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_\mbox{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

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 (t0 == 0.0) {
-        if (strcmp(update->unit_style,"lj") == 0) t0 = 1.0;
-        else t0 = 300.0;
+      if(p_temp_flag) {
+        t0 = p_temp;
+      } else {
+        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