In computational chemistry, a constraint algorithm is a method for satisfying the Newtonian motion of a rigid body which consists of mass points. A restraint algorithm is used to ensure that the distance between mass points is maintained. The general steps involved are: (i) choose novel unconstrained coordinates (internal coordinates), (ii) introduce explicit constraint forces, (iii) minimize constraint forces implicitly by the technique of Lagrange multipliers or projection methods. Constraint algorithms are often applied to molecular dynamics simulations. Although such simulations are sometimes performed using internal coordinates that automatically satisfy the bond-length, bond-angle and torsion-angle constraints, simulations may also be performed using explicit or implicit constraint forces for these three constraints. However, explicit constraint forces give rise to inefficiency; more computational power is required to get a trajectory of a given length. Therefore, internal coordinates and implicit-force constraint solvers are generally preferred. Constraint algorithms achieve computational efficiency by neglecting motion along some degrees of freedom. For instance, in atomistic molecular dynamics, typically the length of covalent bonds to hydrogen are constrained; however, constraint algorithms should not be used if vibrations along these degrees of freedom are important for the phenomenon being studied.
Mathematical background The motion of a set of N particles can be described by a set of second-order ordinary differential equations, Newton's second law, which can be written in matrix form
M ⋅ d 2 q d t 2 = f = − ∂ V ∂ q {\displaystyle \mathbf {M} \cdot {\frac {d^{2}\mathbf {q} }{dt^{2}}}=\mathbf {f} =-{\frac {\partial V}{\partial \mathbf {q} }}}
where M is a mass matrix and q is the vector of generalized coordinates that describe the particles' positions. For example, the vector q may be a 3N Cartesian coordinates of the particle positions rk, where k runs from 1 to N; in the absence of constraints, M would be the 3Nx3N diagonal square matrix of the particle masses. The vector f represents the generalized forces and the scalar V(q) represents the potential energy, both of which are functions of the generalized coordinates q. If M constraints are present, the coordinates must also satisfy M time-independent algebraic equations
g j ( q ) = 0 {\displaystyle g_{j}(\mathbf {q} )=0}
where the index j runs from 1 to M. For brevity, these functions gi are grouped into an M-dimensional vector g below. The task is to solve the combined set of differential-algebraic (DAE) equations, instead of just the ordinary differential equations (ODE) of Newton's second law. This problem was studied in detail by Joseph Louis Lagrange, who laid out most of the methods for solving it. The simplest approach is to define new generalized coordinates that are unconstrained; this approach eliminates the algebraic equations and reduces the problem once again to solving an ordinary differential equation. Such an approach is used, for example, in describing the motion of a rigid body; the position and orientation of a rigid body can be described by six independent, unconstrained coordinates, rather than describing the positions of the particles that make it up and the constraints among them that maintain their relative distances. The drawback of this approach is that the equations may become unwieldy and complex; for example, the mass matrix M may become non-diagonal and depend on the generalized coordinates. A second approach is to introduce explicit forces that work to maintain the constraint; for example, one could introduce strong spring forces that enforce the distances among mass points within a "rigid" body. The two difficulties of this approach are that the constraints are not satisfied exactly, and the strong forces may require very short time-steps, making simulations inefficient computationally. A third approach is to use a method such as Lagrange multipliers or projection to the constraint manifold to determine the coordinate adjustments necessary to satisfy the constraints. Finally, there are various hybrid approaches in which different sets of constraints are satisfied by different methods, e.g., internal coordinates, explicit forces and implicit-force solutions.
Internal coordinate methods The simplest approach to satisfying constraints in energy minimization and molecular dynamics is to represent the mechanical system in so-called internal coordinates corresponding to unconstrained independent degrees of freedom of the system. For example, the dihedral angles of a protein are an independent set of coordinates that specify the positions of all the atoms without requiring any constraints. The difficulty of such internal-coordinate approaches is twofold: the Newtonian equations of motion become much more complex and the internal coordinates may be difficult to define for cyclic systems of constraints, e.g., in ring puckering or when a protein has a disulfide bond. The original methods for efficient recursive energy minimization in internal coordinates were developed by Gō and coworkers. Efficient recursive, internal-coordinate constraint solvers were extended to molecular dynamics. Analogous methods were applied later to other systems.
Lagrange multiplier-based methods
In most of molecular dynamics simulations that use constraint algorithms, constraints are enforced using the method of Lagrange multipliers. Given a set of n linear (holonomic) constraints at the time t,
σ k ( t ) := ‖ x k α ( t ) − x k β ( t ) ‖ 2 − d k 2 = 0 , k = 1 … n {\displaystyle \sigma _{k}(t):=\|\mathbf {x} _{k\alpha }(t)-\mathbf {x} _{k\beta }(t)\|^{2}-d_{k}^{2}=0,\quad k=1\ldots n}
where x k α ( t ) {\displaystyle \scriptstyle \mathbf {x} _{k\alpha }(t)} and x k β ( t ) {\displaystyle \scriptstyle \mathbf {x} _{k\beta }(t)} are the positions of the two particles involved in the kth constraint at the time t and d k {\displaystyle d_{k}} is the prescribed inter-particle distance. The forces due to these constraints are added in the equations of motion, resulting in, for each of the N particles in the system
∂ 2 x i ( t ) ∂ t 2 m i = − ∂ ∂ x i [ V ( x i ( t ) ) − ∑ k = 1 n λ k σ k ( t ) ] , i = 1 … N . {\displaystyle {\frac {\partial ^{2}\mathbf {x} _{i}(t)}{\partial t^{2}}}m_{i}=-{\frac {\partial }{\partial \mathbf {x} _{i}}}\left[V(\mathbf {x} _{i}(t))-\sum _{k=1}^{n}\lambda _{k}\sigma _{k}(t)\right],\quad i=1\ldots N.}
Adding the constraint forces does not change the total energy, as the net work done by the constraint forces (taken over the set of particles that the constraints act on) is zero. Note that the sign on λ k {\displaystyle \lambda _{k}} is arbitrary and some references have an opposite sign. From integrating both sides of the equation with respect to the time, the constrained coordinates of particles at the time, t + Δ t {\displaystyle t+\Delta t} , are given,
x i ( t + Δ t ) = x ^ i ( t + Δ t ) + ∑ k = 1 n λ k ∂ σ k ( t ) ∂ x i ( Δ t ) 2 m i − 1 , i = 1 … N {\displaystyle \mathbf {x} _{i}(t+\Delta t)={\hat {\mathbf {x} }}_{i}(t+\Delta t)+\sum _{k=1}^{n}\lambda _{k}{\frac {\partial \sigma _{k}(t)}{\partial \mathbf {x} _{i}}}\left(\Delta t\right)^{2}m_{i}^{-1},\quad i=1\ldots N}
where x ^ i ( t + Δ t ) {\displaystyle {\hat {\mathbf {x} }}_{i}(t+\Delta t)} is the unconstrained (or uncorrected) position of the ith particle after integrating the unconstrained equations of motion. To satisfy the constraints σ k ( t + Δ t ) {\displaystyle \sigma _{k}(t+\Delta t)} in the next timestep, the Lagrange multipliers should be determined as the following equation,
σ k ( t + Δ t ) := ‖ x k α ( t + Δ t ) − x k β ( t + Δ t ) ‖ 2 − d k 2 = 0. {\displaystyle \sigma _{k}(t+\Delta t):=\left\|\mathbf {x} _{k\alpha }(t+\Delta t)-\mathbf {x} _{k\beta }(t+\Delta t)\right\|^{2}-d_{k}^{2}=0.}
This implies solving a system of n {\displaystyle n} non-linear equations
σ j ( t + Δ t ) := ‖ x ^ j α ( t + Δ t ) − x ^ j β ( t + Δ t ) + ∑ k = 1 n λ k ( Δ t ) 2 [ ∂ σ k ( t ) ∂ x j α m j α − 1 − ∂ σ k ( t ) ∂ x j β m j β − 1 ] ‖ 2 − d j 2 = 0 , j = 1 … n {\displaystyle \sigma _{j}(t+\Delta t):=\left\|{\hat {\mathbf {x} }}_{j\alpha }(t+\Delta t)-{\hat {\mathbf {x} }}_{j\beta }(t+\Delta t)+\sum _{k=1}^{n}\lambda _{k}\left(\Delta t\right)^{2}\left[{\frac {\partial \sigma _{k}(t)}{\partial \mathbf {x} _{j\alpha }}}m_{j\alpha }^{-1}-{\frac {\partial \sigma _{k}(t)}{\partial \mathbf {x} _{j\beta }}}m_{j\beta }^{-1}\right]\right\|^{2}-d_{j}^{2}=0,\quad j=1\ldots n}
simultaneously for the n {\displaystyle n} unknown Lagrange multipliers λ k {\displaystyle \lambda _{k}} . This system of n {\displaystyle n} non-linear equations in n {\displaystyle n} unknowns is commonly solved using Newton–Raphson method where the solution vector λ _ {\displaystyle {\underline {\lambda }}} is updated using
λ _ ( l + 1 ) ← λ _ ( l ) − J σ − 1 σ _ ( t + Δ t ) {\displaystyle {\underline {\lambda }}^{(l+1)}\leftarrow {\underline {\lambda }}^{(l)}-\mathbf {J} _{\sigma }^{-1}{\underline {\sigma }}(t+\Delta t)}
where J σ {\displaystyle \mathbf {J} _{\sigma }} is the Jacobian of the equations σk:
J = ( ∂ σ 1 ∂ λ 1 ∂ σ 1 ∂ λ 2 ⋯ ∂ σ 1 ∂ λ n ∂ σ 2 ∂ λ 1 ∂ σ 2 ∂ λ 2 ⋯ ∂ σ 2 ∂ λ n ⋮ ⋮ ⋱ ⋮ ∂ σ n ∂ λ 1 ∂ σ n ∂ λ 2 ⋯ ∂ σ n ∂ λ n ) . {\displaystyle \mathbf {J} =\left({\begin{array}{cccc}{\frac {\partial \sigma _{1}}{\partial \lambda _{1}}}&{\frac {\partial \sigma _{1}}{\partial \lambda _{2}}}&\cdots &{\frac {\partial \sigma _{1}}{\partial \lambda _{n}}}\\[5pt]{\frac {\partial \sigma _{2}}{\partial \lambda _{1}}}&{\frac {\partial \sigma _{2}}{\partial \lambda _{2}}}&\cdots &{\frac {\partial \sigma _{2}}{\partial \lambda _{n}}}\\[5pt]\vdots &\vdots &\ddots &\vdots \\[5pt]{\frac {\partial \sigma _{n}}{\partial \lambda _{1}}}&{\frac {\partial \sigma _{n}}{\partial \lambda _{2}}}&\cdots &{\frac {\partial \sigma _{n}}{\partial \lambda _{n}}}\end{array}}\right).}
Since not all particles contribute to all of constraints, J σ {\displaystyle \mathbf {J} _{\sigma }} is a block matrix and can be solved individually to block-unit of the matrix. In other words, J σ {\displaystyle \mathbf {J} _{\sigma }} can be solved individually for each molecule. Instead of constantly updating the vector λ _ {\displaystyle {\underline {\lambda }}} , the iteration can be started with λ _ ( 0 ) = 0 {\displaystyle {\underline {\lambda }}^{(0)}=\mathbf {0} } , resulting in simpler expressions for σ k ( t ) {\displaystyle \sigma _{k}(t)} and ∂ σ k ( t ) ∂ λ j {\displaystyle {\frac {\partial \sigma _{k}(t)}{\partial \lambda _{j}}}} . In this case
J i j = ∂ σ j ∂ λ i | λ = 0 = 2 [ x ^ j α − x ^ j β ] [ ∂ σ i ∂ x j α m j α − 1 − ∂ σ i ∂ x j β m j β − 1 ] . {\displaystyle J_{ij}=\left.{\frac {\partial \sigma _{j}}{\partial \lambda _{i}}}\right|_{\mathbf {\lambda } =0}=2\left[{\hat {x}}_{j\alpha }-{\hat {x}}_{j\beta }\right]\left[{\frac {\partial \sigma _{i}}{\partial x_{j\alpha }}}m_{j\alpha }^{-1}-{\frac {\partial \sigma _{i}}{\partial x_{j\beta }}}m_{j\beta }^{-1}\right].}
then λ {\displaystyle \lambda } is updated to
λ j = − J − 1 [ ‖ x ^ j α ( t + Δ t ) − x ^ j β ( t + Δ t ) ‖ 2 − d j 2 ] . {\displaystyle \mathbf {\lambda } _{j}=-\mathbf {J} ^{-1}\left[\left\|{\hat {\mathbf {x} }}_{j\alpha }(t+\Delta t)-{\hat {\mathbf {x} }}_{j\beta }(t+\Delta t)\right\|^{2}-d_{j}^{2}\right].}
After each iteration, the unconstrained particle positions are updated using
x ^ i ( t + Δ t ) ← x ^ i ( t + Δ t ) + ∑ k = 1 n λ k ∂ σ k ∂ x i ( Δ t ) 2 m i − 1 . {\displaystyle {\hat {\mathbf {x} }}_{i}(t+\Delta t)\leftarrow {\hat {\mathbf {x} }}_{i}(t+\Delta t)+\sum _{k=1}^{n}\lambda _{k}{\frac {\partial \sigma _{k}}{\partial \mathbf {x} _{i}}}\left(\Delta t\right)^{2}m_{i}^{-1}.}
The vector is then reset to
λ _ = 0 . {\displaystyle {\underline {\lambda }}=\mathbf {0} .}
The above procedure is repeated until the solution of constraint equations, σ k ( t + Δ t ) {\displaystyle \sigma _{k}(t+\Delta t)} , converges to a prescribed tolerance of a numerical error. Although there are a number of algorithms to compute the Lagrange multipliers, these difference is rely only on the methods to solve the system of equations. For this methods, quasi-Newton methods are commonly used.
The SETTLE algorithm The SETTLE algorithm solves the system of non-linear equations analytically for n = 3 {\displaystyle n=3} constraints in constant time. Although it does not scale to larger numbers of constraints, it is very often used to constrain rigid water molecules, which are present in almost all biological simulations and are usually modelled using three constraints (e.g. SPC/E and TIP3P water models).
The SHAKE algorithm The SHAKE algorithm was first developed for satisfying a bond geometry constraint during molecular dynamics simulations. 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;
λ _ = − J σ − 1 σ _ . {\displaystyle {\underline {\lambda }}=-\mathbf {J} _{\sigma }^{-1}{\underline {\sigma }}.}
This amounts to assuming that J σ {\displaystyle \mathbf {J} _{\sigma }} is diagonally dominant and solving the k {\displaystyle k} th equation only for the k {\displaystyle k} unknown. In practice, we compute
λ k ← σ k ( t ) ∂ σ k ( t ) / ∂ λ k , x k α ← x k α + λ k ∂ σ k ( t ) ∂ x