diff --git a/doc/src/fix_nvt_sllod.rst b/doc/src/fix_nvt_sllod.rst
index 04e60057a1..46e77088a6 100644
--- a/doc/src/fix_nvt_sllod.rst
+++ b/doc/src/fix_nvt_sllod.rst
@@ -17,6 +17,14 @@ Syntax
 
 * ID, group-ID are documented in :doc:`fix <fix>` command
 * nvt/sllod = style name of this fix command
+* zero or more keyword/value pairs may be appended
+* keyword = *psllod*
+
+  .. 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
@@ -67,16 +75,20 @@ equivalent to Newton's equations of motion for shear flow by
 the desired velocity gradient and the correct production of work by
 stresses for all forms of homogeneous flow by :ref:`(Daivis and Todd)
 <Daivis>`.
-The LAMMPS implementation corresponds to the p-SLLOD variant
-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 +171,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 +183,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 +195,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).
diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp
index d0356f10f3..ac582f8db5 100644
--- a/src/fix_nh.cpp
+++ b/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]);
   }
 
diff --git a/src/fix_nvt_sllod.cpp b/src/fix_nvt_sllod.cpp
index bb4b6c4273..6d0b74b859 100644
--- a/src/fix_nvt_sllod.cpp
+++ b/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
@@ -106,10 +118,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];
diff --git a/src/fix_nvt_sllod.h b/src/fix_nvt_sllod.h
index e13eb11311..081f6e657d 100644
--- a/src/fix_nvt_sllod.h
+++ b/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;
 };