diff --git a/doc/src/bond_bpm_rotational.rst b/doc/src/bond_bpm_rotational.rst index ca12d86ccc..cfacfde2cd 100644 --- a/doc/src/bond_bpm_rotational.rst +++ b/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) ` and :ref:`(Wang and Mora) `. -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) `): +dissipative particle dynamics :ref:`(Groot) `: .. 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.) `) along with their associated torques. The rolling and +velocities with a constant prefactor :math:`\gamma_s` :ref:`(Wang et +al.) ` 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 ` 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 ` settings described in the -restrictions. Further details can be found in the `:doc: how to -` page on BPMs. +restrictions. Further details can be found in the :doc:`how to +` 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 ` must be set to bond off. If the -*overlay/pair* option is used, this bond style alternatively requires -setting +and :doc:`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 diff --git a/doc/src/bond_bpm_spring.rst b/doc/src/bond_bpm_spring.rst index d89035dcad..2376b8e077 100644 --- a/doc/src/bond_bpm_spring.rst +++ b/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) `): +dissipative particle dynamics :ref:`(Groot) `: .. 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 ` settings described in the -restrictions. Further details can be found in the `:doc: how to -` page on BPMs. +restrictions. Further details can be found in the :doc:`how to +` 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 ` must be set to bond off. If the -*overlay/pair* option is used, this bond style alternatively requires -setting +and :doc:`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 diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp index da189ccaf5..d965a19067 100644 --- a/src/BPM/bond_bpm.cpp +++ b/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++; diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 0b6bd537e5..30f272791e 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/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];