Merge pull request #2698 from jtclemm/gran_attractive

New option for granular pairstyles and updates to granular walls

doc/src/fix_wall_gran.rst

@@ -29,6 +29,7 @@ Syntax
              gamma_t = damping coefficient for collisions in tangential direction (1/time units or 1/time-distance units - see discussion below)
              xmu = static yield criterion (unitless value between 0.0 and 1.0e4)
              dampflag = 0 or 1 if tangential damping force is excluded or included
+             optional keyword = *limit_damping*, limit damping to prevent attractive interaction
 
   .. parsed-literal::
 
@@ -95,7 +96,8 @@ Specifically, delta = radius - r = overlap of particle with wall, m_eff
 = mass of particle, and the effective radius of contact = RiRj/Ri+Rj is
 set to the radius of the particle.
 
-The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu* and *dampflag*
+The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu*, *dampflag*,
+and the optional keyword *limit_damping*
 have the same meaning and units as those specified with the
 :doc:`pair_style gran/\* <pair_gran>` commands.  This means a NULL can be
 used for either *Kt* or *gamma_t* as described on that page.  If a

doc/src/fix_wall_gran_region.rst

@@ -181,7 +181,8 @@ radius - r = overlap of particle with wall, m_eff = mass of particle,
 and the effective radius of contact is just the radius of the
 particle.
 
-The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu* and *dampflag*
+The parameters *Kn*\ , *Kt*\ , *gamma_n*, *gamma_t*, *xmu*, *dampflag*,
+and the optional keyword *limit_damping*
 have the same meaning and units as those specified with the
 :doc:`pair_style gran/\* <pair_gran>` commands.  This means a NULL can be
 used for either *Kt* or *gamma_t* as described on that page.  If a

doc/src/pair_gran.rst

@@ -26,7 +26,7 @@ Syntax
 
 .. code-block:: LAMMPS
 
-   pair_style style Kn Kt gamma_n gamma_t xmu dampflag
+   pair_style style Kn Kt gamma_n gamma_t xmu dampflag keyword
 
 * style = *gran/hooke* or *gran/hooke/history* or *gran/hertz/history*
 * Kn = elastic constant for normal particle repulsion (force/distance units or pressure units - see discussion below)
@@ -36,6 +36,13 @@ Syntax
 * xmu = static yield criterion (unitless value between 0.0 and 1.0e4)
 * dampflag = 0 or 1 if tangential damping force is excluded or included
 
+* keyword = *limit_damping*
+
+  .. parsed-literal::
+
+      *limit_damping* value = none
+         limit damping to prevent attractive interaction
+
 .. note::
 
    Versions of LAMMPS before 9Jan09 had different style names for
@@ -54,6 +61,8 @@ Examples
 
    pair_style gran/hooke/history 200000.0 NULL 50.0 NULL 0.5 1
    pair_style gran/hooke 200000.0 70000.0 50.0 30.0 0.5 0
+   pair_style gran/hooke 200000.0 70000.0 50.0 30.0 0.5 0 limit_damping
+
 
 Description
 """""""""""
@@ -208,6 +217,12 @@ potential is used as a sub-style of :doc:`pair_style hybrid <pair_hybrid>`, then
 pair_coeff command to determine which atoms interact via a granular
 potential.
 
+If two particles are moving away from each other while in contact, there
+is a possibility that the particles could experience an effective attractive
+force due to damping. If the *limit_damping* keyword is used, this option
+will zero out the normal component of the force if there is an effective
+attractive force.
+
 ----------
 
 .. include:: accel_styles.rst

doc/src/pair_granular.rst

@@ -24,7 +24,7 @@ Examples
    pair_coeff * * hooke 1000.0 50.0 tangential linear_history 500.0 1.0 0.4 damping mass_velocity
 
    pair_style granular
-   pair_coeff * * hertz 1000.0 50.0 tangential mindlin 1000.0 1.0 0.4
+   pair_coeff * * hertz 1000.0 50.0 tangential mindlin 1000.0 1.0 0.4 limit_damping
 
    pair_style granular
    pair_coeff * * hertz/material 1e8 0.3 0.3 tangential mindlin_rescale NULL 1.0 0.4 damping tsuji
@@ -623,6 +623,14 @@ Finally, the twisting torque on each particle is given by:
 
 ----------
 
+If two particles are moving away from each other while in contact, there
+is a possibility that the particles could experience an effective attractive
+force due to damping. If the optional *limit_damping* keyword is used, this option
+will zero out the normal component of the force if there is an effective
+attractive force. This keyword cannot be used with the JKR or DMT models.
+
+----------
+
 The *granular* pair style can reproduce the behavior of the
 *pair gran/\** styles with the appropriate settings (some very
 minor differences can be expected due to corrections in
@@ -657,6 +665,12 @@ then LAMMPS will use that cutoff for the specified atom type
 combination, and automatically set pairwise cutoffs for the remaining
 atom types.
 
+If two particles are moving away from each other while in contact, there
+is a possibility that the particles could experience an effective attractive
+force due to damping. If the *limit_damping* keyword is used, this option
+will zero out the normal component of the force if there is an effective
+attractive force. This keyword cannot be used with the JKR or DMT models.
+
 ----------
 
 .. include:: accel_styles.rst

src/GRANULAR/fix_wall_gran.cpp

@@ -54,7 +54,8 @@ enum {XPLANE=0,YPLANE=1,ZPLANE=2,ZCYLINDER,REGION};
 enum {NONE,CONSTANT,EQUAL};
 enum {DAMPING_NONE, VELOCITY, MASS_VELOCITY, VISCOELASTIC, TSUJI};
 enum {TANGENTIAL_NONE, TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY,
-  TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE};
+  TANGENTIAL_MINDLIN, TANGENTIAL_MINDLIN_RESCALE,
+  TANGENTIAL_MINDLIN_FORCE, TANGENTIAL_MINDLIN_RESCALE_FORCE};
 enum {TWIST_NONE, TWIST_SDS, TWIST_MARSHALL};
 enum {ROLL_NONE, ROLL_SDS};
 
@@ -70,6 +71,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
     error->all(FLERR,"Fix wall/gran requires atom style sphere");
 
   create_attribute = 1;
+  limit_damping = 0;
 
   // set interaction style
   // disable bonded/history option for now
@@ -90,7 +92,6 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
   // wall/particle coefficients
 
   int iarg;
-
   if (pairstyle != GRANULAR) {
     size_history = 3;
     if (narg < 11) error->all(FLERR,"Illegal fix wall/gran command");
@@ -118,6 +119,12 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
       kt /= force->nktv2p;
     }
     iarg = 10;
+
+    if (strcmp(arg[iarg],"limit_damping") == 0) {
+      limit_damping = 1;
+      iarg += 1;
+    }
+
   } else {
     iarg = 4;
     damping_model = VISCOELASTIC;
@@ -211,7 +218,9 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
           iarg += 4;
         } else if ((strcmp(arg[iarg+1], "linear_history") == 0) ||
             (strcmp(arg[iarg+1], "mindlin") == 0) ||
-            (strcmp(arg[iarg+1], "mindlin_rescale") == 0)) {
+            (strcmp(arg[iarg+1], "mindlin_rescale") == 0) ||
+            (strcmp(arg[iarg+1], "mindlin/force") == 0) ||
+            (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)) {
           if (iarg + 4 >= narg)
             error->all(FLERR,"Illegal pair_coeff command, "
                        "not enough parameters provided for tangential model");
@@ -221,8 +230,14 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
             tangential_model = TANGENTIAL_MINDLIN;
           else if (strcmp(arg[iarg+1], "mindlin_rescale") == 0)
             tangential_model = TANGENTIAL_MINDLIN_RESCALE;
+          else if (strcmp(arg[iarg+1], "mindlin/force") == 0)
+            tangential_model = TANGENTIAL_MINDLIN_FORCE;
+          else if (strcmp(arg[iarg+1], "mindlin_rescale/force") == 0)
+            tangential_model = TANGENTIAL_MINDLIN_RESCALE_FORCE;
           if ((tangential_model == TANGENTIAL_MINDLIN ||
-               tangential_model == TANGENTIAL_MINDLIN_RESCALE) &&
+               tangential_model == TANGENTIAL_MINDLIN_RESCALE ||
+               tangential_model == TANGENTIAL_MINDLIN_FORCE ||
+               tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) &&
               (strcmp(arg[iarg+2], "NULL") == 0)) {
             if (normal_model == NORMAL_HERTZ || normal_model == NORMAL_HOOKE) {
               error->all(FLERR, "NULL setting for Mindlin tangential "
@@ -293,13 +308,27 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
           strcmp(arg[iarg], "zcylinder") == 0 ||
           strcmp(arg[iarg], "region") == 0) {
         break;
+      } else if (strcmp(arg[iarg],"limit_damping") == 0) {
+        limit_damping = 1;
+        iarg += 1;
       } else {
         error->all(FLERR, "Illegal fix wall/gran command");
       }
     }
     size_history = 3*tangential_history + 3*roll_history + twist_history;
+    //Unlike the pair style, the wall style does not have a 'touch'
+    //array. Hence, an additional entry in the history is used to
+    //determine if particles previously contacted for JKR cohesion purposes.
     if (normal_model == JKR) size_history += 1;
-    if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) size_history += 1;
+    if (tangential_model == TANGENTIAL_MINDLIN_RESCALE ||
+        tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) size_history += 1;
+
+    if (limit_damping && normal_model == JKR)
+      error->all(FLERR,"Illegal pair_coeff command, "
+          "cannot limit damping with JRK model");
+    if (limit_damping && normal_model == DMT)
+      error->all(FLERR,"Illegal pair_coeff command, "
+          "Cannot limit damping with DMT model");
   }
 
   // wallstyle args
@@ -475,6 +504,7 @@ void FixWallGran::init()
     if (modify->fix[i]->rigid_flag) break;
   if (i < modify->nfix) fix_rigid = modify->fix[i];
 
+  if(pairstyle == GRANULAR) {
     tangential_history_index = 0;
     if (roll_history) {
       if (tangential_history) roll_history_index = 3;
@@ -495,7 +525,8 @@ void FixWallGran::init()
       roll_history_index += 1;
       twist_history_index += 1;
     }
-  if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) {
+    if (tangential_model == TANGENTIAL_MINDLIN_RESCALE ||
+        tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE) {
       roll_history_index += 1;
       twist_history_index += 1;
     }
@@ -507,6 +538,7 @@ void FixWallGran::init()
           4.8218*pow(cor,6);
     }
   }
+}
 
 /* ---------------------------------------------------------------------- */
 
@@ -762,6 +794,7 @@ void FixWallGran::hooke(double rsq, double dx, double dy, double dz,
 
   damp = meff*gamman*vnnr*rsqinv;
   ccel = kn*(radius-r)*rinv - damp;
+  if (limit_damping && (ccel < 0.0)) ccel = 0.0;
 
   // relative velocities
 
@@ -854,6 +887,7 @@ void FixWallGran::hooke_history(double rsq, double dx, double dy, double dz,
 
   damp = meff*gamman*vnnr*rsqinv;
   ccel = kn*(radius-r)*rinv - damp;
+  if (limit_damping && (ccel < 0.0)) ccel = 0.0;
 
   // relative velocities
 
@@ -985,6 +1019,7 @@ void FixWallGran::hertz_history(double rsq, double dx, double dy, double dz,
   if (rwall == 0.0) polyhertz = sqrt((radius-r)*radius);
   else polyhertz = sqrt((radius-r)*radius*rwall/(rwall+radius));
   ccel *= polyhertz;
+  if (limit_damping && (ccel < 0.0)) ccel = 0.0;
 
   // relative velocities
 
@@ -1103,7 +1138,8 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
   double signtwist, magtwist, magtortwist, Mtcrit;
   double tortwist1, tortwist2, tortwist3;
 
-  double shrmag,rsht;
+  double shrmag,rsht,prjmag;
+  bool frameupdate;
 
   r = sqrt(rsq);
   E = normal_coeffs[0];
@@ -1178,11 +1214,18 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
   Fdamp = -damp_normal_prefactor*vnnr;
 
   Fntot = Fne + Fdamp;
+  if (limit_damping && (Fntot < 0.0)) Fntot = 0.0;
 
   //****************************************
   // tangential force, including history effects
   //****************************************
 
+  // For linear, mindlin, mindlin_rescale:
+  // history = cumulative tangential displacement
+  //
+  // For mindlin/force, mindlin_rescale/force:
+  // history = cumulative tangential elastic force
+
   // tangential component
   vt1 = vr1 - vn1;
   vt2 = vr2 - vn2;
@@ -1224,69 +1267,115 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
   int thist2 = thist1 + 1;
 
   if (tangential_history) {
-    if (tangential_model == TANGENTIAL_MINDLIN) {
+    if (tangential_model == TANGENTIAL_MINDLIN ||
+        tangential_model == TANGENTIAL_MINDLIN_FORCE) {
       k_tangential *= a;
     }
-    else if (tangential_model == TANGENTIAL_MINDLIN_RESCALE) {
+    else if (tangential_model ==
+             TANGENTIAL_MINDLIN_RESCALE ||
+             tangential_model ==
+             TANGENTIAL_MINDLIN_RESCALE_FORCE){
       k_tangential *= a;
-      if (a < history[3]) { //On unloading, rescale the shear displacements
+      // on unloading, rescale the shear displacements/force
+      if (a < history[thist2+1]) {
         double factor = a/history[thist2+1];
         history[thist0] *= factor;
         history[thist1] *= factor;
         history[thist2] *= factor;
       }
     }
-    shrmag = sqrt(history[thist0]*history[thist0] +
+
-                  history[thist1]*history[thist1] +
-                  history[thist2]*history[thist2]);
 
     // rotate and update displacements.
     // see e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235
     if (history_update) {
       rsht = history[thist0]*nx + history[thist1]*ny + history[thist2]*nz;
-      if (fabs(rsht) < EPSILON) rsht = 0;
+      if (tangential_model == TANGENTIAL_MINDLIN_FORCE ||
-      if (rsht > 0) {
+          tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE)
-        // if rhst == shrmag, contacting pair has rotated 90 deg in one step,
+        frameupdate = fabs(rsht) > EPSILON*Fscrit;
-        // in which case you deserve a crash!
+      else
-        scalefac = shrmag/(shrmag - rsht);
+        frameupdate = fabs(rsht)*k_tangential > EPSILON*Fscrit;
+      if (frameupdate) {
+        shrmag = sqrt(history[thist0]*history[thist0] +
+            history[thist1]*history[thist1] +
+            history[thist2]*history[thist2]);
+        // projection
         history[thist0] -= rsht*nx;
         history[thist1] -= rsht*ny;
         history[thist2] -= rsht*nz;
+
         // also rescale to preserve magnitude
+        prjmag = sqrt(history[thist0]*history[thist0] +
+            history[thist1]*history[thist1] + history[thist2]*history[thist2]);
+        if (prjmag > 0) scalefac = shrmag/prjmag;
+        else scalefac = 0;
         history[thist0] *= scalefac;
         history[thist1] *= scalefac;
         history[thist2] *= scalefac;
       }
       // update history
+      if (tangential_model == TANGENTIAL_HISTORY ||
+          tangential_model == TANGENTIAL_MINDLIN ||
+          tangential_model == TANGENTIAL_MINDLIN_RESCALE) {
       history[thist0] += vtr1*dt;
       history[thist1] += vtr2*dt;
       history[thist2] += vtr3*dt;
+      } else{
+        // tangential force
+        // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46
+        history[thist0] -= k_tangential*vtr1*dt;
+        history[thist1] -= k_tangential*vtr2*dt;
+        history[thist2] -= k_tangential*vtr3*dt;
+      }
+      if (tangential_model == TANGENTIAL_MINDLIN_RESCALE ||
+          tangential_model == TANGENTIAL_MINDLIN_RESCALE_FORCE)
+        history[thist2+1] = a;
     }
 
     // tangential forces = history + tangential velocity damping
+    if (tangential_model == TANGENTIAL_HISTORY ||
+        tangential_model == TANGENTIAL_MINDLIN ||
+        tangential_model == TANGENTIAL_MINDLIN_RESCALE) {
       fs1 = -k_tangential*history[thist0] - damp_tangential*vtr1;
       fs2 = -k_tangential*history[thist1] - damp_tangential*vtr2;
       fs3 = -k_tangential*history[thist2] - damp_tangential*vtr3;
+    } else {
+      fs1 = history[thist0] - damp_tangential*vtr1;
+      fs2 = history[thist1] - damp_tangential*vtr2;
+      fs3 = history[thist2] - damp_tangential*vtr3;
+    }
 
     // rescale frictional displacements and forces if needed
     Fscrit = tangential_coeffs[2] * Fncrit;
     fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3);
     if (fs > Fscrit) {
+      shrmag = sqrt(history[thist0]*history[thist0] +
+          history[thist1]*history[thist1] +
+          history[thist2]*history[thist2]);
       if (shrmag != 0.0) {
+        if (tangential_model == TANGENTIAL_HISTORY ||
+            tangential_model == TANGENTIAL_MINDLIN ||
+            tangential_model ==
+                TANGENTIAL_MINDLIN_RESCALE) {
           history[thist0] = -1.0/k_tangential*(Fscrit*fs1/fs +
               damp_tangential*vtr1);
           history[thist1] = -1.0/k_tangential*(Fscrit*fs2/fs +
               damp_tangential*vtr2);
           history[thist2] = -1.0/k_tangential*(Fscrit*fs3/fs +
               damp_tangential*vtr3);
+        } else {
+          history[thist0] = Fscrit*fs1/fs + damp_tangential*vtr1;
+          history[thist1] = Fscrit*fs2/fs + damp_tangential*vtr2;
+          history[thist2] = Fscrit*fs3/fs + damp_tangential*vtr3;
+        }
         fs1 *= Fscrit/fs;
         fs2 *= Fscrit/fs;
         fs3 *= Fscrit/fs;
       } else fs1 = fs2 = fs3 = 0.0;
     }
   } else { // classic pair gran/hooke (no history)
-    fs = meff*damp_tangential*vrel;
+    fs = damp_tangential*vrel;
-    if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel;
+    if (vrel != 0.0) Ft = MIN(Fscrit,fs) / vrel;
     else Ft = 0.0;
     fs1 = -Ft*vtr1;
     fs2 = -Ft*vtr2;
@@ -1297,14 +1386,14 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
   // rolling resistance
   //****************************************
 
-  if (roll_model != ROLL_NONE || twist_model != NONE) {
+  if (roll_model != ROLL_NONE || twist_model != TWIST_NONE) {
     relrot1 = omega[0];
     relrot2 = omega[1];
     relrot3 = omega[2];
   }
   if (roll_model != ROLL_NONE) {
-
+    // rolling velocity,
-    // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
+    // see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
     // This is different from the Marshall papers,
     // which use the Bagi/Kuhn formulation
     // for rolling velocity (see Wang et al for why the latter is wrong)
@@ -1316,21 +1405,29 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
     int rhist1 = rhist0 + 1;
     int rhist2 = rhist1 + 1;
 
-    // rolling displacement
+    k_roll = roll_coeffs[0];
+    damp_roll = roll_coeffs[1];
+    Frcrit = roll_coeffs[2] * Fncrit;
+
+    if (history_update) {
+      rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz;
+      frameupdate = fabs(rolldotn)*k_roll > EPSILON*Frcrit;
+      if (frameupdate) { // rotate into tangential plane
         rollmag = sqrt(history[rhist0]*history[rhist0] +
                            history[rhist1]*history[rhist1] +
                            history[rhist2]*history[rhist2]);
-
+        // projection
-    rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz;
-
-    if (history_update) {
-      if (fabs(rolldotn) < EPSILON) rolldotn = 0;
-      if (rolldotn > 0) { // rotate into tangential plane
-        scalefac = rollmag/(rollmag - rolldotn);
         history[rhist0] -= rolldotn*nx;
         history[rhist1] -= rolldotn*ny;
         history[rhist2] -= rolldotn*nz;
+
         // also rescale to preserve magnitude
+        prjmag = sqrt(history[rhist0]*history[rhist0] +
+            history[rhist1]*history[rhist1] +
+            history[rhist2]*history[rhist2]);
+
+        if (prjmag > 0) scalefac = rollmag/prjmag;
+        else scalefac = 0;
         history[rhist0] *= scalefac;
         history[rhist1] *= scalefac;
         history[rhist2] *= scalefac;
@@ -1340,17 +1437,16 @@ void FixWallGran::granular(double rsq, double dx, double dy, double dz,
       history[rhist2] += vrl3*dt;
     }
 
-    k_roll = roll_coeffs[0];
-    damp_roll = roll_coeffs[1];
     fr1 = -k_roll*history[rhist0] - damp_roll*vrl1;
     fr2 = -k_roll*history[rhist1] - damp_roll*vrl2;
     fr3 = -k_roll*history[rhist2] - damp_roll*vrl3;
 
     // rescale frictional displacements and forces if needed
-    Frcrit = roll_coeffs[2] * Fncrit;
-
     fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3);
     if (fr > Frcrit) {
+      rollmag = sqrt(history[rhist0]*history[rhist0] +
+          history[rhist1]*history[rhist1] +
+          history[rhist2]*history[rhist2]);
       if (rollmag != 0.0) {
         history[rhist0] = -1.0/k_roll*(Frcrit*fr1/fr + damp_roll*vrl1);
         history[rhist1] = -1.0/k_roll*(Frcrit*fr2/fr + damp_roll*vrl2);

src/GRANULAR/fix_wall_gran.h

@@ -73,6 +73,7 @@ class FixWallGran : public Fix {
   // for granular model choices
   int normal_model, damping_model;
   int tangential_model, roll_model, twist_model;
+  int limit_damping;
 
   // history flags
   int normal_history, tangential_history, roll_history, twist_history;

src/GRANULAR/pair_gran_hertz_history.cpp

@@ -181,6 +181,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
         ccel = kn*(radsum-r)*rinv - damp;
         polyhertz = sqrt((radsum-r)*radi*radj / radsum);
         ccel *= polyhertz;
+        if (limit_damping && (ccel < 0.0)) ccel = 0.0;
 
         // relative velocities
 
@@ -277,7 +278,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag)
 
 void PairGranHertzHistory::settings(int narg, char **arg)
 {
-  if (narg != 6) error->all(FLERR,"Illegal pair_style command");
+  if (narg != 6 && narg != 7) error->all(FLERR,"Illegal pair_style command");
 
   kn = utils::numeric(FLERR,arg[0],false,lmp);
   if (strcmp(arg[1],"NULL") == 0) kt = kn * 2.0/7.0;
@@ -295,6 +296,12 @@ void PairGranHertzHistory::settings(int narg, char **arg)
       xmu < 0.0 || xmu > 10000.0 || dampflag < 0 || dampflag > 1)
     error->all(FLERR,"Illegal pair_style command");
 
+  limit_damping = 0;
+  if (narg == 7) {
+    if (strcmp(arg[6], "limit_damping") == 0) limit_damping = 1;
+    else error->all(FLERR,"Illegal pair_style command");
+  }
+
   // convert Kn and Kt from pressure units to force/distance^2
 
   kn /= force->nktv2p;
@@ -388,6 +395,7 @@ double PairGranHertzHistory::single(int i, int j, int /*itype*/, int /*jtype*/,
   ccel = kn*(radsum-r)*rinv - damp;
   polyhertz = sqrt((radsum-r)*radi*radj / radsum);
   ccel *= polyhertz;
+  if (limit_damping && (ccel < 0.0)) ccel = 0.0;
 
   // relative velocities
 

src/GRANULAR/pair_gran_hooke.cpp

@@ -158,6 +158,7 @@ void PairGranHooke::compute(int eflag, int vflag)
 
         damp = meff*gamman*vnnr*rsqinv;
         ccel = kn*(radsum-r)*rinv - damp;
+        if (limit_damping && (ccel < 0.0)) ccel = 0.0;
 
         // relative velocities
 
@@ -299,6 +300,7 @@ double PairGranHooke::single(int i, int j, int /*itype*/, int /*jtype*/, double
 
   damp = meff*gamman*vnnr*rsqinv;
   ccel = kn*(radsum-r)*rinv - damp;
+  if (limit_damping && (ccel < 0.0)) ccel = 0.0;
 
   // relative velocities
 

src/GRANULAR/pair_gran_hooke_history.cpp

@@ -237,6 +237,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag)
 
         damp = meff*gamman*vnnr*rsqinv;
         ccel = kn*(radsum-r)*rinv - damp;
+        if (limit_damping && (ccel < 0.0)) ccel = 0.0;
 
         // relative velocities
 
@@ -356,7 +357,7 @@ void PairGranHookeHistory::allocate()
 
 void PairGranHookeHistory::settings(int narg, char **arg)
 {
-  if (narg != 6) error->all(FLERR,"Illegal pair_style command");
+  if (narg != 6 && narg != 7) error->all(FLERR,"Illegal pair_style command");
 
   kn = utils::numeric(FLERR,arg[0],false,lmp);
   if (strcmp(arg[1],"NULL") == 0) kt = kn * 2.0/7.0;
@@ -370,6 +371,12 @@ void PairGranHookeHistory::settings(int narg, char **arg)
   dampflag = utils::inumeric(FLERR,arg[5],false,lmp);
   if (dampflag == 0) gammat = 0.0;
 
+  limit_damping = 0;
+  if (narg == 7) {
+    if (strcmp(arg[6], "limit_damping") == 0) limit_damping = 1;
+    else error->all(FLERR,"Illegal pair_style command");
+  }
+
   if (kn < 0.0 || kt < 0.0 || gamman < 0.0 || gammat < 0.0 ||
       xmu < 0.0 || xmu > 10000.0 || dampflag < 0 || dampflag > 1)
     error->all(FLERR,"Illegal pair_style command");
@@ -686,6 +693,7 @@ double PairGranHookeHistory::single(int i, int j, int /*itype*/, int /*jtype*/,
 
   damp = meff*gamman*vnnr*rsqinv;
   ccel = kn*(radsum-r)*rinv - damp;
+  if(limit_damping && (ccel < 0.0)) ccel = 0.0;
 
   // relative velocities
 

src/GRANULAR/pair_gran_hooke_history.h

@@ -49,6 +49,7 @@ class PairGranHookeHistory : public Pair {
   double dt;
   int freeze_group_bit;
   int history;
+  int limit_damping;
 
   int neighprev;
   double *onerad_dynamic,*onerad_frozen;

src/GRANULAR/pair_granular.cpp

@@ -81,6 +81,8 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp)
   maxrad_dynamic = nullptr;
   maxrad_frozen = nullptr;
 
+  limit_damping = nullptr;
+
   history_transfer_factors = nullptr;
 
   dt = update->dt;
@@ -130,6 +132,7 @@ PairGranular::~PairGranular()
     memory->destroy(tangential_model);
     memory->destroy(roll_model);
     memory->destroy(twist_model);
+    memory->destroy(limit_damping);
 
     delete [] onerad_dynamic;
     delete [] onerad_frozen;
@@ -364,12 +367,13 @@ void PairGranular::compute(int eflag, int vflag)
           damp_normal = a*meff;
         } else if (damping_model[itype][jtype] == TSUJI) {
           damp_normal = sqrt(meff*knfac);
-        }
+        } else damp_normal = 0.0;
 
         damp_normal_prefactor = normal_coeffs[itype][jtype][1]*damp_normal;
         Fdamp = -damp_normal_prefactor*vnnr;
 
         Fntot = Fne + Fdamp;
+        if (limit_damping[itype][jtype] && (Fntot < 0.0)) Fntot = 0.0;
 
         //****************************************
         // tangential force, including history effects
@@ -540,20 +544,17 @@ void PairGranular::compute(int eflag, int vflag)
           relrot1 = omega[i][0] - omega[j][0];
           relrot2 = omega[i][1] - omega[j][1];
           relrot3 = omega[i][2] - omega[j][2];
-          // rolling velocity,
-          // see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
-          // this is different from the Marshall papers,
-          // which use the Bagi/Kuhn formulation
-          // for rolling velocity (see Wang et al for why the latter is wrong)
-          // - 0.5*((radj-radi)/radsum)*vtr1;
-          // - 0.5*((radj-radi)/radsum)*vtr2;
-          // - 0.5*((radj-radi)/radsum)*vtr3;
         }
         //****************************************
         // rolling resistance
         //****************************************
 
         if (roll_model[itype][jtype] != ROLL_NONE) {
+          // rolling velocity,
+          // see eq. 31 of Wang et al, Particuology v 23, p 49 (2015)
+          // this is different from the Marshall papers,
+          // which use the Bagi/Kuhn formulation
+          // for rolling velocity (see Wang et al for why the latter is wrong)
           vrl1 = Reff*(relrot2*nz - relrot3*ny);
           vrl2 = Reff*(relrot3*nx - relrot1*nz);
           vrl3 = Reff*(relrot1*ny - relrot2*nx);
@@ -740,6 +741,7 @@ void PairGranular::allocate()
   memory->create(tangential_model,n+1,n+1,"pair:tangential_model");
   memory->create(roll_model,n+1,n+1,"pair:roll_model");
   memory->create(twist_model,n+1,n+1,"pair:twist_model");
+  memory->create(limit_damping,n+1,n+1,"pair:limit_damping");
 
   onerad_dynamic = new double[n+1];
   onerad_frozen = new double[n+1];
@@ -793,6 +795,7 @@ void PairGranular::coeff(int narg, char **arg)
   roll_model_one = ROLL_NONE;
   twist_model_one = TWIST_NONE;
   damping_model_one = VISCOELASTIC;
+  int ld_flag = 0;
 
   int iarg = 2;
   while (iarg < narg) {
@@ -967,12 +970,15 @@ void PairGranular::coeff(int narg, char **arg)
         error->all(FLERR, "Illegal pair_coeff command, not enough parameters");
       cutoff_one = utils::numeric(FLERR,arg[iarg+1],false,lmp);
       iarg += 2;
-    } else error->all(FLERR, "Illegal pair coeff command");
+    } else if (strcmp(arg[iarg], "limit_damping") == 0) {
+      ld_flag = 1;
+      iarg += 1;
+    } else error->all(FLERR, "Illegal pair_coeff command");
   }
 
   // error not to specify normal or tangential model
   if ((normal_model_one < 0) || (tangential_model_one < 0))
-    error->all(FLERR, "Illegal pair coeff command, "
+    error->all(FLERR, "Illegal pair_coeff command, "
                "must specify normal or tangential contact model");
 
   int count = 0;
@@ -984,6 +990,14 @@ void PairGranular::coeff(int narg, char **arg)
         27.467*powint(cor,4)-18.022*powint(cor,5)+4.8218*powint(cor,6);
   } else damp = normal_coeffs_one[1];
 
+  if (ld_flag && normal_model_one == JKR)
+    error->all(FLERR,"Illegal pair_coeff command, "
+        "Cannot limit damping with JKR model");
+
+  if (ld_flag && normal_model_one == DMT)
+    error->all(FLERR,"Illegal pair_coeff command, "
+        "Cannot limit damping with DMT model");
+
   for (int i = ilo; i <= ihi; i++) {
     for (int j = MAX(jlo,i); j <= jhi; j++) {
       normal_model[i][j] = normal_model[j][i] = normal_model_one;
@@ -1026,11 +1040,14 @@ void PairGranular::coeff(int narg, char **arg)
 
       cutoff_type[i][j] = cutoff_type[j][i] = cutoff_one;
 
+      limit_damping[i][j] = limit_damping[j][i] = ld_flag;
+
       setflag[i][j] = 1;
       count++;
     }
   }
 
+
   if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
 }
 
@@ -1303,6 +1320,7 @@ void PairGranular::write_restart(FILE *fp)
         fwrite(&tangential_model[i][j],sizeof(int),1,fp);
         fwrite(&roll_model[i][j],sizeof(int),1,fp);
         fwrite(&twist_model[i][j],sizeof(int),1,fp);
+        fwrite(&limit_damping[i][j],sizeof(int),1,fp);
         fwrite(normal_coeffs[i][j],sizeof(double),4,fp);
         fwrite(tangential_coeffs[i][j],sizeof(double),3,fp);
         fwrite(roll_coeffs[i][j],sizeof(double),3,fp);
@@ -1333,6 +1351,7 @@ void PairGranular::read_restart(FILE *fp)
           utils::sfread(FLERR,&tangential_model[i][j],sizeof(int),1,fp,nullptr,error);
           utils::sfread(FLERR,&roll_model[i][j],sizeof(int),1,fp,nullptr,error);
           utils::sfread(FLERR,&twist_model[i][j],sizeof(int),1,fp,nullptr,error);
+          utils::sfread(FLERR,&limit_damping[i][j],sizeof(int),1,fp,nullptr,error);
           utils::sfread(FLERR,normal_coeffs[i][j],sizeof(double),4,fp,nullptr,error);
           utils::sfread(FLERR,tangential_coeffs[i][j],sizeof(double),3,fp,nullptr,error);
           utils::sfread(FLERR,roll_coeffs[i][j],sizeof(double),3,fp,nullptr,error);
@@ -1344,6 +1363,7 @@ void PairGranular::read_restart(FILE *fp)
         MPI_Bcast(&tangential_model[i][j],1,MPI_INT,0,world);
         MPI_Bcast(&roll_model[i][j],1,MPI_INT,0,world);
         MPI_Bcast(&twist_model[i][j],1,MPI_INT,0,world);
+        MPI_Bcast(&limit_damping[i][j],1,MPI_INT,0,world);
         MPI_Bcast(normal_coeffs[i][j],4,MPI_DOUBLE,0,world);
         MPI_Bcast(tangential_coeffs[i][j],3,MPI_DOUBLE,0,world);
         MPI_Bcast(roll_coeffs[i][j],3,MPI_DOUBLE,0,world);
@@ -1529,6 +1549,7 @@ double PairGranular::single(int i, int j, int itype, int jtype,
   Fdamp = -damp_normal_prefactor*vnnr;
 
   Fntot = Fne + Fdamp;
+  if (limit_damping[itype][jtype] && (Fntot < 0.0)) Fntot = 0.0;
 
   jnum = list->numneigh[i];
   jlist = list->firstneigh[i];

src/GRANULAR/pair_granular.h

@@ -70,6 +70,7 @@ class PairGranular : public Pair {
   // model choices
   int **normal_model, **damping_model;
   int **tangential_model, **roll_model, **twist_model;
+  int **limit_damping;
 
   // history flags
   int normal_history, tangential_history, roll_history, twist_history;

src/KOKKOS/pair_gran_hooke_history_kokkos.cpp

@@ -392,6 +392,7 @@ void PairGranHookeHistoryKokkos<DeviceType>::operator()(TagPairGranHookeHistoryC
 
     F_FLOAT damp = meff*gamman*vnnr*rsqinv;
     F_FLOAT ccel = kn*(radsum-r)*rinv - damp;
+    if(limit_damping && (ccel < 0.0)) ccel = 0.0;
 
     // relative velocities
 

src/USER-OMP/pair_gran_hertz_history_omp.cpp

@@ -224,6 +224,7 @@ void PairGranHertzHistoryOMP::eval(int iifrom, int iito, ThrData * const thr)
         ccel = kn*(radsum-r)*rinv - damp;
         polyhertz = sqrt((radsum-r)*radi*radj / radsum);
         ccel *= polyhertz;
+        if (limit_damping && (ccel < 0.0)) ccel = 0.0;
 
         // relative velocities
 

src/USER-OMP/pair_gran_hooke_history_omp.cpp

@@ -223,6 +223,7 @@ void PairGranHookeHistoryOMP::eval(int iifrom, int iito, ThrData * const thr)
 
         damp = meff*gamman*vnnr*rsqinv;
         ccel = kn*(radsum-r)*rinv - damp;
+        if (limit_damping && (ccel < 0.0)) ccel = 0.0;
 
         // relative velocities
 

src/USER-OMP/pair_gran_hooke_omp.cpp

@@ -197,6 +197,7 @@ void PairGranHookeOMP::eval(int iifrom, int iito, ThrData * const thr)
 
         damp = meff*gamman*vnnr*rsqinv;
         ccel = kn*(radsum-r)*rinv - damp;
+        if (limit_damping && (ccel < 0.0)) ccel = 0.0;
 
         // relative velocities