where the summation is over geometric constraints, is the Lagrangian for the extended system, and is a Lagrange multiplier associated with a geometric constraint :

(6.21) |

with being a geometric parameter and is the value of fixed during the simulation. In the SHAKE algorithm, Lagrange multipliers are determined in the iterative procedure:

- perform a standard MD step (leap-frog algorithm):
(6.22)

(6.23)

- use the new positions
to compute Lagrange multipliers for all constraints:
(6.24)

- update the velocities and positions by adding a contribution due to restoring forces (proportional
to ):
(6.25)

(6.26)

- repeat steps 2-4 until all are smaller than a predefined tolerance.