[lammps-users] details of SHAKE implementation

I was just wondering if there is any documentation (other than the code itself) on exactly how SHAKE is implemented. I am guessing that it is the implementation outlined by Bruce Palmer (Palmer, B. J., Direct Application of SHAKE to the Velocity Verlet Algorithm. J. Comput. Phys. 1993, 104 (2), 470-472.) or similar, since there is no mention of RATTLE being used and the velocity Verlet algorithm is used. I am trying to figure out exactly how adding another constraint on the velocity of a region will modify the equations of motion when SHAKE is used.


It's conceptually simple. SHAKE adjusts the forces on atoms,
so that when the atom coords are updated (at the beginning
of the next timestep), the bond/angle constraints will be met.
Thus fix shake adds a force to atoms, and does it at the point
in the timestep when other force constraints are typically added,
e.g. wall forces, after pairwise and molecular forces have been