Skip to content

Renaming, adding docs + example #8

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
jtclemm wants to merge 4 commits into from
Closed

Renaming, adding docs + example #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
a02bfaf
Possible example
jtclemm Apr 2, 2026
fff3a9b
Overwriting rotational with corotational
jtclemm Apr 2, 2026
a010ea3
Renaming corotational hclass names
jtclemm Apr 2, 2026
47aa957
Doc updates
jtclemm Apr 2, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
130 changes: 93 additions & 37 deletions doc/src/bond_bpm_rotational.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ Syntax

bond_style bpm/rotational keyword value attribute1 attribute2 ...

* optional keyword = *overlay/pair* or *store/local* or *smooth* or *break*
* optional keyword = *overlay/pair* or *store/local* or *smooth* or *break* or *frame* or *damping*

.. parsed-literal::

Expand All @@ -36,6 +36,17 @@ Syntax
*break* value = *yes* or *no*
indicates whether bonds break during a run

*frame* value = *corotational* or *standard*
what frame is used for calculating rotations

*damping* value = *derivative* or *dem*
damping construction

.. note::

Versions of LAMMPS before .. versionadded:: TBD used *frame* standard
and *dem* damping

Examples
""""""""

Expand All @@ -44,6 +55,9 @@ Examples
bond_style bpm/rotational
bond_coeff 1 1.0 0.2 0.02 0.02 0.20 0.04 0.04 0.04 0.1 0.02 0.002 0.002

bond_style bpm/rotational frame standard damping derivative
bond_coeff 1 1.0 0.2 0.02 0.02 0.20 0.04 0.04 0.04 0.1 0.02 0.002 0.002

bond_style bpm/rotational store/local myfix 1000 time id1 id2
dump 1 all local 1000 dump.broken f_myfix[1] f_myfix[2] f_myfix[3]
dump_modify 1 write_header no
Expand All @@ -65,20 +79,23 @@ has a magnitude of

.. math::

f_r = k_r (r - r_0)
f_\mathrm{radial} = k_\mathrm{radial} (r - r_0)

where :math:`k_r` is a stiffness and :math:`r` is the current distance and
:math:`r_0` is the initial distance between the two particles.
where :math:`k_\mathrm{radial}` is a stiffness and :math:`r` is the
current distance and :math:`r_0` is the initial distance between the
two particles.

A tangential force is applied perpendicular to the normal direction
which is proportional to the tangential shear displacement with a
stiffness of :math:`k_s`. This tangential force also induces a torque.
In addition, bending and twisting torques are also applied to
particles which are proportional to angular bending and twisting
displacements with stiffnesses of :math:`k_b` and :math:`k_t`,
respectively. Details on the calculations of shear displacements and
angular displacements can be found in :ref:`(Wang) <Wang2009>` and
:ref:`(Wang and Mora) <Wang2009b>`.
stiffness of :math:`k_\mathrm{shear}`. This tangential force also
induces a torque. In addition, bending and twisting torques are also
applied to particles which are proportional to angular bending and
twisting displacements with stiffnesses of :math:`k_\mathrm{bend}`
and :math:`k_\mathrm{twist}`, respectively. Details on the calculations
of shear displacements and angular displacements can be found in
:ref:`(Wang) <Wang2009>`, :ref:`(Wang and Mora) <Wang2009b>`, and/or
:ref:`(Alkuino et al.) <Alkuino2026>` depending on the *frame*
(discussed below).

Bonds will break under sufficient stress. A breaking criterion is calculated

Expand All @@ -89,8 +106,9 @@ Bonds will break under sufficient stress. A breaking criterion is calculated

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 torques, respectively. The corresponding variables
:math:`f_{r,c}` :math:`f_{s,c}`, :math:`\tau_{b,c}`, and
bending and twisting torques, respectively, and *r*, *s*, *t* and *b*
are short hand for radial, shear, bend, and twist. 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 the bond type to 0 such that forces and
Expand All @@ -109,35 +127,56 @@ dissipative particle dynamics :ref:`(Groot) <Groot3>`:

.. math::

F_D = - \gamma_n w (\hat{r} \bullet \vec{v})
F_D = - \gamma_\mathrm{radial} w (\hat{r} \bullet \vec{v})

where :math:`\gamma_\mathrm{radial}` 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, additional damping forces/torques
are applied to other modes. These details depend on the *damping*
setting.

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
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.
.. versionadded:: TBD

For *damping* style *derivative* (the default), additional forces/torques
are applied on shear, twisting, and bending modes. These are simply
proportional to the rate in change in the shear, bend, and twisting angle,
respectively with prefactors of :math:`\gamma_\mathrm{shear}`,
:math:`\gamma_\mathrm{twist}`, and :math:`\gamma_\mathrm{bend}`. Details
are described in :ref:`(Alkuino et al.) <Alkuino2026>`.

For the *dem* style, forces are applied to each atom proportional to the
relative differences in sliding velocities with a constant prefactor
:math:`\gamma_\mathrm{slide}` :ref:`(Wang et al.) <Wang20152>` along with
the 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_\mathrm{roll}` and :math:`\gamma_\mathrm{twist}`,
respectively. These modes are commonly used in the discrete element method
(DEM) as in :doc:`pair granular <pair_granular>`.

The following coefficients must be defined for each bond type via the
:doc:`bond_coeff <bond_coeff>` command as in the example above, or in
the data file or restart files read by the :doc:`read_data <read_data>`
or :doc:`read_restart <read_restart>` commands:

* :math:`k_r` (force/distance units)
* :math:`k_s` (force/distance units)
* :math:`k_t` (force*distance/radians units)
* :math:`k_b` (force*distance/radians units)
* :math:`f_{r,c}` (force units)
* :math:`f_{s,c}` (force units)
* :math:`\tau_{t,c}` (force*distance units)
* :math:`\tau_{b,c}` (force*distance units)
* :math:`\gamma_n` (force/velocity units)
* :math:`\gamma_s` (force/velocity units)
* :math:`\gamma_r` (force*distance/velocity units)
* :math:`\gamma_t` (force*distance/velocity units)
* :math:`k_\mathrm{radial}` (force/distance units)
* :math:`k_\mathrm{shear}` (force/distance units)
* :math:`k_\mathrm{twist}` (force*distance/radians units)
* :math:`k_\mathrm{bend}` (force*distance/radians units)
* :math:`f_{\mathrm{radial},c}` (force units)
* :math:`f_{\mathrm{shear},c}` (force units)
* :math:`\tau_{\mathrm{twist},c}` (force*distance units)
* :math:`\tau_{\mathrm{bend},c}` (force*distance units)
* :math:`\gamma_\mathrm{radial}` (force/velocity units)
* :math:`\gamma_\mathrm{shear}` (force/velocity units)
* :math:`\gamma_\mathrm{twist}` (force*distance/velocity units)
* :math:`\gamma_\mathrm{bend}` (force*distance/velocity units)

For damping style *dem*, the last three coefficients are replaced with
(note that the position of the twisting coefficient is swapped)

* :math:`\gamma_\mathrm{slide}` (force/velocity units)
* :math:`\gamma_\mathrm{roll}` (force*distance/velocity units)
* :math:`\gamma_\mathrm{twist}` (force*distance/velocity units)

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`
Expand All @@ -160,6 +199,18 @@ maximum strain used by typical non-bpm bond styles. Similar behavior to *break n
can also be attained by setting arbitrarily high values for all four failure
coefficients. One cannot use *break no* with *smooth yes*.

.. versionadded:: TBD

The *frame* setting determines the reference used to calculate relative
rotations. The *standard* option uses the frame of one particle as
described in :ref:`(Wang) <Wang2009>` and :ref:`(Wang and Mora) <Wang2009b>`.
This determination is based on particle ID in LAMMPS. Alternatively,
the *corotational* option (the default) defines a central frame across
the two particles as described in :ref:`(Alkuino et al.) <Alkuino2026>`.
The latter option implies forces do not depend on particle IDs and can be
more stable, particularly in simulations of thin or high distorted
structures such as the wire example in /examples/bpm.

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
and transferred to an internal fix labeled *fix_ID*. This allows the
Expand Down Expand Up @@ -269,7 +320,7 @@ Related commands
Default
"""""""

The option defaults are *overlay/pair* = *no*, *smooth* = *yes*, *normalize* = *no*, and *break* = *yes*
The option defaults are *overlay/pair* = *no*, *smooth* = *yes*, *normalize* = *no*, *break* = *yes*, *frame* = *corotational*, and *damping* = *derivative*

----------

Expand All @@ -280,9 +331,14 @@ p 117-127 (2009).

.. _Wang2009b:

**(Wang and Mora)** Wang, Mora, Advances in Geocomputing,
**(Wang and Mora)** Wang, and Mora, Advances in Geocomputing,
119, p 183-228 (2009).

.. _Alkuino2026:

**(Alkuino)** Alkuino, Clemmer, Santangelo, and Zhang, arXiv:2603.27279
(2026).

.. _Groot3:

**(Groot)** Groot and Warren, J Chem Phys, 107, 4423-35 (1997).
Expand Down
82 changes: 82 additions & 0 deletions examples/bpm/wires/in.bpm.wires
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,82 @@
units lj
dimension 3
boundary m m m
atom_style hybrid bpm/sphere dipole
special_bonds lj 0.0 1.0 1.0 coul 0.0 1.0 1.0
newton on off
comm_modify vel yes cutoff 2.6
region box block -10 60 -10 30 -10 10
create_box 1 box bond/types 1 extra/bond/per/atom 2 extra/special/per/atom 4

variable Na1 equal 200
variable Na2 equal 480

# Build a straight wire
create_atoms 1 random $(v_Na1) 54435234 box group wire1
variable x1 atom id*0.25
set group wire1 x v_x1 y 0.0 z 0.0


# Build a coiled wire
create_atoms 1 random $(v_Na2) 54435234 box group wire2
variable x2 atom (id-$(v_Na1))*0.1
variable y2 atom 2*sin((id-$(v_Na1))*0.25)+20
variable z2 atom 2*cos((id-$(v_Na1))*0.25)
set group wire2 x v_x2 y v_y2 z v_z2

# Define key groups and properties
group end1 id $(v_Na1)
group end2 id $(v_Na1+v_Na2)
group frozen_ends id 1 $(v_Na1+1)
group loaded_ends union end1 end2

set group all quat 1 0 0 0
set group wire1 dipole 0 0 1
set group wire2 dipole 0 0 0

velocity all create 1e-7 1349832 # random noise to break symmetries
velocity frozen_ends set 0.0 0.0 0.0

# Define interactions
pair_style granular
pair_coeff 1 1 hooke 10.0 0.5 tangential linear_history 5.0 0.5 0.1

neighbor 0.3 bin
create_bonds many wire1 wire1 1 0.0 0.3
create_bonds many wire2 wire2 1 0.0 0.9
special_bonds lj 0.0 0.0 0.0 coul 1.0 1.0 1.0

bond_style bpm/rotational smooth no break no
bond_coeff 1 1.0 0.1 0.1 0.1 -1 -1 -1 -1 0.5 0.5 0.05 0.05

# Control simulation
fix 1 all nve/bpm/sphere update dipole
# dipole only used for tracking/visualizing orientation

fix 2 frozen_ends settorque/atom 0.0 0.0 0.0
fix 3 frozen_ends setforce 0.0 0.0 0.0
fix 4 loaded_ends settorque/atom 0.0 0.0 0.0
fix 5 loaded_ends setforce 0.0 0.0 0.0

fix 6 end1 move transrot -0.001 0 0 0 0 0 1 0 0 1600
fix 7 end2 move linear -0.001 0.0 0.0

fix 8 all viscous 5e-5 #small damping to accelerate decay of global modes

timestep 0.1
thermo 100
thermo_style custom step ke pe pxx pyy pzz

compute quat all property/atom quatw quati quatj quatk
compute tq all property/atom tqx tqy tqz

dump atomDump all custom 400 dump id radius xs ys zs c_quat[*] mux muy muz

run 200000

unfix 4
unfix 5
unfix 6
unfix 7

run 300000
Loading