Misc minor patches/features in BPM/Granular packages

doc/src/bond_bpm_rotational.rst

@@ -24,14 +24,17 @@ Syntax
             *x, y, z* = the center of mass position of the 2 atoms when the bond broke (distance units)
             *x/ref, y/ref, z/ref* = the initial center of mass position of the 2 atoms (distance units)
 
-       *overlay/pair* value = none
+       *overlay/pair* value = *yes* or *no*
           bonded particles will still interact with pair forces
 
        *smooth* value = *yes* or *no*
           smooths bond forces near the breaking point
 
-       *break/no*
-          indicates that bonds should not break during a run
+       *normalize* value = *yes* or *no*
+          normalizes normal and shear forces by the reference length
+
+       *break* value = *yes* or *no*
+          indicates whether bonds break during a run
 
 Examples
 """"""""
@@ -136,16 +139,19 @@ or :doc:`read_restart <read_restart>` commands:
 * :math:`\gamma_r`      (force*distance/velocity units)
 * :math:`\gamma_t`      (force*distance/velocity units)
 
+However, the *normalize* option will normalize the radial and shear forces
+by :math:`r_0` such that :math:`k_r` and :math:`k_s` are unit less.
+
 By default, pair forces are not calculated between bonded particles.
 Pair forces can alternatively be overlaid on top of bond forces using
-the *overlay/pair* keyword. These settings require specific
+the *overlay/pair* option. These settings require specific
 :doc:`special_bonds <special_bonds>` settings described in the
 restrictions.  Further details can be found in the `:doc: how to
 <Howto_BPM>` page on BPMs.
 
 .. versionadded:: 28Mar2023
 
-If the *break/no* keyword is used, then LAMMPS assumes bonds should not break
+If the *break* option is used, then LAMMPS assumes bonds should not break
 during a simulation run. This will prevent some unnecessary calculation.
 However, if a bond does break, it will trigger an error.
 
@@ -251,7 +257,7 @@ Related commands
 Default
 """""""
 
-The option defaults are *smooth* = *yes*
+The option defaults are *overlay/pair* = *no*, *smooth* = *yes*, *normalize* = *no*, and *break* = *yes*
 
 ----------
 

doc/src/bond_bpm_spring.rst

@@ -24,14 +24,17 @@ Syntax
             *x, y, z* = the center of mass position of the 2 atoms when the bond broke (distance units)
             *x/ref, y/ref, z/ref* = the initial center of mass position of the 2 atoms (distance units)
 
-       *overlay/pair* value = none
+       *overlay/pair* value = *yes* or *no*
           bonded particles will still interact with pair forces
 
        *smooth* value = *yes* or *no*
           smooths bond forces near the breaking point
 
-       *break/no*
-          indicates that bonds should not break during a run
+       *normalize* value = *yes* or *no*
+          normalizes bond forces by the reference length
+
+       *break* value = *yes* or *no*
+          indicates whether bonds break during a run
 
 Examples
 """"""""
@@ -66,7 +69,7 @@ particles based on a model described by Clemmer and Robbins
 
    F = k (r - r_0) w
 
-where :math:`k_r` is a stiffness, :math:`r` is the current distance
+where :math:`k` is a stiffness, :math:`r` is the current distance
 and :math:`r_0` is the initial distance between the two particles, and
 :math:`w` is an optional smoothing factor discussed below. Bonds will
 break at a strain of :math:`\epsilon_c`.  This is done by setting by
@@ -102,16 +105,19 @@ the data file or restart files read by the :doc:`read_data
 * :math:`\epsilon_c`    (unit less)
 * :math:`\gamma`        (force/velocity units)
 
+However, the *normalize* option will normalize the elastic bond force by
+:math:`r_0` such that :math:`k` is unit less.
+
 By default, pair forces are not calculated between bonded particles.
 Pair forces can alternatively be overlaid on top of bond forces using
-the *overlay/pair* keyword. These settings require specific
+the *overlay/pair* option. These settings require specific
 :doc:`special_bonds <special_bonds>` settings described in the
 restrictions.  Further details can be found in the `:doc: how to
 <Howto_BPM>` page on BPMs.
 
 .. versionadded:: 28Mar2023
 
-If the *break/no* keyword is used, then LAMMPS assumes bonds should not break
+If the *break* option is used, then LAMMPS assumes bonds should not break
 during a simulation run. This will prevent some unnecessary calculation.
 However, if a bond does break, it will trigger an error.
 
@@ -206,7 +212,7 @@ Related commands
 Default
 """""""
 
-The option defaults are *smooth* = *yes*
+The option defaults are *overlay/pair* = *no*, *smooth* = *yes*, *normalize* = *no*, and *break* = *yes*
 
 ----------
 

doc/src/compute_bond_local.rst

@@ -76,7 +76,10 @@ The value *force* is the magnitude of the force acting between the
 pair of atoms in the bond.
 
 The values *fx*, *fy*, and *fz* are the xyz components of
-*force* between the pair of atoms in the bond.
+*force* between the pair of atoms in the bond. For bond styles that apply
+non-central forces, such as :doc:`bond_style bpm/rotational
+<bond_bpm_rotational>`, these values only include the :math:`(x,y,z)`
+components of the normal force component.
 
 The remaining properties are all computed for motion of the two atoms
 relative to the center of mass (COM) velocity of the 2 atoms in the

doc/src/compute_fabric.rst

@@ -146,13 +146,13 @@ m to :math:`M` (inclusive).  A middle asterisk means all types from m to n
 Output info
 """""""""""
 
-This compute calculates a local vector of doubles and a scalar. The vector
-stores the unique components of the first requested tensor in the order
-:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`
-followed by the same components for all subsequent tensors.
+This compute calculates a global vector of doubles and a global scalar. The
+vector stores the unique components of the first requested tensor in the
+order :math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`,
+:math:`yz` followed by the same components for all subsequent tensors.
 The length of the vector is therefore six times the number of requested
-tensors. The scalar output is the number of pairwise interactions included in
-the calculation of the fabric tensor.
+tensors. The scalar output is the number of pairwise interactions included
+in the calculation of the fabric tensor.
 
 Restrictions
 """"""""""""

doc/src/compute_pair_local.rst

@@ -66,7 +66,9 @@ The value *eng* is the interaction energy for the pair of atoms.
 The value *force* is the force acting between the pair of atoms, which
 is positive for a repulsive force and negative for an attractive
 force.  The values *fx*, *fy*, and *fz* are the :math:`(x,y,z)` components of
-*force* on atom I.
+*force* on atom I. For pair styles that apply non-central forces,
+such as :doc:`granular pair styles <pair_gran>`, these values only include
+the :math:`(x,y,z)` components of the normal force component.
 
 A pair style may define additional pairwise quantities which can be
 accessed as *p1* to *pN*, where :math:`N` is defined by the pair style.

doc/src/pair_granular.rst

@@ -736,7 +736,7 @@ or
 
 .. math::
 
-   E_{eff,ij} = \frac{E_{ij}}{2(1-\nu_{ij})}
+   E_{eff,ij} = \frac{E_{ij}}{2(1-\nu_{ij}^2)}
 
 These pair styles write their information to :doc:`binary restart files <restart>`, so a pair_style command does not need to be
 specified in an input script that reads a restart file.

src/BPM/bond_bpm.cpp

@@ -187,10 +187,12 @@ void BondBPM::settings(int narg, char **arg)
         iarg++;
       }
     } else if (strcmp(arg[iarg], "overlay/pair") == 0) {
-      overlay_flag = 1;
+      if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for overlay/pair");
+      overlay_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
       iarg++;
-    } else if (strcmp(arg[iarg], "break/no") == 0) {
-      break_flag = 0;
+    } else if (strcmp(arg[iarg], "break") == 0) {
+      if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for break");
+      break_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
       iarg++;
     } else {
       leftover_iarg.push_back(iarg);

src/BPM/bond_bpm_rotational.cpp

@@ -50,6 +50,7 @@ BondBPMRotational::BondBPMRotational(LAMMPS *_lmp) :
 {
   partial_flag = 1;
   smooth_flag = 1;
+  normalize_flag = 0;
 
   single_extra = 7;
   svector = new double[7];
@@ -203,6 +204,13 @@ double BondBPMRotational::elastic_forces(int i1, int i2, int type, double r_mag,
   double Ts[3], Tb[3], Tt[3], Tbp[3], Ttp[3], Tsp[3], T_rot[3], Ttmp[3];
 
   double **quat = atom->quat;
+  double r0_mag_inv = 1.0 / r0_mag;
+  double Kr_type = Kr[type];
+  double Ks_type = Ks[type];
+  if (normalize_flag) {
+    Kr_type *= r0_mag_inv;
+    Ks_type *= r0_mag_inv;
+  }
 
   q1[0] = quat[i1][0];
   q1[1] = quat[i1][1];
@@ -217,26 +225,26 @@ double BondBPMRotational::elastic_forces(int i1, int i2, int type, double r_mag,
   // Calculate normal forces, rb = bond vector in particle 1's frame
   MathExtra::qconjugate(q2, q2inv);
   MathExtra::quatrotvec(q2inv, r, rb);
-  Fr = Kr[type] * (r_mag - r0_mag);
+  Fr = Kr_type * (r_mag - r0_mag);
 
   MathExtra::scale3(Fr * r_mag_inv, rb, F_rot);
 
   // Calculate forces due to tangential displacements (no rotation)
   r0_dot_rb = MathExtra::dot3(r0, rb);
-  c = r0_dot_rb * r_mag_inv / r0_mag;
+  c = r0_dot_rb * r_mag_inv * r0_mag_inv;
   gamma = acos_limit(c);
 
   MathExtra::cross3(rb, r0, rb_x_r0);
   MathExtra::cross3(rb, rb_x_r0, s);
   MathExtra::norm3(s);
 
-  MathExtra::scale3(Ks[type] * r_mag * gamma, s, Fs);
+  MathExtra::scale3(Ks_type * r_mag * gamma, s, Fs);
 
   // Calculate torque due to tangential displacements
   MathExtra::cross3(r0, rb, t);
   MathExtra::norm3(t);
 
-  MathExtra::scale3(0.5 * r_mag * Ks[type] * r_mag * gamma, t, Ts);
+  MathExtra::scale3(0.5 * r_mag * Ks_type * r_mag * gamma, t, Ts);
 
   // Relative rotation force/torque
   // Use representation of X'Y'Z' rotations from Wang, Mora 2009
@@ -316,12 +324,12 @@ double BondBPMRotational::elastic_forces(int i1, int i2, int type, double r_mag,
   Ttp[1] = 0.0;
   Ttp[2] = Kt[type] * psi;
 
-  Fsp[0] = -0.5 * Ks[type] * r_mag * theta * cos_phi;
-  Fsp[1] = -0.5 * Ks[type] * r_mag * theta * sin_phi;
+  Fsp[0] = -0.5 * Ks_type * r_mag * theta * cos_phi;
+  Fsp[1] = -0.5 * Ks_type * r_mag * theta * sin_phi;
   Fsp[2] = 0.0;
 
-  Tsp[0] = 0.25 * Ks[type] * r_mag * r_mag * theta * sin_phi;
-  Tsp[1] = -0.25 * Ks[type] * r_mag * r_mag * theta * cos_phi;
+  Tsp[0] = 0.25 * Ks_type * r_mag * r_mag * theta * sin_phi;
+  Tsp[1] = -0.25 * Ks_type * r_mag * r_mag * theta * cos_phi;
   Tsp[2] = 0.0;
 
   // Rotate forces/torques back to 1st particle's frame
@@ -667,6 +675,10 @@ void BondBPMRotational::settings(int narg, char **arg)
       if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for smooth");
       smooth_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
       i += 1;
+    } else if (strcmp(arg[iarg], "normalize") == 0) {
+      if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for normalize");
+      normalize_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
+      i += 1;
     } else {
       error->all(FLERR, "Illegal bond bpm command, invalid argument {}", arg[iarg]);
     }
@@ -793,7 +805,7 @@ double BondBPMRotational::single(int type, double rsq, int i, int j, double &ffo
   double breaking = elastic_forces(i, j, type, r_mag, r0_mag, r_mag_inv, rhat, r, r0, force1on2,
                                    torque1on2, torque2on1);
   damping_forces(i, j, type, rhat, r, force1on2, torque1on2, torque2on1);
-  fforce = MathExtra::dot3(force1on2, r);
+  fforce = MathExtra::dot3(force1on2, rhat);
   fforce *= -1;
 
   double smooth = 1.0;

src/BPM/bond_bpm_rotational.h

@@ -41,7 +41,7 @@ class BondBPMRotational : public BondBPM {
  protected:
   double *Kr, *Ks, *Kt, *Kb, *gnorm, *gslide, *groll, *gtwist;
   double *Fcr, *Fcs, *Tct, *Tcb;
-  int smooth_flag;
+  int smooth_flag, normalize_flag;
 
   double elastic_forces(int, int, int, double, double, double, double *, double *, double *,
                         double *, double *, double *);

src/BPM/bond_bpm_spring.cpp

@@ -37,6 +37,7 @@ BondBPMSpring::BondBPMSpring(LAMMPS *_lmp) :
 {
   partial_flag = 1;
   smooth_flag = 1;
+  normalize_flag = 0;
 
   single_extra = 1;
   svector = new double[1];
@@ -190,7 +191,10 @@ void BondBPMSpring::compute(int eflag, int vflag)
     }
 
     rinv = 1.0 / r;
-    fbond = k[type] * (r0 - r);
+    if (normalize_flag)
+      fbond = -k[type] * e;
+    else
+      fbond = k[type] * (r0 - r);
 
     delvx = v[i1][0] - v[i2][0];
     delvy = v[i1][1] - v[i2][1];
@@ -302,6 +306,10 @@ void BondBPMSpring::settings(int narg, char **arg)
       if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for smooth");
       smooth_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
       i += 1;
+    } else if (strcmp(arg[iarg], "normalize") == 0) {
+      if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command, missing option for normalize");
+      normalize_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
+      i += 1;
     } else {
       error->all(FLERR, "Illegal bond bpm command, invalid argument {}", arg[iarg]);
     }
@@ -376,7 +384,11 @@ double BondBPMSpring::single(int type, double rsq, int i, int j, double &fforce)
 
   double r = sqrt(rsq);
   double rinv = 1.0 / r;
-  fforce = k[type] * (r0 - r);
+
+  if (normalize_flag)
+    fforce = k[type] * (r0 - r) / r0;
+  else
+    fforce = k[type] * (r0 - r);
 
   double **x = atom->x;
   double **v = atom->v;

src/BPM/bond_bpm_spring.h

@@ -40,7 +40,7 @@ class BondBPMSpring : public BondBPM {
 
  protected:
   double *k, *ecrit, *gamma;
-  int smooth_flag;
+  int smooth_flag, normalize_flag;
 
   void allocate();
   void store_data();