(add numerical scheme here to show how lambda and mu are related)
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\tau^2\) is the normalized Lagrange multiplier.
Note that since \(\delta\mathbf{\upsilon}=\frac{\delta\mathbf{r}}{\tau}\), the scaled Lagrange multiplier for the velocities is \(\mu_c^\upsilon=\frac{\mu_c}{\tau}\).
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 notations. First, let us us the \(r_i^0\) for the coordinates of the atom before the integration and \(r_i\) as the coordinates for the same particle after the integration and constraining is done. Second, let us assume that the constrain connecting atoms \(m\) and \(l\) has an index \(c\), and we will use the constraint index \(c\) and indices of the connected particles \(m\) and \(l\) interchangeably (we already did this for \(d_c=d_{ml}\) above). We will be assigning these indexes for each case before deriving the formulas. This should allow us to have less noise in the formulas and the generalization for the \(n\)-th step and for arbitrary indices of particles should be straightforward.
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{r}_{01}=\mathbf{r}_0-\mathbf{r}_1\), 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:
As expected for quadratic equation, we have two solutions. However only one is valid for our problem: the one that gives smaller absolute value of \(\mu\). This is because we assume that the bond length after the numerical integration without constraints is close to the target length, i.e. the displacements of the atoms is small. Indeed, we can move the atoms so they satisfy the constraint two ways, one of which employs flipping the bond.
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 we are interested in is the one with the plus sign before the square root. So the final solution for \(\mu_1\) is:
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.
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 5 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 5 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 (i.e. we are looking for \(mu\)’s with smallest absolute value). 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 Jakobian matrix for the system. Using Equation 5, we can compute \(J_\mathbf{F}\left(\mu_1, \mu_2\right)\):
With everything in the Equation 7 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). We can get slightly better initial approximation by assuming that the constraints are not coupled, i.e. by using Equation 3 or Equation 4 for two constraints we have in our case:
The numerical procedure to evaluate \(\mu_1\) and \(\mu_2\) is then contains following steps. (1) Compute coefficients in Equation 6. (2) Use Equation 8 to get initial approximations \(\mu_1^0\) and \(\mu_2^0\) for \(\mu_1\) and \(\mu_2\). (3) Compute \(\mathbf{F}(\mu_1, \mu_2)\) using Equation 5 and \(J_\mathbf{F}\left(\mu_1, \mu_2\right)\), invert the matrix \(J_\mathbf{F}\) (one can derive formulas to compute the inverse Jakobian directly, without computing the Jakobian itself). (4) Apply the Equation 7 to get the next iteration values for \(\mu_1\) and \(\mu_2\). (5) Repeat steps (3) and (4) until convergence.
Three constraints coupled through the central atom
Fig. 9 A schematic representation of three atoms constrained to one central atom.
Note that \(\left(\mathbf{r}_1^0\cdot\mathbf{r}_2^0\right)=\left(\mathbf{r}_2^0\cdot\mathbf{r}_3^0\right)=\left(\mathbf{r}_3^0\cdot\mathbf{r}_1^0\right)=-2S_0\), where \(S_0\) is the area of the triangle.
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.
We can keep expanding our formulas to the cases with more coupled constraints, and they will get more and more complex. However, there are only a few cases that is useful to solve this way: four constraints coupled to the central (e.g. methane molecule) atom and various rings that should be kept planar (e.g. benzene ring in phenylalanine or double ring in the tryptophane side-chain). These are not very common, the rings are usually kept planar with improper dihedral potentials, so we will not derive them here. Much more important is the case where all covalent bonds are constraints, where universal formalism is needed. The two most common algorithms for this case are SHAKE [J.P. Ryckaert and Berendsen, 1977] (and its “velocity version” RATTLE [Andersen, 1983]) and LINCS [Hess, 2008, Hess et al., 1997]. Note that 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.
