236 Free energy, entropy and potential of mean force

Therefore,

u (r)

V mf (r ) = constant ’ kT ln wu (r )e+βV

= constant ’ kT ln[wu (r )] ’ V u (r), (7.69)

which says that the potential of mean force can be reconstructed in the

neighborhood of the restrained con¬guration by keeping track of the dis-

tribution over the primed coordinates and correcting the bias caused by

the umbrella potential. This reconstruction is only accurate in a region

where su¬cient statistics is obtained to determine wu (r ) accurately. The

full potential of mean force can again be reconstructed by ¬tting adjacent

umbrella distributions to each other.

An alternative to reconstructing the local V mf from the distribution

function is averaging of the umbrella force, which is easily monitored in a

MD simulation.11 In the case of a harmonic umbrella, the average force

is also equal to the mean displacement of the coordinate(s) on which the

umbrella is imposed (with respect to the umbrella center), divided by the

harmonic force constant. The average umbrella force is approximately,

but not exactly, equal to the derivative of the potential of mean force at

the umbrella center. In fact, it is given exactly by a weighted average of

the derivative over the umbrella distribution wu (r ):

dr wu (r )[‚V mf (r )/‚r i ]

u

= , (7.70)

F u

dr wu (r )

which is accurate to second order (i.e., including the second derivative

of the potential) to the derivative at the average position of the primed

coordinate in the umbrella ensemble. This average can also be used to

reconstruct the potential of mean force.

Proof The average umbrella force cancels the average internal force acting

on the primed coordinate, and thus is equal to the average derivative of

the total potential V = V (r , r ) in the umbrella ensemble:

dr [‚V /‚r i ] exp[’βV ’ βV u ]

dr

‚V

u

= = . (7.71)

F u

dr dr exp[’βV ’ βV u ]

‚r i u

In the nominator of (7.71) the term exp(βV u ) can be taken out of the

integration over r and the remainder can be replaced by the derivative

11 K¨stner and Thiel (2005) describe a method for the determination of the derivative of V mf .

a

7.7 Methods to derive the potential of mean force 237

of V mf (8.11):

‚V mf

‚V ’βV

e’βV dr .

e dr = (7.72)

‚r i ‚r i

We now obtain

dr exp(’βV u )[‚V mf /‚r i ] dr exp(’βV )

u

= , (7.73)

F u

dr dr exp(’βV ’ βV u )

from which (7.70) follows.

• By particle insertion by the method of Widom (1963) in the special case

that the potential of mean force is a function of the position of a speci¬c

particle type and in addition the density of the medium is low enough to

allow successful insertions. Let us consider the case that the potential

of mean force of a system that is inhomogeneous in the z-direction but

homogeneous in the x, y-plane (e.g., containing a planar phase bound-

ary situated at z = 0 between two immiscible ¬‚uids, or between a ¬‚uid

and a polymer, or containing a membrane between two liquid phases).

Assume an equilibrium simulation is available. One may place a ghost

particle (not exerting any forces on its environment) at a given z- but

randomly chosen x, y-position in a con¬guration of the generated ensem-

ble, and registers the total interaction energy V between the particle and

the surrounding molecules. This insertion is repeated many times, and

μex (z) = ’kT ln exp(’βV ) z is determined. The standard chemical po-

tential μ0 (z) of the particle type at position z in the inhomogeneous sys-

tem is equal to the standard chemical potential μ0 of that particle type as

id

an ideal gas, plus the measured excess μex (z); the potential of mean force

as a function of z is, but for an additive constant, equal to μex (z). One

may choose the number of insertions per value of z to satisfy statistical

requirements.

• By directly measuring the particle concentration c(z) (as a number density

per unit volume) in an equilibrium simulation in the special case (as above)

that the potential of mean force is a function of the position of a speci¬c

particle type, in regions where that concentration is high enough to allow

its determination with su¬cient accuracy. In those regions the potential

of mean force can be constructed from

V mf (z) = const ’ kT ln c(z). (7.74)

• By enforcing the system to move from one part of con¬gurational space

to another. In such pulling simulations, also called steering molecular

238 Free energy, entropy and potential of mean force

dynamics (SMD),12 an extra external force is exerted on the system such

that it will move much more quickly over intervening barriers than it would

do in an equilibrium simulation. The advantage is that the steering can be

executed in a global way without the need for a detailed description of the

reaction coordinate; e.g., one may choose to distribute the pulling force

over many atoms and let the system decide to ¬nd the easiest pathway

under such a global force. The disadvantage is that the work exerted by

the external force contains a frictional contribution and is therefore not

related to the potential of mean force in a straightforward way. Only in

the uninteresting limit of zero pulling rate the force acts in a reversible

manner and equals the derivative of a potential of mean force.

However, in the case of small, but ¬nite forces or pulling rates, it is pos-

sible to derive the potential of mean force V mf from the work W exerted

by the external force. It is necessary to use a thermostat that prevents

local heating by frictional forces. Although W will always exceed the

reversible work ”V mf , it can be shown (Jarzynsky, 1997a, 1997b) that

e’β”V = e’βW ,

mf

(7.75)

which is valid if a su¬ciently large set of pulling simulations, starting from

samples from an initial equilibrium ensemble, is performed. It is quite

di¬cult to obtain su¬cient statistics for the evaluation of exp(’βW ) ;

Park et al. (2003) found that more reliable results are obtained with a

second-order expansion:

β2

’βW

= ’β W + W2 ’ W 2

ln e , (7.76)

2

using a time-dependent external potential of the form

k

[ξ(r) ’ »]2 ,

V ext = »(t) = »0 + vt, (7.77)

2

with a large force constant k. The method is very similar to the imposi-

tion of a slowly changing constraint. Except for near-equilibrium SMD,

the method is not preferred above thermodynamic integration using con-

strained or umbrella-restrained intermediate simulations. In the next sec-

12 Simulations involving pulling were ¬rst performed with the explicit purpose to mimic exper-

iments with the atomic force microscope (AFM), e.g., pulling a ligand bound to a protein

away from the protein. In such simulations the external force is due to a spring that is slowly

displaced, like the lever of an AFM, but for the computation of a potential of mean force the

spring is not required. The ¬rst simulation of this kind was on the “unbinding” of biotin from

the protein streptavidin by Grubm¨ ller et al. (1996), soon followed by Izrailev et al. (1997).

u