The ¬rst term is the thermodynamic potential of the i-th species in an ideal

gas of non-interacting particles at a density (Ni + 1)/V . For large Ni this

density is to a good approximation equal to Ni /V . Equation (7.26) is also

valid for small Ni ; for example, when a water molecule is inserted into a

¬‚uid consisting of another species (e.g., hexane), Ni = 0 and the excess

thermodynamic potential is obtained with respect to the thermodynamic

potential of ideal-gas water at a density of one molecule in the entire volume.

For solutions we are interested in the standard thermodynamic potential

and activity coe¬cient of the solute and “ in some cases “ the solvent.

How are these obtained from particle insertion? On the basis of molar

concentration c, the thermodynamic potential (16.55) is expressed as

γc c

μ(c) = μ0 + RT ln , (7.27)

c

c0

where c0 is an agreed standard concentration (e.g., 1 molar) and μ0 is de¬ned

c

by (16.58):

c

def

μ0 = lim μ(c) ’ RT ln 0 . (7.28)

c

c

c’0

This de¬nition guarantees that the activity coe¬cient γc approaches the

value 1 for in¬nitely dilute solutions. A single measurement of μ(c) at one

concentration can never determine both μ0 and γc . The standard thermo-

c

dynamic potential requires a measurement at “in¬nite” dilution. A single

solute molecule in a pure solvent can be considered as in¬nitely dilute since

there is no interaction between solute particles. Thus, for inserting a single

solute particle, the activity coe¬cient will be unity and the thermodynamic

potential is given by

c

μ = μ0 (solution) + RT ln 0 , (7.29)

c

c

where

1

c= . (7.30)

NA V

But since μ is also given by

μ = μ(id.gas, c) + μexc (7.31)

c

= μ(id.gas, c0 ) + RT ln + μexc , (7.32)

0

c

7.4 Free energy by perturbation and integration 221

it follows that

μ0 (solution) = μ0 (id.gas) + μexc . (7.33)

c c

Here, μexc is “measured” by particle insertion according to (7.23). The

standard concentration for the solution and the ideal gas must be the same.

Note Thus far we have treated particle insertion into a canonical (N, V, T ) en-

semble, which yields a relation between the averaged Boltzmann factor and the

Helmholtz free energy A, based on (7.22). This equation is not exact and not valid

for small numbers Ni . However, (7.21) is exact, as G = ni μi under conditions of

constant pressure and temperature (see (16.12) on page 428); this relation is valid

because both p and T are intensive quantities. Such a relation does not exist for

A. It is possible to make corrections to A, using the compressibility, but it is more

elegant to use a (N, p, T ) ensemble. The N, p, T average of the Boltzmann factor

yields the thermodynamic potential exactly (see (17.33) and (17.31) on page 461).

The problem with the particle-insertion method is that in realistic dense

¬‚uids the insertion in a random position nearly always results in a high, re-

pulsive, interaction energy and hence in a negligible Boltzmann factor. Even

with computational tricks that avoid the full computation of all interaction

energies it is very di¬cult to obtain good statistics on the ensemble aver-

age. A way out is to insert a smaller particle and “ in a second step “ let

the particle grow to its full size and determine the change in free energy by

thermodynamic integration (see next section). In cases where the di¬erence

in thermodynamic standard potential between two coexisting phases is re-

quired “ as for the determination of partition coe¬cients “ a suitable method

is to determine the potential of mean force over a path that leads from one

phase into the other. The di¬erence between the two plateau levels of the

PMF in the two phases is also the di¬erence in standard thermodynamic

potential.

7.4 Free energy by perturbation and integration

Consider a potential function V (r, ») with a parametric dependence on a

coupling parameter 0 ¤ » ¤ 1 that modi¬es the interaction. The two

extremes, » = 0 and » = 1, correspond to two di¬erent systems, A and B,

respectively, with interaction functions VA (r) and VB (r):

VA (r) = V (r, » = 0), (7.34)

VB (r) = V (r, » = 1). (7.35)

For example, system A may consist of two neon atoms dissolved in 1000

water molecules, while system B consists of one sodium ion and one ¬‚uoride

222 Free energy, entropy and potential of mean force

ion dissolved in 1000 water molecules. The parameter » changes the neon“

water interaction into an ion“water interaction, essentially switching on the

Coulomb interactions. Or, system A may correspond to a protein with a

bound ligand LA in solution, while system B corresponds to the same protein

in solution, but with a slightly modi¬ed ligand LB . The end states A and

B represent real physical systems, but intermediate states with a » unequal

to either zero or one are arti¬cial constructs. The dependence on » is not

prescribed; it can be a simple linear relation like

V (») = (1 ’ »)VA + »VB , (7.36)

or have a complex non-linear form. The essential features are that the

potential is a continuous function of » that satis¬es (7.34) and (7.35).

Now consider the Helmholtz free energy of the system at a given value of

»:

e’βV (r,») dr .

A(») = ’kB T ln c (7.37)

It is impossible to compute this multidimensional integral from simulations.

But it is possible to compute A(» + ”») as a perturbation from an ensemble

average:4

exp[’βV (r, » + ”»)] dr

A(» + ”») ’ A(») = ’kB T ln (7.38)

exp[’βV (r, »)] dr

= ’kB T ln e’β[V (»+”»)’V (»)] . (7.39)

»

It is also possible to compute dA/d» from an ensemble average:

‚V

‚» (r, ») exp[’βV (r, »)], dr

dA ‚V

= = . (7.40)

d» ‚»

exp[’βV (r, »)], dr »

The averages must be taken over an equilibrium ensemble using V (»). Thus,

if the »-path from 0 to 1 is constructed from a number of intermediate points,

then the total ”A = AB ’ AA = A(» = 1) ’ A(» = 0) can be reconstructed

from the ensemble averages at the intermediate points. In general the most

convenient and accurate reconstruction is from the derivatives at intermedi-

ate points by integration with an appropriate numerical procedure (Press et

al., 1993), e.g., by computing a cubic spline with the given derivatives (see

4 This equation was ¬rst given by Torrie and Valleau (1974) in the context of Monte Carlo

simulations of

Lennard“Jones ¬‚uids. Pearlman and Kollman (1989a) have re¬ned windowing techniques for

thermodynamic integration.

7.4 Free energy by perturbation and integration 223

Chapter 19).5 This procedure to ¬nd di¬erences in free energies is called

thermodynamic integration. Integration can also be accomplished from a

series of perturbations, using (7.39); in that case there may be systematic

deviations if the interval is not very small and it is recommended to ap-

ply both positive and negative perturbations from each point and check the

closure.