2 2 2 2

the constraint q0 + q1 + q2 + q3 = 1, which is usually compensated by regular

scaling of the q™s.

15.8 Holonomic constraints

Holonomic constraints depend only on coordinates and can be described by a

constraint equation σ(r) = 0 that should be satis¬ed at all times. For every

constraint there is such an equation. Examples are (we use the notation

r ij = r i ’ r j ):

• distance constraint between two particles: |r 12 |’d12 = 0, or, alternatively,

(r 12 )2 ’ d2 = 0;

12

• angle 1-2-3 constraint between two constrained bonds: r 12 · r 32 ’ c = 0,

where c = d12 d32 cos φ, or, alternatively, r 2 ’ d2 = 0.

13 13

The way to introduce holonomic constraints into the equations of motion is

by minimizing the action while preserving the constraints, using Lagrange

multipliers. Thus we add each of the m constraints σs (q), s = 1, . . . , m with

their undetermined Lagrange multipliers »s to the Lagrangian and minimize

the action. As can be easily veri¬ed, this results in the modi¬ed Lagrange

equations (compare (15.2))

d ‚L ‚L

’ =0 i = 1, . . . , n, (15.91)

dt ‚ qk

™ ‚qk

where

m

def

L = L+ »s σs (q), (15.92)

s=1

while for all q along the path

σs (q) = 0, s = 1, . . . , m (15.93)

Equations (15.91) and (15.93) fully determine the path, i.e., both q(t) and

»(t). In fact the path is restricted to a hypersurface determined by the con-

straint equations. There are n+m variables (the n q™s and the m »™s) and an

418 Lagrangian and Hamiltonian mechanics

equal number of equations (n Lagrange equations and m constraint equa-

tions). Note that the generalized momenta are not modi¬ed by holonomic

constraints17 because the constraints are not functions of q, but the forces

™

are. The total generalized force is built up from an unconstrained force and

a constraint force:

m

‚L ‚L ‚σs

= + »s . (15.94)

‚qk ‚qk ‚qk

s=1

If the constraints are not eliminated by the use of generalized coordinates,

the »™s must be solved from the constraint equations. We can distinguish

two ways to obtain the solution, both of which will be worked out in more

detail in the following subsections. The ¬rst method, which is historically

the oldest and in practice the most popular one, was devised by Ryckaert

et al. (1977). It resets the coordinates after an unconstrained time step,

so as to satisfy the constraints to within a given numerical precision, and

therefore prevents the propagation of errors. The method is most suitable

in conjunction with integration algorithms that do not contain explicit ve-

locities, although a velocity variant is also available (Andersen, 1983). The

second method rewrites the equations of motion to include the constraints

by solving the »™s from the fact that the σs ™s are zero at all times: hence

all time derivatives of σs are also zero. This allows an explicit solution for

the Lagrange multipliers, but the solutions contain the velocities, and since

only derivatives of the constraints appear, errors may propagate. We shall

call this class of solutions projection methods because in fact the accelera-

tions are projected onto the hypersurface of constraint. The ¬rst algorithm

using this method for molecules was published by Edberg et al. (1986);

the solution was cast in more general terms by de Leeuw et al. (1990) and

discussed in the context of various types of di¬erential equations by Bekker

(1996). A similar method was devised by Yoneya et al. (1994), and an e¬-

cient algorithm (LINCS) was published by Hess et al. (1997). We note that

matrix methods to solve for holonomic as well as non-holonomic (velocity-

dependent) constraints were already known in the ¬eld of macroscopic rigid

body dynamics.18

17 That is true for our de¬nition of constraints; de Leeuw et al. (1990) also consider a de¬nition

based on the constancy of the time derivative of σ, in which case the generalized momenta are

modi¬ed, but the ¬nal equations are the same.

18 See, e.g., Section 5.3 of Wittenburg (1977).

15.8 Holonomic constraints 419

15.8.1 Generalized coordinates

A straightforward, but often not feasible, method to implement constraints

is the use of generalized coordinates in such a way that the constraints are

themselves equivalent to generalized coordinates. Suppose that we transform

(r 1 . . . r N ) ’ (q , q ), (15.95)

where q = q1 , . . . , qn=3N ’m and q = qn+1 , . . . qn+m=3N are the free and

constrained coordinates, respectively. The constrained coordinates ful¬ll

the m constraint equations:

σs = qn+s ’ cs = 0. (15.96)

Because q = c, q = 0 and the kinetic energy does not contain q . The

™ ™

Lagrangian does also not depend on q as variables, but contains the c™s

only as ¬xed parameters.Thus the q do not ¬gure at all in the equations

of motion and the Lagrangian or Hamiltonian mechanics simply preserves

the constraints. The dynamics is equivalent to that produced by the use

of Lagrange multipliers. The di¬culty with this method is that only in

simple cases the equations of motion can be conveniently written in such

generalized coordinates.

15.8.2 Coordinate resetting

One popular method of solving the constraint equations (Ryckaert et al.,

1977) can be used in conjunction with the Verlet algorithm, usually applied

with Cartesian coordinates:

(”t)2 u

r i (t + ”t) = 2r i (t) ’ r i (t ’ ”t) + [F i (t) + F c (t)], (15.97)

i

mi

where F u are the forces disregarding the constraints, and the constraint

force on particle i at time t is given by

‚σs

F c (t) = »s (t) . (15.98)

i

‚r i

s

The e¬ect of the constraint force is to add a second contribution to the dis-

placement of the particles. The algorithm ¬rst computes the new positions

r i disregarding the constraints:

(”t)2 u

r i = 2r i (t) ’ r i (t ’ ”t) + F i (t), (15.99)

mi

420 Lagrangian and Hamiltonian mechanics

” rj

” ri

r™

rj (t + ” t) j

ri (t + ” t )

r™

i

d

d

ri ( t ) r tj ( )

Figure 15.4 Coordinate resetting to realize a bond length constraint.