= 2mγkB T (”t). (6.102)

In fact, it is not required to draw the change in p from a normal distribution,

as long as the distribution has zero mean and ¬nite variance. We now show

that this procedure yields extra terms in ‚f /‚t that force the distribution

functions of the momenta to which the noise and friction are applied to a

normal distribution. The random force on p causes a di¬usion of p with

a di¬usion constant D given by the mean-square displacement of p (one

dimension) in the time interval ”t:

(”p)2 = 2D”t. (6.103)

Consequently the di¬usion constant is given by

D = mγkB T. (6.104)

Di¬usion leads to Fick™s equation for the distribution function:

‚ 2 f (p, t)

‚f (p, t)

=D (6.105)

‚p2

‚t

and the friction term leads to an extra ¬‚ux f p and therefore to

™

‚f ‚(f p)

™ ‚f

=’ = γf + γp . (6.106)

‚t ‚p ‚p

If these two contributions cancel each other, the distribution function will

be stationary. The solution of the equation

‚2f ‚f

D 2 + γf + γp =0 (6.107)

‚p ‚p

is

p2

f (p) ∝ exp ’ . (6.108)

2mkB T

The total distribution function must be proportional to this normal distri-

bution of pi , and that applies to all momenta to

which random and friction forces are applied. Since H contains p2 /2m

as an additive contribution, and the Hamiltonian terms will force the sta-

tionary distribution function to be a function of H, it follows that the total

distribution will be:

H(z)

f (z) ∝ exp ’ . (6.109)

kB T

198 Molecular dynamics

Thus the distribution will be canonical and the temperature of the ensem-

ble is set by the relation between applied friction and noise. In a sense

the friction and noise driving the canonical distribution compete with ac-

cumulating disturbing errors in the numerical solution of the Hamiltonian

equations. The applied damping enforces an extra ¬rst-order decay on ve-

locity correlation functions and thereby disturb the dynamic behavior on

time scales comparable to 1/γ. To minimize dynamic disturbance, the γ™s

should be taken as small as possible, and need not be applied to all degrees

of freedom. But with large errors a strong damping on many particles is

required. The Langevin thermostat provides a smooth decay to the desired

temperature with ¬rst-order kinetics. Note that we have given the equa-

tions for cartesian coordinates, for which the mass tensor is diagonal and

constant; for generalized coordinates the equations become complicated and

di¬cult to handle because the mass tensor is a function of the con¬guration.

The velocity rescaling thermostat of Andersen (1980) has a similar ef-

fect, but does not distribute the external disturbance as smoothly as the

Langevin thermostat. Andersen™s method consist of reassigning the velocity

of a randomly selected molecule at certain time intervals from the appropri-

ate Maxwellian distribution. As expected, a canonical ensemble is generated.

6.5.2 Strong-coupling methods

The strong-coupling methods arti¬cially constrain a property to the desired

value, e.g., the total kinetic energy to a prescribed value determined by the

desired temperature. This is accomplished by scaling the velocities with a

multiplicative factor that preserves the shape of the distribution function.

This amounts to adding an acceleration v = ’±v t every degree of freedom.

™

This isokinetic or Gauss thermostat was introduced by Hoover et al. (1982)

and Evans (1983). Following Tuckerman et al. (1999), who use the notion

of phase-space compressibility (see Chapter 17, Section 17.8, page 495), we

shall show that the isokinetic thermostat produces a canonical distribution

in con¬guration space.

Start with the de¬nition of phase-space compressibility κ for non-Hamilton-

ian ¬‚ow (17.165):

κ = ∇ · z.

™ (6.110)

Tuckerman et al. showed (see page 495) that, if a function w(z) can be

de¬ned whose time derivative equals κ, then exp[’w(z)]dz1 , . . . , dz2n is an

invariant volume element along the trajectory, meaning that exp[’w(z)] is

6.5 Controlling the system 199

the equilibrium weight function in phase space. Now κ follows from ±:

n

‚

pi = ’n±,

κ= ™ (6.111)

‚pi

i=1

where n is the number of degrees of freedom (= 3N ’ nc ). We can express

± in phase-space functions by realizing that (in cartesian coordinates)

n n

p2

d pi pi

™

i

= = 0, (6.112)

dt 2mi 2mi

i=1 i=1

and therefore

n n

p2

pi Fi

’± i

= 0. (6.113)

mi mi

i=1 i=1

Hence

i pi Fi /mi pi Fi dV (q)

κ = ’n = ’β =β , (6.114)

2 mi dt

i pi /mi i

since

dV (q) ‚V Fi pi

qi = ’

= ™ .

dt ‚qi mi

i i

Thus the w = βV and the weight function equals exp[’βV ]. We conclude

that ensemble averages of a variable A(z) over the isokinetic thermostat are

given by

A(z) exp[’βV (q)]δ( i p2 /2mi ’ nkB T ) dz

i

A= . (6.115)

2 /2m ’ nk T ) dz

exp[’βV (q)]δ( i pi i B

Strong coupling can also be applied to constrain the pressure to a preset

value (Evans and Morriss, 1983a; 1984). Combined with the isokinetic ther-

mostat Evans and Morriss (1983b) have shown that a NPT ensemble (on the

hypersurface of constrained pressure and kinetic energy) is obtained with a

weight factor equal to exp[’βH], where H is the enthalpy U + pV . We

do not pursue this type of pressure control further as it seems not to have

become common practice, but refer instead to Allen and Tildesley (1987).

6.5.3 Weak-coupling methods

Weak-coupling methods (Berendsen et al., 1984) are not stochastic in na-