 The method was then generalised to handle any holonomic constraint, such as those required to maintain constant bond angles, or molecular rigidity.In SHAKE algorithm, the system of non-linear constraint equations is solved using the Gauss–Seidel method which approximates the solution of the linear system of equations using the Newton–Raphson method;

This amounts to assuming that formula_21 is diagonally dominant and solving the formula_38th equation only for the formula_38 unknown.In practice, we compute

for all formula_41 iteratively until the constraint equations formula_34 are solved to a given tolerance.