diff --git a/doc/src/fix_deform.rst b/doc/src/fix_deform.rst
index 9146b987c8..2c76463369 100644
--- a/doc/src/fix_deform.rst
+++ b/doc/src/fix_deform.rst
@@ -64,6 +64,8 @@ Syntax
                  effectively an engineering shear strain rate
            *erate* value = R
              R = engineering shear strain rate (1/time units)
+           *erate/rescale* value = R (ONLY available in :doc:`fix deform/pressure <fix_deform_pressure>` command)
+             R = engineering shear strain rate (1/time units)
            *trate* value = R
              R = true shear strain rate (1/time units)
            *wiggle* values = A Tp
diff --git a/doc/src/fix_wall_gran.rst b/doc/src/fix_wall_gran.rst
index f6465d1159..25d659241c 100644
--- a/doc/src/fix_wall_gran.rst
+++ b/doc/src/fix_wall_gran.rst
@@ -115,6 +115,18 @@ friction and twisting friction supported by the :doc:`pair_style granular <pair_
 supported for walls. These are discussed in greater detail on the doc
 page for :doc:`pair_style granular <pair_granular>`.
 
+.. note::
+   When *fstyle* *granular* is specified, the associated *fstyle_params* are taken as
+   those for a wall/particle interaction. For example, for the *hertz/material* normal
+   contact model with :math:`E = 960` and :math:`\nu = 0.2`, the effective Young's
+   modulus for a wall/particle interaction is computed as
+   :math:`E_{eff} = \frac{960}{2(1-0.2^2)} = 500`. Any pair coefficients defined by
+   :doc:`pair_style granular <pair_granular>` are not taken into consideration. To
+   model different wall/particle interactions for particles of different material
+   types, the user may define multiple fix wall/gran commands operating on separate
+   groups (e.g. based on particle type) each with a different wall/particle effective
+   Young's modulus.
+
 Note that you can choose a different force styles and/or different
 values for the wall/particle coefficients than for particle/particle
 interactions.  E.g. if you wish to model the wall as a different
diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst
index 597869a05b..dca85d8be1 100644
--- a/doc/src/pair_granular.rst
+++ b/doc/src/pair_granular.rst
@@ -111,7 +111,7 @@ For the *hertz* model, the normal component of force is given by:
 
    \mathbf{F}_{ne, Hertz} = k_n R_{eff}^{1/2}\delta_{ij}^{3/2} \mathbf{n}
 
-Here, :math:`R_{eff} = \frac{R_i R_j}{R_i + R_j}` is the effective
+Here, :math:`R_{eff} = R = \frac{R_i R_j}{R_i + R_j}` is the effective
 radius, denoted for simplicity as *R* from here on.  For *hertz*, the
 units of the spring constant :math:`k_n` are *force*\ /\ *length*\ \^2, or
 equivalently *pressure*\ .
@@ -120,13 +120,14 @@ For the *hertz/material* model, the force is given by:
 
 .. math::
 
-   \mathbf{F}_{ne, Hertz/material} = \frac{4}{3} E_{eff} R_{eff}^{1/2}\delta_{ij}^{3/2} \mathbf{n}
+   \mathbf{F}_{ne, Hertz/material} = \frac{4}{3} E_{eff} R^{1/2}\delta_{ij}^{3/2} \mathbf{n}
 
-Here, :math:`E_{eff} = E = \left(\frac{1-\nu_i^2}{E_i} + \frac{1-\nu_j^2}{E_j}\right)^{-1}` is the effective Young's
-modulus, with :math:`\nu_i, \nu_j` the Poisson ratios of the particles of
-types *i* and *j*\ . Note that if the elastic modulus and the shear
-modulus of the two particles are the same, the *hertz/material* model
-is equivalent to the *hertz* model with :math:`k_n = 4/3 E_{eff}`
+Here, :math:`E_{eff} = E = \left(\frac{1-\nu_i^2}{E_i} + \frac{1-\nu_j^2}{E_j}\right)^{-1}`
+is the effective Young's modulus, with :math:`\nu_i, \nu_j` the Poisson ratios
+of the particles of types *i* and *j*. :math:`E_{eff}` is denoted as *E* from here on.
+Note that if the elastic modulus and the shear modulus of the two particles are the
+same, the *hertz/material* model is equivalent to the *hertz* model with
+:math:`k_n = 4/3 E`
 
 The *dmt* model corresponds to the
 :ref:`(Derjaguin-Muller-Toporov) <DMT1975>` cohesive model, where the force
@@ -270,7 +271,8 @@ where :math:`k_n = \frac{4}{3} E_{eff}` for the *hertz/material* model. Since
 *coeff_restitution* accounts for the effective mass, effective radius, and
 pairwise overlaps (except when used with the *hooke* normal model) when calculating
 the damping coefficient, it accurately reproduces the specified coefficient of
-restitution for both monodisperse and polydisperse particle pairs.
+restitution for both monodisperse and polydisperse particle pairs.  This damping
+model is not compatible with cohesive normal models such as *JKR* or *DMT*.
 
 The total normal force is computed as the sum of the elastic and
 damping components:
@@ -441,11 +443,11 @@ discussion above. To match the Mindlin solution, one should set
 
    G_{eff} = \left(\frac{2-\nu_i}{G_i} + \frac{2-\nu_j}{G_j}\right)^{-1}
 
-where :math:`G` is the shear modulus, related to Young's modulus :math:`E`
-and Poisson's ratio :math:`\nu` by :math:`G = E/(2(1+\nu))`. This can also be
-achieved by specifying *NULL* for :math:`k_t`, in which case a
-normal contact model that specifies material parameters :math:`E` and
-:math:`\nu` is required (e.g. *hertz/material*, *dmt* or *jkr*\ ). In this
+where :math:`G_i` is the shear modulus of a particle of type :math:`i`, related to Young's
+modulus :math:`E_i` and Poisson's ratio :math:`\nu_i` by :math:`G_i = E_i/(2(1+\nu_i))`.
+This can also be achieved by specifying *NULL* for :math:`k_t`, in which case a
+normal contact model that specifies material parameters :math:`E_i` and
+:math:`\nu_i` is required (e.g. *hertz/material*, *dmt* or *jkr*\ ). In this
 case, mixing of the shear modulus for different particle types *i* and
 *j* is done according to the formula above.
 
@@ -575,7 +577,7 @@ opposite torque on each particle, according to:
 
 .. math::
 
-   \tau_{roll,i} =  R_{eff} \mathbf{n} \times \mathbf{F}_{roll}
+   \tau_{roll,i} =  R \mathbf{n} \times \mathbf{F}_{roll}
 
 .. math::
 
diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp
index 9c2c680cc5..e8521df32a 100644
--- a/src/BPM/bond_bpm.cpp
+++ b/src/BPM/bond_bpm.cpp
@@ -352,19 +352,32 @@ void BondBPM::process_broken(int i, int j)
 {
   if (!break_flag)
     error->one(FLERR, "BPM bond broke with break no option");
-  if (fix_store_local) {
-    for (int n = 0; n < nvalues; n++) (this->*pack_choice[n])(n, i, j);
 
-    fix_store_local->add_data(output_data, i, j);
+  int nlocal = atom->nlocal;
+  if (fix_store_local) {
+    // If newton off, bond can break on two procs so only record if proc owns lower tag
+    //    (BPM bond styles should sort so i -> atom with lower tag)
+    if (force->newton_bond || (i < nlocal)) {
+      for (int n = 0; n < nvalues; n++) (this->*pack_choice[n])(n, i, j);
+      fix_store_local->add_data(output_data, i, j);
+    }
   }
 
-  if (fix_update_special_bonds) fix_update_special_bonds->add_broken_bond(i, j);
+  if (fix_update_special_bonds) {
+    // If this processor owns two copies of the bond (i.e. if the domain is periodic and 1 proc thick),
+    //   skip instance where larger tag (j) owned
+    int check = 1;
+    if (i >= nlocal) {
+      int imap = atom->map(atom->tag[i]);
+      if (imap < nlocal) check = 0;
+    }
+    if (check) fix_update_special_bonds->add_broken_bond(i, j);
+  }
 
   // Manually search and remove from atom arrays
   // need to remove in case special bonds arrays rebuilt
-  int m, n;
-  int nlocal = atom->nlocal;
 
+  int m, n;
   tagint *tag = atom->tag;
   tagint **bond_atom = atom->bond_atom;
   int **bond_type = atom->bond_type;
diff --git a/src/BPM/compute_nbond_atom.cpp b/src/BPM/compute_nbond_atom.cpp
index 4f0fc4c3f0..85ef6f3fce 100644
--- a/src/BPM/compute_nbond_atom.cpp
+++ b/src/BPM/compute_nbond_atom.cpp
@@ -65,14 +65,13 @@ void ComputeNBondAtom::compute_peratom()
   tagint **bond_atom = atom->bond_atom;
   int **bond_type = atom->bond_type;
 
-  int ntotal = nlocal;
-  if (force->newton) ntotal += atom->nghost;
-
   // set local nbond array
   int i, j, k;
   int *num_bond = atom->num_bond;
   int newton_bond = force->newton_bond;
 
+  int ntotal = nlocal;
+  if (newton_bond) ntotal += atom->nghost;
   for (i = 0; i < ntotal; i++) nbond[i] = 0;
 
   for (i = 0; i < nlocal; i++) {
@@ -88,7 +87,7 @@ void ComputeNBondAtom::compute_peratom()
   }
 
   // communicate ghost nbond between neighbor procs
-  if (force->newton) comm->reverse_comm(this);
+  if (newton_bond) comm->reverse_comm(this);
 
   // zero nbond of atoms not in group
   // only do this after comm since ghost contributions must be included
diff --git a/src/BPM/fix_update_special_bonds.cpp b/src/BPM/fix_update_special_bonds.cpp
index cdc72ee987..1b408d6d4b 100644
--- a/src/BPM/fix_update_special_bonds.cpp
+++ b/src/BPM/fix_update_special_bonds.cpp
@@ -100,6 +100,7 @@ void FixUpdateSpecialBonds::pre_exchange()
       n1 = nspecial[i][0];
       for (m = 0; m < n1; m++)
         if (slist[m] == tagj) break;
+      if (m == n1) error->one(FLERR, "Special bond {} {} not found", tagi, tagj);
       for (; m < n1 - 1; m++) slist[m] = slist[m + 1];
       nspecial[i][0]--;
       nspecial[i][1] = nspecial[i][2] = nspecial[i][0];
@@ -110,6 +111,7 @@ void FixUpdateSpecialBonds::pre_exchange()
       n1 = nspecial[j][0];
       for (m = 0; m < n1; m++)
         if (slist[m] == tagi) break;
+      if (m == n1) error->one(FLERR, "Special bond {} {} not found", tagi, tagj);
       for (; m < n1 - 1; m++) slist[m] = slist[m + 1];
       nspecial[j][0]--;
       nspecial[j][1] = nspecial[j][2] = nspecial[j][0];
diff --git a/src/GRANULAR/gran_sub_mod_damping.cpp b/src/GRANULAR/gran_sub_mod_damping.cpp
index df7146a619..1745f6addf 100644
--- a/src/GRANULAR/gran_sub_mod_damping.cpp
+++ b/src/GRANULAR/gran_sub_mod_damping.cpp
@@ -156,6 +156,7 @@ double GranSubModDampingTsuji::calculate_forces()
 GranSubModDampingCoeffRestitution::GranSubModDampingCoeffRestitution(GranularModel *gm, LAMMPS *lmp) :
     GranSubModDamping(gm, lmp)
 {
+  allow_cohesion = 0;
 }
 
 void GranSubModDampingCoeffRestitution::init()
@@ -173,6 +174,12 @@ void GranSubModDampingCoeffRestitution::init()
 
 double GranSubModDampingCoeffRestitution::calculate_forces()
 {
-  damp_prefactor = damp * sqrt(gm->meff * gm->Fnormal / gm->delta);
+  // in case argument <= 0 due to precision issues
+  double sqrt1;
+  if (gm->delta > 0.0)
+    sqrt1 = MAX(0.0, gm->meff * gm->Fnormal / gm->delta);
+  else
+    sqrt1 = 0.0;
+  damp_prefactor = damp * sqrt(sqrt1);
   return -damp_prefactor * gm->vnnr;
 }
diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp
index 6e14964a35..6f887a33e6 100644
--- a/src/delete_atoms.cpp
+++ b/src/delete_atoms.cpp
@@ -22,6 +22,7 @@
 #include "comm.h"
 #include "domain.h"
 #include "error.h"
+#include "fix_bond_history.h"
 #include "force.h"
 #include "group.h"
 #include "input.h"
@@ -754,6 +755,10 @@ void DeleteAtoms::bondring(int nbuf, char *cbuf, void *ptr)
 
   int nlocal = daptr->atom->nlocal;
 
+  // find instances of bond history to delete data
+  auto histories = daptr->modify->get_fix_by_style("BOND_HISTORY");
+  int n_histories = histories.size();
+
   // cbuf = list of N deleted atom IDs from other proc, put them in hash
 
   hash->clear();
@@ -771,6 +776,11 @@ void DeleteAtoms::bondring(int nbuf, char *cbuf, void *ptr)
         if (hash->find(bond_atom[i][m]) != hash->end()) {
           bond_type[i][m] = bond_type[i][n - 1];
           bond_atom[i][m] = bond_atom[i][n - 1];
+          if (n_histories > 0)
+            for (auto &ihistory: histories) {
+              dynamic_cast<FixBondHistory *>(ihistory)->shift_history(i,m,n-1);
+              dynamic_cast<FixBondHistory *>(ihistory)->delete_history(i,n-1);
+            }
           n--;
         } else
           m++;
diff --git a/src/fix_bond_history.h b/src/fix_bond_history.h
index fafcf52bd9..e19deee82f 100644
--- a/src/fix_bond_history.h
+++ b/src/fix_bond_history.h
@@ -44,6 +44,7 @@ class FixBondHistory : public Fix {
 
   void update_atom_value(int, int, int, double);
   double get_atom_value(int, int, int);
+  int get_ndata() const { return ndata; }
 
   // methods to reorder/delete elements of atom->bond_atom
   void delete_history(int, int);