Fix shake/rattle first timestep is unconstrained

Good Morning,

I’ve been reading through the SHAKE/RATTLE source and the forum threads from its original implementation and was wondering if the first position and velocity update are left unconstrained? It seems these algorithms act in-between timesteps so the step from time 0 to time 1 would be unconstrained. I doubt this is a big problem just trying to make sure I understand what’s going on.

Also is there any interest in SETTLE being added in place of shake3() and rattle3()?


It is not completely unconstrained, but it is not fully working for algorithmic reasons either.
The issue is that the SHAKE algorithm is derived for regular Verlet time stepping where you have the positions at the current and the previous step and velocities are “hidden”. LAMMPS, however, uses velocity Verlet time stepping using positions and velocities. To make it work regardless, you have to store the positions for the previous step, but in the initial setup in LAMMPS you can only recover unconstrained old positions, thus the first step is rather inaccurate and constraints are not so well maintained.

To have more choices is always a good thing, for as long as it is done in a sufficiently clean way and does not negatively impact the implementation of the existing features by imposing limits or making the code harder to read than it already is.

1 Like