Constraint algorithm
Encyclopedia
In mechanics
, a constraint algorithm is a method for satisfying constraints for bodies that obey Newton's equations of motion. There are three basic approaches to satisfying such constraints: choosing novel unconstrained coordinates ("internal coordinates"), introducing explicit constraint forces, and minimizing 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 carried out in internal coordinates that automatically satisfy the bond-length and bond-angle constraints, they may also be carried out with explicit or implicit constraint forces for the bond lengths and bond angles. Explicit constraint forces typically shorten the time-step significantly, making the simulation less efficient computationally; in other words, more computer power is required to compute a trajectory of a given length. Therefore, internal coordinates and implicit-force constraint solvers are generally preferred.
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
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.
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.
In most molecular dynamics
simulation, constraints are enforced using the method of Lagrange multipliers
. Given a set of n linear (holonomic
) constraints at the time t,
where
and 
are the positions of the two particles involved in the kth constraint at the time t and 
is the prescribed inter-particle distance.
These constraint equations, are added to the potential energy function in the equations of motion, resulting in, for each of the N particles in the system
Adding the constraint equations to the potential does not change it, since all
should, ideally, be zero.
Integrating both sides of the equations of motion twice in time yields the constrained particle positions at the time
where
is the unconstrained (or uncorrected) position of the ith particle after integrating the unconstrained equations of motion.
To satisfy the constraints
in the next timestep, the Lagrange multipliers
must be chosen such that
This implies solving a system of
non-linear equations
simultaneously for the
unknown Lagrange multipliers 
.
This system of
non-linear equations in 
unknowns is best solved using Newton's method
where the solution vector
is updated using
where
is the Jacobian
of the equations σk:
Since not all particles are involved in all constraints,
is blockwise-diagonal and can be solved blockwise, i.e. molecule for molecule.
Furthermore, instead of constantly updating the vector
, the iteration is startedwith 
, resulting in simpler expressions for 
and 
. After each iteration, the unconstrained particle positions are updated using
.
The vector is then reset to
This is repeated until the constraint equations
are satisfied up to a prescribed tolerance.
Although there are a number of algorithms to compute the Lagrange multipliers, they differ only in how they solve the system of equations, usually using Quasi-Newton methods.
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 model
s).
It solves the system of non-linear constraint equations using the Gauss-Seidel method
to approximate the solution of the linear system of equations
in the Newton iteration. This amounts to assuming that
is diagonally dominant and solving the 
th equation only for the 
unknown. In practice, we compute
for all
iteratively until the constraint equations 
are solved to a given tolerance.
Each iteration of the SHAKE algorithm costs
operations and the iterations themselves converge linearly.
A noniterative form of SHAKE was developed later.
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 limited to mechanical systems with a tree structure, i.e., no closed loops of constraints. A later extension of the method, QSHAKE (Quaternion
SHAKE) was developed to amend this. It works satisfactorily for rigid loops such as aromatic ring systems but fails for flexible loops, such as when a protein has a disulfide bond.
Further extensions include RATTLE, WIGGLE and MSHAKE.
RATTLE works the same way as SHAKE, yet using the Velocity Verlet time integration scheme. WIGGLE extends SHAKE and RATTLE by using an initial estimate for the Lagrange multipliers
based on the particle velocities. Finally, MSHAKE computes corrections on the constraint forces, achieving better convergence.
A final modification is the P-SHAKE algorithm for rigid or semi-rigid molecules. P-SHAKE computes and updates a pre-conditioner which is applied to the constraint gradients before the SHAKE iteration, causing the Jacobian
to become diagonal or strongly diagonally dominant. The thus de-coupled constraints converge much faster (quadratically as opposed to linearly) at a cost of 
.
directly. In each iteration, the linear system of equations
is solved exactly using an LU decomposition
. Each iteration costs
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
directly, yet doing so only once, in the first iteration. The first iteration then costs 
operations, whereas the following iterations cost only 
operations (for the matrix-vector multiplication). This improvement comes at a cost though, since the Jacobian is no longer updated, convergence is only linear, 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..
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
:
in each step of the Newton iteration. This approximation only works for matrices with Eigenvalues 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.
, and result in Lagrange equations of the mixed type.
Mechanics
Mechanics is the branch of physics concerned with the behavior of physical bodies when subjected to forces or displacements, and the subsequent effects of the bodies on their environment....
, a constraint algorithm is a method for satisfying constraints for bodies that obey Newton's equations of motion. There are three basic approaches to satisfying such constraints: choosing novel unconstrained coordinates ("internal coordinates"), introducing explicit constraint forces, and minimizing constraint forces implicitly by the technique of Lagrange multipliers
Lagrange multipliers
In mathematical optimization, the method of Lagrange multipliers provides a strategy for finding the maxima and minima of a function subject to constraints.For instance , consider the optimization problem...
or projection methods.
Constraint algorithms are often applied to molecular dynamics
Molecular dynamics
Molecular dynamics is a computer simulation of physical movements of atoms and molecules. The atoms and molecules are allowed to interact for a period of time, giving a view of the motion of the atoms...
simulations. Although such simulations are sometimes carried out in internal coordinates that automatically satisfy the bond-length and bond-angle constraints, they may also be carried out with explicit or implicit constraint forces for the bond lengths and bond angles. Explicit constraint forces typically shorten the time-step significantly, making the simulation less efficient computationally; in other words, more computer power is required to compute a trajectory of a given length. Therefore, internal coordinates and implicit-force constraint solvers are generally preferred.
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 formwhere 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
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
Joseph Louis Lagrange
Joseph-Louis Lagrange , born Giuseppe Lodovico Lagrangia, was a mathematician and astronomer, who was born in Turin, Piedmont, lived part of his life in Prussia and part in France, making significant contributions to all fields of analysis, to number theory, and to classical and celestial mechanics...
, 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
Lagrange multipliers
In mathematical optimization, the method of Lagrange multipliers provides a strategy for finding the maxima and minima of a function subject to constraints.For instance , consider the optimization problem...
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
Molecular dynamics
Molecular dynamics is a computer simulation of physical movements of atoms and molecules. The atoms and molecules are allowed to interact for a period of time, giving a view of the motion of the atoms...
simulation, constraints are enforced using the method of Lagrange multipliers
Lagrange multipliers
In mathematical optimization, the method of Lagrange multipliers provides a strategy for finding the maxima and minima of a function subject to constraints.For instance , consider the optimization problem...
. Given a set of n linear (holonomic
Holonomic
In mathematics and physics, the term holonomic may occur with several different meanings.-Holonomic basis:A holonomic basis for a manifold is a set of basis vectors ek for which all Lie derivatives vanish:[e_j,e_k]=0 \,...
) constraints at the time t,
where
and are the positions of the two particles involved in the kth constraint at the time t and is the prescribed inter-particle distance.These constraint equations, are added to the potential energy function in the equations of motion, resulting in, for each of the N particles in the system
Adding the constraint equations to the potential does not change it, since all
should, ideally, be zero.Integrating both sides of the equations of motion twice in time yields the constrained particle positions at the time
where
is the unconstrained (or uncorrected) position of the ith particle after integrating the unconstrained equations of motion.To satisfy the constraints
in the next timestep, the Lagrange multipliersLagrange multipliers
In mathematical optimization, the method of Lagrange multipliers provides a strategy for finding the maxima and minima of a function subject to constraints.For instance , consider the optimization problem...
must be chosen such that
This implies solving a system of
non-linear equationssimultaneously for the
unknown Lagrange multipliers .This system of
non-linear equations in unknowns is best solved using Newton's methodNewton's method
In numerical analysis, Newton's method , named after Isaac Newton and Joseph Raphson, is a method for finding successively better approximations to the roots of a real-valued function. The algorithm is first in the class of Householder's methods, succeeded by Halley's method...
where the solution vector
is updated usingwhere
is the JacobianJacobian
In vector calculus, the Jacobian matrix is the matrix of all first-order partial derivatives of a vector- or scalar-valued function with respect to another vector. Suppose F : Rn → Rm is a function from Euclidean n-space to Euclidean m-space...
of the equations σk:
Since not all particles are involved in all constraints,
is blockwise-diagonal and can be solved blockwise, i.e. molecule for molecule.Furthermore, instead of constantly updating the vector
, the iteration is startedwith , resulting in simpler expressions for and . After each iteration, the unconstrained particle positions are updated using.The vector is then reset to
This is repeated until the constraint equations
are satisfied up to a prescribed tolerance.Although there are a number of algorithms to compute the Lagrange multipliers, they differ only in how they solve the system of equations, usually using Quasi-Newton methods.
The SETTLE algorithm
The SETTLE algorithm solves the system of non-linear equations analytically for 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 modelWater model
In computational chemistry, classical water models are used for the simulation of water clusters, liquid water, and aqueous solutions with explicit solvent. These models use the approximations of molecular mechanics...
s).
The SHAKE algorithm
The SHAKE algorithm was the first algorithm developed to satisfy bond geometry constraints during molecular dynamics simulations.It solves the system of non-linear constraint equations using the Gauss-Seidel method
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 linear system of equations. It is named after the German mathematicians Carl Friedrich Gauss and Philipp Ludwig von Seidel, and...
to approximate the solution of the linear system of equations
in the Newton iteration. This amounts to assuming that
is diagonally dominant and solving the th equation only for the unknown. In practice, we compute |
 |
 |
 |
 |
 |
 |
 |
 |
for all
iteratively until the constraint equations are solved to a given tolerance.Each iteration of the SHAKE algorithm costs
operations and the iterations themselves converge linearly.A noniterative form of SHAKE was developed later.
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
Lagrange multipliers
In mathematical optimization, the method of Lagrange multipliers provides a strategy for finding the maxima and minima of a function subject to constraints.For instance , consider the optimization problem...
which are computed using the Gauss-Seidel method
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 linear system of equations. It is named after the German mathematicians Carl Friedrich Gauss and Philipp Ludwig von Seidel, and...
.
The original SHAKE algorithm is limited to mechanical systems with a tree structure, i.e., no closed loops of constraints. A later extension of the method, QSHAKE (Quaternion
Quaternion
In mathematics, the quaternions are a number system that extends the complex numbers. They were first described by Irish mathematician Sir William Rowan Hamilton in 1843 and applied to mechanics in three-dimensional space...
SHAKE) was developed to amend this. It works satisfactorily for rigid loops such as aromatic ring systems but fails for flexible loops, such as when a protein has a disulfide bond.
Further extensions include RATTLE, WIGGLE and MSHAKE.
RATTLE works the same way as SHAKE, yet using the Velocity Verlet time integration scheme. WIGGLE extends SHAKE and RATTLE by using an initial estimate for the Lagrange multipliers
Lagrange multipliers
In mathematical optimization, the method of Lagrange multipliers provides a strategy for finding the maxima and minima of a function subject to constraints.For instance , consider the optimization problem...
based on the particle velocities. Finally, MSHAKE computes corrections on the constraint forces, achieving better convergence.A final modification is the P-SHAKE algorithm for rigid or semi-rigid molecules. P-SHAKE computes and updates a pre-conditioner which is applied to the constraint gradients before the SHAKE iteration, causing the Jacobian
to become diagonal or strongly diagonally dominant. The thus de-coupled constraints converge much faster (quadratically as opposed to linearly) at a cost of .The M-SHAKE algorithm
The M-SHAKE algorithm solves the non-linear system of equations using Newton's methodNewton's method
In numerical analysis, Newton's method , named after Isaac Newton and Joseph Raphson, is a method for finding successively better approximations to the roots of a real-valued function. The algorithm is first in the class of Householder's methods, succeeded by Halley's method...
directly. In each iteration, the linear system of equations
is solved exactly using an LU decomposition
LU decomposition
In linear algebra, LU decomposition is a matrix decomposition which writes a matrix as the product of a lower triangular matrix and an upper triangular matrix. The product sometimes includes a permutation matrix as well. This decomposition is used in numerical analysis to solve systems of linear...
. Each iteration costs
operations, yet the solution converges quadratically, requiring fewer iterations than SHAKE.This solution was first proposed in 1986 by Ciccotti
Giovanni Ciccotti
Giovanni Ciccotti is an Italian physicist.Ciccotti holds the position of Professor of the Structure of Matter at the University of Rome La Sapienza. He's the author of more than a hundred articles on molecular dynamics and statistical mechanics. He worked with J.P. Ryckaert on new methods for...
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
directly, yet doing so only once, in the first iteration. The first iteration then costs operations, whereas the following iterations cost only operations (for the matrix-vector multiplication). This improvement comes at a cost though, since the Jacobian is no longer updated, convergence is only linear, 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..
The 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
Lagrange multipliers
In mathematical optimization, the method of Lagrange multipliers provides a strategy for finding the maxima and minima of a function subject to constraints.For instance , consider the optimization problem...
to the constraint forces and solves for the multipliers by using a series expansion to approximate the inverse of the Jacobian
:in each step of the Newton iteration. This approximation only works for matrices with Eigenvalues 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 LagrangeJoseph Louis Lagrange
Joseph-Louis Lagrange , born Giuseppe Lodovico Lagrangia, was a mathematician and astronomer, who was born in Turin, Piedmont, lived part of his life in Prussia and part in France, making significant contributions to all fields of analysis, to number theory, and to classical and celestial mechanics...
, and result in Lagrange equations of the mixed type.
See also
- Molecular dynamicsMolecular dynamicsMolecular dynamics is a computer simulation of physical movements of atoms and molecules. The atoms and molecules are allowed to interact for a period of time, giving a view of the motion of the atoms...
- Software for molecular mechanics modeling