Lesson 4.4 · 15 min

Constraints and SHAKE

The fastest motions in a molecule, bond vibrations involving hydrogen, force a tiny time step. Constraints freeze those bond lengths. SHAKE, RATTLE, SETTLE and LINCS enforce them every step, which lets simulations use 2 fs instead of 0.5 fs.

Why constraining X–H bonds allows larger time steps
How SHAKE corrects positions and RATTLE corrects velocities (Lagrange multipliers)
SETTLE and LINCS, the effect on degrees of freedom, and hydrogen mass repartitioning

The idea

Lesson 4.1 showed that an O–H or C–H bond vibrates with a period of about 10 fs, which limits flexible simulations to Δt ≈ 0.5–1 fs. These vibrations are also not very realistic classically: at room temperature they are almost entirely in their quantum ground state. So most simulations simply freeze them: the bond length is held fixed at its equilibrium value.

How. A constraint is a condition, |r_i − r_j| = d_ij, that must hold at every step. It is enforced by a constraint force along the bond, whose size is a Lagrange multiplier chosen exactly so the condition holds. In practice:

SHAKE (Ryckaert, Ciccotti and Berendsen, 1977): take the normal, unconstrained step, then move each pair of atoms back along their old bond direction until the length is right. Constraints share atoms, so this is repeated (sweeps) until all lengths agree within a tolerance.
RATTLE (Andersen, 1983): the velocity Verlet version. It also removes the velocity components along each constrained bond, so the bonds do not try to stretch.
SETTLE (Miyamoto and Kollman, 1992): an exact, non-iterative solution for rigid three-site water. Used for water in essentially all codes.
LINCS (Hess et al., 1997): solves all constraints at once with a matrix expansion. Fast and parallel; the GROMACS default for bonds.

Consequences. With all bonds to hydrogen constrained (and rigid water), Δt = 2 fs is standard. Each constraint removes one degree of freedom, which must be subtracted when computing the temperature. Making hydrogens heavier by moving mass from the atoms they are bonded to (hydrogen mass repartitioning, HMR) slows the remaining fast motions enough for 4 fs.

Key idea. Constrain X–H bonds (and water) to remove the fastest vibrations. SHAKE/RATTLE iterate, SETTLE is exact for water, LINCS handles the rest in parallel. Count N_f = 3N − N_c − 3.

The algorithm

Velocity Verlet with RATTLE:

v += (dt/2) * F / m
x_new = x + dt * v                         # unconstrained step
repeat until all |x_i - x_j| = d_ij within tol:      # SHAKE
    for each constraint (i, j, d):
        s = x_new_i - x_new_j;  r = x_i - x_j      # r = old bond vector
        g = (s·s - d²) / (2 (1/m_i + 1/m_j) s·r)
        x_new_i -= g r / m_i;  x_new_j += g r / m_j
v = (x_new - x) / dt;  x = x_new
F = forces(x)
v += (dt/2) * F / m
repeat until all r_ij·v_ij = 0:                       # RATTLE
    for each constraint: k = r·v_ij / ((1/m_i + 1/m_j) d²)
        v_i -= k r / m_i;  v_j += k r / m_j

Use a tight tolerance (10⁻⁶ to 10⁻⁸ relative). A loose tolerance lets the bond lengths drift and the energy is no longer conserved.
Constrain only bonds (to hydrogen, or all bonds). Constraining angles as well makes the system too rigid and is avoided, except in rigid water.
Thermostats and the temperature must use the reduced N_f; MD codes do this automatically from the topology.
GROMACS: constraints = h-bonds, constraint-algorithm = lincs (water uses SETTLE). LAMMPS: fix shake or fix rattle. AMBER: ntc = 2, ntf = 2 (SHAKE on bonds to H). HMR: parmed HMassRepartition, or the force field's HMR topology.