computational chemistry Computational chemistry is a branch of chemistry that uses computer simulation to assist in solving chemical problems. It uses methods of theoretical chemistry, incorporated into computer programs, to calculate the structures and properties of m ...

, 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 In mathematical optimization, the method of Lagrange multipliers is a strategy for finding the local maxima and minima of a function subject to equality constraints (i.e., subject to the condition that one or more equations have to be satisfied ex ...

or projection methods. Constraint algorithms are often applied to

molecular dynamics Molecular dynamics (MD) is a computer simulation method for analyzing the physical movements of atoms and molecules. The atoms and molecules are allowed to interact for a fixed period of time, giving a view of the dynamic "evolution" of the ...

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 :

\mathbf \cdot \frac = \mathbf = -\frac

where M is a ''mass matrix'' and q is the

vector Vector most often refers to: *Euclidean vector, a quantity with a magnitude and a direction *Vector (epidemiology), an agent that carries and transmits an infectious pathogen into another living organism Vector may also refer to: Mathematic ...

generalized coordinate In analytical mechanics, generalized coordinates are a set of parameters used to represent the state of a system in a configuration space. These parameters must uniquely define the configuration of the system relative to a reference state.,p. 39 ...

s that describe the particles' positions. For example, the vector q may be a ''3N'' Cartesian coordinates of the particle positions r_''k'', where ''k'' runs from 1 to ''N''; in the absence of constraints, M would be the ''3N''x''3N'' 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_(\mathbf) = 0

where the index ''j'' runs from 1 to ''M''. For brevity, these functions ''g''_''i'' 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 Joseph-Louis Lagrange (born Giuseppe Luigi Lagrangia 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'', :

\sigma_k(t) := \,  \mathbf x_(t) - \mathbf x_(t) \, ^2 - d_k^2 = 0, \quad k=1 \ldots n

where

\scriptstyle \mathbf x_(t)

and

\scriptstyle\mathbf x_(t)

are the positions of the two particles involved in the ''k''th constraint at the time ''t'' and

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 :

\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

\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 + \Delta t

, are given, :

\mathbf x_i(t + \Delta t) = \hat_i(t + \Delta t) + \sum_^n \lambda_k \frac\left(\Delta t\right)^2m_i^, \quad i=1 \ldots N

where

\hat_i(t + \Delta t)

is the unconstrained (or uncorrected) position of the ''i''th particle after integrating the unconstrained equations of motion. To satisfy the constraints

\sigma_k(t + \Delta t)

in the next timestep, the

should be determined as the following equation, :

\sigma_k(t + \Delta t) := \left\,  \mathbf x_(t+\Delta t) - \mathbf x_(t+\Delta t)\right\, ^2 - d_k^2 = 0.

This implies solving a system of

n

non-linear equations :

\right\, ^2 - d_j^2 = 0, \quad j = 1 \ldots n

simultaneously for the

n

unknown Lagrange multipliers

\lambda_k

. This system of

n

non-linear equations in

n

unknowns is commonly solved using

Newton–Raphson method In numerical analysis, Newton's method, also known as the Newton–Raphson method, named after Isaac Newton and Joseph Raphson, is a root-finding algorithm which produces successively better approximations to the roots (or zeroes) of a real-valu ...

where the solution vector

\underline

is updated using :

\underline^ \leftarrow \underline^ - \mathbf J_\sigma^ \underline(t+\Delta t)

where

\mathbf J_\sigma

is the Jacobian of the equations σ_''k'': :

\frac & \frac & \cdots & \frac \end\right).

Since not all particles contribute to all of constraints,

\mathbf J_\sigma

is a

block matrix In mathematics, a block matrix or a partitioned matrix is a matrix that is '' interpreted'' as having been broken into sections called blocks or submatrices. Intuitively, a matrix interpreted as a block matrix can be visualized as the original mat ...

and can be solved individually to block-unit of the matrix. In other words,

\mathbf J_\sigma

can be solved individually for each molecule. Instead of constantly updating the vector

\underline

, the iteration can be started with

\underline^ = \mathbf 0

, resulting in simpler expressions for

\sigma_k(t)

and

\frac

. In this case :

J_ = \left.\frac\_ =  2\left hat_ - \hat_\right left frac m_^ - \frac m_^ \right

then

\lambda

is updated to :

\mathbf \lambda_j = - \mathbf J^\left \hat_(t+\Delta t) - \hat_(t+\Delta t)\right\, ^2 - d_j^2\right

After each iteration, the unconstrained particle positions are updated using :

\hat_i(t+\Delta t) \leftarrow \hat_i(t+\Delta t) + \sum_^n \lambda_k\frac \left(\Delta t\right)^2m_i^.

The vector is then reset to :

\underline = \mathbf 0.

The above procedure is repeated until the solution of constraint equations,

\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 method Quasi-Newton methods are methods used to either find zeroes or local maxima and minima of functions, as an alternative to Newton's method. They can be used if the Jacobian or Hessian is unavailable or is too expensive to compute at every iteration. ...

s are commonly used.

The SETTLE algorithm

The SETTLE algorithm solves the system of non-linear equations analytically for

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 In computational chemistry, a water model is used to simulate and thermodynamically calculate water clusters, liquid water, and aqueous solutions with explicit solvent. The models are determined from quantum mechanics, molecular mechanics, experim ...

water model In computational chemistry, a water model is used to simulate and thermodynamically calculate water clusters, liquid water, and aqueous solutions with explicit solvent. The models are determined from quantum mechanics, molecular mechanics, experime ...

s).

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 In numerical linear algebra, the Gauss–Seidel method, also known as the Liebmann method or the method of successive displacement, is an iterative method used to solve a system of linear equations. It is named after the German mathematicians Carl ...

which approximates the solution of the linear system of equations using the

; :

\underline = -\mathbf J_\sigma^ \underline.

This amounts to assuming that

\mathbf J_\sigma

is diagonally dominant and solving the

k

th equation only for the

k

unknown. In practice, we compute :

\mathbf x_ & \leftarrow \mathbf x_ + \lambda_k \frac, \end

for all

k=1\ldots n

iteratively until the constraint equations

\sigma_k(t+\Delta t)

are solved to a given tolerance. The calculation cost of each iteration is

\mathcal O(n)

, and the iterations themselves converge linearly. A non

iterative Iteration is the repetition of a process in order to generate a (possibly unbounded) sequence of outcomes. Each repetition of the process is a single iteration, and the outcome of each iteration is then the starting point of the next iteration. ...

form of SHAKE was developed later on. Several variants of the SHAKE algorithm exist. Although they differ in how they compute or apply the constraints themselves, the constraints are still modelled using Lagrange multipliers which are computed using the Gauss–Seidel method. The original SHAKE algorithm is capable of constraining both rigid and flexible molecules (eg. water, benzene and

biphenyl Biphenyl (also known as diphenyl, phenylbenzene, 1,1′-biphenyl, lemonene or BP) is an organic compound that forms colorless crystals. Particularly in older literature, compounds containing the functional group consisting of biphenyl less one ...

) and introduces negligible error or

energy drift In computer simulations of mechanical systems, energy drift is the gradual change in the total energy of a closed system over time. According to the laws of mechanics, the energy should be a constant of motion and should not change. However, in s ...

into a molecular dynamics simulation. One issue with SHAKE is that the number of iterations required to reach a certain level of convergence does rise as molecular geometry becomes more complex. To reach 64 bit computer accuracy (a relative tolerance of

\approx 10^

) in a typical molecular dynamics simulation at a temperature of 310K, a 3-site water model having 3 constraints to maintain molecular geometry requires an average of 9 iterations (which is 3 per site per time-step). A 4-site butane model with 5 constraints needs 17 iterations (22 per site), a 6-site benzene model with 12 constraints needs 36 iterations (72 per site), while a 12-site

model with 29 constraints requires 92 iterations (229 per site per time-step). Hence the CPU requirements of the SHAKE algorithm can become significant, particularly if a molecular model has a high degree of rigidity. A later extension of the method, QSHAKE (

Quaternion In mathematics, the quaternion number system extends the complex numbers. Quaternions were first described by the Irish mathematician William Rowan Hamilton in 1843 and applied to mechanics in three-dimensional space. Hamilton defined a quatern ...

SHAKE) was developed as a faster alternative for molecules composed of rigid units, but it is not as general purpose. It works satisfactorily for ''rigid'' loops such as

aromatic ring In chemistry, aromaticity is a chemical property of cyclic ( ring-shaped), ''typically'' planar (flat) molecular structures with pi bonds in resonance (those containing delocalized electrons) that gives increased stability compared to saturat ...

systems but QSHAKE fails for flexible loops, such as when a

protein Proteins are large biomolecules and macromolecules that comprise one or more long chains of amino acid residues. Proteins perform a vast array of functions within organisms, including catalysing metabolic reactions, DNA replication, respo ...

has a disulfide bond. Further extensions include RATTLE, WIGGLE, and MSHAKE. While RATTLE works the same way as SHAKE, yet using the

Velocity Verlet Verlet integration () is a numerical method used to integrate Newton's equations of motion. It is frequently used to calculate trajectories of particles in molecular dynamics simulations and computer graphics. The algorithm was first used in 179 ...

time integration scheme, WIGGLE extends SHAKE and RATTLE by using an initial estimate for the Lagrange multipliers

\lambda_k

based on the particle velocities. It is worth mentioning that MSHAKE computes corrections on the constraint ''forces'', achieving better convergence. A final modification to the SHAKE algorithm is the P-SHAKE algorithm that is applied to very rigid or semi-rigid

molecule A molecule is a group of two or more atoms held together by attractive forces known as chemical bonds; depending on context, the term may or may not include ions which satisfy this criterion. In quantum physics, organic chemistry, and bioch ...

s. P-SHAKE computes and updates a pre-conditioner which is applied to the constraint gradients before the SHAKE iteration, causing the Jacobian

\mathbf J_\sigma

to become diagonal or strongly diagonally dominant. The thus de-coupled constraints converge much faster (quadratically as opposed to linearly) at a cost of

\mathcal O(n^2)

The M-SHAKE algorithm

The M-SHAKE algorithm solves the non-linear system of equations using

Newton's method In numerical analysis, Newton's method, also known as the Newton–Raphson method, named after Isaac Newton and Joseph Raphson, is a root-finding algorithm which produces successively better approximations to the roots (or zeroes) of a real-valu ...

directly. In each iteration, the linear system of equations :

\underline = -\mathbf J_\sigma^ \underline

is solved exactly using an

LU decomposition In numerical analysis and linear algebra, lower–upper (LU) decomposition or factorization factors a matrix as the product of a lower triangular matrix and an upper triangular matrix (see matrix decomposition). The product sometimes includes a pe ...

. Each iteration costs

\mathcal O(n^3)

operations, yet the solution converges quadratically, requiring fewer iterations than SHAKE. This solution was first proposed in 1986 by Ciccotti and Ryckaert under the title "the matrix method", yet differed in the solution of the linear system of equations. Ciccotti and Ryckaert suggest inverting the matrix

\mathbf J_\sigma

directly, yet doing so only once, in the first iteration. The first iteration then costs

\mathcal O(n^3)

operations, whereas the following iterations cost only

\mathcal O(n^2)

operations (for the matrix-vector multiplication). This improvement comes at a cost though, since the Jacobian is no longer updated, convergence is only

linear Linearity is the property of a mathematical relationship (''function'') that can be graphically represented as a straight line. Linearity is closely related to '' proportionality''. Examples in physics include rectilinear motion, the linear r ...

, albeit at a much faster rate than for the SHAKE algorithm. Several variants of this approach based on sparse matrix techniques were studied by Barth ''et al.''.

SHAPE algorithm

The SHAPE algorithm is a multicenter analog of SHAKE for constraining rigid bodies of three or more centers. Like SHAKE, an unconstrained step is taken and then corrected by directly calculating and applying the rigid body rotation matrix that satisfies: :

L^\text \left( t + \frac 2 \right) = L^\text \left( t + \frac 2 \right)

This approach involves a single 3×3 matrix diagonalization followed by three or four rapid Newton iterations to determine the rotation matrix. SHAPE provides the identical trajectory that is provided with fully converged iterative SHAKE, yet it is found to be more efficient and more accurate than SHAKE when applied to systems involving three or more centers. It extends the ability of SHAKE like constraints to linear systems with three or more atoms, planar systems with four or more atoms, and to significantly larger rigid structures where SHAKE is intractable. It also allows rigid bodies to be linked with one or two common centers (e.g. peptide planes) by solving rigid body constraints iteratively in the same basic manner that SHAKE is used for atoms involving more than one SHAKE constraint.

LINCS algorithm

An alternative constraint method, LINCS (Linear Constraint Solver) was developed in 1997 by Hess, Bekker, Berendsen and Fraaije, and was based on the 1986 method of Edberg, Evans and Morriss (EEM), and a modification thereof by Baranyai and Evans (BE). LINCS applies Lagrange multipliers to the constraint forces and solves for the multipliers by using a series expansion to approximate the inverse of the Jacobian

\mathbf J_\sigma

: :

(\mathbf I - \mathbf J_\sigma)^ = \mathbf I + \mathbf J_\sigma + \mathbf J_\sigma^2 + \mathbf J_\sigma^3 + \cdots

in each step of the Newton iteration. This approximation only works for matrices with

eigenvalues In linear algebra, an eigenvector () or characteristic vector of a linear transformation is a nonzero vector that changes at most by a scalar factor when that linear transformation is applied to it. The corresponding eigenvalue, often denoted b ...

smaller than 1, making the LINCS algorithm suitable only for molecules with low connectivity. LINCS has been reported to be 3–4 times faster than SHAKE.

Hybrid methods

Hybrid methods have also been introduced in which the constraints are divided into two groups; the constraints of the first group are solved using internal coordinates whereas those of the second group are solved using constraint forces, e.g., by a Lagrange multiplier or projection method. This approach was pioneered by Lagrange, and result in ''Lagrange equations of the mixed type''.

References and footnotes

{{DEFAULTSORT:Constraint Algorithm Molecular dynamics Computational chemistry Molecular physics Computational physics Numerical differential equations

Mathematical background

Internal coordinate methods

Lagrange multiplier-based methods

The SETTLE algorithm

The SHAKE algorithm

The M-SHAKE algorithm

SHAPE algorithm

LINCS algorithm

Hybrid methods

See also

References and footnotes