for large frequencies. Expressed in the variance x2 the classical entropy is

1 kB T

ho

m x2

Scl = kB + kB ln , (7.12)

2

2

which becomes negative for small ¬‚uctuations. This is entirely due to the

neglect of quantum behavior. The e¬ect is unphysical, as even a constraint

with no freedom to move “ which should physically have a zero contribution

to the entropy “ has a negative, in¬nite entropy. Equation (7.8) must be

considered wrong: if any eigenvalue of the correlation matrix is zero, the

determinant vanishes and no entropy calculation is possible.

The correct quantum expression for the h.o. is

kB ξ

Squ = ’kB ln 1 ’ e’ξ +

ho

, (7.13)

eξ ’ 1

where ξ = ω/kB T . This can be expressed in terms of the classical variance,

using mω 2 x2 = kB T ,

ξ= . (7.14)

kB T m x2

The entropy now behaves as expected; it goes properly to zero when the

¬‚uctuation goes to zero. In the multidimensional case one can use (7.13)

after diagonalizing the correlation matrix of mass-weighted positional ¬‚uc-

tuations:

Cij = (”xi )(”xj ) , (7.15)

where

√

xi = xi mi . (7.16)

Let us denote the eigenvalues of the C matrix by »k . Each eigenvalue

corresponds to an independent harmonic mode with

ξk = √ (7.17)

kB T »k

and each mode contributes the the total entropy according to (7.13). Now

zero eigenvalues do not contribute to the entropy and do not cause the total

entropy to diverge.

For the exact quantum case one can no longer express the entropy in

218 Free energy, entropy and potential of mean force

terms of the determinant of the correlation matrix as in (7.8). However, by

a clever invention of Schlitter (1993), there is an approximate, but good and

e¬cient, solution. Equation (7.13) is well approximated by S :

e2

¤ S = 0.5kB ln 1 + 2

ho

Squ , (7.18)

ξ

yielding for the multidimensional case

e2 kB T

S = 0.5kB ln Πk 1 + »k . (7.19)

2

Since the diagonal matrix » is obtained from the mass-weighted correla-

tion matrix C by an orthogonal transformation, leaving the determinant

invariant, (7.19) can be rewritten as

e2 kB T

S = 0.5kB ln det 1 + . (7.20)

C

2

This equation is a convenient and more accurate alternative to (7.8). Note

that the mass-weighted positional ¬‚uctuations are needed.

In an interesting study, Sch¨fer et al. (2000) have applied the Schlit-

a

ter version of the entropy calculation to systems where the validity of a

maximum-entropy approach based on covariances is doubtful, such as an

ideal gas, a Lennard“Jones ¬‚uid and a peptide in solution. As expected,

the ideal gas results deviate appreciably from the exact value, but for the

Lennard“Jones ¬‚uid the computed entropy comes close (within ≈ 5%) to

the real entropy. For a β-heptapeptide in methanol solution, which shows

reversible folding in the simulations, the con¬gurational entropy of the pep-

tide itself can be calculated on the basis of the positional ¬‚uctuations of the

solute atoms. However, the contribution of the solvent that is missing in

such calculations appears to be essential for a correct determination of free

energy di¬erences between conformational clusters.

7.3 Thermodynamic potentials and particle insertion

Thermodynamic potentials μi of molecular species i are very important to

relate microscopic to macroscopic quantities. Equilibrium constants of reac-

tions, including phase equilibria, partition coe¬cients of molecules between

phases and binding and dissociation constants, are all expressed as changes

in standard Gibbs free energies, which are composed of standard thermo-

dynamic potentials of the participating molecules (see Section 16.7). How

do thermodynamic potentials relate to free energies and potentials of mean

force and how can they be computed from simulations?

7.3 Thermodynamic potentials and particle insertion 219

The thermodynamic potential of a molecular species is de¬ned as the

derivative of the total free energy of a system of particles with respect to

the number of moles ni of the considered species (see (16.8) on page 428):

‚G

= NA {G(p, T, Ni + 1, Nj ) ’ G(p, T, Ni , Nj )},

μi = (7.21)

‚ni p,T,nj=i

or, equivalently (see (16.27) on page 430),

‚A

≈ NA {A(V, T, Ni + 1, Nj ) ’ A(V, T, Ni , Nj )}, (7.22)

μi =

‚ni V,T,nj=i

where Ni is the number of particles of the i-th species, assumed to be large.

The latter equation is the basis of the particle insertion method of Widom

(1963) to determine the thermodynamic potential from simulations. The

method is as ingenious as it is simple: place a ghost particle of species i in

a random position and orientation in a simulated equilibrium con¬guration.

It is immaterial how the con¬guration is obtained, but it is assumed that

an equilibrium ensemble at a given temperature “ and hence β “ of con-

¬gurations is available. Compute the ensemble average of the Boltzmann

factor exp(’βVint ), where Vint is the interaction energy of the ghost particle

with the real particles in the system. For a homogeneous system ensem-

ble averaging includes averaging over space. The ghost particle senses the

real particles, but does not interact with them and does not in¬‚uence the

ensemble. Now the thermodynamic potential is given by

μexc = ’RT ln e’βVint , (7.23)

where μexc is the excess thermodynamic potential, i.e., the di¬erence be-

tween the thermodynamic potential of species i in the simulated ensemble

and the thermodynamic potential of species i in the ideal gas phase at the

same density. This results follows from (7.22) and A = ’kB T ln Q, yielding

Q(Ni + 1)

μ = ’RT ln . (7.24)

Q(Ni )

With (7.48) we can write

3

dr exp[’β{V (r) + Vint }]

dr ghost

Q(Ni + 1) 2πmi kB T 1

= ,

h2

Q(Ni ) Ni + 1 dr exp[’βV (r)]

(7.25)

where r stands for the coordinates of all real particles and r ghost for the co-

ordinates of the ghost particle. The ratio of integrals is the ensemble average

of exp(’βVint ), integrated over the volume. But since “ in a homogeneous

220 Free energy, entropy and potential of mean force

system “ the ensemble average does not depend on the position of the ghost

particle, this integration simply yields the volume V . Therefore:

3

2πmi kB T V

’ RT ln e’βVint .

μ = ’RT ln (7.26)