represent pure Hamiltonian potentials in which particles move, and since

the reaction coordinates are generalized coordinates which are (in general)

non-linear functions of the particle coordinates, the e¬ective masses (or a

mass tensor in multidimensional cases) are complicated functions of the

coordinates. Instead of computing such e¬ective masses, the frequencies of

the normal modes can much more easily be determined from a relatively

short MD run with full detail in the potential well. Monitoring and Fourier-

™

transforming the velocities ξ(t) of the reaction coordinates will reveal the

eigenfrequencies of the motion of the reaction coordinates in the well of the

7.6 Reconstruction of free energy from PMF 233

βA

8

βAqu

6

βAqu - Acl

4

2

βAcl

0

-2

-4

0.5 1 1.5 2 2.5 3 3.5 4

βhν

Figure 7.4 Helmholtz free energies divided by kB T for a single harmonic oscillator,

as a function of hν/kB T , for both classical and quantum-mechanical statistics. The

drawn line gives the quantum correction to a classical free energy.

PMF without bothering about the e¬ective masses of the resulting motion.

There are as many eigenfrequencies (but they may be degenerate) as there

are independent reaction coordinates. According to quantum statistics (see

Chapter 17), each eigenfrequency ν leads to a contribution to the free energy

of

1 1

Aqu = kB T ln sinh βhν , (7.64)

ν

2 2

which is to be compared to the classical contribution

Acl = kB T ln(βhν). (7.65)

ν

One may use the di¬erence to correct the classical free energy determination,

and “ from temperature derivatives “ the enthalpy

and entropy (see Fig. 7.4).

234 Free energy, entropy and potential of mean force

7.7 Methods to derive the potential of mean force

In general a potential of mean force V mf (r ) describes the e¬ective potential

that determines the motion of coordinates r in a reduced system, averaged

over an equilibrium ensemble of the other coordinates r . In Chapter 8

the use of potentials of mean force in reduced systems is treated in detail.

For simplicity we write r , r as a subdivision of cartesian space, but often

the reduced system is described by a set of generalized coordinates. In

this section we look at methods to derive potentials of mean force, which

may then be useful for implementation in the reduced systems dynamics of

Chapter 8.

In most cases of interest (containing improbable areas of the primed space)

it is impossible to determine V mf directly from an equilibrium simulation of

the whole system. If it were, there would not be much point in reducing the

number of degrees of freedom in the ¬rst place. The following possibilities

are open to derive a suitable potential of mean force:

• From a macroscopic (generally a mean-¬eld) theory. For example, if we

wish to treat a solvent as “irrelevant,” its in¬‚uence on the electrostatic

interactions of charges within the “relevant” particles and on the elec-

trostatic contribution to the solvation free energy of (partially) charged

particles, can be computed from electrostatic continuum theory (see Sec-

tion 13.7). This requires solving the Poisson equation (or the Poisson“

Boltzmann equation) with a ¬nite-di¬erence or Fourier method on a grid

or with a boundary-element method on a triangulated surface. A com-

putationally less demanding approximation is the generalized Born model

(see Section 13.7.5 on page 351). Since such a treatment cannot be ac-

curate on the atomic scale and misses non-electrostatic contributions, the

electrostatic potential of mean force must be augmented by local interac-

tion terms depending on the chemical nature and the surface accessibility

of the primed particles. Another example is the treatment of all parti-

cles outside a de¬ned boundary as “irrelevant”. If the boundary of the

primed system is taken to be spherical, the electrostatic terms may be

represented by a reaction ¬eld that is much simpler to compute than the

Poisson equation for an irregular surface (see Section 13.7.4).

• By thermodynamic integration. It is possible to obtain the derivative(s)

of the potential of mean force at a given con¬guration of r by performing

a constrained equilibrium simulation of the full system and averaging the

constraint forces (which are easily obtained from the simulation) over

the double-primed ensemble. By performing a su¬cient number of such

simulations at strategically chosen con¬gurations of the primed particles,

7.7 Methods to derive the potential of mean force 235

the potential of mean force can be obtained from numerical integration

of the average constraint forces. This method is only feasible for a few

(one to three) dimensions of the primed degrees of freedom, because the

number of points, and hence full simulations, that is needed to reconstruct

a V mf surface in n dimensions increases with the number of points in one

dimension (say, 10) to the power n. By taking the gradient of (8.11), we

¬nd that

‚V (r ,r ) ’βV (r ,r )

e dr

‚V mf (r ) ‚r i

=

e’βV (r r ) dr

‚r i

‚V (r , r )

=

‚r i

= Fc . (7.66)

i

The second line in the above equation gives “ except for the sign “ the

average over the constrained ensemble of the internal force acting on the

i-th primed particle; in a constrained simulation this force is exactly bal-

anced by the constraint force F c on that particle. These equations are

i

modi¬ed for generalized coordinates (see den Otter and Briels, 1998; Sprik

and Ciccotti, 1998; den Otter, 2000).

• By thermodynamic perturbation. Instead of averaging the derivatives of

the potential, we may also average the Boltzmann factor of a (small but

¬nite) perturbation:

e’βV (r +”r ,r ) dr

(r + ”r ) ’ V (r ) = ’kB T ln

mf mf

V

e’βV (r ,r ) dr

= ’kB T ln e’β[V (r +”r ,r )’V (r ,r )]

. (7.67)

This equation is exact, but statistically only accurate for small displace-

ments. By choosing a su¬ciently dense net of con¬gurations to generate

the ensembles, the potentials of mean force can be reconstructed by ¬tting

perturbations of one point to those of a nearby point.

• By umbrella sampling. This method, pioneered by Torrie and Valleau

(1977), restrains, rather than constrains, the primed coordinates around a

given con¬guration by adding a restraining potential V u (r ) to the poten-

tial V (r , r ). This umbrella potential could, for example, be harmonic in

shape. The resulting canonical umbrella distribution wu (r ) in the primed

coordinates will in equilibrium be given by

dr e’βV (r ,r )’βV U (r )