Timestep with atom_style body and fix nve/body

[He_Huang]

Dear LAMMPS Users,

Does anyone has experience using the BODY module for running rigid body dynamics? I am aware that the fix rigid command is often the more mature and goto method in LAMMPS for rigid body dynamics. I am interested in using the atom_style body because I’m implementing a new pair style which treats a rigid body as a single particle.

I have run a couple simple tests of this module with a system consisting of a few hundred water shaped rigid-body particles with the prototypic body pair_style and the nve/body fix. It seems that the integrator is not quite stable at 1 fs timestep. Does anyone know if rigid body dynamics in general requires a smaller timestep or is this specific to the integration algorithm implemented for the BODY module? I checked the manual and code and do know that the ‘fix nve/body’ uses a different integration algorithm as compared to ‘fix rigid/nve'. The time cost of the planned simulation is high and therefore a longer timestep means a lot to me. If the allowed timestep lengths between the two integrators do differ by a significant factor, I probably should try transplanting the more stable one into the BODY framework.

Any advices would be appreciated.

Best regards,
He

[akohlmey]

Dear LAMMPS Users,

Does anyone has experience using the BODY module for running rigid body
dynamics? I am aware that the fix rigid command is often the more mature
and goto method in LAMMPS for rigid body dynamics. I am interested in using
the atom_style body because I’m implementing a new pair style which treats
a rigid body as a single particle.

I have run a couple simple tests of this module with a system consisting
of a few hundred water shaped rigid-body particles with the prototypic body
pair_style and the nve/body fix. It seems that the integrator is not quite
stable at 1 fs timestep. Does anyone know if rigid body dynamics in general
requires a smaller timestep or is this specific to the integration
algorithm implemented for the BODY module?

​the only way to reliably integrate rigid water molecules with a time step
of 1fs or more is via SHAKE (2fs is commonly used at room ​temperature).
but that isn't really a rigid body integration. instead you are moving
point particles with holonomic contraints applied, which remove the fast
O-H vibrations, that would otherwise limit the possible time step. with
rigid body integrators like fix rigid (and i do mean *any* variant of
them), you more likely need to stick timesteps of under 0.5fs, possibly
even under 0.25fs.

the reason for this is fundamental: the stability of numerical integration
of the equation of motions depends on time step and the velocities, the
latter of which in turn depends on temperature and mass: if you increase
velocity, either through raising the temperature or decreasing the mass,
you have to reduce your time step. correspondingly, for integrating
rotations you need to look at the corresponding moments of inertia instead
of masses. if you look at the near-prolate geometry of a water molecule, it
should be obvious that this will curb the allowed time step. just think of
a spinning figure skater: through pulling in the arms, the moment of
intertia is reduced and due to energy conservation, the rotational speed
increases.

I checked the manual and code and do know that the ‘fix nve/body’ uses a
different integration algorithm as compared to ‘fix rigid/nve'.

​not really. ​fix rigid/nve use a slightly different formulation, since it
is implemented to support a nose-hoover chain thermostat on the rigid
bodies and a nose-hoover barostat. but the crucial steps and time step
related properties are the same as in fix rigid and fix rigid/body. the
actual time integration algorithm employed, velocity verlet, is determined
by the selected integrate style, not the rigid fix.

The time cost of the planned simulation is high and therefore a longer

timestep means a lot to me. If the allowed timestep lengths between the two
integrators do differ by a significant factor, I probably should try
transplanting the more stable one into the BODY framework.

​you could easily test this yourself.

axel.

p.s.: thanks for taking the effort to have your inquiry formulated in a
much more comprehensible way. it is much appreciated.

[He_Huang]

Hi Alex,

Thanks for taking the time to write the detailed answer. The point that I couldn’t quite get is what makes constrained dynamics which is numerically approximated by SHAKE/RATTLE/LINCS differ from the real rigid body equations of motion. The fundamental reason you gave suggests that these are different physical processes. However, in my understanding, they are mathematically equivalent in the continuous scenario. In other words, integrating EOM for individual particles plus applying constrains or directly integrating EOM for the body are two numerically different ways to approach the same physical trajectory. If this understanding is correct, starting from the same point in phase space, exact integration of both schemes should land on the same place after the same time interval, meaning all DOFs in the system should make the same amount of displacement. Any discrepancy should be due to numerical instead of physical reasons.

As to the integration algorithm with rigid and body, what I found in the code is that fix nve/body uses a Richardson iteration to update the quaternion after the angular momentum is updated by the torque. According to the documentation, fix rigid/nve uses Miller’s nosquish algorithm for integrating quaternion dynamics. That algorithm separates the rotation into two parts, one caused by torque and one free of torque. The latter part is integrated with a finer discretization within each timestep. I haven’t fully grasped the details of the algorithm but am wondering if that is sufficient for a significant boost of timestep. You’re absolutely right that I should simply run a test. I haven’t checked any fix with thermostat but NVE is sufficient for us for now.

Thanks again for your kind answer!

He

P.S., you nailed it that the previous question was also from us. It was sent by a student in the group who just started working on the project. Thanks for being patient.

[sjplimp]

The stable timestep for a rigid body depends on its mass and moments

of inertia (for rotation), and the forces you are computing, etc. There

is no single good answer to that Q. In general the rigid/nve formalism

is better than just the NVE rigid and nve/body. As I recall, the latter

two are essentially the same.

Steve

[akohlmey]

Hi Alex,

Thanks for taking the time to write the detailed answer. The point that I
couldn’t quite get is what makes constrained dynamics which is numerically
approximated by SHAKE/RATTLE/LINCS differ from the real rigid body
equations of motion. The fundamental reason you gave suggests that these
are different physical processes. However, in my understanding, they are
mathematically equivalent in the continuous scenario. In other words,
integrating EOM for individual particles plus applying constrains or
directly integrating EOM for the body are two numerically different ways to
approach the same physical trajectory. If this understanding is correct,
starting from the same point in phase space, exact integration of both
schemes should land on the same place after the same time interval, meaning
all DOFs in the system should make the same amount of displacement. Any
discrepancy should be due to numerical instead of physical reasons.

​SHAKE and related algorithms work well for "local" rigidity, e.g. keep all
bonds to hydrogens in a (large) molecule at fixed length, or keep small
molecules rigid, as that keeps it relatively easy and efficient to
(iteratively and approximately) solve the system of constraint equations,
especially when running in parallel. so it is only an approximation to a
rigid body integrator.

with SHAKE (or SHAKE+RATTLE) you first do an unconstrained position (and
velocity) update and then apply the constraints, i.e. modify the
position/velocity update, so that the constrains are satisfied to the
requested degree. that basically makes the time scale for updates, the time
scale of non-bonded position updates. with typical MD systems, you have
quite soft potentials, which can tolerate larger time steps without too
much of an error. however, it is crucial, that the constraint equations
converge well, which requires them to have near diagonal coefficient
matrix. also, please keep in mind that you approximate a system non-linear
equations with linear equations.

in rigid body integrators, you do not do position updates for individual
atoms, but separate the center of mass and rotation updates, so now the
time scales for each of those matter.

As to the integration algorithm with rigid and body, what I found in the
code is that fix nve/body uses a Richardson iteration to update the
quaternion after the angular momentum is updated by the torque. According
to the documentation, fix rigid/nve uses Miller’s nosquish algorithm for
integrating quaternion dynamics. That algorithm separates the rotation into
two parts, one caused by torque and one free of torque. The latter part is
integrated with a finer discretization within each timestep. I haven’t
fully grasped the details of the algorithm but am wondering if that is
sufficient for a significant boost of timestep. You’re absolutely right
that I should simply run a test. I haven’t checked any fix with thermostat
but NVE is sufficient for us for now.

yeah, some kind of multi-timestepping might be helpful. not sure, though,
how far you can get with typical molecular systems assuming the torque
doesn't change (much). it may be worth considering/testing distance based
multi-timestepping via run-style respa. depending on your system, it may be
quite reasonable to assume that the interactions affecting torque are short
ranged and thus only a subset of interactions will need to be recomputed at
every time step.

keep in mind, that the computational efficiency of your simulations is the
product of how long a timestep you can afford and how costly the
computation of interactions is. there are some examples in the literature,
that may not be intuitive, e.g. it may be computationally more efficient to
represent a complex "shaped" object by a superposition of - extremely fast
to compute - point particles rather than having to compute relative
orientations and surface shape dependent interactions.

axel.