The constraint force is evaluated along the bond directions at the old positions \(\mathbf{r}^n\); its multipliers are chosen so that the new positions satisfy the constraints.
The numerical integration can be done in two steps: (1) integration as if there are no constraints and (2) adding the displacement due to ‘’constrain force’’ \(\delta\mathbf{f}_k\), so that the constraints are satisfied:
Here, \(\mu_c=-2\lambda_c^n\tau^2\) is the displacement multiplier for this leapfrog scheme and the sign convention used above. The relation between force and displacement multipliers depends on the integration scheme.
For the leapfrog half-step velocities, \(\delta\mathbf{v}_k^{n+1/2}=\delta\mathbf{r}_k/\tau\), so the corresponding velocity multiplier is \(\mu_c^v=\mu_c/\tau\). This update does not by itself impose the end-of-step velocity constraints in a velocity-Verlet scheme; those require a separate velocity projection, as in RATTLE.
For each constrain c, after the numerical integration, the distance between constrained atoms should be equal to the target distance of that constraint, \(d_c\). Hence, for each constraint we have \(|\mathbf{r}_{ij}^{n+1}|=d_{ij}\) or \({\mathbf{r}_{ij}^{n+1}}^2=d_{ij}^2\), where \(\mathbf{r}^{n+1}_{ij}=\mathbf{r}^{n+1}_{i}-\mathbf{r}^{n+1}_{j}\).
These conditions should be satisfied by finding appropriate values for \(\mu_c\) (or equally \(\lambda_c\)). The general form for these quadratic equations is:
Note that the sum over all constraints \(c=\{m,l\}\) should include the pair \(\{i,j\}\), in which case only two of the kronecker’s deltas will survive. For a coupled constraint, when either \(i\) or \(j\) is equal to one index in the \(\{m,l\}\) pair, only one delta will not be zero. Since \(i\ne j\) and \(m\ne l\), there is never going to be the case when more than two deltas are not zero.
Before applying the above formula to specific cases, let us simplify the notation. We write \(\mathbf{x}_i^0\), \(\mathbf{x}_i^*\), and \(\mathbf{x}_i\) for the old, trial, and corrected positions of atom \(i\). For an oriented constraint \(c=(i,j)\), we use \(\mathbf{r}_c^0=\mathbf{x}_i^0-\mathbf{x}_j^0\), \(\mathbf{r}_c^*=\mathbf{x}_i^*-\mathbf{x}_j^*\), and \(\mathbf{r}_c=\mathbf{x}_i-\mathbf{x}_j\) for its bond vectors. Thus atom indices on \(\mathbf{x}\) and constraint indices on \(\mathbf{r}\) are distinct. On position and bond vectors, the superscript \(0\) denotes the old integration step, not a Newton iteration. Constraint orientations and indices are specified separately for each case. All bond vectors must use consistent periodic images.
The simplest case by far is when we have a single constraint, which is not coupled to any other. An example should be the constrained C-H bond in the protein backbone in the case when only bonds with hydrogens are constrained (with the only exception of glycine, where there is another hydrogen instead of a side-chain). Indeed, this case is very simple: we have the direction of the bond in which atoms should be moved. And the sum of the displacements for two moved atoms should be so that the final distance is equal to the target. The displacements are reversely proportional to masses of atoms. Let us derive this formally, using the Equation 1 above.
Fig. 6 A schematic representation of a constrained bond. Index on top of the bond represents the constraint index \(c\). Arrow indicates the direction of the vector \(\mathbf{r}_1\). Masses subscripts correspond to aton indices.
We have one constraint \(c=1=\{0,1\}\) that connects two atoms \(i=0\) and \(j=1\) with masses \(m_0\) and \(m_1\). We will also denote all the connecting vectors with constrain index, i.e. \(\mathbf{r}_1=\mathbf{x}_1-\mathbf{x}_0\), etc. so there is less noise in formulas. Note that we could remove the index altogether since there is only one constrained bond, but we will keep the index here for consistency with the rest of the cases. The sole equation derived from Equation 1 is:
For a positive discriminant, there are two real solutions. We assume that the trial bond vector is close to the old one, not merely that its length is close to the target. In this regime, \(\mathbf{r}_1^*\cdot\mathbf{r}_1^0>0\), and the solution of smaller absolute magnitude gives the small correction connected to the previous configuration. The other root reverses the component of the bond along its old direction. A zero discriminant gives a repeated root; a negative discriminant means that no real correction along the old bond direction can satisfy this isolated constraint.
Fig. 7 Two roots give two different directions of the bond after constraining. Given that the timestep is small, the displacement should be small as well. Hence the correct solution for \(\mu_1\) is the smallest by absolute value.
The solution with the plus sign before the square root is the small root under the positive-dot-product assumption above. The resulting expression for \(\mu_1\) is:
For numerical evaluation near an already satisfied constraint, the subtraction in Equation 3 can lose significant digits. Define \(a={\mathbf{r}_1^0}^2\), \(b=\mathbf{r}_1^*\cdot\mathbf{r}_1^0\), and \(c={\mathbf{r}_1^*}^2-d_1^2\). An algebraically equivalent, cancellation-resistant expression for \(b>0\) is
The discriminant must be checked before taking the square root. If \(b\leq0\), the small-step assumption has failed and the plus-root choice is no longer justified by this argument.
Though this equation is quite simple, we need to prepare ourselves to the following cases. So let us introduce the notations for all the coefficients in the Equation 2, re-writing it as:
\[k_1^{11}\mu_1^2+k_1^1\mu_1+k_1^0=0\]
Here, the indices for parameters \(k_x^y\) are chosen the following way. The bottom index is the number of equation (we have one so far). The top indexes indicate what combination of variables \(\mu\) this coefficient is multiplied by: \(k_1^{11}\) is for \(\mu_1\mu_1=\mu^2\), \(k_1^1\) is for \(\mu_1\) and \(0\) is for free coefficient. Later we will encounter other variables (\(\mu_2\), \(\mu_3\), etc.), for which we will have \(k_1^{22}\), \(k_1^{33}\) as well as some mixed products, e.g. \(\mu_1\mu_2\) with \(k_1^{12}\), etc. As follows from Equation 2, in the case of one uncoupled constraint, there are three coefficient:
Fig. 8 A schematic representation of two coupled constrained bonds.
Orient both bond vectors from the central atom to the outer atoms: \(\mathbf{r}_1=\mathbf{x}_1-\mathbf{x}_0\) and \(\mathbf{r}_2=\mathbf{x}_2-\mathbf{x}_0\).
Two coupled constraints will give us a much more complicated case of two coupled quadratic equations. Assuming that the central atom has index \(0\), with atoms \(1\) and \(2\) connected to it with two constrained bonds \(c=1=\{0,1\}\) and \(c=2=\{0,2\}\), and using the reduced masses \(M_1=\frac{m_0m_1}{m_0+m_1}\) and \(M_2=\frac{m_0m_2}{m_0+m_2}\), from Equation 1 we get:
Equation 6 are two coupled quadratic equations of a general form. Unfortunately, there is no simple analytical solution to it. Hence the numerical method should be used. Before describing the methods used currently in most MD software packages, let us try to introduce a method for the specific cases we are considering here.
If we introduce two dimensional vector-function \(F(\mu_1, \mu_2)\) in Equation 6 it becomes:
\[\mathbf{F}(\mu_1, \mu_2) = \mathbf{0}\]
This is a system of non-linear equations, which are notoriously hard to solve if we (1) need to find all the solutions and (2) don’t have good initial guess for these solutions. However, this is not the case here. Firstly, we assume that atoms did satisfy the constraints on the previous step (so we seek the solution near the zero multiplier vector). Secondly, we only need to find this solution (see how we dropped one of the roots when solving for uncoupled constraint). Hence we can apply the Newtons method (also known as Newton-Raphson method). The general form of in in multi-dimensional case is:
Here, \(\mathbf{\mu}^n=(\mu_1^n,\mu_2^n)^T\) is a vector with the solutions on \(n\)-th iteration, \(J_\mathbf{F}\left(\mathbf{\mu}\right)=J_\mathbf{F}\left(\mu_1, \mu_2\right)\) is the Jacobian matrix for the system. Using Equation 6, we can compute \(J_\mathbf{F}\left(\mu_1, \mu_2\right)\):
With everything in the Equation 8 derived, we need good initial approximation for \(\mu_1\) and \(\mu_2\). One can use zeroes, assuming that atoms did not moved much during the integration time step atoms did not moved far away from their previous positions (which did satisfy the constraints). An alternative initial approximation is obtained by assuming that the constraints are not coupled, i.e. by using Equation 3 or Equation 5 for two constraints we have in our case:
The numerical procedure is as follows. (1) Compute the coefficients in Equation 7. (2) Use the independent-constraint roots in Equation 9 as initial estimates, evaluating them with the stable form Equation 4 when applicable. (3) Evaluate \(\mathbf{F}\) and \(J_\mathbf{F}\). (4) Solve \(J_\mathbf{F}\Delta\boldsymbol{\mu}=-\mathbf{F}\) and update \(\boldsymbol{\mu}\leftarrow\boldsymbol{\mu}+\Delta\boldsymbol{\mu}\). Solving this small linear system avoids explicitly forming the inverse. (5) Repeat until \(\max_c |F_c|/d_c^2\) is below the chosen tolerance, then apply the atomic displacements and verify the resulting bond lengths. This is a relative squared-distance tolerance; for small errors it is approximately twice the relative length error. Newton convergence is local, so invalid initial roots, singular Jacobians, and failure to converge within an iteration limit must be handled explicitly.
Three constraints coupled through the central atom
Fig. 9 A schematic representation of three atoms constrained to one central atom.
As in the two-bond case, atom \(0\) is the center. For \(c=1,2,3\), define \(\mathbf{r}_c=\mathbf{x}_c-\mathbf{x}_0\) and \(M_c=m_0m_c/(m_0+m_c)\). Use these orientations for the old and trial bond vectors as well.
Four constraints coupled through the central atom
Consider a central atom \(0\) connected to four outer atoms \(1,2,3,4\), with no additional constraints between them. A typical example is methane with its four C-H bond lengths constrained. This fixes the four distances to the center, but does not fix the bond angles or make the whole group rigid.
As above, orient each bond from the central atom to an outer atom:
For each \(c=1,2,3,4\), define \(F_c=\left\lVert\mathbf{r}_c^*+\sum_{a=1}^4H_{ca}\mu_a\mathbf{r}_a^0\right\rVert^2-d_c^2\). Expanding gives four quadratic equations, each with fifteen coefficients:
For example, \(k_1^{11}=\lVert\mathbf{r}_1^0\rVert^2/M_1^2\), \(k_1^{44}=\lVert\mathbf{r}_4^0\rVert^2/m_0^2\), \(k_1^{14}=2(\mathbf{r}_1^0\cdot\mathbf{r}_4^0)/(M_1m_0)\), and \(k_1^{24}=2(\mathbf{r}_2^0\cdot\mathbf{r}_4^0)/m_0^2\). These coefficients are computed once for a given trial step and reused during the Newton iterations.
The Jacobian is a \(4\times4\) matrix. With the shorthand \(k_c^{ba}=k_c^{ab}\) for mixed coefficients, its entries are
As in the two- and three-bond cases, initialize each multiplier with its isolated-bond solution. Let \(a_c=\lVert\mathbf{r}_c^0\rVert^2\), \(b_c=\mathbf{r}_c^*\cdot\mathbf{r}_c^0\), and \(g_c=\lVert\mathbf{r}_c^*\rVert^2-d_c^2\). Under the same conditions as Equation 4, use
Stop when \(\max_{c=1,2,3,4}|F_c|/d_c^2\) satisfies the requested tolerance. Check invalid discriminants, singular or poorly conditioned linear systems, and the iteration limit as in the two-bond case. Finally, apply Equation 12 and verify all four corrected distances. The masses and bond lengths need not be equal, and the reference geometry need not be tetrahedral.
Use the same orientations for the old and trial vectors. The reduced masses are \(M_1=m_1m_2/(m_1+m_2)\), \(M_2=m_2m_3/(m_2+m_3)\), and \(M_3=m_3m_1/(m_3+m_1)\). Then Equation 1 gives
The triangle of constraints case is very common for the biomolecular simulations since most of water models constrain both bonds between hydrogens and oxygen and H-O-H angle. And there is an analytical solution for this specific case, called SETTLE [Miyamoto and Kollman, 1992]. This algorithm is in fact used in most of the MD software. In brief, the idea of SETTLE is first to reduce the problem into two dimensions, and than compute the rotation angle of the constraints triangle using two-dimensional Euler angles.
The same approach can be extended to stars with more than four bonds and to other small connected groups, including rings. Additional constraints between outer atoms change the coupling matrix and must be included in the same system. We do not derive those cases here. When all covalent bonds are constrained, the connected groups can become large and a general constraint solver is needed. Two common algorithms are SHAKE [J.P. Ryckaert and Berendsen, 1977] (and its velocity counterpart RATTLE [Andersen, 1983]) and LINCS [Hess, 2008, Hess et al., 1997]. The latter is the default option in GROMACS.
H.C. Andersen. RATTLE: A “velocity” version of the shake algorithm for molecular dynamics calculations. J. Comput. Phys., 52(1):24–34, 1983. doi:10.1016/0021-9991(83)90014-1.
J.P. Ryckaert, G. Ciccotti and H.J.C. Berendsen. Numerical integration of the cartesian equations of motion of a system with constraints: Molecular Dynamics of n-alkanes. J. Comp. Phys., 23(3):327–341, 1977. doi:10.1016/0021-9991(77)90098-5.
S. Miyamoto and P.A. Kollman. SETTLE: An analytical version of the SHAKE and RATTLE algorithm for rigid water models. J. Comput. Chem., 13(8):952–962, 1992. doi:10.1002/jcc.540130805.
