ular dynamics (or Monte Carlo) simulations at M di¬erent temperatures

T1 , . . . , TM or inverse temperatures divided by kB , β1 , . . . , βM . Initially,

system Si is at temperature Ti , but we allow exchange between the tem-

peratures of two systems, so that in general system Si has temperature

Tm . We order the temperatures always sequentially: m = 1, 2, . . . , M ,

but the sequence {i} = i(1), i(2), . . . , i(M ) of the systems is a permuta-

39 We use the canonical probability in con¬guration space only, not in the full phase space.

206 Molecular dynamics

tion of the sequence 1, 2, . . . , M . A state X in the generalized ensemble

consists of M con¬gurations r i(1) , r i(2) , . . . , r i(M ) , with potential energies

Ei(1) , Ei(2) , . . . , Ei(M ) . Because there are no interactions between the sys-

tems and each temperature occurs exactly once, the probability of this state

is

w ∝ exp[’(β1 Ei(1) + β2 Ei(2) + · · · + βM Ei(M ) )]. (6.136)

Now consider a possible exchange between two systems, e.g., systems Si at

βm and Sj at βn . After the exchange, system Si will be at βn and Sj will be

at βm . The probabilities before and after the exchange must be, respectively,

wbefore ∝ exp[’(βm Ei + βn Ej )], (6.137)

and

wafter ∝ exp[’(βn Ei + βm Ej )]. (6.138)

Their ratio is

wafter

= e’” , (6.139)

wbefore

where

” = (βn ’ βm )(Ei ’ Ej ). (6.140)

The transition probabilities W’ (meaning from “˜before” to “after”) and

W← (meaning from “˜after” to “before”) must ful¬ll the detailed balance

condition:

wbefore W’ = wafter W← . (6.141)

Thus it follows that

W’

= e’” . (6.142)

W←

This is accomplished by the Metropolis acceptance criterion:

” ¤ 0 : W’ = 1

for (6.143)

” > 0 : W’ = e’” ,

for (6.144)

as is easily seen by considering the backward transition probability:

” < 0 : W← = e”

for (6.145)

” ≥ 0 : W← = 1,

for (6.146)

which ful¬lls (6.142).

Although exchange between any pair may be attempted, in practice only

neighboring temperatures yield non-zero acceptance ratios and the exchange

6.7 Applications of molecular dynamics 207

attempt can be limited to neighbors. An acceptance ratio of 20% is consid-

ered reasonable. One should choose the set of temperatures such that the

acceptance ratio is more or less uniform over the full temperature range.

This will depend on the system; Sugita and Okamoto (1999) ¬nd for a pep-

tide that an exponential distribution (equal ratios) is satisfactory. They use

eight temperatures between 200 and 700 K, but a higher number (10 to 20)

is recommended.

The exchange, once accepted, can be accomplished in various ways. In

Monte Carlo simulations, one exchanges both the random step sizes and the

β™s in the acceptance criterion. In dynamics, the least disturbing implemen-

tation is to use a weak-coupling or a Langevin thermostat and switch the

reference temperatures of the thermostats. The simulation should then ex-

tend over several time constants of the thermostat before another exchange

is attempted. Using an isokinetic thermostat, the velocities should be scaled

proportional to the square root of the temperature ratio upon exchange. In

principle, it is also possible to switch not only the thermostats, but also

all velocities of the particles between the two systems. This will drastically

break up the time correlation of concerted motions; it is not clear whether

this is advantageous or disadvantageous for the sampling e¬ciency.

6.7 Applications of molecular dynamics

Molecular dynamics simulations with atomic detail can be routinely applied

to systems containing up to a million particles over several nanoseconds.

The time step for stable simulations is determined by the highest frequen-

cies in the system; as a rule of thumb one may assume that at least ten,

but preferably 50 time steps should be taken within the shortest period of

oscillation.40 If the system contains mobile hydrogen atoms, bond oscilla-

tion periods may be as short as 10 fs; bond vibration involving heavy atoms

typically exceed 20 fs. When covalent bonds are constrained, the highest

frequencies are rotational and librational modes that involve hydrogen; di-

hedral angle rotations involving hydroxyl groups have similar periods. For

example, in liquid water librational frequencies up to 800 cm’1 (40 fs) occur.

A usual time step for simulations of systems containing liquid water is 1 to

2 fs when internal bond constraints are imposed; with further restrictions

of hydrogen motions, “hydrogen-rich” systems as hydrated proteins remain

stable with time steps up to 7 fs (Feenstra et al., 1999). In 2005, simulation

of a typical medium-sized protein (lysozyme) in water, totalling some 30 000

40 See Berendsen and van Gunsteren (1981). Mazur (1997) concludes that even less than 10 steps

per period su¬ce for the leap-frog algorithm.

208 Molecular dynamics

atoms, reached about 1 ns per day on a single state-of-the-art processor (van

der Spoel et al., 2005), and this performance is expected to increase by a

factor of ten every ¬ve years, according to Murphy™s law, which has been

followed quite closely over the last decades. The availability of massively

parallel clusters of processors allows simulation of much larger system sizes

and much longer time scales. Proteins can be followed over microseconds,

which is not yet su¬cient to simulate realistic folding processes and reli-

ably predict protein structures from sequence data. With the near future

peta¬‚op computers, the protein folding problem, which has been called the

Holy Grail of biophysics (Berendsen, 1998), is likely to be solved. In ma-

terial science, micron-size solids simulated for microseconds will become a

reality.

Figure 6.13 shows several snapshots of a simulation involving more than

a billion (109 ) particles by Abraham et al. (2002).41 The simulated system

is a crystal of Lennard“Jones particles, modeled to mimic a copper crystal,

which is subjected to external tension forces that cause a crack to increase

in size. The purpose of this simulation is to investigate the formation and

propagation of dislocations that characterize the crack and model the process

of work-hardening of metals. The system is a slab with 1008 atoms along

the three orthogonal sides. Two notches are centered midway along the x-

direction, at y = 0 and y = Ly , with a y-extension of 90 atomic layers which

extends through the entire thickness Lz . The exposed notch faces are in

the y ’ z planes with (110) faces, and the notch is pointed in the (1, ’1, 0)

direction. Periodic boundary conditions are imposed between the x’y faces

at z = 0 and z = Lz . This notched slab geometry has a total of 1 023 103 872

atoms. The total simulation time for this study is 200 000 time-steps or 2 ns.

The slab is initialized at zero temperature, and an outward strain of 4% is

imposed on the outermost columns of atoms de¬ning the opposing vertical

yz faces of the slab. The ¬gures show only atoms with a potential energy

less than 97 % of the bulk value magnitude.

In the ¬gures, we see a spaghetti-like network of atomic strings ¬‚ying

from the vertices of the two opposing crack edges. This is simply a large

number of mobile dislocations being created at each crack edge, rapidly

¬‚owing through the stretched solid in an erratic manner, and eventually col-

liding with intersecting dislocations from the opposite edge. For the simple

face-centered-cubic solid, dislocations are easily created at the apex of the

two microcracks where the stress is at a maximum and easily ¬‚ow through

41 The author is indebted to Dr Farid Abraham for providing the pictures in Fig. 6.13 and

the accompanying description. The interested reader is referred to Abraham (2003) for an

introductory text on cracks and defects, and ductile and brittle behavior of metals.

Exercises 209

the solid giving rise to the ductility of the solid. The simulation supports

the prevailing view that even though there may not be enough dislocations

originally present in a crystal to account for the extensive slip in a ductile

material (in this simulation there are initially no dislocations), their cre-

ation in vast amounts can occur at severe stress concentrations, such as at

crack tips, enabling a stressed solid to be rapidly ¬lled with dislocations and

giving rise to material deformation under a steady load. The ¬gures show

snapshots of the propagating dislocations and rigid junctions evolving into

a complex topology of the defect-solid landscape.

Colliding dislocations can cause permanent atomic relocation along a

line, called a rigid junction or sessile dislocation. A coarse-grain three-

dimensional skeleton of such sessile dislocations becomes apparent from a

distant view. They are obstacles to further dislocation mobility. If their

density is su¬ciently high, dislocation mobility becomes insigni¬cant, and