tion motion: the ductile solid becomes brittle through this work-hardening

process. Thus the simulation, with an impressive billion atoms still repre-

senting only a very small solid of 0.3 μm size, gives detailed insight into the

mechanisms of work-hardening. The dynamical time span is on the order

of a few nanoseconds, enough time for the phenomenon to achieve a ¬nal

structure state for the small size solid cube.

Exercises

Write out the operator product U(”t) ≈ Uf (”t/2)Uv (”t)Uf (”t/2)

6.1

to obtain the velocity-Verlet algorithm.

6.2 Obtain another algorithm by interchanging Uv and Uf .

6.3 Solve (6.107).

6.4 Compute (in reduced units) the period of oscillation of two Lennard“

Jones particles at the distance where their interaction energy is min-

imal. What would be an appropriate time step (in reduced units)

for a leap-frog simulation of a Lennard“Jones ¬‚uid? How many fs is

that time step for argon?

210 Molecular dynamics

Figure 6.13 Five snapshots from a one-billion particle MD simulation of the propa-

gation of a crack in a copper crystal under tension, showing the massive formation of

dislocations. See text for details. Figures were kindly provided by Dr Farid Abra-

ham of IBM Research and Lawrence Livermore National Laboratory, Livermore,

CA, USA.

7

Free energy, entropy and potential of mean force

7.1 Introduction

As we know from the applications of thermodynamics, free energy is much

more important than energy, since it determines phase equilibria, such as

melting and boiling points and the pressure of saturated vapors, and chemi-

cal equilibria such as solubilities, binding or dissociation constants and con-

formational changes. Unfortunately, it is generally much more di¬cult to

derive free energy di¬erences from simulations than it is to derive energy

di¬erences. The reason for this is that free energy incorporates an entropic

term ’T S; entropy is given by an integral over phase space, while energy

is an ensemble average. Only when the system is well localized in space

(as a vibrating solid or a macromolecule with a well-de¬ned average struc-

ture) is it possible to approximate the multidimensional integral for a direct

determination of entropy. This case will be considered in Section 7.2.

Free energies of substates can be evaluated directly from completely equi-

librated trajectories or ensembles that contain all accessible regions of con-

¬gurational space. In practice it is hard to generate such complete ensembles

when there are many low-lying states separated by barriers, but the ideal

distribution may be approached by the replica exchange method (see Section

6.6). Once the con¬gurational space has been subdivided into substates or

conformations (possibly based on a cluster analysis of structures), the free

energy of each substate is determined by the number of con¬gurations ob-

served in each substate. One may also observe the density of con¬gurations

along a de¬ned parameter (often called an order parameter or a reaction

coordinate, which is a function of the coordinates) and derive the potential

of mean force along that parameter. When a replica-exchange method has

been used, the free energies of substates are obtained simultaneously for a

range of temperatures, providing energies, entropies and speci¬c heats.

211

212 Free energy, entropy and potential of mean force

(2)

S´

S

(2)

N

(3) (1)

(3) (1)

N

(4) S´

S

E (4) E

(a) (b)

Figure 7.1 Two thermodynamic cycles, allowing the replacement of one path by

the sum of three other paths (see text). (a) Hydration of an ion (+) through the

intermediate of a neutral atom (N), (b) binding of a substrate S to an enzyme E

(dark gray), compared to the binding of another, similar, substrate S. Light gray

squares represent an aqueous environment.

In general we cannot evaluate highly multidimensional integrals and all we

have is a trajectory or ensemble in phase space or in con¬gurational space.

Then we must apply tricks to derive free energy di¬erences from ensemble

averages.

This chapter is mostly about such tricks. In Section 7.3 Widom™s particle

insertion method is treated, which relates the ensemble average of the Boltz-

mann factor of an inserted particle to its thermodynamic potential. Section

7.4 is about perturbation and integration methods that relate the ensemble

average of the Boltzmann factor of a small perturbation to the di¬erence

in free energy, or the ensemble average of a Hamiltonian derivative to the

derivative of the free energy. In Section 7.5 we ¬rst make a clear distinction

between two kinds of free energy: free energy of a thermodynamic state, and

free energy in a restricted space as a function of one or more reaction coor-

dinates. The latter is called a “potential of mean force,” or PMF, because

its derivative is the ensemble-averaged force. The PMF is indispensable for

simulations in reduced dimensionality (Chapter 8), as it provides the sys-

tematic forces in such cases. In Section 7.6 the connection between PMF and

free energy is made and Section 7.7 lists a number of methods to determine

potentials of mean force. Finally, Section 7.8 considers the determination of

free energy di¬erences from the work done in non-equilibrium pathways.

Practical questions always concern free energy di¬erences between ther-

modynamic states. So we need to obtain free energy di¬erences from simula-

tions. But, by the use of thermodynamic cycles, the computed pathways and

7.2 Free energy determination by spatial integration 213

intermediates may di¬er from the experimental ones as long as the end re-

sults match. Thus one may choose even physically unrealistic or impossible

intermediates to arrive at the required free energy di¬erences. For example,

if we wish to compute the free energy of hydration of a sodium ion, the

interest lies in the di¬erence in standard free energy ”G0 of process (1) be-

1

1 Here, (g) means ideal gas referred to standard concentration and (aq)

low.

means in¬nitely dilute aqueous solution referred to standard concentration.

For all processes, constant pressure and temperature are assumed. Now,

process (1) can be built up from processes (2), (3) and (4), with a neutral,

sodium-like atom N as intermediate (see Fig. 7.1a). Here N may be any

convenient intermediate, e.g., a repulsive cavity with a size close to that of

a sodium ion.

(1) Na+ (g) ’ Na+ (aq) ”G0 1

’ Na+ (g) ”G0

(2) N(g) 2

’ ”G0

(3) N(g) N(aq) 3

(4) N(aq) ’ Na + (aq) ”G0

4

Going clockwise around the thermodynamic cycle, the total free energy

change must be zero:

”G0 ’ ”G0 ’ ”G0 + ”G0 = 0. (7.1)

1 4 3 2

Therefore, ”G0 can be determined from three di¬erent processes, which are

1

each simpler and more e¬cient to compute. The intermediate N can be

chosen to optimize computational e¬ciency, but it can also be chosen to

provide a single intermediate for the hydration of many di¬erent ions.