Constraint algorithm

From Wikipedia, the free 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.

Contents

[edit] Mathematical background

The motion of a set of N particles can be described by a set of second-order ordinary differential equations, Newton's third law, which can be written in matrix form

\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}(\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 third law.

This problem was studied in detail by Joseph Louis Lagrange, who laid out most of the methods for solving it.[1] 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 a 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.

[edit] 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 two-fold: 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.[2][3]

Efficient recursive, internal-coordinate constraint solvers were extended to molecular dynamics.[4][5] Analogous methods were applied later to other systems.[6][7][8]

[edit] Implicit constraint force methods

With constraints, the system cannot wander as it chooses through its coordinate space; rather, its coordinates q must lie on the constraint manifold of allowed conformations. The implicit-force methods enforce this by computing the coordinate changes necessary to keep the system from straying; the methods are implicit because the constraint forces themselves are not explicitly computed, nor is there functional form given in advance. The standard approach for calculating such changes is the Lagrange-multiplier method; however, other methods return a system to its constraint manifold by projecting it there, i.e., by finding the point on the constraint manifold closest to the errant system's present coordinates.

[edit] SETTLE

For three atoms satisfying three constraints, it is trivial to solve the constraint equations analytically. The SETTLE algorithm[9] is commonly used to satisfy such sets of constraints, which commonly occur in simulations of three-atom molecules such as water.

[edit] SHAKE

The SHAKE algortihm was the first algorithm developed to satisfy bond geometry constraints during molecular dynamics simulations.[10]

A noniterative form of SHAKE was developed later.[11]

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.[12] 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.[13]

[edit] LINCS

An alternative constraint method, LINCS (Linear Constraint Solver) was developed in 1997,[14] but was based on earlier method, EEM,[15] and a modification thereof.[16] LINCS has been reported to be 3-4 times faster than SHAKE.[14]

The LINCS algorithm adopts a Lagrange-multiplier approach. Denoting the vector of Lagrange multipliers as \boldsymbol\lambda, the modified differential equations can be written

\mathbf{M} \cdot \frac{d^{2}\mathbf{q}}{dt^{2}} = \frac{\partial}{\partial \mathbf{q}} \left( -V + \boldsymbol\lambda \cdot \mathbf{g} \right) = \mathbf{f} + \mathbf{B^{T}} \cdot \boldsymbol\lambda

where the matrix B contains the derivatives of the constraint equations

\mathbf{B}_{mn} = \frac{\partial g_{m}}{\partial q_{n}}

Since the constraint equations are time-independent and zero, two additional equations may be derived

\frac{d\mathbf{g}}{dt} = \mathbf{B} \cdot \frac{d\mathbf{q}}{dt} = 0
\frac{d^{2}\mathbf{g}}{dt^{2}} = \mathbf{B} \cdot \frac{d^{2}\mathbf{q}}{dt^{2}} + \frac{d\mathbf{B}}{dt} \cdot \frac{d\mathbf{q}}{dt} = 0

The latter equation may be exploited by multiplying the above Lagrange-multiplier equation by B M-1, yielding

\mathbf{B} \cdot \frac{d^{2}\mathbf{q}}{dt^{2}} = - \frac{d\mathbf{B}}{dt} \cdot \frac{d\mathbf{q}}{dt} = \mathbf{B} \mathbf{M}^{-1} \mathbf{f} + \mathbf{B} \mathbf{M}^{-1} \mathbf{B^{T}} \cdot \boldsymbol\lambda

Solving this equation for the Lagrange multipliers \boldsymbol\lambda and substitution into the original Lagrange-multiplier equation yields the final equation

\frac{d^{2}\mathbf{q}}{dt^{2}} = \left(\mathbf{I} - \mathbf{TB} \right) \mathbf{M}^{-1} \mathbf{f} - \mathbf{T} \frac{d\mathbf{B}}{dt} \cdot \frac{d\mathbf{q}}{dt}

where I is the identity matrix and the matrix T is defined

\mathbf{T} = \mathbf{M}^{-1} \mathbf{B^{T}} \left( \mathbf{B} \mathbf{M}^{-1} \mathbf{B^{T}} \right)^{-1}

[edit] 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.[17][18] This approach was pioneered by Lagrange,[1] and result in Lagrange equations of the mixed type.[19]

[edit] See also

[edit] References and footnotes

  1. ^ a b Laplace, PS (1788). Mécanique analytique. 
  2. ^ Noguti T; Gō N (1983). "A Method of Rapid Calculation of a 2nd Derivative Matrix of Conformational Energy for Large Molecules". Journal of the Physical Society of Japan 52: 3685–3690. 
  3. ^ Abe, H; Braun W, Noguti T, Gō N (1984). "Rapid Calculation of 1st and 2nd Derivatives of Conformational Energy with respect to Dihedral Angles for Proteins: General Recurrent Equations". Computers and Chemistry 8: 239–247. 
  4. ^ Bae, D-S; Haug EJ (1988). "A Recursive Formulation for Constrained Mechanical System Dynamics: Part I. Open Loop Systems". Mechanics of Structures and Machines 15: 359–382. 
  5. ^ Jain, A; Vaidehi N, Rodriguez G (1993). "A Fast Recursive Algorithm for Molecular Dynamics Simulation". Journal of Computational Physics 106: 258–268. 
  6. ^ Rice, LM; Brünger AT (1994). "Torsion Angle Dynamics: Reduced Variable Conformational Sampling Enhances Crystallographic Structure Refinement". Proteins: Structure, Function, and Genetics 19: 277–290. 
  7. ^ Mathiowetz, AM; Jain A, Karasawa N, Goddard III, WA (1994). "Protein Simulations Using Techniques Suitable for Very Large Systems: The Cell Multipole Method for Nonbond Interactions and the Newton-Euler Inverse Mass Operator Method for Internal Coordinate Dynamics". Proteins: Structure, Function, and Genetics 20: 227–247. 
  8. ^ Mazur, AK (1997). "Quasi-Hamiltonian Equations of Motion for Internal Coordinate Molecular Dynamics of Polymers". Journal of Computational Chemistry 18: 1354–1364. 
  9. ^ Miyamoto, S; Kollman PA (1992). "SETTLE: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models". Journal of Computational Chemistry 13: 952–962. 
  10. ^ Ryckaert, J-P; Ciccotti G, Berendsen HJC (1977). "Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes". Journal of Computational Physics 23: 327–341. 
  11. ^ Yoneya, M; Berendsen HJC, Hirasawa K. "A Noniterative Matrix Method for Constraint Molecular-Dynamics Simulations". Molecular Simulations 13: 395–405. 
  12. ^ Forester, TR; Smith W (1998). "SHAKE, Rattle, and Roll: Efficient Constraint Algorithms for Linked Rigid Bodies". Journal of Computational Chemistry 19: 102–111. 
  13. ^ McBride, C; Wilson MR, Howard JAK (1998). "Molecular dynamics simulations of liquid crystal phases using atomistic potentials". Molecular Physics 93: 955–964. 
  14. ^ a b Hess, B; Bekker H, Berendsen HJC, Fraaije JGEM (1997). "LINCS: A Linear Constraint Solver for Molecular Simulations". Journal of Computational Chemistry 18: 1463–1472. 
  15. ^ Edberg, R; Evans DJ, Morriss GP (1986). "Constrained Molecular-Dynamics Simulations of Liquid Alkanes with a New Algorithm". Journal of Chemical Physics 84: 6933–6939. 
  16. ^ Baranyai, A; Evans DJ (1990). "New Algorithm for Constrained Molecular-Dynamics Simulation of Liquid Benzene and Naphthalene". Molecular Physics 70: 53–63. 
  17. ^ Bae, D-S; Haug EJ (1988). "A Recursive Formulation for Constrained Mechanical System Dynamics: Part II. Closed Loop Systems". Mechanics of Structures and Machines 15: 481–506. 
  18. ^ Rodriguez, G; Jain A, Kreutz-Delgado K (1991). "A Spatial Operator Algebra for Manipulator Modeling and Control". The International Journal for Robotics Research 10: 371–381. 
  19. ^ Sommerfeld, Arnold (1952). Lectures on Theoretical Physics, Vol. I: Mechanics. New York: Academic Press. ISBN-10 0126546703. 

Bold text


CONTRAINTS ARE SOME THING THAT CAN LIMIT A MOTION. IT CAN BE FRAME OF REFRERENCE,OR ANY OTHER THING. THEY ARE MAINLY 4 TYPES OF CONTRAINTS

                1) HOLONOMIC CONTRAINTS
                2) NON-HOLONOMIC CONTRAINTS
                3) RHEONOMOUS CONTRAINTS
                4) SCLERONOMOUS CONTRAINTS
    IN THIS HOLONOMIC CONTRAINTS ARE THOSE ONE THAT DOES NOT DEPEND                                      ON VELOCITY. Eg:RIGID BODY 
       IT CAN BE RESPRESENTED BY THE EQUATION
                   f(r1,r2,........rn,t)=0

IN THE CASE OF RIGID BIDIES WE CAN RESPRESENT IT AS r < math > Insertformulahere

     IN NON HOLO                 IT DEPEND ON VELOCITY