Merge pull request #3582 from athomps/sllod_variants

Implement sllod variants

doc/src/fix_nvt_sllod.rst

@@ -17,7 +17,15 @@ Syntax
 
 * ID, group-ID are documented in :doc:`fix <fix>` command
 * nvt/sllod = style name of this fix command
-* additional thermostat related keyword/value pairs from the :doc:`fix nvt <fix_nh>` command can be appended
+* zero or more keyword/value pairs may be appended
+
+  .. parsed-literal::
+
+    keyword = *psllod*
+      *psllod* value = *no* or *yes* = use SLLOD or p-SLLOD variant, respectively
+
+* additional thermostat related keyword/value pairs from the :doc:`fix nvt <fix_nh>`
+  command can be appended
 
 Examples
 """"""""
@@ -61,22 +69,27 @@ large x velocity.  These velocities do not contribute to the thermal
    consistent.
 
 The SLLOD equations of motion, originally proposed by Hoover and Ladd
-(see :ref:`(Evans and Morriss) <Evans3>`), were proven to be
+(see :ref:`(Evans and Morriss) <Evans3>`), were proven to be equivalent
-equivalent to Newton's equations of motion for shear flow by
+to Newton's equations of motion for shear flow by :ref:`(Evans and
-:ref:`(Evans and Morriss) <Evans3>`. They were later shown to generate
+Morriss) <Evans3>`. They were later shown to generate the desired
-the desired velocity gradient and the correct production of work by
+velocity gradient and the correct production of work by stresses for all
-stresses for all forms of homogeneous flow by :ref:`(Daivis and Todd)
+forms of homogeneous flow by :ref:`(Daivis and Todd) <Daivis>`.
-<Daivis>`.
+
-The LAMMPS implementation corresponds to the p-SLLOD variant
+.. versionchanged:: TBD
-of SLLOD. :ref:`(Edwards) <Edwards>`.
+
-The equations of motion are coupled to a
+For the default (*psllod* = *no*), the LAMMPS implementation adheres to
+the standard SLLOD equations of motion, as defined by :ref:`(Evans and
+Morriss) <Evans3>`.  The option *psllod* = *yes* invokes the slightly
+different SLLOD variant first introduced by :ref:`(Tuckerman et al.)
+<Tuckerman>` as g-SLLOD and later by :ref:`(Edwards) <Edwards>` as
+p-SLLOD.  In all cases, the equations of motion are coupled to a
 Nose/Hoover chain thermostat in a velocity Verlet formulation, closely
 following the implementation used for the :doc:`fix nvt <fix_nh>`
 command.
 
 .. note::
 
-   A recent (2017) book by :ref:`(Daivis and Todd) <Daivis-sllod>`
+   A recent (2017) book by :ref:`(Todd and Daivis) <Todd-sllod>`
    discusses use of the SLLOD method and non-equilibrium MD (NEMD)
    thermostatting generally, for both simple and complex fluids,
    e.g. molecular systems.  The latter can be tricky to do correctly.
@@ -159,7 +172,7 @@ Restrictions
 """"""""""""
 
 This fix works best without Nose-Hoover chain thermostats, i.e. using
-tchain = 1.  Setting tchain to larger values can result in poor
+*tchain* = 1.  Setting *tchain* to larger values can result in poor
 equilibration.
 
 Related commands
@@ -171,7 +184,7 @@ Related commands
 Default
 """""""
 
-Same as :doc:`fix nvt <fix_nh>`, except tchain = 1.
+Same as :doc:`fix nvt <fix_nh>`, except *tchain* = 1, psllod = *no*.
 
 ----------
 
@@ -183,11 +196,16 @@ Same as :doc:`fix nvt <fix_nh>`, except tchain = 1.
 
 **(Daivis and Todd)** Daivis and Todd, J Chem Phys, 124, 194103 (2006).
 
+.. _Todd-sllod:
+
+**(Todd and Daivis)** Todd and Daivis, Nonequilibrium Molecular Dynamics (book),
+Cambridge University Press, (2017) https://doi.org/10.1017/9781139017848.
+
+.. _Tuckerman:
+
+**(Tuckerman et al.)** Tuckerman, Mundy, Balasubramanian, and Klein, J Chem Phys 106, 5615 (1997).
+
 .. _Edwards:
 
 **(Edwards)** Edwards, Baig, and Keffer, J Chem Phys 124, 194104 (2006).
 
-.. _Daivis-sllod:
-
-**(Daivis and Todd)** Daivis and Todd, Nonequilibrium Molecular Dynamics (book),
-Cambridge University Press, https://doi.org/10.1017/9781139017848, (2017).

doc/src/fix_nvt_sllod_eff.rst

@@ -12,7 +12,15 @@ Syntax
 
 * ID, group-ID are documented in :doc:`fix <fix>` command
 * nvt/sllod/eff = style name of this fix command
-* additional thermostat related keyword/value pairs from the :doc:`fix nvt/eff <fix_nh_eff>` command can be appended
+* zero or more keyword/value pairs may be appended
+
+  .. parsed-literal::
+
+     keyword = *psllod*
+       *psllod* value = *no* or *yes* = use SLLOD or p-SLLOD variant, respectively
+
+* additional thermostat related keyword/value pairs from the :doc:`fix
+  nvt/eff <fix_nh_eff>` command may be appended, too.
 
 Examples
 """"""""
@@ -25,18 +33,20 @@ Examples
 Description
 """""""""""
 
-Perform constant NVT integration to update positions and velocities
+Perform constant NVT integration to update positions and velocities each
-each timestep for nuclei and electrons in the group for the :doc:`electron force field <pair_eff>` model, using a Nose/Hoover temperature
+timestep for nuclei and electrons in the group for the :doc:`electron
+force field <pair_eff>` model, using a Nose/Hoover temperature
 thermostat.  V is volume; T is temperature.  This creates a system
 trajectory consistent with the canonical ensemble.
 
-The operation of this fix is exactly like that described by the :doc:`fix nvt/sllod <fix_nvt_sllod>` command, except that the radius and
+The operation of this fix is exactly like that described by the
+:doc:`fix nvt/sllod <fix_nvt_sllod>` command, except that the radius and
 radial velocity of electrons are also updated and thermostatted.
 Likewise the temperature calculated by the fix, using the compute it
 creates (as discussed in the :doc:`fix nvt, npt, and nph <fix_nh>` doc
-page), is performed with a :doc:`compute temp/deform/eff <compute_temp_deform_eff>` command that includes
+page), is performed with a :doc:`compute temp/deform/eff
-the eFF contribution to the temperature from the electron radial
+<compute_temp_deform_eff>` command that includes the eFF contribution to
-velocity.
+the temperature from the electron radial velocity.
 
 Restart, fix_modify, output, run start/stop, minimize info
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""

doc/utils/sphinx-config/false_positives.txt

@@ -2854,6 +2854,7 @@ pscrozi
 Pseudocode
 pseudodynamics
 pseudopotential
+psllod
 pSp
 pSPICA
 Pstart

src/EFF/fix_nh_eff.cpp

@@ -33,7 +33,7 @@ enum{NOBIAS,BIAS};
 FixNHEff::FixNHEff(LAMMPS *lmp, int narg, char **arg) : FixNH(lmp, narg, arg)
 {
   if (!atom->electron_flag)
-    error->all(FLERR,"Fix nvt/nph/npt/eff requires atom style electron");
+    error->all(FLERR,"Fix {} requires atom style electron", style);
 }
 
 /* ----------------------------------------------------------------------

src/EFF/fix_nvt_sllod_eff.cpp

@@ -42,8 +42,21 @@ FixNVTSllodEff::FixNVTSllodEff(LAMMPS *lmp, int narg, char **arg) :
 
   // default values
 
+  psllod_flag = 0;
   if (mtchain_default_flag) mtchain = 1;
 
+  // select SLLOD/p-SLLOD/g-SLLOD variant
+
+  int iarg = 3;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"psllod") == 0) {
+      if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix nvt/sllod/eff psllod", error);
+      psllod_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
+      iarg += 2;
+    } else iarg++;
+  }
+
   // create a new compute temp style
   // id = fix-ID + temp
 
@@ -108,10 +121,11 @@ void FixNVTSllodEff::nh_v_temp()
 
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
-      temperature->remove_bias(i,v[i]);
+      if (!psllod_flag) temperature->remove_bias(i,v[i]);
       vdelu[0] = h_two[0]*v[i][0] + h_two[5]*v[i][1] + h_two[4]*v[i][2];
       vdelu[1] = h_two[1]*v[i][1] + h_two[3]*v[i][2];
       vdelu[2] = h_two[2]*v[i][2];
+      if (psllod_flag) temperature->remove_bias(i,v[i]);
       v[i][0] = v[i][0]*factor_eta - dthalf*vdelu[0];
       v[i][1] = v[i][1]*factor_eta - dthalf*vdelu[1];
       v[i][2] = v[i][2]*factor_eta - dthalf*vdelu[2];

src/EFF/fix_nvt_sllod_eff.h

@@ -31,6 +31,7 @@ class FixNVTSllodEff : public FixNHEff {
 
  private:
   int nondeformbias;
+  int psllod_flag;
 
   void nh_v_temp() override;
 };

src/INTEL/fix_nvt_sllod_intel.cpp

@@ -35,14 +35,26 @@ FixNVTSllodIntel::FixNVTSllodIntel(LAMMPS *lmp, int narg, char **arg) :
   FixNHIntel(lmp, narg, arg)
 {
   if (!tstat_flag)
-    error->all(FLERR,"Temperature control must be used with fix nvt/sllod");
+    error->all(FLERR,"Temperature control must be used with fix nvt/sllod/intel");
   if (pstat_flag)
-    error->all(FLERR,"Pressure control can not be used with fix nvt/sllod");
+    error->all(FLERR,"Pressure control can not be used with fix nvt/sllod/intel");
 
   // default values
 
+  psllod_flag = 0;
   if (mtchain_default_flag) mtchain = 1;
 
+  // select SLLOD/p-SLLOD/g-SLLOD variant
+
+  int iarg = 3;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"psllod") == 0) {
+      if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix nvt/sllod/intel psllod", error);
+      psllod_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
+      iarg += 2;
+    } else iarg++;
+  }
 
   // create a new compute temp style
   // id = fix-ID + temp
@@ -61,23 +73,21 @@ void FixNVTSllodIntel::init()
   FixNHIntel::init();
 
   if (!temperature->tempbias)
-    error->all(FLERR,"Temperature for fix nvt/sllod does not have a bias");
+    error->all(FLERR,"Temperature for fix nvt/sllod/intel does not have a bias");
 
   nondeformbias = 0;
   if (strcmp(temperature->style,"temp/deform") != 0) nondeformbias = 1;
 
   // check fix deform remap settings
 
-  int i;
+  auto deform = modify->get_fix_by_style("^deform");
-  for (i = 0; i < modify->nfix; i++)
+  if (deform.size() < 1) error->all(FLERR,"Using fix nvt/sllod/intel with no fix deform defined");
-    if (strncmp(modify->fix[i]->style,"deform",6) == 0) {
+
-      if ((dynamic_cast<FixDeform *>(modify->fix[i]))->remapflag != Domain::V_REMAP)
+  for (auto ifix : deform) {
-        error->all(FLERR,"Using fix nvt/sllod with inconsistent fix deform "
+    auto f = dynamic_cast<FixDeform *>(ifix);
-                   "remap option");
+    if (f && (f->remapflag != Domain::V_REMAP))
-      break;
+      error->all(FLERR,"Using fix nvt/sllod/intel with inconsistent fix deform remap option");
-    }
+  }
-  if (i == modify->nfix)
-    error->all(FLERR,"Using fix nvt/sllod with no fix deform defined");
 }
 
 /* ----------------------------------------------------------------------
@@ -105,10 +115,11 @@ void FixNVTSllodIntel::nh_v_temp()
 
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
+      if (!psllod_flag) temperature->remove_bias(i,v[i]);
       vdelu[0] = h_two[0]*v[i][0] + h_two[5]*v[i][1] + h_two[4]*v[i][2];
       vdelu[1] = h_two[1]*v[i][1] + h_two[3]*v[i][2];
       vdelu[2] = h_two[2]*v[i][2];
-      temperature->remove_bias(i,v[i]);
+      if (psllod_flag) temperature->remove_bias(i,v[i]);
       v[i][0] = v[i][0]*factor_eta - dthalf*vdelu[0];
       v[i][1] = v[i][1]*factor_eta - dthalf*vdelu[1];
       v[i][2] = v[i][2]*factor_eta - dthalf*vdelu[2];

src/INTEL/fix_nvt_sllod_intel.h

@@ -36,6 +36,7 @@ class FixNVTSllodIntel : public FixNHIntel {
 
  private:
   int nondeformbias;
+  int psllod_flag;     // 0 for SLLOD, 1 for p-SLLOD
 
   void nh_v_temp() override;
 };

src/KOKKOS/fix_nvt_sllod_kokkos.cpp

@@ -47,12 +47,25 @@ FixNVTSllodKokkos<DeviceType>::FixNVTSllodKokkos(LAMMPS *lmp, int narg, char **a
   this->domainKK = (DomainKokkos *) this->domain;
 
   if (!this->tstat_flag)
-    this->error->all(FLERR,"Temperature control must be used with fix nvt/kk");
+    this->error->all(FLERR,"Temperature control must be used with fix nvt/sllod/kk");
   if (this->pstat_flag)
-    this->error->all(FLERR,"Pressure control can not be used with fix nvt/kk");
+    this->error->all(FLERR,"Pressure control can not be used with fix nvt/sllod/kk");
 
+  this->psllod_flag = 0;
   if (this->mtchain_default_flag) this->mtchain = 1;
 
+  // select SLLOD/p-SLLOD/g-SLLOD variant
+
+  int iarg = 3;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"psllod") == 0) {
+      if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix nvt/sllod/kk psllod", this->error);
+      this->psllod_flag = utils::logical(FLERR,arg[iarg+1],false,this->lmp);
+      iarg += 2;
+    } else iarg++;
+  }
+
   this->id_temp = utils::strdup(std::string(this->id)+"_temp");
   this->modify->add_compute(fmt::format("{} all temp/deform/kk",this->id_temp));
   this->tcomputeflag = 1;
@@ -67,22 +80,22 @@ void FixNVTSllodKokkos<DeviceType>::init()
   FixNHKokkos<DeviceType>::init();
 
   if (!this->temperature->tempbias)
-    this->error->all(FLERR,"Temperature for fix nvt/sllod does not have a bias");
+    this->error->all(FLERR,"Temperature for fix nvt/sllod/kk does not have a bias");
 
   nondeformbias = 0;
   if (utils::strmatch(this->temperature->style,"^temp/deform")) nondeformbias = 1;
 
   // check fix deform remap settings
 
-  int i;
+  auto deform = this->modify->get_fix_by_style("^deform");
-  for (i = 0; i < this->modify->nfix; i++)
+  if (deform.size() < 1)
-    if (utils::strmatch(this->modify->fix[i]->style,"^deform")) {
+    this->error->all(FLERR,"Using fix nvt/sllod/kk with no fix deform defined");
-      if (((FixDeform *) this->modify->fix[i])->remapflag != Domain::V_REMAP)
+
-        this->error->all(FLERR,"Using fix nvt/sllod with inconsistent fix deform remap option");
+  for (auto ifix : deform) {
-      break;
+    auto f = dynamic_cast<FixDeform *>(ifix);
-    }
+    if (f && (f->remapflag != Domain::V_REMAP))
-  if (i == this->modify->nfix)
+      this->error->all(FLERR,"Using fix ntv/sllod/kk with inconsistent fix deform remap option");
-    this->error->all(FLERR,"Using fix nvt/sllod with no fix deform defined");
+  }
 }
 
 /* ----------------------------------------------------------------------
@@ -104,6 +117,7 @@ void FixNVTSllodKokkos<DeviceType>::nh_v_temp()
     this->temperature->compute_scalar();
     atomKK->modified(this->temperature->execution_space,this->temperature->datamask_modify);
   }
+
   v = atomKK->k_v.view<DeviceType>();
   mask = atomKK->k_mask.view<DeviceType>();
   int nlocal = atomKK->nlocal;
@@ -117,17 +131,35 @@ void FixNVTSllodKokkos<DeviceType>::nh_v_temp()
   if (vdelu.extent(0) < atomKK->nmax)
     vdelu = typename AT::t_v_array(Kokkos::NoInit("nvt/sllod/kk:vdelu"), atomKK->nmax);
 
+  if (!this->psllod_flag) {
+    atomKK->sync(this->temperature->execution_space,this->temperature->datamask_read);
+    this->temperature->remove_bias_all();
+    atomKK->modified(this->temperature->execution_space,this->temperature->datamask_modify);
+  }
+
+  atomKK->sync(this->execution_space,V_MASK | MASK_MASK);
+
   this->copymode = 1;
   Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagFixNVTSllod_temp1>(0,nlocal),*this);
   this->copymode = 0;
 
-  this->temperature->remove_bias_all();
+  if (this->psllod_flag) {
+    atomKK->sync(this->temperature->execution_space,this->temperature->datamask_read);
+    this->temperature->remove_bias_all();
+    atomKK->modified(this->temperature->execution_space,this->temperature->datamask_modify);
+  }
+
+  atomKK->sync(this->execution_space,V_MASK | MASK_MASK);
 
   this->copymode = 1;
   Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagFixNVTSllod_temp2>(0,nlocal),*this);
   this->copymode = 0;
 
+  atomKK->modified(this->execution_space,V_MASK);
+
+  atomKK->sync(this->temperature->execution_space,this->temperature->datamask_read);
   this->temperature->restore_bias_all();
+  atomKK->modified(this->temperature->execution_space,this->temperature->datamask_modify);
 }
 
 template<class DeviceType>
@@ -157,4 +189,3 @@ template class FixNVTSllodKokkos<LMPHostType>;
 #endif
 }
 
-

src/KOKKOS/fix_nvt_sllod_kokkos.h

@@ -50,6 +50,7 @@ class FixNVTSllodKokkos : public FixNHKokkos<DeviceType> {
 
  private:
   int nondeformbias;
+  int psllod_flag;
 
   void nh_v_temp() override;
 

src/OPENMP/fix_nvt_sllod_omp.cpp

@@ -43,14 +43,26 @@ FixNVTSllodOMP::FixNVTSllodOMP(LAMMPS *lmp, int narg, char **arg) :
   FixNHOMP(lmp, narg, arg)
 {
   if (!tstat_flag)
-    error->all(FLERR,"Temperature control must be used with fix nvt/sllod");
+    error->all(FLERR,"Temperature control must be used with fix nvt/sllod/omp");
   if (pstat_flag)
-    error->all(FLERR,"Pressure control can not be used with fix nvt/sllod");
+    error->all(FLERR,"Pressure control can not be used with fix nvt/sllod/omp");
 
   // default values
 
+  psllod_flag = 0;
   if (mtchain_default_flag) mtchain = 1;
 
+  // select SLLOD/p-SLLOD/g-SLLOD variant
+
+  int iarg = 3;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"psllod") == 0) {
+      if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix nvt/sllod/omp psllod", error);
+      psllod_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
+      iarg += 2;
+    } else iarg++;
+  }
 
   // create a new compute temp style
   // id = fix-ID + temp
@@ -75,15 +87,14 @@ void FixNVTSllodOMP::init()
 
   // check fix deform remap settings
 
-  int i;
+  auto deform = modify->get_fix_by_style("^deform");
-  for (i = 0; i < modify->nfix; i++)
+  if (deform.size() < 1) error->all(FLERR,"Using fix nvt/sllod/omp with no fix deform defined");
-    if (utils::strmatch(modify->fix[i]->style,"^deform")) {
+
-      if ((dynamic_cast<FixDeform *>(modify->fix[i]))->remapflag != Domain::V_REMAP)
+  for (auto ifix : deform) {
-        error->all(FLERR,"Using fix nvt/sllod/omp with inconsistent fix deform remap option");
+    auto f = dynamic_cast<FixDeform *>(ifix);
-      break;
+    if (f && (f->remapflag != Domain::V_REMAP))
-    }
+      error->all(FLERR,"Using fix nvt/sllod/omp with inconsistent fix deform remap option");
-  if (i == modify->nfix)
+  }
-    error->all(FLERR,"Using fix nvt/sllod/omp with no fix deform defined");
 }
 
 /* ----------------------------------------------------------------------
@@ -114,10 +125,11 @@ void FixNVTSllodOMP::nh_v_temp()
   for (int i = 0; i < nlocal; i++) {
     double vdelu0,vdelu1,vdelu2,buf[3];
     if (mask[i] & groupbit) {
+      if (!psllod_flag) temperature->remove_bias_thr(i,&v[i].x,buf);
       vdelu0 = h_two[0]*v[i].x + h_two[5]*v[i].y + h_two[4]*v[i].z;
       vdelu1 = h_two[1]*v[i].y + h_two[3]*v[i].z;
       vdelu2 = h_two[2]*v[i].z;
-      temperature->remove_bias_thr(i,&v[i].x,buf);
+      if (psllod_flag) temperature->remove_bias_thr(i,&v[i].x,buf);
       v[i].x = v[i].x*factor_eta - dthalf*vdelu0;
       v[i].y = v[i].y*factor_eta - dthalf*vdelu1;
       v[i].z = v[i].z*factor_eta - dthalf*vdelu2;

src/OPENMP/fix_nvt_sllod_omp.h

@@ -31,6 +31,7 @@ class FixNVTSllodOMP : public FixNHOMP {
 
  private:
   int nondeformbias;
+  int psllod_flag;     // 0 for SLLOD, 1 for p-SLLOD
 
   void nh_v_temp() override;
 };

src/fix_nh.cpp

@@ -359,6 +359,11 @@ FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) :
     } else if (strcmp(arg[iarg],"ext") == 0) {
       iarg += 2;
 
+    // keyword psllod is parsed in fix/nvt/sllod
+
+    } else if (strcmp(arg[iarg],"psllod") == 0) {
+      iarg += 2;
+
     } else error->all(FLERR,"Unknown fix nvt/npt/nph keyword: {}", arg[iarg]);
   }
 

src/fix_nvt_sllod.cpp

@@ -44,8 +44,20 @@ FixNVTSllod::FixNVTSllod(LAMMPS *lmp, int narg, char **arg) :
 
   // default values
 
+  psllod_flag = 0;
   if (mtchain_default_flag) mtchain = 1;
 
+  // select SLLOD/p-SLLOD/g-SLLOD variant
+
+  int iarg = 3;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"psllod") == 0) {
+      if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix nvt/sllod psllod", error);
+      psllod_flag = utils::logical(FLERR,arg[iarg+1],false,lmp);
+      iarg += 2;
+    } else iarg++;
+  }
 
   // create a new compute temp style
   // id = fix-ID + temp
@@ -63,22 +75,21 @@ void FixNVTSllod::init()
   FixNH::init();
 
   if (!temperature->tempbias)
-    error->all(FLERR,"Temperature for fix nvt/sllod does not have a bias");
+    error->all(FLERR,"Temperature for fix {} does not have a bias", style);
 
   nondeformbias = 0;
   if (strcmp(temperature->style,"temp/deform") != 0) nondeformbias = 1;
 
   // check fix deform remap settings
 
-  int i;
+  auto deform = modify->get_fix_by_style("^deform");
-  for (i = 0; i < modify->nfix; i++)
+  if (deform.size() < 1) error->all(FLERR,"Using fix {} with no fix deform defined", style);
-    if (strncmp(modify->fix[i]->style,"deform",6) == 0) {
+
-      if ((dynamic_cast<FixDeform *>(modify->fix[i]))->remapflag != Domain::V_REMAP)
+  for (auto ifix : deform) {
-        error->all(FLERR,"Using fix nvt/sllod with inconsistent fix deform remap option");
+    auto f = dynamic_cast<FixDeform *>(ifix);
-      break;
+    if (f && (f->remapflag != Domain::V_REMAP))
-    }
+      error->all(FLERR,"Using fix {} with inconsistent fix deform remap option", style);
-  if (i == modify->nfix)
+  }
-    error->all(FLERR,"Using fix nvt/sllod with no fix deform defined");
 }
 
 /* ----------------------------------------------------------------------
@@ -106,10 +117,11 @@ void FixNVTSllod::nh_v_temp()
 
   for (int i = 0; i < nlocal; i++) {
     if (mask[i] & groupbit) {
+      if (!psllod_flag) temperature->remove_bias(i,v[i]);
       vdelu[0] = h_two[0]*v[i][0] + h_two[5]*v[i][1] + h_two[4]*v[i][2];
       vdelu[1] = h_two[1]*v[i][1] + h_two[3]*v[i][2];
       vdelu[2] = h_two[2]*v[i][2];
-      temperature->remove_bias(i,v[i]);
+      if (psllod_flag) temperature->remove_bias(i,v[i]);
       v[i][0] = v[i][0]*factor_eta - dthalf*vdelu[0];
       v[i][1] = v[i][1]*factor_eta - dthalf*vdelu[1];
       v[i][2] = v[i][2]*factor_eta - dthalf*vdelu[2];

src/fix_nvt_sllod.h

@@ -32,6 +32,7 @@ class FixNVTSllod : public FixNH {
 
  private:
   int nondeformbias;
+  int psllod_flag;     // 0 for SLLOD, 1 for p-SLLOD
 
   void nh_v_temp() override;
 };

unittest/force-styles/tests/fix-timestep-nvt-psllod.yaml

@@ -0,0 +1,82 @@
+---
+lammps_version: 22 Dec 2022
+date_generated: Tue Jan  3 17:19:09 2023
+epsilon: 5e-12
+skip_tests:
+prerequisites: ! |
+  atom full
+  fix nvt/sllod
+pre_commands: ! ""
+post_commands: ! |
+  change_box all triclinic
+  variable gamma index 0.05
+  variable vx atom y*${gamma}
+  if $(step)==0 then "velocity all set v_vx NULL NULL sum yes"
+  fix all all deform 2 xy erate ${gamma} remap v
+  fix test all nvt/sllod temp 50.0 ${t_target} 1.0 psllod yes
+input_file: in.fourmol
+natoms: 29
+global_scalar: 2854867.1697558644
+global_vector: ! |-
+  12 1167.0668964303359 582.8541838571305 19481.306509850383 2835385.863246014 0 0 0 0 0 0 0 0
+run_pos: ! |2
+    1  7.3766923105571358e-01  2.4745425115631381e+00 -1.7150321559212145e-01
+    2  1.3664612528682669e+00  2.9513079914030929e+00 -8.5728427260886253e-01
+    3  1.9844034438342817e-01  1.2423847783028510e+00 -6.2275616693380798e-01
+    4 -6.5919858368017614e-01  1.4910537181348680e+00 -1.2483715460591922e+00
+    5 -3.3421547180145339e-02  9.3441475316631362e-01  4.0189034091250875e-01
+    6  1.0791403982726226e+00  2.3587900895167735e-01 -1.2773831741505397e+00
+    7  1.1086458909800943e+00 -1.0907385860336127e-02 -2.4728565776108500e+00
+    8  1.8864879896626023e+00 -4.8393570832815591e-01 -6.6998265853493333e-01
+    9  2.1216439487920260e+00 -2.5289625793938608e-01  2.7015879660912390e-01
+   10  2.6502404066034719e+00 -1.4343374038132586e+00 -9.6995371157599219e-01
+   11  2.3630358464762002e+00 -1.9872831599030887e+00 -1.8838953155606841e+00
+   12  3.7234311952442782e+00 -4.8922907163933527e-01 -1.6189953186835542e+00
+   13  4.7239138043397295e+00 -9.0130311982189415e-01 -1.6384795118705959e+00
+   14  3.3322391969848182e+00 -4.0806672377347741e-01 -2.6554845853739399e+00
+   15  3.8013299252875621e+00  5.6299116229155821e-01 -1.2434058969031581e+00
+   16  3.1848210365073433e+00 -2.3981796396728949e+00  2.3939136626907498e-02
+   17  2.7876270656144042e+00 -2.1004939222590853e+00  1.1560660679674388e+00
+   18  3.2095148146891033e+00  3.0177491615332581e+00 -3.5162947400390170e+00
+   19  2.5678305107909827e+00  2.6318700691921251e+00 -4.2467067539085264e+00
+   20  3.9013627348415882e+00  3.6829761983330425e+00 -3.9378172327571779e+00
+   21  5.2682659583371061e+00 -4.0751064692764150e+00 -3.6215029530747600e+00
+   22  4.7170760198580988e+00 -4.2055811767551088e+00 -4.4654878517895158e+00
+   23  6.1493858803934662e+00 -3.5762372876530470e+00 -3.8819088926442529e+00
+   24  3.1529793979324401e+00  3.1517692477958761e+00  3.1555262612569699e+00
+   25  2.3961773755311784e+00  3.2756195279484714e+00  2.5088714725751604e+00
+   26  3.7513477858704301e+00  4.0121604598214713e+00  3.2133597352292220e+00
+   27 -1.6273511722770202e+00 -4.3557333906459954e+00  2.1103388391541440e+00
+   28 -2.3744938524038011e+00 -4.0202427421428917e+00  1.5821223930194463e+00
+   29 -9.0154250196648267e-01 -3.5993683848967324e+00  2.2680257402204091e+00
+run_vel: ! |2
+    1  5.0836177409811223e-01  3.1370612048802763e-35  5.6842516004313395e-36
+    2  5.3221189005127478e-01 -1.1568811919842230e-35 -3.5425288652190980e-36
+    3  4.4676818218962833e-01 -2.0462738591065681e-35 -1.9816457116709468e-36
+    4  4.5919768180460213e-01 -6.1493933436381889e-36  6.5001492988802745e-36
+    5  4.3137141867770246e-01 -2.6763870980544101e-35 -9.6990952860385141e-36
+    6  3.9637638800124680e-01  2.2068966032648931e-34  1.0145991376231060e-34
+    7  3.8410804818903493e-01 -3.6686393855848731e-35 -1.9564741985042386e-34
+    8  3.6044445720876200e-01  7.4601232742643210e-36  1.8489721871948319e-34
+    9  3.7199364832455928e-01  1.7545603923683604e-35  3.1265484809293659e-35
+   10  3.1297495052774027e-01 -1.7593806209949413e-34 -7.7046944954625270e-35
+   11  2.8527982378990069e-01 -2.4721798297694001e-36  4.9188732460764729e-36
+   12  3.6018219953618469e-01 -1.2804205152436534e-36 -9.3769006282735508e-37
+   13  3.3957996710649124e-01 -6.6106246330704225e-36 -1.7427644528182118e-36
+   14  3.6424039896494054e-01 -1.5778843791068475e-36 -1.2101058571618634e-36
+   15  4.1279504160528241e-01 -7.9174244692607290e-36  6.7664714295642654e-36
+   16  2.6474997690126106e-01 -5.6779741842920302e-35 -1.9293097988877318e-34
+   17  2.7960772444465343e-01  3.9449292471834482e-35  1.4916612742698992e-34
+   18  5.3553032772200015e-01  1.5604831891093189e-36 -6.1723357450219555e-36
+   19  5.1623422200556024e-01  9.3604554340506648e-36  1.5850254034795202e-35
+   20  5.6879641471194298e-01 -1.5609359591659303e-35  8.8070498729597009e-36
+   21  1.8088744775485388e-01  1.9173821972353103e-36  6.7827061363565996e-38
+   22  1.7436495263537166e-01 -2.5798299640495988e-36 -5.7654851673571996e-36
+   23  2.0583283295383026e-01 -5.0667978738978151e-36  5.5016847005387114e-36
+   24  5.4223134481326019e-01  1.5141491823979144e-36  3.5248099786960212e-36
+   25  5.4842320602394556e-01  3.8812169659023934e-36 -1.2261583985187966e-35
+   26  5.8525406175525085e-01 -9.9308682506483793e-36 -1.8201170311975618e-36
+   27  1.6685679719496241e-01 -6.0455970911467563e-37  6.6604942772249726e-36
+   28  1.8362680298442527e-01  1.5809202035279489e-35 -2.2926494732590911e-35
+   29  2.0467857431672881e-01 -1.3393001700153206e-35 -3.6636010064204463e-36
+...

unittest/force-styles/tests/fix-timestep-nvt-sllod.yaml

@@ -0,0 +1,82 @@
+---
+lammps_version: 22 Dec 2022
+date_generated: Tue Jan  3 17:19:09 2023
+epsilon: 5e-12
+skip_tests:
+prerequisites: ! |
+  atom full
+  fix nvt/sllod
+pre_commands: ! ""
+post_commands: ! |
+  change_box all triclinic
+  variable gamma index 0.05
+  variable vx atom y*${gamma}
+  if $(step)==0 then "velocity all set v_vx NULL NULL sum yes"
+  fix all all deform 2 xy erate ${gamma} remap v
+  fix test all nvt/sllod temp 50.0 ${t_target} 1.0 psllod no
+input_file: in.fourmol
+natoms: 29
+global_scalar: 2854867.1697558644
+global_vector: ! |-
+  12 1167.0668964303359 582.8541838571305 19481.306509850383 2835385.863246014 0 0 0 0 0 0 0 0
+run_pos: ! |2
+    1  7.3766923105571358e-01  2.4745425115631381e+00 -1.7150321559212145e-01
+    2  1.3664612528682669e+00  2.9513079914030929e+00 -8.5728427260886253e-01
+    3  1.9844034438342817e-01  1.2423847783028510e+00 -6.2275616693380798e-01
+    4 -6.5919858368017614e-01  1.4910537181348680e+00 -1.2483715460591922e+00
+    5 -3.3421547180145339e-02  9.3441475316631362e-01  4.0189034091250875e-01
+    6  1.0791403982726226e+00  2.3587900895167735e-01 -1.2773831741505397e+00
+    7  1.1086458909800943e+00 -1.0907385860336127e-02 -2.4728565776108500e+00
+    8  1.8864879896626023e+00 -4.8393570832815591e-01 -6.6998265853493333e-01
+    9  2.1216439487920260e+00 -2.5289625793938608e-01  2.7015879660912390e-01
+   10  2.6502404066034719e+00 -1.4343374038132586e+00 -9.6995371157599219e-01
+   11  2.3630358464762002e+00 -1.9872831599030887e+00 -1.8838953155606841e+00
+   12  3.7234311952442782e+00 -4.8922907163933527e-01 -1.6189953186835542e+00
+   13  4.7239138043397295e+00 -9.0130311982189415e-01 -1.6384795118705959e+00
+   14  3.3322391969848182e+00 -4.0806672377347741e-01 -2.6554845853739399e+00
+   15  3.8013299252875621e+00  5.6299116229155821e-01 -1.2434058969031581e+00
+   16  3.1848210365073433e+00 -2.3981796396728949e+00  2.3939136626907498e-02
+   17  2.7876270656144042e+00 -2.1004939222590853e+00  1.1560660679674388e+00
+   18  3.2095148146891033e+00  3.0177491615332581e+00 -3.5162947400390170e+00
+   19  2.5678305107909827e+00  2.6318700691921251e+00 -4.2467067539085264e+00
+   20  3.9013627348415882e+00  3.6829761983330425e+00 -3.9378172327571779e+00
+   21  5.2682659583371061e+00 -4.0751064692764150e+00 -3.6215029530747600e+00
+   22  4.7170760198580988e+00 -4.2055811767551088e+00 -4.4654878517895158e+00
+   23  6.1493858803934662e+00 -3.5762372876530470e+00 -3.8819088926442529e+00
+   24  3.1529793979324401e+00  3.1517692477958761e+00  3.1555262612569699e+00
+   25  2.3961773755311784e+00  3.2756195279484714e+00  2.5088714725751604e+00
+   26  3.7513477858704301e+00  4.0121604598214713e+00  3.2133597352292220e+00
+   27 -1.6273511722770202e+00 -4.3557333906459954e+00  2.1103388391541440e+00
+   28 -2.3744938524038011e+00 -4.0202427421428917e+00  1.5821223930194463e+00
+   29 -9.0154250196648267e-01 -3.5993683848967324e+00  2.2680257402204091e+00
+run_vel: ! |2
+    1  5.0836177409811223e-01  3.1370612048802763e-35  5.6842516004313395e-36
+    2  5.3221189005127478e-01 -1.1568811919842230e-35 -3.5425288652190980e-36
+    3  4.4676818218962833e-01 -2.0462738591065681e-35 -1.9816457116709468e-36
+    4  4.5919768180460213e-01 -6.1493933436381889e-36  6.5001492988802745e-36
+    5  4.3137141867770246e-01 -2.6763870980544101e-35 -9.6990952860385141e-36
+    6  3.9637638800124680e-01  2.2068966032648931e-34  1.0145991376231060e-34
+    7  3.8410804818903493e-01 -3.6686393855848731e-35 -1.9564741985042386e-34
+    8  3.6044445720876200e-01  7.4601232742643210e-36  1.8489721871948319e-34
+    9  3.7199364832455928e-01  1.7545603923683604e-35  3.1265484809293659e-35
+   10  3.1297495052774027e-01 -1.7593806209949413e-34 -7.7046944954625270e-35
+   11  2.8527982378990069e-01 -2.4721798297694001e-36  4.9188732460764729e-36
+   12  3.6018219953618469e-01 -1.2804205152436534e-36 -9.3769006282735508e-37
+   13  3.3957996710649124e-01 -6.6106246330704225e-36 -1.7427644528182118e-36
+   14  3.6424039896494054e-01 -1.5778843791068475e-36 -1.2101058571618634e-36
+   15  4.1279504160528241e-01 -7.9174244692607290e-36  6.7664714295642654e-36
+   16  2.6474997690126106e-01 -5.6779741842920302e-35 -1.9293097988877318e-34
+   17  2.7960772444465343e-01  3.9449292471834482e-35  1.4916612742698992e-34
+   18  5.3553032772200015e-01  1.5604831891093189e-36 -6.1723357450219555e-36
+   19  5.1623422200556024e-01  9.3604554340506648e-36  1.5850254034795202e-35
+   20  5.6879641471194298e-01 -1.5609359591659303e-35  8.8070498729597009e-36
+   21  1.8088744775485388e-01  1.9173821972353103e-36  6.7827061363565996e-38
+   22  1.7436495263537166e-01 -2.5798299640495988e-36 -5.7654851673571996e-36
+   23  2.0583283295383026e-01 -5.0667978738978151e-36  5.5016847005387114e-36
+   24  5.4223134481326019e-01  1.5141491823979144e-36  3.5248099786960212e-36
+   25  5.4842320602394556e-01  3.8812169659023934e-36 -1.2261583985187966e-35
+   26  5.8525406175525085e-01 -9.9308682506483793e-36 -1.8201170311975618e-36
+   27  1.6685679719496241e-01 -6.0455970911467563e-37  6.6604942772249726e-36
+   28  1.8362680298442527e-01  1.5809202035279489e-35 -2.2926494732590911e-35
+   29  2.0467857431672881e-01 -1.3393001700153206e-35 -3.6636010064204463e-36
+...