Bios 203

BIOS 203:

Free Energy
Methods

Tom Markland

1
Motivation
Lots of interesting properties require
knowledge of free energies:
• Phase diagrams/coexistence.
• Drug binding affinities.
• Rates of reactions.
• Equilibrium constants.
• Solvation properties.
• Acid-base equilibria.
• Isotope effects.

2
Statistical Mechanics
Recap
Partition function in classical mechanics,

Momentum part Position part

Position and momentum integrals are separable,

Today we will
Yielding the configuration integral, primarily deal with
this.
Observables obtained performing trajectories which sample points in the ensemble and averaging
the observables over them (assuming ergodic hypothesis),

If observable depends
Expectation value of observable O only on position.
3
Simulations
Basic simulation scheme

Get energy/forces
Initial Structure Pair/Empirical
From known Ab initio
crystal structure or
starting from
lattice and
equilibrating.

Evolve system: Repeat until

Monte Carlo average of
Molecular Dynamics observables are
converged.

Calculate observables
energy, pressure,
structure etc.

4
Free energy from statistical mechanics
How to obtain free energy?

For simplicity everything will be presented in terms of Helmholtz energy A (i.e. working in
the NVT ensemble) but extension to Gibbs (NPT) is straightforward.

• Not good news: to get an absolute free energy of a system we need to know the
entire partition function i.e. to count all the states for a given N,V,T!

 In practice calculating the partition function is impossible for a system of

more than a few degree of freedom (with the exception of special cases).

• However, calculating differences in free energies is possible.

5
Free Energy Perturbation

• The above expression is not in a nice form to extract from a simulation.

Rearrange Inserting 1

Now take ratio,

Of the form:

Free energy perturbation (FEP)

(Zwanzig 1954)

6
Free Energy Perturbation

What does this translate to?

• Do a simulation of system 0.

• During the simulation accumulate the average:

• At the end take the natural logarithm of the observable

and multiply by –kBT and you obtain the free energy
difference between system 1 and 2.
Example:

I want to develop a new drug that binds into an active site more strongly
than my current drug.
Can I do this by replacing a Hydrogen by a Methyl Group?
 Potential V0 is force field with H and V1 is that with Me. Hence can use
FEP to calculate the difference in free energy of binding.

7
Free Energy Perturbation
Unfortunately FEP is not perfect
• Suppose the two potentials (V0 and V1) are very different.
• Consider in 1D two systems with different potentials:
Remember in FEP we do a
V(x) simulation in potential V0
Potential and accumulate the average
V0 of
V1

x
P(x)
Probability

In this region In this region

is small but potential V0 samples it a lot. is big but potential V0 samples it very little.

8
Free Energy Perturbation
V(x) V0 V1
Potential

P(x) x

Probability

So what happens when we do a simulation? x

Do a simulation in potential V1 and accumulate the average of
Massive jumps in expectation value correspond to the rare events of entering a region where
potential V2-V1 is small giving rise to a massive contribution to

Simulation time/ MC steps

FEP fails when the difference potential shows a
large standard deviation compared to .

9
Free Energy Perturbation
How to avoid this problem?
0
• Good news: free energies are state functions.
i.e. they only depend on the current state of
the system and not how the system acquired 1
that state.
• This allows us to take any path between the two states
we are interested in whether or not the intermediates
represent anything “real”.
Hence we can break the journey from transforming from potential V0 to V2 by splitting it
into M sections:
V0 V1 Here M=5
V0.2 V0.4 V0.6 V0.8
We can then calculate the free energy change for each section and add them together to generate
the total free energy change on going from 1 to 2:

Analysis: Proc. Roy. Soc. A 468, 2-17 (2012). 10

Free energy from statistical mechanics
General expression for making M steps between start and end-point

How does one pick intermediate potentials (V0.2, V0.4 etc.)?

• In principles could pick anything.

• However, if intermediates are widely spaced or take a long path between the two end
potentials then evaluation will be inefficient.
• Hence a wise choice might be a linear interpolation between V0 and V1.

When λ=1 all potential

V0.5 1 and when λ=0 all
V(x) V0 V1 potential 0. In between
Potential a mixture of them.

Inserting a middle state which has overlap with both V0 and V1 allows
for much more efficient estimates of the free energy difference.

11
Thermodynamic Integration
Another approach is to use thermodynamic integration
As before we want to calculate the free energy change between two potentials V0 and V1.

Define a potential:

Where f(λ) and g(λ) switch between the potentials V0 and V1 and satisfy:
• f(0)=1, f(1)=0, g(0)=0 and g(1)=1

Note: Although f(λ) = (1- λ) and g(λ) = λ are a valid choice one can pick other forms that
satisfy the above constraint.

• The free energy is now a function of the phase point (N,V,T) and the parameter λ

• Taking derivative of A w.r.t. λ yields:

12
Thermodynamic Integration

What is ?

Combining gives β cancels kBT

Recognize this is of form:

Relates change
in potential
with change in
free energy.

13
Thermodynamic Integration
Integrate

Thermodynamic integration:
Kirkwood (1935)
If we choose:

Implementation
• Perform a series of simulations at different values of λ to get the gradients:
• Perform the integral of the gradients by midpoint rule or other numerical
integration schemes.
More advanced techniques (beyond this course): Adiabatic free energy dynamics
• Allow λ to evolve continuously in the simulation by including it as an additional coordinate
which is dynamically, As long as this coordinate moves much slower than the system relaxes
(adiabatic separation) one can obtain the free energy from one simulation.
Adiabatic free energy dynamics: J. Chem. Phys., 116, 4389 (2002) 14
Free energy along collective variables
• In many cases in chemistry one is interested in the probability of finding a particle at some particular
point along a single coordinate of the system called a collective variable, R(r).
• The collective variable can be defined as a general function of any of the positions in the system r.
• For example in transition state theory a key input is knowing the probability of reaching the top of the
barrier along the “reaction coordinate”.
Energy
• Although we will refer to R(r) as a
collective variable (CV) it can also be
called an order parameter or reaction
coordinate.
Collective variable

Example: to calculate a transition

stat theory prediction of the rate of
hydrogen diffusion through ice one
can compute the free energy
associated with moving the H
between the 6 O’s defining ice’s
hexagonal lattice.

15
Free energy along collective variables

Probability of obtaining a given configuration:

Consider a system where one coordinate R (which can be a complex function of all coordinates
in the system i.e. R(r)) is constrained to hold a particular value R* but all other coordinates are
distributed according to Boltzmann,

It follows the probability distribution along the coordinate is given by,

16
Free energy along collective variables
Constrained
partition function.

Recall, and hence

Potential of
mean force

In practice we are (like before) more interested in differences in free energy in going from one
position along the coordinate R to another e.g. R0 to R1.

Really useful equation: allows one to convert a histogram of probabilities along a particular
coordinate to the free energy change associated with moving along that coordinate.

Figure: Ann. Rep. in Comput. Chem. 6, 280–296 (2010) 17

Free energy along collective variables
Example 1:
Radial distribution gives probability of being found a given distance away from a
particle at the center.
Can convert probability
Using: to a free energy change
as a function of particle
separation.
g(R)
A(R)
V(R)
Note: The free energy as a function of
particle separation, A(R) is not the
Energy

same as pair potential, V(R) unless

only two particles are present.

18
Free energy along collective variables
Example 2:

Can also extend to more constrained dimensions e.g. the Ramachandran plot of
alanine dipeptide shows the free energy as a function of two torsions angles of the
molecule.

Red regions are

low free energy,
blue is high free
energy.

In both this and the previous case the problem is sufficiently simple that one can run
for long enough to see point in all regions of interest  suppose there is a high barrier?

19
Free energy along order parameters
Methods of obtaining free energy when high barriers are present
Umbrella
Umbrella sampling Energy
potential
• Add an additional force to the potential to keep
the system close to a particular value of the CV,
R* e.g. a harmonic force.

R* R = Collective
• The biased simulation can then be converted variable
into an unbiased probability by using techniques
such as WHAM1 or umbrella integration2. 1.) Kumar et. al. J. Comp. Chem. 13, 1011 (1992)
2.) Kästner and Thiel, J. Chem. Phys. 123, 144104
Hard constraints (2005)

• One can also use hard constraints to perform molecular dynamics

with the collective variable held at one position. Caveat: Can also get an
additional term due to the
• The free energy is then obtained by integrating the average force in Jacobian of the change into
the direction of the collective variable. generalized coordinates
• This gives rise to the name “Potential of mean force” for . (depends on form of R).

Umbrella sampling: Torrie and Valleau, J. Comput. Phys. 2 187-199 (1977) 20

Metadynamics
A more recent approach to obtaining the free energy along an order parameter is
metadynamics.

Basic idea
• Each time a position along the order parameter is visited place down a Gaussian to “fill up”
the potential.
• The Gaussians placed bias against revisiting the same region over and over.
• At the end of the simulation the probability density can be obtained simply by looking at
the density of the Gaussians placed at each position along the collective variable.
Place Gaussians on places
previously visited Note: to obtain the free energy
Energy accurately the rate at which the
Gaussians are dropped must be
slow compared to the rate at which
the other coordinates move.

Collective variable
Video by Giovanni Bussi showing metadynamics in action:
https://www.youtube.com/watch?feature=player_embedded&v=IzEBpQ0c8TA

Metadynamics: Laio and Parrinello PNAS 99 (20): 12562–12566 (2002) 21

Phase coexistence
• Suppose we know a coexistence point of two phases.
TIP4P/05

Model

Experiment
• Can we trace the rest of the phase diagram from one
coexistence point?
 Yes, by Gibbs-Duhem integration TIP3P
• Integrate the Clapeyron equation numerically:
Model

Run NPT simulations of each phase at the

Freezes to Ice II.
current point and calculate difference in Hexagonal ice Experiment
enthalpy and volume. Use that to obtain only stable at
new coexistence point and repeat. negative pressure!

Gibbs-Duhem: Kofke Mol. Phys. 78, 6 1331 (1993) Phase Diagrams: C. Vega et. al., Faraday Discuss., 141, 251-276 (2009) 22
Summary

Thanks for listening.

Books:

D. Frenkel and B. Smit , Understanding molecular simulation: from algorithms to

applications

M. E. Tuckerman, Statistical mechanics: theory and molecular simulation

M.P. Allen and D.J. Tildesley, Computer simulation of liquids

23