Constraint Algorithm

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.

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{M} \cdot rac{d^{2}\mathbf{q}}{dt^{2}} = \mathbf{f} = -rac{\partial V}{\partial \mathbf{q}}

where M is a ''mass matrix'' and '''q''' is the Vector of Generalized Coordinate 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_{j}(\mathbf{q}) = 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 , 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.

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.² ³

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 Molecular Dynamics simulation, constraints are enforced using the method of Lagrange Multipliers . Given a set of

n

linear (holonomic) constraints at the time

t

Adding the constraint equations to the potential does not change it, since all

\sigma_k^{(t)}

should, ideally, be zero.

Integrating both sides of the equations of motion twice in time yields the constrained particle positions at the time

t + \Delta t

\mathbf x^{(t + \Delta t)}_i = \hat{\mathbf x}^{(t + \Delta t)}_i + \sum_{k=1}^n \lambda_k rac{\partial\sigma_k^{(t)}}{\partial \mathbf x_i}\left(\Delta t
ight)^2m_i^{-1}, \quad i=1 \dots N

where

\hat{\mathbf x}^{(t + \Delta t)}_i

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 Lagrange Multipliers must be chosen such that

:

\hat{\mathbf x}_i^{(t+\Delta t)} \leftarrow \hat{\mathbf x}_i^{(t+\Delta t)} + \sum_{k=1}^n \lambda_k^{(l)}rac{\partial \sigma_k}{\partial \mathbf x_i}

.

The vector is then reset to

:

\underline{\lambda}^{(l)} = \mathbf 0

.

This is repeated until the constraint equations

\sigma_k^{(t+\Delta t)}

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

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 Water Models ).

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 to approximate the solution of the linear system of equations

:

\lambda = -\mathbf J_\sigma^{-1} \underline{\sigma}

in the Newton Iteration . 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

	CATEGORIES ABOUT CONSTRAINT ALGORITHM

	molecular dynamics

	computational chemistry

	molecular physics

	computational physics

	numerical differential equations

	APPAREL
	BABY
	BEAUTY
	BOOKS
	CAR TOYS
	CELL PHONES
	DVD'S
	ELECTRONICS
	GOURMET FOOD
	GROCERIES
	HEALTH & PERSONAL
	HOME & GARDEN
	JEWELRY
	MUSIC
	MUSIC INSTRUMENTS
	OFFICE PRODUCTS
	SOFTWARE
	SPORTING GOODS
	TOOLS & HARDWARE
	TOYS
	VIDEO GAMES
	SHOPPING HOME

	MORE SHOPPING...

Information About

Constraint Algorithm