SHAKE Algorithm

SHAKE operates by iteratively resetting atomic positions after an unconstrained integration step. It applies a correction vector along each constrained bond to satisfy a fixed distance criterion. The algorithm loops through all constraints until all bond lengths converge to within a specified tolerance, typically 10⁻⁴ Å. This iterative Gauss-Seidel approach is simple but can struggle with highly interconnected constraint networks, such as those involving three-center bonds or ring systems.

The primary purpose of SHAKE is to eliminate the fastest vibrational modes from a simulation, specifically bonds involving hydrogen atoms (e.g., O-H, N-H, C-H). These bonds vibrate with a period of ~10 fs, which would otherwise dictate a maximum time step of ~1 fs. By constraining them, the integration time step can be safely increased to 2 fs in atomistic simulations, effectively doubling computational throughput without sacrificing accuracy for slower, more relevant conformational dynamics.

SHAKE only corrects atomic positions, leaving velocities unconstrained. RATTLE is the complementary algorithm that applies a second set of constraints to the velocities, ensuring they are orthogonal to the bond direction. This is essential for generating correct kinetic energy distributions and is required for proper NVT and NPT ensemble sampling. Most modern MD engines implement both SHAKE (for positions) and RATTLE (for velocities) as a unified constraint solver.

For the specific case of rigid water models like TIP3P or SPC, the iterative SHAKE procedure is computationally wasteful. SETTLE provides an exact, analytical solution to reset the positions and velocities of a three-site water molecule in a single step. By solving the constraint equations directly from the molecular geometry, SETTLE eliminates iteration overhead and is the standard method for maintaining water rigidity in virtually all biomolecular simulations.

The LINCS (LINear Constraint Solver) algorithm is a non-iterative alternative to SHAKE that resets bonds in a single step using a matrix inversion. After an unconstrained update, LINCS projects the new bond lengths back onto the constraint manifold using a power-series expansion of the inverse constraint coupling matrix. It is significantly faster than SHAKE for large, interconnected constraint networks and is the default constraint algorithm in the GROMACS simulation package.

SHAKE's iterative nature means it can fail to converge if the constraint network is over-specified or if atoms move too far in a single step. A relative tolerance (typically 10⁻⁴) defines when a bond is considered constrained. If the maximum number of iterations is exceeded, the time step is often rejected and retried with a smaller step. This is a common source of simulation crashes, particularly during energy minimization or when using hydrogen mass repartitioning schemes with aggressive time steps.

The foundational symplectic integrator for Newton's equations of motion. The Velocity Verlet variant calculates positions and velocities at the same time step, providing excellent energy conservation. SHAKE operates directly on the positions generated by this integrator, resetting bond lengths before the velocity update.

• Standard time step: 1-2 fs without constraints
• Error per step: O(Δt⁴) for positions
• Symplectic nature preserves phase space volume

The mathematical formalism underlying holonomic constraint enforcement. SHAKE iteratively solves for the undetermined multipliers that represent the constraint forces required to maintain fixed bond lengths. Each multiplier corresponds to the magnitude of the force acting along a specific rigid bond vector.

• Converts constrained dynamics to an unconstrained optimization problem
• RATTLE extends this to velocity constraints for kinetic energy consistency
• Direct matrix inversion is avoided via iterative Gauss-Seidel relaxation

A non-iterative constraint solver that resets bonds in a single pass using a power series expansion of the inverse constraint matrix. Significantly faster than SHAKE for large molecules and parallelizes efficiently on GPUs. The default constraint algorithm in GROMACS for bonds involving hydrogen atoms.

• Expansion order: typically 4-8 terms
• Handles bond-angle coupling via matrix inversion
• Preferred for GPU-accelerated MD due to deterministic execution time

An analytical constraint solver specifically for rigid water models (TIP3P, SPC/E). Solves the constraint equations exactly without iteration by treating the water molecule as a rigid triangle. Vastly more efficient than applying SHAKE to each water bond and angle individually.

• Constrains two O-H bonds and one H-O-H angle simultaneously
• Closed-form solution using quaternion rotation
• Essential for simulations where water constitutes >80% of atoms