Extreme Twisting & Plectoneme Formation

Based on Lazarus et al. (2013)

A heavy elastic rod is axially compressed until it buckles under gravity, then twisted up to 9 full rotations. The simulation captures writhing instabilities and plectoneme formation with self-contact via Hertz–Mindlin interactions. Two cases are considered: an initially straight rod and a naturally curved rod.

Parameters

Parameter	Value
Rod length \(L\)	300 mm
Rod diameter \(d\)	3.1 mm
Density \(\rho\)	1.2 g/cm\(^3\)
Young's modulus \(E\)	1.3 MPa
Shear modulus \(G\)	433 kPa
Number of bonds \(N\)	300 (bond length \(\ell_0 \approx 1\) mm)
Curvature (curved rod)	\(\kappa = 44.84\) m\(^{-1}\) (2.14 coils)

Setup

Both cases share the same material properties and boundary conditions. The left end is clamped (position and orientation fixed); the right end is displaced and rotated using fix move.

Straight rod (twist-straight.in)

1. Compress: the right end is displaced by 80 mm at a constant velocity.
2. Relax: the rod equilibrates under Langevin dynamics with gravity.
3. Twist: the right end is rotated up to 9 full turns at a constant rate using fix move rotate.

Curved rod (twist-curved.in)

The rod is initialized in its stress-free coiled configuration and bonds are created in this state. An additional stage precedes the compression:

1. Straighten: a pure bending moment \(M = \kappa EI\) is applied at the free end via fix addtorque/atom to uncoil the rod. Motion is constrained to 2D using fix setforce NULL NULL 0. Contact and gravity are disabled during this stage.
2. Compress: similar to the straight rod.
3. Relax: similar to the straight rod.
4. Twist: similar to the straight rod, with an additional frictionless wall at the rotating end (see below).

Atom types and pair style. The left and right end atoms are assigned different atom types (set group left type 2, set group right type 3). A pair_style hybrid zero ... granular is used so that during straightening (Stage 0), no granular contact is active — the coiled rod overlaps itself and contact would interfere with the uncoiling. The pair_coeff 2 3 granular ... assignment exists only to satisfy LAMMPS's requirement that every sub-style in a hybrid pair style is referenced by at least one pair_coeff; types 2 and 3 are at opposite ends of the rod and never contact each other. After straightening, type 1–1 granular contact is enabled for the compression and twist stages.

Wall at the rotating end. During twisting, the plectoneme loop can migrate toward the rotating end and slip off, causing the rod to untwist. A frictionless wall (fix wall/gran ... xplane 0 220) is placed near the rotating end to prevent this. The straight rod does not need this wall because its plectoneme forms farther from the end.

Key LAMMPS settings

• Units: units lj with milli units (see Unit System)
• Atom style: atom_style hybrid bpm/sphere dipole — a "dipole" is assigned to each atom to visualize the twist. Note that the dipole vector does not capture the full 3D rotation.
• Contact: pair_style granular with Hertz–Mindlin; atom diameter set to \(\approx 3 \ell_0\) for smooth contact surface; special_bonds lj 0 0 0 disables pair forces for bonded neighbors up to 3 bonds away. This is necessary because the atom diameter causes bonded neighbors to overlap.
• Time integration: fix nve/bpm/sphere update dipole
• Damping: fix langevin for quasi-static equilibration

Note

The atom diameter is independent of the rod diameter \(d\). The latter is only used to compute the BPM parameters.

Files

Data file generator

create_straight_rod.py ⬇ Download

import numpy as np

N = 300  # number of bonds

# Units: (mm, mg, ms)

# Physical parameters
L = 300.
d = 3.1
rho = 1.2

# Discretization
x = np.linspace(0, L, N+1)

dx = L/N
A = np.pi * d**2 / 4
mass = rho * dx * A

diam = 3 * dx # atom diameter for contact != rod diameter
m = mass * np.ones(N+1)
m[0] = m[-1] = 0.5 * mass
dens = m / (np.pi * diam**3 / 6)

with open('rod-N{}.lam'.format(int(N)), 'w') as file:
    file.write('LAMMPS DATA FILE \n\n')

    file.write(str(N+1) + ' atoms \n')
    file.write(str(N) + ' bonds \n\n')

    file.write('1 atom types \n')
    file.write('1 bond types \n\n')

    # Create box
    file.write('%f %f xlo xhi \n'%(-0.1*L, 1.1*L))
    file.write('%f %f ylo yhi \n'%(-0.6*L, 0.6*L))
    file.write('%f %f zlo zhi \n\n'%(-0.6*L, 0.6*L))

    # atom_style hybrid bpm/sphere dipole
    file.write('Atoms \n\n')
    for k in range(N+1):
        # atom-id atom-type x y z molid diameter density q mux muy muz
        file.write('%i 1 %f 0 0 1 %f %f 0 0 1 0\n'%(k+1, x[k], diam, dens[k]))
    file.write('\n')

    file.write('Bonds \n\n')
    for k in range(N):
        file.write('%i %i %i %i \n'%(k+1, 1, k+1, k+2))
    file.write('\n')

create_curved_rod.py ⬇ Download

import numpy as np

N = 300  # number of bonds

# Units: (mm, mg, ms)

# Physical parameters
L = 300.
d = 3.1
rho = 1.2

# Curvature
R = 1/44.84 * 1000
theta = L/R

# Discretization
t = np.linspace(0, theta, N+1)
x = R * np.cos(t-np.pi/2)
y = R * (1 + np.sin(t-np.pi/2))

dx = L/N
A = np.pi * d**2 / 4
mass = rho * dx * A

diam = 3 * dx # atom diameter for contact != rod diameter
m = mass * np.ones(N+1)
m[0] = m[-1] = 0.5 * mass
dens = m / (np.pi * diam**3 / 6)

with open('curved-rod-N{}.lam'.format(int(N)), 'w') as file:
    file.write('LAMMPS DATA FILE \n\n')

    file.write(str(N+1) + ' atoms \n')
    file.write(str(N) + ' bonds \n\n')

    file.write('3 atom types \n')
    file.write('1 bond types \n\n')

    # Create box
    file.write('%f %f xlo xhi \n'%(-1.1*L, 1.1*L))
    file.write('%f %f ylo yhi \n'%(-1.1*L, 1.1*L))
    file.write('%f %f zlo zhi \n\n'%(-0.6*L, 0.6*L))

    # atom_style hybrid bpm/sphere dipole
    file.write('Atoms \n\n')
    for k in range(N+1):
        # atom-id atom-type x y z molid diameter density q mux muy muz
        file.write('%i 1 %f %f 0 1 %f %f 0 0 1 0\n'%(k+1, x[k], y[k], diam, dens[k]))
    file.write('\n')

    file.write('Bonds \n\n')
    for k in range(N):
        file.write('%i %i %i %i \n'%(k+1, 1, k+1, k+2))
    file.write('\n')

LAMMPS input scripts

twist-straight.in ⬇ Download

# =========================================================================
# Twisting a pre-buckled straight rod
# =========================================================================

variable    N       equal 300

# Beam geometry (mm)
variable    L       equal 300
variable    l0      equal $(v_L/v_N)
variable    d       equal 3.1

# Cross-section properties
variable    A       equal (PI/4)*v_d^2
variable    I       equal (PI/64)*v_d^4
variable    J       equal (PI/32)*v_d^4

# Elastic moduli (kPa)
variable    E       equal 1300
variable    G       equal 433

# =========================================================================
# Initialization
# =========================================================================

units           lj
dimension       3
boundary        f f f
processors      * 1 1

atom_style      hybrid bpm/sphere dipole
special_bonds   lj 0.0 0.0 0.0 coul 1.0 1.0 1.0
newton          on off

neighbor        $(5*v_l0) bin
comm_modify     mode single cutoff $(5*v_l0) vel yes

read_data       rod-N${N}.lam

# =========================================================================
# Groups
# =========================================================================

group           l1 id 1
group           l2 id 2
group           r2 id $(v_N)
group           r1 id $(v_N+1)

group           left  union l1 l2
group           right union r1 r2
group           ends  union l1 r1
group           free  subtract all ends
group           mobile union free right

# =========================================================================
# BPM bond style
# =========================================================================

variable        Kr equal v_E*v_A/v_l0
variable        Ks equal 12*v_E*v_I/v_l0^3
variable        Kt equal v_G*v_J/v_l0
variable        Kb equal v_E*v_I/v_l0

bond_style      bpm/rotational smooth no break no overlay/pair no
bond_coeff      1 ${Kr} ${Ks} ${Kt} ${Kb} &
                -1 -1 -1 -1 &
                1e-4 1e-4 1e-4 1e-4

# =========================================================================
# Granular contact (Hertz-Mindlin)
# =========================================================================

# Made up values
variable        YoungsModulus    equal 10*v_E
variable        PoissonsRatio    equal 0.4
variable        FrictionCoeff    equal 0.1
variable        CoeffRestitution equal 0.7

pair_style      granular limit_damping yes
pair_coeff      * * hertz/material $(v_YoungsModulus) $(v_CoeffRestitution) $(v_PoissonsRatio) &
                tangential mindlin/force NULL 0.5 $(v_FrictionCoeff) &
                damping coeff_restitution

# =========================================================================
# Computes and output
# =========================================================================

compute         mytemp all temp/sphere
compute         kerot  all erotate/sphere

# Note: variable torque and compute torque/chunk only works for group of point particles
compute         ft l1 property/atom fx fy fz tqx tqy tqz
compute         f  l1 reduce sum c_ft[1] c_ft[2] c_ft[3]
compute         t  l1 reduce sum c_ft[4] c_ft[5] c_ft[6]

variable        fx equal c_f[1]
variable        tx equal c_t[1]

compute         xmax    all reduce max x
compute         xmin    all reduce min x
variable        dx      equal "c_xmax-c_xmin"
compute         quat    all property/atom quatw quati quatj quatk

compute         ymaxabs all reduce maxabs y
compute         zmaxabs all reduce maxabs z

variable        cycles equal 0

thermo          10000
thermo_style    custom step ke c_kerot pe v_cycles v_dx c_ymaxabs c_zmaxabs v_fx v_tx
thermo_modify   temp mytemp norm no

# =========================================================================
# Stage 1: Axial compression (80 mm)
# =========================================================================

dump            1 all custom 100000 ./compress-N${N}.dump &
                id diameter x y z mux muy muz fx fy fz tqx tqy tqz c_quat[*]

variable        timestep equal 1e-3
timestep        ${timestep}

variable        nsteps equal 5e6
variable        period equal v_timestep*v_nsteps
variable        vx     equal -80/(v_timestep*v_nsteps)

fix             1 free nve/bpm/sphere update dipole
fix             2 r1 move linear ${vx} 0 0
fix             3 free gravity 9.81e-3 vector 0 -1 0
fix             4 free langevin 0 0 1e3 12345 omega yes

run             $(v_nsteps)

# =========================================================================
# Stage 2: Relaxation
# =========================================================================

unfix           2
fix             4 free langevin 0 0 1e3 12345 omega yes

variable        nsteps equal 1e7
run             $(v_nsteps)

# write_restart   compressed-beam.restart

undump          1
unfix           1
unfix           3
unfix           4

# =========================================================================
# Stage 3: Twist (9 full rotations)
# =========================================================================

variable        timestep equal 1e-5
timestep        ${timestep}
reset_timestep  0

balance         1.2 shift x 10 1.1

dump            1 all custom 10000000 ./twist-N${N}.dump &
                id diameter x y z mux muy muz fx fy fz tqx tqy tqz c_quat[*]

variable        nturns       equal 9
variable        stepsperturn equal 1e9
variable        period       equal v_timestep*v_stepsperturn
variable        cycles       equal step/v_stepsperturn
variable        theta        equal 2*PI*v_cycles

bond_coeff      1 ${Kr} ${Ks} ${Kt} ${Kb} &
                -1 -1 -1 -1 &
                1e-2 1e-2 1e-2 1e-2

fix             1 free nve/bpm/sphere update dipole
fix             2 r1 move rotate 0 0 0 -1 0 0 $(v_period) update dipole
fix             3 free gravity 9.81e-3 vector 0 -1 0
fix             4 free langevin 0 0 1e2 12345 omega yes

fix             5 all print 1000000 &
                "${theta} $(c_ymaxabs) $(c_zmaxabs) ${tx}" &
                file twist-N${N}-theta-yzmaxabs-Tx.txt screen no &
                title "theta [rad], y_maxabs [mm], z_maxabs [mm], Tx [uJ]"

variable        i loop ${nturns}
label           loopstart
run             $(v_stepsperturn)
next            i
jump            SELF loopstart

twist-curved.in ⬇ Download

# =========================================================================
# Twisting a pre-buckled naturally curved rod
# =========================================================================

variable    N       equal 300

# Beam geometry (mm)
variable    R       equal 1000/44.84
variable    L       equal 300
variable    l0      equal $(v_L/v_N)
variable    d       equal 3.1

# Cross-section properties
variable    A       equal (PI/4)*v_d^2
variable    I       equal (PI/64)*v_d^4
variable    J       equal (PI/32)*v_d^4

# Elastic moduli (kPa)
variable    E       equal 1300
variable    G       equal 433

# =========================================================================
# Initialization
# =========================================================================

units           lj
dimension       3
boundary        f f f
processors      * 1 1

atom_style      hybrid bpm/sphere dipole
special_bonds   lj 0.0 0.0 0.0 coul 1.0 1.0 1.0
newton          on off

neighbor        $(5*v_l0) bin
comm_modify     mode single cutoff $(5*v_l0) vel yes

read_data       curved-rod-N${N}.lam

# =========================================================================
# Groups
# =========================================================================

group           l1 id 1
group           l2 id 2
group           r2 id $(v_N)
group           r1 id $(v_N+1)

group           left  union l1 l2
group           right union r1 r2
group           ends  union l1 r1
group           free  subtract all ends
group           mobile union free right

set             group left type 2
set             group right type 3

# =========================================================================
# BPM bond style
# =========================================================================

variable        Kr equal v_E*v_A/v_l0
variable        Ks equal 12*v_E*v_I/v_l0^3
variable        Kt equal v_G*v_J/v_l0
variable        Kb equal v_E*v_I/v_l0

bond_style      bpm/rotational smooth no break no overlay/pair no
bond_coeff      1 ${Kr} ${Ks} ${Kt} ${Kb} &
                -1 -1 -1 -1 &
                1e-2 1e-2 1e-2 1e-2

# =========================================================================
# Granular contact (Hertz-Mindlin)
# =========================================================================

# Made up values
variable        YoungsModulus    equal 10*v_E
variable        PoissonsRatio    equal 0.4
variable        FrictionCoeff    equal 0.1
variable        CoeffRestitution equal 0.7

pair_style      hybrid zero $(5*v_l0) granular limit_damping yes
pair_coeff      * * zero
pair_coeff      2 3 granular hertz/material $(v_YoungsModulus) $(v_CoeffRestitution) $(v_PoissonsRatio) &
                tangential mindlin/force NULL 0.5 $(v_FrictionCoeff) &
                damping coeff_restitution

# =========================================================================
# Computes and output
# =========================================================================

compute         mytemp  all temp/sphere
compute         kerot   all erotate/sphere

# Note: variable torque and compute torque/chunk only works for group of point particles
compute         ft l1 property/atom fx fy fz tqx tqy tqz
compute         f  l1 reduce sum c_ft[1] c_ft[2] c_ft[3]
compute         t  l1 reduce sum c_ft[4] c_ft[5] c_ft[6]

variable        fx equal c_f[1]
variable        tx equal c_t[1]

compute         xmax    all reduce max x
compute         xmin    all reduce min x
variable        dx      equal "c_xmax-c_xmin"
compute         quat    all property/atom quatw quati quatj quatk

compute         ymaxabs all reduce maxabs y
compute         zmaxabs all reduce maxabs z

variable        cycles equal 0

thermo          10000
thermo_style    custom step ke c_kerot pe v_cycles v_dx c_ymaxabs c_zmaxabs v_fx v_tx
thermo_modify   temp mytemp norm no

# =========================================================================
# Stage 0: Straighten by applying end moment M = kappa * EI
#          (no gravity, no contact)
# =========================================================================

dump            1 all custom 100000 ./bend-N${N}.dump &
                id diameter x y z mux muy muz fx fy fz tqx tqy tqz c_quat[*]

variable        timestep equal 1e-3
timestep        ${timestep}

variable        nsteps equal 5e6
variable        Mz     equal v_E*v_I/v_R
variable        mz     equal -v_Mz*ramp(0,1)

fix             1 mobile nve/bpm/sphere update dipole
fix             2 r1 addtorque/atom 0 0 v_mz
fix             3 free gravity 0 vector 0 -1 0
fix             4 free langevin 0 0 1e3 12345 omega yes
fix             5 all setforce NULL NULL 0

run             $(v_nsteps)

variable        mz equal -v_Mz
fix             4 free langevin 0 0 1e2 12345 omega yes

run             $(14*v_nsteps)

unfix           1
unfix           2
unfix           3
unfix           4
unfix           5
undump          1

set             group all dipole 0 1 0

# write_restart   straightened-beam.restart

# =========================================================================
# Stage 1: Axial compression (80 mm)
# =========================================================================

balance         1.2 shift x 10 1.1

pair_coeff      1 1 granular hertz/material $(v_YoungsModulus) $(v_CoeffRestitution) $(v_PoissonsRatio) &
                tangential mindlin/force NULL 0.5 $(v_FrictionCoeff) &
                damping coeff_restitution

dump            1 all custom 100000 ./compress-N${N}.dump &
                id diameter x y z mux muy muz fx fy fz tqx tqy tqz c_quat[*]

variable        timestep equal 1e-5
timestep        ${timestep}

variable        nsteps equal 5e7
variable        period equal v_timestep*v_nsteps
variable        vx     equal -80/(v_timestep*v_nsteps)

fix             1 free nve/bpm/sphere update dipole
fix             2 r1 move linear ${vx} 0 0
fix             3 free gravity 9.81e-3 vector 0 -1 0
fix             4 free langevin 0 0 1e2 12345 omega yes

run             $(v_nsteps)

# =========================================================================
# Stage 2: Relaxation
# =========================================================================

unfix           2
fix             4 free langevin 0 0 1e3 12345 omega yes

run             $(4*v_nsteps)

fix             4 free langevin 0 0 1e2 12345 omega yes

run             $(5*v_nsteps)

# write_restart   compressed-beam.restart

undump          1
unfix           1
unfix           3
unfix           4

# =========================================================================
# Stage 3: Twist (9 full rotations, with wall)
# =========================================================================

variable        timestep equal 1e-4
timestep        ${timestep}
reset_timestep  0

balance         1.2 shift x 10 1.1

dump            1 all custom 1000000 ./twist-N${N}.dump &
                id diameter x y z mux muy muz fx fy fz tqx tqy tqz c_quat[*]

variable        nturns       equal 9
variable        stepsperturn equal 1e8
variable        period       equal v_timestep*v_stepsperturn
variable        cycles       equal step/v_stepsperturn
variable        theta        equal 2*PI*v_cycles

fix             1 free nve/bpm/sphere update dipole
fix             2 r1 move rotate 0 0 0 -1 0 0 $(v_period) update dipole
fix             3 free gravity 9.81e-3 vector 0 -1 0
fix             4 free langevin 0 0 1e2 12345 omega yes
fix             5 free wall/gran granular &
                hertz/material $(v_YoungsModulus) $(v_CoeffRestitution) $(v_PoissonsRatio) &
                tangential mindlin/force NULL 0.5 $(v_FrictionCoeff) &
                damping coeff_restitution &
                xplane 0 220

fix             6 all print 100000 &
                "${theta} $(c_ymaxabs) $(c_zmaxabs) ${tx}" &
                file twist-N${N}-theta-yzmaxabs-Tx.txt screen no &
                title "theta [rad], y_maxabs [mm], z_maxabs [mm], Tx [uJ]"

variable        i loop ${nturns}
label           loopstart
run             $(v_stepsperturn)
next            i
jump            SELF loopstart

References

• Lazarus, A., Miller, J., Reis, P.M. (2013). Continuation of equilibria and stability of slender elastic rods using an asymptotic numerical method. J. Mech. Phys. Solids, 61(8), 1712–1736. doi:10.1016/j.jmps.2013.04.002
• Lazarus, A., Miller, J.T., Metlitz, M.M., Reis, P.M. (2013). Contorting a heavy and naturally curved elastic rod. Soft Matter, 9(34), 8274–8281. doi:10.1039/c3sm50873k
• Grange, S., Bertrand, D. (2024). Co-rotational 3D beam element using quaternion algebra. Int. J. Solids Struct., 293, 112746. doi:10.1016/j.ijsolstr.2024.112746

---

← Back to Examples