ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #3827 from jtclemm/BPM

Browse Source 
Bugfix and documentation corrections for BPM+Granular packages
This commit is contained in:
Axel Kohlmeyer
2023-07-06 00:35:58 -04:00
committed by GitHub
parent c5d9f901d9 838be24af4
commit ad14862cd0
4 changed files with 46 additions and 44 deletions
Download Patch File Download Diff File Expand all files Collapse all files

45
doc/src/bond_bpm_rotational.rst
View File

@ -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

37
doc/src/bond_bpm_spring.rst
View File

@ -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

4
src/BPM/bond_bpm.cpp
View File

@ -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++;

4
src/GRANULAR/pair_granular.cpp
View File

@ -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];