which the derivative or perturbation is averaged, is a proper equilibrium en-

semble. In slowly relaxing systems there may be remains of history from the

previous »-point and the integration may develop a systematic deviation.

It is recommended to perform the thermodynamic integration from both

directions, i.e., changing » from 0 to 1 as well as from 1 to 0. Systematic de-

viations due to insu¬cient equilibration are expected to have opposite signs

in both cases. So the obtained hysteresis is an indication of the equilibration

error.

A limiting case of thermodynamic integration is the slow-growth method,

in which » is changed with a small increment (”»)i at the i-th step in

a molecular dynamics simulation, starting at 0 and ending at 1 (or vice

versa).6 This increment may or may not be taken as a constant. Then the

total change in free energy is approximated by

‚V (»)

”A = (”»)i . (7.41)

‚» i

i

Only in the limit of an in¬nitely slow change of » a true free energy di¬erence

will be obtained; if the growth is too rapid, the ensemble will “lag behind”

the proper equilibrium at the actual value of ». This will lead to an easily

detected hysteresis when slow-growth simulations in forward and backward

directions are compared. The average between the results of a forward and

backward integration is always more accurate than either value. Figure 7.2

gives an early example of hysteresis in a simulation that changes a model

neon atom into a sodium ion in aqueous solution by charging the atom

proportional to time (Straatsma and Berendsen, 1988). In this case the free

energy appears to change quadratically with time and the ensemble appears

to relax quickly. Pearlman and Kollman (1989b) and Wood (1991) have

5 The standard error in the integral can be evaluated if the standard error σi in each ensemble

average Ai is known. Numerical integration yields an integral that can be expressed as ”A =

wi (Ai ± σi ), where wi are weights depending on the interval and the procedure used. The

22

standard error in ”A equals wi σ i .

6 The slow-growth method was pioneered by Postma (1985), see Berendsen et al. (1985), and

applied “ among others “ by Straatsma and Berendsen (1988) and by Pearlman and Kollman

(1989b).

224 Free energy, entropy and potential of mean force

± ”G (kJ/mol)

440

420

400

20 40 60 80

Total simulation time T (ps)

Figure 7.2 Free energy change resulting from transforming a model neon atom into

a sodium ion in a bath of 216 water molecules, by linearly increasing the charge on

the atom in a total simulation time T . Upward triangles are for growth from Ne

to Na+ , yielding a negative ”G; downward triangles are for the opposite change.

Dashed curves are predictions from the theory of Wood (1991) assuming a relaxation

time for the ensemble of 0.15 ps. Data are from Straatsma and Berendsen (1988).

analyzed the e¬ects of the rate of slow-growth free energy determinations.

At present slow growth is used less frequently than integration based on a

number of well-equilibrated intermediate points because the latter allows a

better evaluation and optimization of the overall accuracy.

Note A similar remark as was made in connection with particle insertion can be

made here as well. In most applications one is interested in constant pressure

rather than constant volume conditions. For example, in order to ¬nd equilibrium

constants, one needs ”G0 . Using N, V, T ensembles one may need corrections to

connect the end points at constant pressure rather than constant volume. It is

much more elegant to use N, p, T ensembles, with partition function ” (see (17.31)

on page 461) that relate directly to Gibbs free energies. Ensemble averages of

Hamiltonian derivatives now yield derivatives of G rather than A.

There are several practical considerations concerning the method used for

the integration of free energy from initial to ¬nal state. As computational

integration is not limited to physically realistic systems (i.e., as long as the

initial and ¬nal states are realistic), there is almost no bound to the phantasy

that can be introduced into the methodology. The word “computational

7.4 Free energy by perturbation and integration 225

alchemy” is not misplaced, as one may choose to change lead into gold, be

it that the gold must “ unfortunately “ be returned to lead before a realistic

result is obtained. We list a few tricks and warnings.

• The free energy as a function of » may not be well-behaved, so that nu-

merical integration from a limited number of points becomes inaccurate.

The density of points in the range 0 ¤ » ¤ 1 can be chosen to optimize the

integration accuracy, but ideally A(») should be a smooth function with-

out steep derivatives, being well-represented by a polynomial of low order.

One can manipulate the function by changing the functional dependence

of the Hamiltonian on ».

• Ideally the free energy curve should be monotonous; if a large interme-

diate maximum or minimum occurs, computational e¬ort must be spent

to compute compensating free-energy changes. A maximum may eas-

ily occur when there are highly repulsive con¬gurations at intermediate

values of »™s. Such repulsive intermediates can be avoided by choosing

appropriately smoothed potentials.

• Replacing diverging functions as the repulsive r’12 or dispersion and Cou-

lomb interactions by soft-core interactions for intermediate »™s removes

the singularities and allows particles to move through each other rather

than having to avoid each other.7 The GROMACS software (van der

Spoel et al., 2005) uses a modi¬cation of the distance between particles

of the form

V (r) = (1 ’ »)VA (rA ) + »VB (rB ), (7.42)

rA = (c»2 + r6 )1/6 , (7.43)

rB = [c(1 ’ »)2 + r6 ]1/6 , (7.44)

while Tappura et al. (2000) switch to a function ar6 + b below a speci¬ed

(short) distance, with a and b such that the potential function and its

derivative are continuous at the switch distance.

• Another non-physical intervention is to allow particles to move into a

fourth spatial dimension for intermediate »™s (van Gunsteren et al., 1993).

Since there is much more space in four than in three dimensions, parti-

cles can easily avoid repulsive con¬gurations. But they also loose their

structural coherence and the motion in the fourth dimension must be

carefully restrained. The method is more aesthetically appealing than it

is practical.

7 Soft-core potentials were originally developed for structure optimization and protein folding

(Levitt, 1983; Huber et al., 1997; Tappura et al., 2000).

226 Free energy, entropy and potential of mean force

• There can be problems when particles “vanish” at either end point of

the integration path. When the interactions of an atom with its environ-

ment are made to vanish, the particle is still there in the simulation as

an ideal gas atom. It has mass and velocity, but is uncoupled to other

degrees of freedom and therefore does not equilibrate properly. Prob-

lems are avoided by constraining the vanishing particle to a ¬xed position

where it has neither a kinetic energy nor a con¬gurational entropy and

does not contribute to the free energy. One should take care that the

»-dependence is not diverging near the value of » where the particle van-

ishes. Strictly speaking, one should also correct for the vanishing kinetic

term (2πmkB T /h2 )’1/2 in the free energy, but that term will always be

compensated when a complete, physically realistic, cycle is completed.

• When particles are changed into other particles with di¬erent mass, the

free energy change has a di¬erent kinetic term. It is possible to change

the masses also with a coupling parameter, but there is no need to do

that, as “ just as in the case of a vanishing particle “ the kinetic e¬ect

will always be compensated when a complete, physically realistic, cycle is

completed. Real free energy di¬erences always concern the same number

and type of particles on both sides of the reaction.

• Be careful when the coupling parameter involves a constraint.8 For ex-