Bugfix and documentation corrections for BPM+Granular packages

doc/src/bond_bpm_rotational.rst

@@ -10,7 +10,7 @@ Syntax
 
    bond_style bpm/rotational keyword value attribute1 attribute2 ...
 
-* optional keyword = *overlay/pair* or *store/local* or *smooth* or *break/no*
+* optional keyword = *overlay/pair* or *store/local* or *smooth* or *break*
 
   .. parsed-literal::
 
@@ -80,32 +80,32 @@ respectively.  Details on the calculations of shear displacements and
 angular displacements can be found in :ref:`(Wang) <Wang2009>` and
 :ref:`(Wang and Mora) <Wang2009b>`.
 
-Bonds will break under sufficient stress. A breaking criteria is calculated
+Bonds will break under sufficient stress. A breaking criterion is calculated
 
 .. math::
 
-   B = \mathrm{max}\{0, \frac{f_r}{f_{r,c}} + \frac{|f_s|}{f_{s,c}} +
-       \frac{|\tau_b|}{\tau_{b,c}} + \frac{|\tau_t|}{\tau_{t,c}} \}
+   B = \mathrm{max}\left\{0, \frac{f_r}{f_{r,c}} + \frac{|f_s|}{f_{s,c}} +
+       \frac{|\tau_b|}{\tau_{b,c}} + \frac{|\tau_t|}{\tau_{t,c}} \right\}
 
 where :math:`|f_s|` is the magnitude of the shear force and
 :math:`|\tau_b|` and :math:`|\tau_t|` are the magnitudes of the
-bending and twisting forces, respectively. The corresponding variables
+bending and twisting torques, respectively. The corresponding variables
 :math:`f_{r,c}` :math:`f_{s,c}`, :math:`\tau_{b,c}`, and
 :math:`\tau_{t,c}` are critical limits to each force or torque.  If
 :math:`B` is ever equal to or exceeds one, the bond will break.  This
-is done by setting by setting its type to 0 such that forces and
+is done by setting the bond type to 0 such that forces and
 torques are no longer computed.
 
 After computing the base magnitudes of the forces and torques, they
 can be optionally multiplied by an extra factor :math:`w` to smoothly
 interpolate forces and torques to zero as the bond breaks. This term
-is calculated as :math:`w = (1.0 - B^4)`. This smoothing factor can be
-added or removed using the *smooth* keyword.
+is calculated as :math:`w = (1.0 - B^4)`. This smoothing factor can be added
+or removed by setting the *smooth* keyword to *yes* or *no*, respectively.
 
 Finally, additional damping forces and torques are applied to the two
 particles. A force is applied proportional to the difference in the
 normal velocity of particles using a similar construction as
-dissipative particle dynamics (:ref:`(Groot) <Groot3>`):
+dissipative particle dynamics :ref:`(Groot) <Groot3>`:
 
 .. math::
 
@@ -115,8 +115,8 @@ where :math:`\gamma_n` is the damping strength, :math:`\hat{r}` is the
 radial normal vector, and :math:`\vec{v}` is the velocity difference
 between the two particles. Similarly, tangential forces are applied to
 each atom proportional to the relative differences in sliding
-velocities with a constant prefactor :math:`\gamma_s` (:ref:`(Wang et
-al.) <Wang20152>`) along with their associated torques. The rolling and
+velocities with a constant prefactor :math:`\gamma_s` :ref:`(Wang et
+al.) <Wang20152>` along with their associated torques. The rolling and
 twisting components of the relative angular velocities of the two
 atoms are also damped by applying torques with prefactors of
 :math:`\gamma_r` and :math:`\gamma_t`, respectively.
@@ -139,21 +139,23 @@ 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.
+If the *normalize* keyword is set to *yes*, the radial and shear forces
+will be normalized by :math:`r_0` such that :math:`k_r` and :math:`k_s`
+must be given in force units.
 
 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* option. These settings require specific
+Pair forces can alternatively be overlaid on top of bond forces by setting
+the *overlay/pair* keyword to *yes*. 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.
+restrictions.  Further details can be found in the :doc:`how to
+<Howto_bpm>` page on BPMs.
 
 .. versionadded:: 28Mar2023
 
-If the *break* option is used, then LAMMPS assumes bonds should not break
+If the *break* keyword is set to *no*, 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.
+However, if a bond reaches a damage criterion greater than one,
+it will trigger an error.
 
 If the *store/local* keyword is used, an internal fix will track bonds that
 break during the simulation. Whenever a bond breaks, data is processed
@@ -239,9 +241,8 @@ requires setting
 
    special_bonds lj 0 1 1 coul 1 1 1
 
-and :doc:`newton <newton>` must be set to bond off.  If the
-*overlay/pair* option is used, this bond style alternatively requires
-setting
+and :doc:`newton <newton>` must be set to bond off.  If the *overlay/pair*
+keyword is set to *yes*, this bond style alternatively requires setting
 
 .. code-block:: LAMMPS
 

doc/src/bond_bpm_spring.rst

@@ -10,7 +10,7 @@ Syntax
 
    bond_style bpm/spring keyword value attribute1 attribute2 ...
 
-* optional keyword = *overlay/pair* or *store/local* or *smooth* or *break/no*
+* optional keyword = *overlay/pair* or *store/local* or *smooth* or *break*
 
   .. parsed-literal::
 
@@ -72,13 +72,13 @@ particles based on a model described by Clemmer and Robbins
 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
-setting its type to 0 such that forces are no longer computed.
+break at a strain of :math:`\epsilon_c`.  This is done by setting
+the bond type to 0 such that forces are no longer computed.
 
 An additional damping force is applied to the bonded
 particles.  This forces is proportional to the difference in the
 normal velocity of particles using a similar construction as
-dissipative particle dynamics (:ref:`(Groot) <Groot4>`):
+dissipative particle dynamics :ref:`(Groot) <Groot4>`:
 
 .. math::
 
@@ -88,9 +88,10 @@ where :math:`\gamma` is the damping strength, :math:`\hat{r}` is the
 radial normal vector, and :math:`\vec{v}` is the velocity difference
 between the two particles.
 
-The smoothing factor :math:`w` can be added or removed using the
-*smooth* keyword. It is constructed such that forces smoothly go
-to zero, avoiding discontinuities, as bonds approach the critical strain
+The smoothing factor :math:`w` can be added or removed by setting the
+*smooth* keyword to *yes* or *no*, respectively. It is constructed such
+that forces smoothly go to zero, avoiding discontinuities, as bonds
+approach the critical strain
 
 .. math::
 
@@ -105,21 +106,22 @@ 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.
+If the *normalize* keyword is set to *yes*, the elastic bond force will be
+normalized by :math:`r_0` such that :math:`k` must be given in force units.
 
 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* option. These settings require specific
+Pair forces can alternatively be overlaid on top of bond forces by setting
+the *overlay/pair* keyword to *yes*. 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.
+restrictions.  Further details can be found in the :doc:`how to
+<Howto_bpm>` page on BPMs.
 
 .. versionadded:: 28Mar2023
 
-If the *break* option is used, then LAMMPS assumes bonds should not break
+If the *break* keyword is set to *no*, 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.
+However, if a bond reaches a strain greater than :math:`\epsilon_c`,
+it will trigger an error.
 
 If the *store/local* keyword is used, an internal fix will track bonds that
 break during the simulation. Whenever a bond breaks, data is processed
@@ -196,9 +198,8 @@ requires setting
 
    special_bonds lj 0 1 1 coul 1 1 1
 
-and :doc:`newton <newton>` must be set to bond off.  If the
-*overlay/pair* option is used, this bond style alternatively requires
-setting
+and :doc:`newton <newton>` must be set to bond off.  If the *overlay/pair*
+keyword is set to *yes*, this bond style alternatively requires setting
 
 .. code-block:: LAMMPS
 

src/BPM/bond_bpm.cpp

@@ -189,11 +189,11 @@ void BondBPM::settings(int narg, char **arg)
     } else if (strcmp(arg[iarg], "overlay/pair") == 0) {
       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++;
+      iarg += 2;
     } 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++;
+      iarg += 2;
     } else {
       leftover_iarg.push_back(iarg);
       iarg++;

src/GRANULAR/pair_granular.cpp

@@ -767,10 +767,10 @@ double PairGranular::single(int i, int j, int itype, int jtype,
   model->history = history;
 
   model->calculate_forces();
-  double *forces = model->forces;
 
   // apply forces & torques
-  fforce = MathExtra::len3(forces);
+  // Calculate normal component, normalized by r
+  fforce = model->Fnormal * model->rinv;
 
   // set single_extra quantities
   svector[0] = model->fs[0];