From 9b99ad271f6f7d00de092e2575633ad55dc153b8 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Mon, 11 May 2020 14:34:00 -0600 Subject: [PATCH 1/5] Adding new ptemp keyword to fix nh --- doc/src/fix_nh.rst | 14 +++++++++++++- src/fix_nh.cpp | 25 ++++++++++++++++++++----- src/fix_nh.h | 11 +++++++++++ 3 files changed, 44 insertions(+), 6 deletions(-) mode change 100644 => 100755 doc/src/fix_nh.rst mode change 100644 => 100755 src/fix_nh.cpp mode change 100644 => 100755 src/fix_nh.h diff --git a/doc/src/fix_nh.rst b/doc/src/fix_nh.rst old mode 100644 new mode 100755 index 6927b47fcb..ee1c0e86b0 --- 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,16 @@ the :doc:`units ` command). 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_\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) `. 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 diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp old mode 100644 new mode 100755 index 339407f6d5..de066403df --- 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 From 4c6c4ac405b27757ea2a46976a7f408648a51c19 Mon Sep 17 00:00:00 2001 From: jtclemm <63308469+jtclemm@users.noreply.github.com> Date: Tue, 9 Jun 2020 12:20:21 -0600 Subject: [PATCH 2/5] Update fix_nh.rst Adding linebreaks for legibility --- doc/src/fix_nh.rst | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/doc/src/fix_nh.rst b/doc/src/fix_nh.rst index ee1c0e86b0..e4b2b6294a 100755 --- a/doc/src/fix_nh.rst +++ b/doc/src/fix_nh.rst @@ -216,9 +216,21 @@ The relaxation rate of the barostat is set by its inertia :math:`W`: 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) `. 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. +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) `. +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. +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 From 833ebeb0a0db34f4b25d6089282efb887c6bba27 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 9 Jun 2020 16:42:26 -0400 Subject: [PATCH 3/5] remove \mbox{} macros which are incompatible with `make pdf` --- doc/src/fix_nh.rst | 32 ++++++++++----------- doc/utils/sphinx-config/false_positives.txt | 1 + 2 files changed, 17 insertions(+), 16 deletions(-) diff --git a/doc/src/fix_nh.rst b/doc/src/fix_nh.rst index e4b2b6294a..a02ff3ae99 100755 --- a/doc/src/fix_nh.rst +++ b/doc/src/fix_nh.rst @@ -214,23 +214,23 @@ 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) `. -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. + W = (N + 1) k T_{target} P_{damp}^2 -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. +where :math:`N` is the number of atoms, :math:`k` is the Boltzmann constant, +and :math:`T_{target}` is the target temperature of the barostat :ref:`(Martyna) `. +If a thermostat is defined, :math:`T_{target}` is the target temperature +of the thermostat. If a thermostat is not defined, :math:`T_{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_{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_{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 diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index ae8eee9b2a..82c095d213 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 From 3ea39f5a23be82b6de0268155c3df69f5ce1705a Mon Sep 17 00:00:00 2001 From: jtclemm <63308469+jtclemm@users.noreply.github.com> Date: Tue, 9 Jun 2020 14:58:25 -0600 Subject: [PATCH 4/5] Update fix_nh.rst Removed mbox commands --- doc/src/fix_nh.rst | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/doc/src/fix_nh.rst b/doc/src/fix_nh.rst index a02ff3ae99..6f394dfa8d 100755 --- a/doc/src/fix_nh.rst +++ b/doc/src/fix_nh.rst @@ -214,23 +214,23 @@ The relaxation rate of the barostat is set by its inertia :math:`W`: .. math:: - W = (N + 1) k T_{target} P_{damp}^2 + 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. -where :math:`N` is the number of atoms, :math:`k` is the Boltzmann constant, -and :math:`T_{target}` is the target temperature of the barostat :ref:`(Martyna) `. -If a thermostat is defined, :math:`T_{target}` is the target temperature -of the thermostat. If a thermostat is not defined, :math:`T_{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_{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_{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. +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 From 2f5a8d093712cafad378969e3c541966de76771b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 9 Jun 2020 17:28:26 -0400 Subject: [PATCH 5/5] remove trailing whitespace --- doc/src/fix_nh.rst | 30 +++++++++++++++--------------- src/fix_nh.cpp | 2 +- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/doc/src/fix_nh.rst b/doc/src/fix_nh.rst index 6f394dfa8d..53eb1dfd25 100755 --- a/doc/src/fix_nh.rst +++ b/doc/src/fix_nh.rst @@ -215,22 +215,22 @@ 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. +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 diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index de066403df..014a2b8855 100755 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -775,7 +775,7 @@ void FixNH::setup(int /*vflag*/) // 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 - // error if T less than 1e-6 + // error if T less than 1e-6 // if it was read in from a restart file, leave it be if (t0 == 0.0) {