and then corrects the positions with ”r i such that

σs (r + (”r) = 0, s = 1, . . . , m, (15.100)

where

(”t)2 ‚σs (r(t))

”r i = »s (t) . (15.101)

mi ‚r i

s

These equations represent a set of m (generally non-linear) coupled equa-

tions for the m »™s, which can be solved in several ways, but as a result

of the nonlinear character always requiring iteration. They can be either

linearized and then solved as a set of linear equations, or the constraints

can be solved sequentially and the whole procedure iterated to convergence.

The latter method is easy to implement and is used in the routine SHAKE

(Ryckaert et al., 1977).

Let us illustrate how one distance constraint between particles i and j:

σ = r 2 ’ d2 = 0 (15.102)

ij

will be reset in a partial iteration step of SHAKE. See Fig. 15.4. The particle

positions are ¬rst displaced to r . Because ‚σ/‚r i = ’‚σ/‚r j = 2r ij , the

displacements must be in the direction of r ij (t) and proportional to the

inverse mass of each particle:

2(”t)2

”r i = » r ij (t),

mi

2(”t)2

=’

”r j »r ij (t). (15.103)

mj

The variable » is determined such that the distance between r i + ”r i and

r j + ”r j is equal to d. This procedure is repeated for all constraints until

all constraints have converged within a given tolerance.

It is illustrative to consider the motion of a particle that is constrained to

15.8 Holonomic constraints 421

r(t’”t)

r(t)

r´

r(t+”t)

Figure 15.5 The action of SHAKE for a rotating particle constrained to move on a

circle.

move on a circle by a distance constraint to a ¬xed origin (Fig. 15.5). There

is no external force acting on the particle, so it should move with constant

angular velocity. The positions at times t ’ ”t and t are given. According

to (15.99), the position is ¬rst linearly extrapolated to r . Subsequently the

position is reset in the direction of the constraint r(t) until the constraint

is satis¬ed. Thus r(t + ”t) is obtained. It is easily seen that this algorithm

gives exact results up to an angular displacement per step of close to 90—¦

(four steps per period), beyond which the algorithm is unstable and fails to

¬nd a solution. It is also seen that resetting in the direction of the constraint

at time t is correct.

15.8.3 Projection methods

Consider a cartesian system of point masses with equations of motion

nc

‚σs

Fu

¨

mi r i = + , i = 1, . . . , N, (15.104)

i

‚r i

s=1

which we shall write in matrix notation as

M¨ = f + CT », (15.105)

x

where we use x for the 3N — 1 column matrix (x1 , y1 , z1 , x2 , . . . , yN , zN )T ,

similarly f for F u , M is the 3N — 3N diagonal matrix of masses, and the

constraint matrix C is de¬ned by

‚σs

Csi = . (15.106)

‚xi

422 Lagrangian and Hamiltonian mechanics

By taking the time derivative of

™

σs = (Cx)s = 0

™ (15.107)

the following relation is found:

™™

C¨ = ’Cx. (15.108)

x

By left-multiplying (15.105) ¬rst by M’1 and then by C, and substituting

(15.108), » can be solved and we obtain

» = ’(CM’1 CT )’1 (CM’1 f + Cx).

™™ (15.109)

The matrix CM’1 CT is non-singular and can be inverted if the constraints

are independent.19 Substituting (15.109) into (15.105) we obtain the follow-

ing equation of motion for x:

x = (1 ’ TC)M’1 f ’ TCx,

™™

¨ (15.110)

where

T = M’1 CT(CM’1 CT )’1 .

def

(15.111)

The matrix 1 ’ TC projects the accelerations due to the unconstrained

forces onto the constraint hypersurface. The ¬rst term in (15.110) gives

the constrained accelerations due to the systematic forces (derivatives of the

potential) and the second term gives the constrained accelerations due to

centripetal forces.

Equation (15.110) contains the velocities at the same time as when the

forces are evaluated, but these velocities are not known at that time. There-

fore this equation is in principle implicit and needs an iterative solution. In

practice this is not done and not necessary. Hess et al. (1997) show that

the velocities at the previous half step can be used in conjunction with ap-

propriate corrections. The corrections are made in such a way that a stable

algorithm results without drift.

Although the SHAKE algorithm of Ryckaert et al. (1977) is easier to

implement, the LINCS (LINear Constraint Solver) algorithm of Hess et al.

(1997) is faster, more robust, more accurate and more suitable for parallel

computers. For special cases it can be advantageous to solve the equations

analytically, as in the SETTLE algorithm of Miyamoto and Kollman (1992)

for water molecules.

19 De Leeuw et al. (1990) show that the matrix is not only non-singular but also positive de¬nite.

16

Review of thermodynamics

16.1 Introduction and history

This book is not a textbook on thermodynamics or statistical mechanics.

The reason to incorporate these topics nevertheless is to establish a common

frame of reference for the readers of this book, including a common nomen-

clature and notation. For details, commentaries, proofs and discussions, the

reader is referred to any of the numerous textbooks on these topics.