temperature control they do have the same e¬ect as a Langevin thermostat

on the variance of velocities (i.e., on the temperature). The idea is to rescale

200 Molecular dynamics

velocities per step in such a way that the total temperature T of the system

will decay with a ¬rst-order process to the desired temperature T0 :

T0 ’ T

dT

= . (6.116)

dt „

This rate equation would cause a temperature deviation from T0 to decay

exponentially to zero with a time constant „ . The implementation in terms

of a scaling factor » for the velocities is given by

”t T0

’1 ,

»2 = 1 + (6.117)

„ T

where T is given by the kinetic energy found after updating the velocities

in a normal dynamic step. For the smallest possible value of time con-

stant „ = ”t the scaling is complete and the temperature is exactly con-

served. This corresponds to the Gauss isokinetic thermostat which produces

a canonical ensemble in con¬guration space. For „ much longer than the

intrinsic correlation times for internal exchange of energy, the scaling has no

e¬ect and a microcanonical ensemble is obtained. This is borne out by the

¬‚uctuations in kinetic and potential energy: for small „ the kinetic energy

does not ¬‚uctuate but the potential energy does; as „ increases, ¬‚uctuations

in kinetic energy appear at the expense of potential energy ¬‚uctuations, to

become equal and opposite at large „ . The cross-over occurs roughly from a

factor of 10 below to a factor of 10 above the intrinsic relaxation time for the

exchange between kinetic and potential energy, which is system-dependent.

Morishita (2000), in an attempt to characterize the weak-coupling ensem-

ble, derived the following equation for the compressibility κ:

22 dV (q)

β δK + O(N ’1 )

κ= β’ , (6.118)

n dt

plus some unknown function of p. Here δK is the ¬‚uctuation of K, which

depends on „ . For small „ , when δK = 0, this form reduces to βdV /dt,

yielding the canonical con¬guration space distribution derived above for the

isokinetic thermostat. For large „ , when δK = ’δV , the con¬guration space

distribution tends to

12

f (q) = exp ’βV (q) ’ β (δV )2 , (6.119)

n

which equals the microcanonical con¬guration space distribution already

derived earlier by Nos´ (1984a). For intermediate values of „ Morishita

e

made the assumption that the ¬‚uctuation of K is related to that of V by

6.5 Controlling the system 201

δK = ’±δV , with ± depending on „ , and arrives at a con¬guration space

distribution

±

f (q) = exp ’βV (q) ’ β 2 (δV )2 . (6.120)

n

The distribution in momentum space remains unknown, but is less important

as the integration over canonical momenta can be carried out separately.

Note again that this is valid for cartesian coordinates only. Also note that

the weak-coupling algorithm is no longer time-reversible, unless in the two

extreme cases „ = ”t and „ ’ ∞.

Pressure control by weak coupling is possible by scaling coordinates. In

the spirit of weak coupling one attempts to regulate the pressure P according

to the ¬rst-order equation

dP 1

= (P0 ’ P ). (6.121)

dt „p

Assume the isothermal compressibility βT = ’(1/V )‚V /‚P is known, then

scaling coordinates and volume,

r = χr, (6.122)

V = χ3 V, (6.123)

every time step with a scaling factor χ, given by

”t

χ3 = 1 ’ βT (P0 ’ P ), (6.124)

„p

will accomplish that task. As the compressibility only enters the algorithm

in conjunction with the time constant, its value need not be precisely known.

The weak pressure coupling has the advantage of smooth response, but the

disadvantages are that it does not generate a known ensemble and ¬‚uctua-

tions cannot be used.

6.5.4 Extended system dynamics

The idea to extend the system with an extra degree of freedom that can be

used to control a variable in the system, was introduced by Nos´ (1984a,

e

1984b) for temperature control. The method involved a somewhat inconve-

nient time scaling and was modi¬ed by Hoover (1985) into a scheme known

as the Nos´“Hoover thermostat, which has been widely used since. We treat

e

the thermostat case, but the same principle can be applied to control pres-

sure in a barostat. An extra variable · is introduced, which is a factor

scaling the velocities. It has an associated “momentum” p· = Q·, where Q

™

202 Molecular dynamics

is the “mass” of the extra degree of freedom. The equations of motion are

(p, q, m stand for all pi , qi , mi ):

p

q=

™ , (6.125)

m

p·

p = F (q) ’ p ,

™ (6.126)

Q

2

pi

’ nkB T.

p· =

™ (6.127)

2mi

The temperature deviation from the bath temperature drives the time deriv-

ative of the velocity scaling factor, rather than the scaling factor itself, as is

the case in the weak-coupling method. This makes the equations of motion

time reversible again, and allows to compute the phase space distribution.

On the other hand, it has the practical disadvantage that the temperature

control is now a second-order di¬erential equation in time, which leads to

oscillatory approaches to equilibrium. Hoover shows that the equilibrium

phase space distribution is given by

p2 1

f (q, p, ·, p· ) ∝ exp ’β + Q· 2

i

V (q) + , (6.128)

2mi 2

i

which is canonical. The extra variable is statistically independent of posi-

tions and velocities.

The Nos´“Hoover thermostat has been criticized because its behavior is

e

non-ergodic (Toxvaerd and Olsen, 1990), which led Martyna et al. (1992) to

formulation of the Nos´“Hoover chain thermostat. In this thermostat there

e

is a sequence of M additional variables ·1 , ·2 , . . . , ·M with their masses and

conjugate momenta, each scaling its predecessor in the chain:

p

q=

™ , (6.129)

m

p·

F (q) ’ p 1 ,

p=

™ (6.130)

Q1

p·1

·1 =

™ , (6.131)

Q1