zyme E. ”G0 of this process yields the equilibrium binding constant KES .

1

The direct determination of ”G0 by simulation is di¬cult, but the ther-

1

modynamic cycle allows to determine the binding constant of S relative to

another similar substrate S, for which the binding constant KES of process

(3) is known, by two simple processes (2) and (4).

7.2 Free energy determination by spatial integration

Consider a system of N particles that are on the average situated at positions

r i , i = 1, . . . , N , and ¬‚uctuate randomly about those positions. Such sys-

tems are called non-di¬usive when the mean-squared ¬‚uctuations are ¬nite

1 This is a simpli¬ed model; even process (1) is not a physically realistic process. It must be

considered as part of a larger cycle involving a negative ion as well, as realistic thermodynamic

experiments require electroneutral species.

214 Free energy, entropy and potential of mean force

and stationary.2 Without loss of generality we may assume that the system

as a whole has no translational and rotational motion and only consider the

intra-system free energy. If the system is freely translating and rotating,

as a macromolecule in dilute solution, there will be ideal-gas translational

and rotational contributions to the total free energy (see Section 17.5.3 on

page 468); the partition function (p.f.) is the product of a translational p.f

(17.75), a rotational p.f. (17.82) and an internal partition function Qint . If

the system is a non-translating, non-rotating crystal, the total free energy

will consist of the free energy due to lattice vibrations and the internal free

energy, assuming that the coupling between the two is negligible. What

interests us here is the computation of the internal free energy, based on an

ensemble of con¬gurations in a system-¬xed cartesian coordinate system,

obtained by a proper molecular dynamics or Monte Carlo procedure. More-

over, as the interest lies in the relative stability of di¬erent conformations

(or clusters of con¬gurations), we are not interested in the absolute value

of free energies or entropies. Kinetic contributions will cancel out in the

di¬erence between conformations.

The determination of entropy directly from simulations was ¬rst discussed

by Karplus and Kushick (1981). Applications were published “ among oth-

ers “ by Edholm et al. (1983) on the entropy of bilayer membranes and

by DiNola et al.(1984) on the entropy di¬erences between macromolecular

conformations. The method has been evaluated with an error analysis by

Edholm and Berendsen (1984) and the systematic errors on incomplete equi-

libration were discussed by Berendsen (1991a). A major improvement has

been proposed by Schlitter (1993). The method has not become a standard:

it is not applicable to di¬usive systems, it cannot be easily applied to macro-

molecules in solution and considerable computational e¬ort is required for

slowly relaxing systems.

The essential contribution to the free energy that depends on a multi-

dimensional integral is the con¬gurational entropy:

Sconf = const ’ kB w(q) ln w(q) dq, (7.2)

where q is the set of generalized coordinates that exclude translational and

rotational degrees of freedom (and also constrained degrees of freedom if

applicable) and w(q) is the joint probability distribution for all the remaining

n = 3N ’ 6 ’ nc coordinates. The constant in this equation comes from

integration over the kinetic degrees of freedom (the conjugated momenta),

2 This is in contrast to di¬usive systems, such as liquids, in which the mean-squared ¬‚uctuations

increase with “ usually proportional to “ time.

7.2 Free energy determination by spatial integration 215

but also contains a term due to the mass tensor (see page 401) that may

depend on q. The in¬‚uence of the mass tensor on conformational di¬erences

is often negligible and usually neglected or disregarded.

Equation 7.2 is still an unsolvable multidimensional integral. But in non-

di¬usive systems it is possible to derive the most relevant information on

the multidimensional distribution. For example, we can construct the cor-

relation matrix of the coordinate ¬‚uctuations:

Cij = (qi ’ qi )(qj ’ qj ) , (7.3)

or

C = (”q)(”q)T , (7.4)

where i, j = 1, . . . , n run over the generalized coordinates and ”q = q ’ q .

It is generally not possible to assess higher correlations with any accuracy. If

only the matrix of ¬‚uctuations C is known, one can estimate the maximum

entropy of any multidimensional distribution with a given C. By maximizing

Sconf with respect to w under the conditions:

w(q) dq = 1, (7.5)

w(q)”qi ”qj dq = Cij , (7.6)

using Lagrange multipliers (see page 456), it appears that w(q) must be a

multivariate Gaussian distribution in q:

1

w(q) = (2π)’n/2 (det C)’1/2 exp[’ ”qT C’1 ”q]. (7.7)

2

The entropy of this distribution is

1 1

Smax = kB T n(1 + ln 2π) + kB T ln(det C). (7.8)

2 2

Thus, if the distribution really is a multivariate Gaussian, Smax is the con-

¬gurational entropy; for any other distribution Smax is an upper bound for

the entropy of the distribution:

Sconf ¤ Smax . (7.9)

The constants in (7.8) are irrelevant and all we need is the determinant of the

correlation matrix of the positional ¬‚uctuations. Generally it is possible to

determine this entropy accurately; when equilibration is slow, the computed

entropy tends to increase with the length of the simulation and approach

the limit with a di¬erence that is inversely proportional to the length of the

simulation (DiNola et al., 1984; Berendsen, 1991).

216 Free energy, entropy and potential of mean force

It is also possible to derive the entropy from a principal component anal-

ysis of the positional ¬‚uctuations in cartesian coordinates. In that case

translational and rotational degrees of freedom must have been constrained,

which is most easily done by subjecting all con¬gurations to a standard

translational“rotational ¬t. The principal component analysis, often re-

ferred to as “essential dynamics” (Amadei et al., 1993), diagonalizes the

correlation matrix of the positional ¬‚uctuations, thus producing a new set

of collective coordinates with uncorrelated ¬‚uctuations. Each eigenvalue

is proportional to the contribution of its corresponding collective degree of

freedom to the total ¬‚uctuation. There are 6 + nc zero (or very small) eigen-

values corresponding to the translation, rotation and internal constraints;

these should be omitted from the entropy calculation. The determinant of

C in (7.8) is now the product of all remaining eigenvalues.

When there is more knowledge on the ¬‚uctuations than the positional

correlation matrix, the value of the entropy can be re¬ned. Each re¬nement

on the basis of additional information will decrease the computed entropy.

For example, the marginal distributions wi (qi ) over single coordinates3 can

usually be evaluated in more detail than just its variance. This is partic-

ularly important for dihedral angle distributions that have more than one

maximum and that deviate signi¬cantly from a Gaussian distribution. Di-

hedral angles in alkane chains have three populated ranges corresponding to

trans, gauche ’ and gauche + con¬gurations. It is then possible to compute

the con¬gurational entropy for each degree of freedom from

Smarg = ’kB wi (qi ) ln wi (qi ) dqi (7.10)

and use that - after subtraction of the entropy of the marginal distribution

had the latter been Gaussian with the same variance “ as a correction to

the entropy computed from the correlation matrix (Edholm and Berendsen,

1984). Another re¬nement on the basis of extra knowledge is to exploit an

observed clustering of con¬gurations in con¬gurational space: each cluster

may be considered as a di¬erent species and its entropy determined; the

total entropy consists of a weighted average with an additional mixing term

(see (7.55)).

This determination of the classical entropy on the basis of positional ¬‚uc-

tuations has an important drawback: it computes in fact the classical en-

tropy of a harmonic oscillator (h.o.), which is very wrong for high frequen-

3 The marginal distribution of qi is the full distribution integrated over all coordinates except

qi .

7.2 Free energy determination by spatial integration 217

cies. For a one-dimensional classical h.o., the entropy is given by

kB T

ho

Scl = kB + kB ln , (7.11)

ω