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.