diff --git a/doc/src/Howto_bpm.rst b/doc/src/Howto_bpm.rst
index 0ca4e85fbb..5aa277275d 100644
--- a/doc/src/Howto_bpm.rst
+++ b/doc/src/Howto_bpm.rst
@@ -15,7 +15,8 @@ orientation for rotational models. This produces a stress-free initial
 state. Furthermore, bonds are allowed to break under large strains,
 producing fracture. The examples/bpm directory has sample input scripts
 for simulations of the fragmentation of an impacted plate and the
-pouring of extended, elastic bodies.
+pouring of extended, elastic bodies. See :ref:`(Clemmer) <howto-Clemmer>`
+for more general information on the approach and the LAMMPS implementation.
 
 ----------
 
@@ -150,3 +151,9 @@ the following are currently compatible with BPM bond styles:
    interactions, one will need to switch between different *special_bonds*
    settings in the input script. An example is found in
    ``examples/bpm/impact``.
+
+----------
+
+.. _howto-Clemmer:
+
+**(Clemmer)** Clemmer, Monti, Lechman, Soft Matter, 20, 1702 (2024).
diff --git a/doc/src/fix_deform_pressure.rst b/doc/src/fix_deform_pressure.rst
index 8e848b3969..1490390988 100644
--- a/doc/src/fix_deform_pressure.rst
+++ b/doc/src/fix_deform_pressure.rst
@@ -29,10 +29,12 @@ Syntax
            NOTE: All other styles are documented by the :doc:`fix deform <fix_deform>` command
 
        *xy*, *xz*, *yz* args = style value
-         style = *final* or *delta* or *vel* or *erate* or *trate* or *wiggle* or *variable* or *pressure*
+         style = *final* or *delta* or *vel* or *erate* or *trate* or *wiggle* or *variable* or *pressure* or *erate/rescale*
            *pressure* values = target gain
              target = target pressure (pressure units)
              gain = proportional gain constant (1/(time * pressure) or 1/time units)
+           *erate/rescale* value = R
+             R = engineering shear strain rate (1/time units)
            NOTE: All other styles are documented by the :doc:`fix deform <fix_deform>` command
 
        *box* = style value
@@ -159,6 +161,21 @@ details of a simulation and testing different values is
 recommended. One can also apply a maximum limit to the magnitude of
 the applied strain using the :ref:`max/rate <deform_max_rate>` option.
 
+The *erate/rescale* style operates similarly to the *erate* style with
+a specified strain rate in units of 1/time. The difference is that
+the change in the tilt factor will depend on the current length of
+the box perpendicular to the shear direction, L, instead of the
+original length, L0. The tilt factor T as a function of time will
+change as
+
+.. parsed-literal::
+
+   T(t) = T(t-1) + L\*erate\* \Delta t
+
+where T(t-1) is the tilt factor on the previous timestep and :math:`\Delta t`
+is the timestep size. This option may be useful in scenarios where
+L changes in time.
+
 ----------
 
 The *box* parameter provides an additional control over the *x*, *y*,
diff --git a/doc/src/fix_nonaffine_displacement.rst b/doc/src/fix_nonaffine_displacement.rst
index 0a271ebc32..fd9830cc48 100644
--- a/doc/src/fix_nonaffine_displacement.rst
+++ b/doc/src/fix_nonaffine_displacement.rst
@@ -8,7 +8,7 @@ Syntax
 
 .. parsed-literal::
 
-   fix ID group nonaffine/displacement style args reference/style nstep
+   fix ID group nonaffine/displacement style args reference/style nstep keyword values
 
 * ID, group are documented in :doc:`fix <fix>` command
 * nonaffine/displacement = style name of this fix command
@@ -32,6 +32,13 @@ Syntax
        *update* = update the reference frame every *nstep* timesteps
        *offset* = update the reference frame *nstep* timesteps before calculating the nonaffine displacement
 
+* zero or more keyword/value pairs may be appended
+
+  .. parsed-literal::
+
+       *z/min* values = zmin
+         zmin = minimum coordination number to calculate D2min
+
 Examples
 """"""""
 
@@ -76,6 +83,12 @@ is the identity tensor. This calculation is only performed on timesteps that
 are a multiple of *nevery* (including timestep zero). Data accessed before
 this occurs will simply be zeroed.
 
+For particles with low coordination numbers, calculations of :math:`D^2_\mathrm{min}`
+may not be accurate. An optional minimum coordination number can be defined using
+the *z/min* keyword. If any particle has fewer than the specified number of particles
+in the cutoff distance or in contact, the above calculations will be skipped and the
+corresponding peratom array entries will be zero.
+
 The *integrated* style simply integrates the velocity of particles
 every timestep to calculate a displacement. This style only works if
 used in conjunction with another fix that deforms the box and displaces
diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp
index f1482d4203..9c2c680cc5 100644
--- a/src/BPM/bond_bpm.cpp
+++ b/src/BPM/bond_bpm.cpp
@@ -14,6 +14,7 @@
 #include "bond_bpm.h"
 
 #include "atom.h"
+#include "citeme.h"
 #include "comm.h"
 #include "domain.h"
 #include "error.h"
@@ -30,6 +31,19 @@
 
 using namespace LAMMPS_NS;
 
+static const char cite_bpm[] =
+  "BPM bond style: doi:10.1039/D3SM01373A\n\n"
+  "@Article{Clemmer2024,\n"
+  " author =  {Clemmer, Joel T. and Monti, Joseph M. and Lechman, Jeremy B.},\n"
+  " title =   {A soft departure from jamming: the compaction of deformable\n"
+  "            granular matter under high pressures},\n"
+  " journal = {Soft Matter},\n"
+  " year =    2024,\n"
+  " volume =  20,\n"
+  " number =  8,\n"
+  " pages =   {1702--1718}\n"
+  "}\n\n";
+
 /* ---------------------------------------------------------------------- */
 
 BondBPM::BondBPM(LAMMPS *_lmp) :
@@ -55,6 +69,8 @@ BondBPM::BondBPM(LAMMPS *_lmp) :
 
   id_fix_dummy2 = utils::strdup("BPM_DUMMY2");
   modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy2));
+
+  if (lmp->citeme) lmp->citeme->add(cite_bpm);
 }
 
 /* ---------------------------------------------------------------------- */
diff --git a/src/EXTRA-FIX/fix_deform_pressure.cpp b/src/EXTRA-FIX/fix_deform_pressure.cpp
index 51ea75cfed..e2bcdd7f4e 100644
--- a/src/EXTRA-FIX/fix_deform_pressure.cpp
+++ b/src/EXTRA-FIX/fix_deform_pressure.cpp
@@ -110,6 +110,12 @@ FixDeformPressure::FixDeformPressure(LAMMPS *lmp, int narg, char **arg) :
         }
         set_extra[index].pgain = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
         i += 4;
+      } else if (strcmp(arg[iarg + 1], "erate/rescale") == 0) {
+        if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure erate/rescale", error);
+        set[index].style = ERATERS;
+        set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
+        iarg += 3;
+        i += 3;
       } else error->all(FLERR, "Illegal fix deform/pressure command: {}", arg[iarg + 1]);
 
     } else if (strcmp(arg[iarg], "box") == 0) {
@@ -424,16 +430,31 @@ void FixDeformPressure::init()
     if (!pressure)
       error->all(FLERR, "Pressure ID {} for fix deform/pressure does not exist", id_press);
   }
+
+  // if yz [3] changes and will cause box flip, then xy [5] cannot be changing
+  // this is b/c the flips would induce continuous changes in xz
+  //   in order to keep the edge vectors of the flipped shape matrix
+  //   an integer combination of the edge vectors of the unflipped shape matrix
+  // error if style PRESSURE/ERATEER for yz, can't calculate if box flip occurs
+
+  if (set[3].style && set[5].style) {
+    int flag = 0;
+    double lo,hi;
+    if (flipflag && set[3].style == PRESSURE)
+      error->all(FLERR, "Fix {} cannot use yz pressure with xy", style);
+    if (flipflag && set[3].style == ERATERS)
+      error->all(FLERR, "Fix {} cannot use yz erate/rescale with xy", style);
+  }
 }
 
 /* ----------------------------------------------------------------------
-   compute T,P if needed before integrator starts
+   compute T,P before integrator starts
 ------------------------------------------------------------------------- */
 
 void FixDeformPressure::setup(int /*vflag*/)
 {
-  // trigger virial computation on next timestep
-  if (pressure_flag) pressure->addstep(update->ntimestep+1);
+  // trigger virial computation, if needed, on next timestep
+  if (pressure_flag) pressure->addstep(update->ntimestep + 1);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -446,7 +467,20 @@ void FixDeformPressure::end_of_step()
 
   // set new box size for strain-based dims
 
-  if (strain_flag) FixDeform::apply_strain();
+  if (strain_flag) {
+    FixDeform::apply_strain();
+
+    for (int i = 3; i < 6; i++) {
+      if (set[i].style == ERATERS) {
+        double L = domain->zprd;
+        if (i == 5) L = domain->yprd;
+
+        h_rate[i] = set[i].rate * L;
+        set_extra[i].cumulative_shift += update->dt * h_rate[i];
+        set[i].tilt_target = set[i].tilt_start + set_extra[i].cumulative_shift;
+      }
+    }
+  }
 
   // set new box size for pressure-based dims
 
@@ -479,12 +513,33 @@ void FixDeformPressure::end_of_step()
     for (int i = 0; i < 3; i++) {
       set_extra[i].prior_pressure = pressure->vector[i];
       set_extra[i].prior_rate = ((set[i].hi_target - set[i].lo_target) /
-                           (domain->boxhi[i] - domain->boxlo[i]) - 1.0)  / update->dt;
+                           domain->prd[i] - 1.0)  / update->dt;
     }
   }
 
   if (varflag) modify->addstep_compute(update->ntimestep + nevery);
 
+  // If tilting while evolving linear dimension, sum remapping effects
+  // otherwise, update_domain() will inaccurately use the current
+  // linear dimension to apply prior remappings
+
+  for (int i = 3; i < 6; i++) {
+    int idenom = 0;
+    if (i == 3) idenom = 1;
+    if (set[i].style && (set_box.style || set[idenom].style) && domain->periodicity[idenom]) {
+      // Add prior remappings. If the box remaps this timestep, don't
+      // add it yet so update_domain() will first detect the remapping
+      set[i].tilt_target += set_extra[i].cumulative_remap;
+
+      // Update remapping for next timestep
+      double prd = set[idenom].hi_target - set[idenom].lo_target;
+      double prdinv = 1.0 / prd;
+      if (set[i].tilt_target * prdinv < -0.5)
+        set_extra[i].cumulative_remap += prd;
+      if (set[i].tilt_target * prdinv > 0.5)
+        set_extra[i].cumulative_remap -= prd;
+    }
+  }
 
   FixDeform::update_domain();
 
@@ -502,7 +557,7 @@ void FixDeformPressure::apply_pressure()
 {
   // If variable pressure, calculate current target
   for (int i = 0; i < 6; i++)
-    if (set[i].style == PRESSURE)
+    if (set[i].style == PRESSURE || set[i].style == PMEAN)
       if (set_extra[i].pvar_flag)
         set_extra[i].ptarget = input->variable->compute_equal(set_extra[i].pvar);
 
@@ -556,26 +611,24 @@ void FixDeformPressure::apply_pressure()
 
     h_ratelo[i] = -0.5 * h_rate[i];
 
-    double offset = 0.5 * (domain->boxhi[i] - domain->boxlo[i]) * (1.0 + update->dt * h_rate[i]);
-    set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - offset;
-    set[i].hi_target = 0.5 * (set[i].lo_start + set[i].hi_start) + offset;
+    double shift = domain->prd[i] * update->dt * h_rate[i];
+    set_extra[i].cumulative_shift += shift;
+    set[i].lo_target = set[i].lo_start - 0.5 * set_extra[i].cumulative_shift;
+    set[i].hi_target = set[i].hi_start + 0.5 * set_extra[i].cumulative_shift;
   }
 
   for (int i = 3; i < 6; i++) {
     if (set[i].style != PRESSURE) continue;
 
-    double L, tilt, pcurrent;
+    double L, pcurrent;
     if (i == 3) {
       L = domain->zprd;
-      tilt = domain->yz;
       pcurrent = tensor[5];
     } else if (i == 4) {
       L = domain->zprd;
-      tilt = domain->xz + update->dt;
       pcurrent = tensor[4];
     } else {
       L = domain->yprd;
-      tilt = domain->xy;
       pcurrent = tensor[3];
     }
 
@@ -592,7 +645,8 @@ void FixDeformPressure::apply_pressure()
       if (fabs(h_rate[i]) > max_h_rate)
         h_rate[i] = max_h_rate * h_rate[i] / fabs(h_rate[i]);
 
-    set[i].tilt_target = tilt + update->dt * h_rate[i];
+    set_extra[i].cumulative_shift += update->dt * h_rate[i];
+    set[i].tilt_target = set[i].tilt_start + set_extra[i].cumulative_shift;
   }
 }
 
@@ -629,9 +683,9 @@ void FixDeformPressure::apply_volume()
         double dt = update->dt;
         double e1i = set_extra[i].prior_rate;
         double e2i = set_extra[fixed].prior_rate;
-        double L1i = domain->boxhi[i] - domain->boxlo[i];
-        double L2i = domain->boxhi[fixed] - domain->boxlo[fixed];
-        double L3i = domain->boxhi[dynamic1] - domain->boxlo[dynamic1];
+        double L1i = domain->prd[i];
+        double L2i = domain->prd[fixed];
+        double L3i = domain->prd[dynamic1];
         double L3 = (set[dynamic1].hi_target - set[dynamic1].lo_target);
         double Vi = L1i * L2i * L3i;
         double V = L3 * L1i * L2i;
@@ -680,7 +734,7 @@ void FixDeformPressure::apply_volume()
       }
     }
 
-    h_rate[i] = (2.0 * shift / (domain->boxhi[i] - domain->boxlo[i]) - 1.0) / update->dt;
+    h_rate[i] = (2.0 * shift / domain->prd[i] - 1.0) / update->dt;
     h_ratelo[i] = -0.5 * h_rate[i];
 
     set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - shift;
@@ -742,7 +796,7 @@ void FixDeformPressure::apply_box()
       set[i].hi_target = 0.5 * (set[i].lo_start + set[i].hi_start) + shift;
 
       // Recalculate h_rate
-      h_rate[i] = (set[i].hi_target - set[i].lo_target) / (domain->boxhi[i] - domain->boxlo[i]) - 1.0;
+      h_rate[i] = (set[i].hi_target - set[i].lo_target) / domain->prd[i] - 1.0;
       h_rate[i] /= update->dt;
       h_ratelo[i] = -0.5 * h_rate[i];
     }
@@ -767,14 +821,21 @@ void FixDeformPressure::apply_box()
       if (fabs(v_rate) > max_h_rate)
         v_rate = max_h_rate * v_rate / fabs(v_rate);
 
-    scale = (1.0 + update->dt * v_rate);
     for (i = 0; i < 3; i++) {
-      shift = 0.5 * (set[i].hi_target - set[i].lo_target) * scale;
-      set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - shift;
-      set[i].hi_target = 0.5 * (set[i].lo_start + set[i].hi_start) + shift;
+      shift = (set[i].hi_target - set[i].lo_target) * update->dt * v_rate;
+      set_extra[6].cumulative_vshift[i] += shift;
+
+      if (set[i].style == NONE) {
+        // Overwrite default targets of current length
+        set[i].lo_target = set[i].lo_start;
+        set[i].hi_target = set[i].hi_start;
+      }
+
+      set[i].lo_target -= 0.5 * set_extra[6].cumulative_vshift[i];
+      set[i].hi_target += 0.5 * set_extra[6].cumulative_vshift[i];
 
       // Recalculate h_rate
-      h_rate[i] = (set[i].hi_target - set[i].lo_target) / (domain->boxhi[i] - domain->boxlo[i]) - 1.0;
+      h_rate[i] = (set[i].hi_target - set[i].lo_target) / domain->prd[i] - 1.0;
       h_rate[i] /= update->dt;
       h_ratelo[i] = -0.5 * h_rate[i];
     }
@@ -788,7 +849,7 @@ void FixDeformPressure::apply_box()
 void FixDeformPressure::write_restart(FILE *fp)
 {
   if (comm->me == 0) {
-    int size = 9 * sizeof(double) + 7 * sizeof(Set) + 7 * sizeof(SetExtra);
+    int size = 7 * sizeof(Set) + 7 * sizeof(SetExtra);
     fwrite(&size, sizeof(int), 1, fp);
     fwrite(set, sizeof(Set), 6, fp);
     fwrite(&set_box, sizeof(Set), 1, fp);
@@ -803,22 +864,16 @@ void FixDeformPressure::write_restart(FILE *fp)
 void FixDeformPressure::restart(char *buf)
 {
   int n = 0;
-  auto list = (double *) buf;
-  for (int i = 0; i < 6; i++)
-    h_rate[i] = list[n++];
-  for (int i = 0; i < 3; i++)
-    h_ratelo[i] = list[n++];
-
-  n = n * sizeof(double);
   int samestyle = 1;
-  Set *set_restart = (Set *) &buf[n];
+  Set *set_restart = (Set *) buf;
   for (int i = 0; i < 6; ++i) {
     // restore data from initial state
     set[i].lo_initial = set_restart[i].lo_initial;
     set[i].hi_initial = set_restart[i].hi_initial;
     set[i].vol_initial = set_restart[i].vol_initial;
     set[i].tilt_initial = set_restart[i].tilt_initial;
-    // check if style settings are consistent (should do the whole set?)
+
+    // check if style settings are consistent
     if (set[i].style != set_restart[i].style)
       samestyle = 0;
     if (set[i].substyle != set_restart[i].substyle)
diff --git a/src/EXTRA-FIX/fix_deform_pressure.h b/src/EXTRA-FIX/fix_deform_pressure.h
index 10af1e5ba3..7ce69b9bc5 100644
--- a/src/EXTRA-FIX/fix_deform_pressure.h
+++ b/src/EXTRA-FIX/fix_deform_pressure.h
@@ -51,6 +51,9 @@ class FixDeformPressure : public FixDeform {
   struct SetExtra {
     double ptarget, pgain;
     double prior_pressure, prior_rate;
+    double cumulative_shift;
+    double cumulative_vshift[3];
+    double cumulative_remap;
     int saved;
     char *pstr;
     int pvar, pvar_flag;
diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp
index eaf45f4e59..b651d5dc5e 100644
--- a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp
+++ b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp
@@ -46,6 +46,8 @@ enum { TYPE, RADIUS, CUSTOM };
 enum { INTEGRATED, D2MIN };
 enum { FIXED, OFFSET, UPDATE };
 
+static constexpr double EPSILON = 1.0e-15;
+
 static const char cite_nonaffine_d2min[] =
   "@article{PhysRevE.57.7192,\n"
   " title = {Dynamics of viscoplastic deformation in amorphous solids},\n"
@@ -66,7 +68,7 @@ static const char cite_nonaffine_d2min[] =
 
 FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg), id_fix(nullptr), fix(nullptr), D2min(nullptr), X(nullptr), Y(nullptr),
-  F(nullptr), norm(nullptr)
+  F(nullptr), norm(nullptr), singular(nullptr)
 {
   if (narg < 4) utils::missing_cmd_args(FLERR,"fix nonaffine/displacement", error);
 
@@ -74,6 +76,8 @@ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char *
   if (nevery <= 0) error->all(FLERR,"Illegal nevery value {} in fix nonaffine/displacement", nevery);
 
   reference_timestep = update_timestep = offset_timestep = -1;
+  z_min = 0;
+
   int iarg = 4;
   if (strcmp(arg[iarg], "integrated") == 0) {
     nad_style = INTEGRATED;
@@ -117,6 +121,16 @@ FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char *
     if ((offset_timestep <= 0) || (offset_timestep > nevery))
       error->all(FLERR, "Illegal offset timestep {} in fix nonaffine/displacement", arg[iarg + 1]);
   } else error->all(FLERR,"Illegal reference style {} in fix nonaffine/displacement", arg[iarg]);
+  iarg += 2;
+
+  while (iarg < narg) {
+    if (strcmp(arg[iarg], "z/min") == 0) {
+      if (iarg + 2 > narg) utils::missing_cmd_args(FLERR,"fix nonaffine/displacement", error);
+      z_min = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
+      if (z_min < 0) error->all(FLERR, "Minimum coordination must be positive");
+      iarg += 2;
+    } else error->all(FLERR,"Illegal keyword {} in fix nonaffine/displacement", arg[iarg]);
+  }
 
   if (nad_style == D2MIN)
     if (cut_style == RADIUS && (!atom->radius_flag))
@@ -151,6 +165,7 @@ FixNonaffineDisplacement::~FixNonaffineDisplacement()
     memory->destroy(Y);
     memory->destroy(F);
     memory->destroy(norm);
+    memory->destroy(singular);
     memory->destroy(D2min);
   }
 
@@ -395,6 +410,7 @@ void FixNonaffineDisplacement::calculate_D2Min()
       }
     }
     norm[i] = 0;
+    singular[i] = 0;
     D2min[i] = 0;
   }
 
@@ -471,14 +487,29 @@ void FixNonaffineDisplacement::calculate_D2Min()
     }
 
     if (dim == 3) {
-      invert3(Y_tmp, Y_inv);
+      denom = det3(Y_tmp);
+      if (fabs(denom) < EPSILON) {
+        singular[i] = 1;
+        for (j = 0; j < 3; j++)
+          for (k = 0; k < 3; k++)
+            Y_inv[j][k] = 0.0;
+      } else {
+        invert3(Y_tmp, Y_inv);
+      }
     } else {
       denom = Y_tmp[0][0] * Y_tmp[1][1] - Y_tmp[0][1] * Y_tmp[1][0];
-      if (denom != 0.0) denom = 1.0 / denom;
-      Y_inv[0][0] = Y_tmp[1][1] * denom;
-      Y_inv[0][1] = -Y_tmp[0][1] * denom;
-      Y_inv[1][0] = -Y_tmp[1][0] * denom;
-      Y_inv[1][1] = Y_tmp[0][0] * denom;
+      if (fabs(denom) < EPSILON) {
+        singular[i] = 1;
+        for (j = 0; j < 2; j++)
+          for (k = 0; k < 2; k++)
+            Y_inv[j][k] = 0.0;
+      } else {
+        denom = 1.0 / denom;
+        Y_inv[0][0] = Y_tmp[1][1] * denom;
+        Y_inv[0][1] = -Y_tmp[0][1] * denom;
+        Y_inv[1][0] = -Y_tmp[1][0] * denom;
+        Y_inv[1][1] = Y_tmp[0][0] * denom;
+      }
     }
 
     times3(X_tmp, Y_inv, F_tmp);
@@ -559,10 +590,16 @@ void FixNonaffineDisplacement::calculate_D2Min()
   for (i = 0; i < nlocal; i++) {
     if (!(mask[i] & groupbit)) continue;
 
-    if (norm[i] != 0)
-      D2min[i] /= norm[i];
-    else
-      D2min[i] = 0.0;
+    if (norm[i] < z_min || singular[i] == 1) {
+      if (norm[i] >= z_min)
+        error->warning(FLERR, "Singular matrix detected for atom {}, defaulting output to zero", atom->tag[i]);
+      array_atom[i][0] = 0.0;
+      array_atom[i][1] = 0.0;
+      array_atom[i][2] = 0.0;
+      continue;
+    }
+
+    D2min[i] /= norm[i];
 
     for (j = 0; j < 3; j++)
       for (k = 0; k < 3; k++)
@@ -743,10 +780,12 @@ void FixNonaffineDisplacement::grow_arrays(int nmax_new)
     memory->destroy(F);
     memory->destroy(D2min);
     memory->destroy(norm);
+    memory->destroy(singular);
     memory->create(X, nmax, 3, 3, "fix_nonaffine_displacement:X");
     memory->create(Y, nmax, 3, 3, "fix_nonaffine_displacement:Y");
     memory->create(F, nmax, 3, 3, "fix_nonaffine_displacement:F");
     memory->create(D2min, nmax, "fix_nonaffine_displacement:D2min");
     memory->create(norm, nmax, "fix_nonaffine_displacement:norm");
+    memory->create(singular, nmax, "fix_nonaffine_displacement:singular");
   }
 }
diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.h b/src/EXTRA-FIX/fix_nonaffine_displacement.h
index c7177bd3d9..b0e9c464ca 100644
--- a/src/EXTRA-FIX/fix_nonaffine_displacement.h
+++ b/src/EXTRA-FIX/fix_nonaffine_displacement.h
@@ -48,12 +48,12 @@ class FixNonaffineDisplacement : public Fix {
   int nmax, comm_flag;
   int nad_style, cut_style;
   int reference_style, offset_timestep, reference_timestep, update_timestep;
-  int reference_saved;
+  int reference_saved, z_min;
   double cutoff_custom, cutsq_custom, mycutneigh;
   double xprd0, yprd0, zprd0, xprd0_half, yprd0_half, zprd0_half, xy0, xz0, yz0;
 
   double *D2min, ***X, ***Y, ***F;
-  int *norm;
+  int *norm, *singular;
 
   class NeighList *list;    // half neighbor list
 
diff --git a/src/GRANULAR/fix_heat_flow.cpp b/src/GRANULAR/fix_heat_flow.cpp
index b7643c2c24..be8d93839f 100644
--- a/src/GRANULAR/fix_heat_flow.cpp
+++ b/src/GRANULAR/fix_heat_flow.cpp
@@ -16,6 +16,7 @@
 #include "atom.h"
 #include "comm.h"
 #include "error.h"
+#include "force.h"
 #include "memory.h"
 #include "modify.h"
 #include "update.h"
@@ -127,7 +128,7 @@ void FixHeatFlow::final_integrate()
   if (igroup == atom->firstgroup) nlocal = atom->nfirst;
 
   // add ghost contributions to heatflow if first instance of fix
-  if (first_flag) comm->reverse_comm(this);
+  if (force->newton_pair && first_flag) comm->reverse_comm(this);
 
   if (rmass) {
     for (int i = 0; i < nlocal; i++)
diff --git a/src/GRANULAR/granular_model.cpp b/src/GRANULAR/granular_model.cpp
index 14431f41b4..9764ec42e9 100644
--- a/src/GRANULAR/granular_model.cpp
+++ b/src/GRANULAR/granular_model.cpp
@@ -216,9 +216,15 @@ int GranularModel::define_classic_model(char **arg, int iarg, int narg)
   // manually parse coeffs
   normal_model->coeffs[0] = kn;
   normal_model->coeffs[1] = gamman;
-  tangential_model->coeffs[0] = kt;
-  tangential_model->coeffs[1] = gammat / gamman;
-  tangential_model->coeffs[2] = xmu;
+
+  if (tangential_model->num_coeffs == 2) {
+    tangential_model->coeffs[0] = gammat / gamman;
+    tangential_model->coeffs[1] = xmu;
+  } else {
+    tangential_model->coeffs[0] = kt;
+    tangential_model->coeffs[1] = gammat / gamman;
+    tangential_model->coeffs[2] = xmu;
+  }
 
   normal_model->coeffs_to_local();
   tangential_model->coeffs_to_local();
diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp
index bb27faeaa8..135d7176e6 100644
--- a/src/fix_deform.cpp
+++ b/src/fix_deform.cpp
@@ -63,7 +63,7 @@ irregular(nullptr), set(nullptr)
   int nskip;
   if (utils::strmatch(style, "^deform/pressure")) {
     child_parameters.insert("box");
-    child_styles.insert({{"pressure", 4}, {"pressure/mean", 4}, {"volume", 2}});
+    child_styles.insert({{"pressure", 4}, {"pressure/mean", 4}, {"erate/rescale", 3}, {"volume", 2}});
   }
 
   // set defaults
diff --git a/src/fix_deform.h b/src/fix_deform.h
index b133729444..c524c2fe6c 100644
--- a/src/fix_deform.h
+++ b/src/fix_deform.h
@@ -29,7 +29,7 @@ class FixDeform : public Fix {
   int remapflag;     // whether x,v are remapped across PBC
   int dimflag[6];    // which dims are deformed
 
-  enum { NONE, FINAL, DELTA, SCALE, VEL, ERATE, TRATE, VOLUME, WIGGLE, VARIABLE, PRESSURE, PMEAN };
+  enum { NONE, FINAL, DELTA, SCALE, VEL, ERATE, TRATE, VOLUME, WIGGLE, VARIABLE, PRESSURE, PMEAN, ERATERS };
   enum { ONE_FROM_ONE, ONE_FROM_TWO, TWO_FROM_ONE };
 
   FixDeform(class LAMMPS *, int, char **);