diff --git a/doc/src/fix_nh.rst b/doc/src/fix_nh.rst old mode 100644 new mode 100755 index 636d70140e..9b73a8b4ed --- a/doc/src/fix_nh.rst +++ b/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 ` 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) `. +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 diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index ffe18f5cb1..a517ce290d 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -135,6 +135,7 @@ atan atc AtC ATC +athermal athomps atm atomeye diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp old mode 100644 new mode 100755 index 339407f6d5..014a2b8855 --- a/src/fix_nh.cpp +++ b/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; diff --git a/src/fix_nh.h b/src/fix_nh.h old mode 100644 new mode 100755 index 807ea2a71f..624b0c8bb4 --- a/src/fix_nh.h +++ b/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