Advanced sampling ChE210D

Today's lecture: methods for facilitating equilibration and sampling in complex,

frustrated, or slow-evolving systems

Difficult-to-simulate systems
Practically speaking, one is always limited in the length and time scales accessible from simula-
tion. The former relates to the maximum system size that can be used (in terms of number of
atoms), and the latter to the length of the run that can be performed (in terms of number of
MD integration or MC steps).

If the intrinsic, physical time scale for some system of interest is beyond that accessible from
simulation, we cannot investigate it using straightforward dynamical methods like molecular
dynamics because it will be challenging to equilibrate the system. However, if we are only
interested in thermodynamic and not kinetic properties, then we can construct artificial dynam-
ics in the system that accelerate convergence towards equilibrium. Such methods inevitably
involve a Monte Carlo component.

What makes a system difficult to sample in an equilibrium sense? Typically, these systems
have frustrated or rugged energy landscapes. This means that the underlying potential energy
surface contains many deep minima separated by high energy barriers. At low temperatures, it
is challenging for the system to move between the relevant low-energy regions because it must
surmount these barriers—a rare event process that takes time.

A wide range of methods have been designed to accelerate equilibration in simulations of

systems with rugged energy landscapes. Here, we will only review a small selection of these.
Some kinds of systems that might fall into this category include:

• systems with strong electrostatic interactions

• fluids of strongly orientation-dependent interactions (e.g., dipoles, hydrogen bonds)

• macromolecular and polymeric systems

• biomolecular systems (proteins, lipids)

• self-organizing or self-structuring systems (micelles, bilayers)

• supercooled liquid and amorphous systems at very low temperature

© M. S. Shell 2009 1/10 last modified 12/4/2019

Simulated tempering
The simulated tempering algorithm is equivalent to a canonical Monte Carlo simulation in
which the temperature changes randomly. That is, there are two kinds of moves:

• energy fluctuation moves – single-particle displacements, orientational displacements,

etc. The usual canonical acceptance criteria, at the current temperature, are used.

• temperature fluctuation moves – the temperature is periodically increased or de-

creased by fixed amounts so that the temperature can acquire one of a fixed set of val-
ues 𝑇min , 𝑇min + Δ𝑇, 𝑇min + 2Δ𝑇, … , 𝑇max .

These temperature fluctuations help the system equilibrate by allowing it to fluctuate to higher
temperatures where it can more readily cross energy barriers.

In these simulations, the microstate is now a function of both the configurations and the
temperature:
𝑁)
℘(𝐫 𝑁 , 𝑇) ∝ 𝑒 −𝛽𝑈(𝐫

What is the acceptance criterion for a temperature change? We want all temperatures to
appear with equal probability. That is, if we measured the distribution of temperature over the
simulation run, we would want it to be roughly uniform. To do this, we introduce a weighting
function in temperature:
𝑁 )+𝜂(𝑇)
℘(𝐫 𝑁 , 𝑇) ∝ 𝑒 −𝛽𝑈(𝐫

This means that the acceptance criterion for a random temperature perturbation 𝑇2 = 𝑇1 + Δ𝑇
is given by
acc
𝑃12 = min[1, 𝑒 −(𝛽2 −𝛽1 )𝑈+𝜂(𝑇2 )−𝜂(𝑇1 ) ]
= min[1, 𝑒 −𝑈Δ𝛽+Δ𝜂 ]

The distribution of temperatures is found by integrating this joint probability over configura-
tions:

℘(𝑇) ∝ ∫ ℘(𝐫 𝑁 , 𝑇)𝑑𝐫 𝑁

= 𝑍(𝑇, 𝑉, 𝑁)𝑒 𝜂(𝑇)

= 𝑒 𝜂(𝑇)−𝛽𝐴(𝑇)

To achieve a uniform distribution in temperature, we demand that this expression is constant.

This means that we should choose:

© M. S. Shell 2009 2/10 last modified 12/4/2019

 𝜂(𝑇) = 𝛽𝐴(𝑇)

We do not know the free energies at the different temperatures a priori. Instead, we must
iteratively determine them over the course of a simulation. Many methods can be used to do
this, including:

• flat histogram methods –the flat histogram coordinate is the temperature

• multiple histogram-reweighting – we maintain a histogram of energies at each temper-

ature and periodically use the reweighting equations to determine the relative free en-
ergies for each temperature

• free energy perturbation – we use, for example, Bennett’s method between adjacent
temperatures to compute free energy differences at periodic intervals in the simulation

How do we compute properties at a temperature of interest from this approach? We can use
the WHAM-based reweighting equation applied to all temperatures and all trajectories. We
first compute a configurational weight according to

𝑒 −𝛽𝑈𝑖
𝑤𝑖 =
∑𝐽𝑗=1 𝑒 𝛽𝑗𝐴𝑗−𝛽𝑗𝑈𝑖

Here, 𝑖 is an index over simulation observations. The inverse temperature 𝛽 corresponds to an

arbitrary reweighting temperature 𝑇. One should choose 𝑇min ≤ 𝑇 ≤ 𝑇max for good statistical
accuracy.

For an arbitrary observable 𝑋, we can then compute averages and distributions using

∑𝑖 𝑤𝑖 𝑋𝑖
⟨𝑋⟩ =
∑𝑖 𝑤𝑖
℘(𝑋) ∝ ∑ 𝑤𝑖 𝛿𝑋𝑖 ,𝑋
𝑖

Replica exchange
The simulated tempering approach is useful for facilitating equilibrium because it allows the
system to explore multiple temperatures. It has two drawbacks, however:

• We need to compute the free energies 𝐴(𝑇𝑗 ) at every temperature in order to properly
sample all temperatures.

© M. S. Shell 2009 3/10 last modified 12/4/2019

 • We need to wait for the system to traverse the entire temperature range many times in
order to accumulate good statistics at each. This can be a very long time if we have
many temperatures.

An alternative approach, which has become the method of choice for sampling challenging
systems is the replica exchange method. It overcomes these limitations by having the following
general construction:

• 𝐽 simulations (“replicas”) of the same system are performed simultaneously at different

temperatures 𝑇𝑗 .

• Each simulation is evolved independently, either through MD or MC methods, at the

corresponding temperature.

• At set intervals, replica swap moves are performed between adjacent temperature rep-
licas. In a swap move, the instantaneous configurations are exchanged between the
two temperatures.

Swap move acceptance criterion

The replica exchange simulation performs a Markov chain in the entire 𝐽-system ensemble. We
need to determine the acceptance criterion for performing swap moves. To do that, we need
to determine the total probability of one microstate in the entire ensemble. Here, a microstate
is the list of all of the positions in each of the replicas, 𝐑 = (𝐫1𝑁 , 𝐫2𝑁 , … , 𝐫𝐽𝑁 ). Since the replicas
do not interact (there are no energy terms between atoms in different replicas), we can write

℘(𝐑) = ∏ ℘𝑗 (𝐫𝑗𝑁 )
𝑗

Using canonical probabilities in each replica,

𝑁
𝑒 −𝛽𝑗𝑈(𝐫𝑗 )
℘(𝐑) = ∏
𝑍𝑗
𝑗

Consider a swap move between two temperatures 1 and 2. Initially the configuration in tem-
perature 1 is 𝐫1𝑁 and in temperature 2 is 𝐫2𝑁 . The acceptance criterion for the move stems from
the detailed balance equation,

acc
℘(𝐑 2 )
𝑃12,swap = min [1, ]
℘(𝐑1 )

© M. S. Shell 2009 4/10 last modified 12/4/2019

Here, 𝐑 2 corresponds to the set of configurations in which 𝐫1𝑁 and 𝐫2𝑁 are transposed to differ-
ent temperatures. We have:

𝑁 𝑁 𝐽 𝑁
𝑒 −𝛽𝑗𝑈(𝐫𝑗 )
𝑒 −𝛽1 𝑈(𝐫1 ) 𝑒 −𝛽2𝑈(𝐫2 )
℘(𝐑1 ) = × ×∏
𝑍1 𝑍2 𝑍𝑗
𝑗=3

𝑁 𝑁 𝐽 𝑁
𝑒 −𝛽𝑗𝑈(𝐫𝑗 )
𝑒 −𝛽1 𝑈(𝐫2 ) 𝑒 −𝛽2 𝑈(𝐫1 )
℘(𝐑 2 ) = × ×∏
𝑍1 𝑍2 𝑍𝑗
𝑗=3

Plugging these probabilities into the acceptance criterion,

acc 𝑁 𝑁 𝑁 𝑁
𝑃12,swap = min [1, 𝑒 𝛽1 𝑈(𝐫1 )+𝛽2𝑈(𝐫2 )−𝛽1 𝑈(𝐫2 )−𝛽2𝑈(𝐫1 ) ]

Notice that the partition functions in the denominator cancel. This means that we do not need
to know the free energies in each temperature when evaluating the acceptance criterion. It is
because we swap systems at two temperatures, rather than perturb a single system’s tempera-
ture, that we no longer need to know the free energy (as compared to simulated tempering).

Simplifying the above expression,

acc
𝑃12,swap = min[1, 𝑒 Δ𝛽Δ𝑈 ] where Δ𝛽 = 𝛽2 − 𝛽1 , Δ𝑈 = 𝑈(𝐫2𝑁 ) − 𝑈(𝐫1𝑁 )

Notice that we are required to compute the potential energy difference between the instanta-
neous configurations in the two temperature replicas.

Swap moves are performed like any other MC moves: the move is proposed, the acceptance
criterion is computed, a random number is drawn, and it is decided whether or not to perform
the move. Unlike many other MC moves, swap moves are very inexpensive to perform: they
only require the current energies in each temperature, which is typically maintained through-
out the simulation anyways.

Considerations for acceptance ratios and the temperature schedule

Consider a swap between two temperatures 𝑇2 > 𝑇1. We expect that the configuration drawn
from the higher temperature will have a higher energy. That is, 𝑈(𝐫2𝑁 ) > 𝑈(𝐫1𝑁 ). Thus the
quantity Δ𝛽Δ𝑈 will most often be negative, which will make the acceptance probability small.
For very large temperature differences, this quantity becomes even more negative. The result
is that the temperatures must be spaced close enough together in order to achieve a good rate
of accepted swaps.

© M. S. Shell 2009 5/10 last modified 12/4/2019

We can think about this problem in terms of the distributions ℘1 (𝑈), ℘2 (𝑈) between neigh-
boring temperatures. If there is a substantial overlap, we will have a high frequency of swap
move acceptance:

𝑇1 𝑇1
℘(𝑈) ℘(𝑈)
𝑇2 > 𝑇1 𝑇2 > 𝑇1

swaps frequently accepted 𝑈 swaps infrequently accepted 𝑈

These considerations play into the way in which we pick the temperatures in our replica ex-
change simulation?:

• 𝑇min - typically we pick a minimum temperature to be the temperature of interest that is

difficult to simulate

• 𝑇max - one chooses a high temperature where free energy barriers can be crossed, but
not so high as to require many intermediate temperatures

• 𝐽 - the number of temperatures is usually chosen so as to achieve ~50% acceptance of

swap moves between adjacent replicas

• 𝑇𝑗 - for a system with a constant heat capacity, it can be shown that a constant rate of
acceptance between adjacent temperatures corresponds to an exponential distribution
in temperature, equivalently, a power law in replica number. Thus, one normally picks
𝑗
𝑇max 𝐽−1
𝑇𝑗 = 𝑇min ( )
𝑇min

Scaling with system size

As the system size grows, the distribution in energy at a given temperature becomes increasing-
ly narrow with respect to the average energy, as 1/√𝑁. Without overlap between the energy
fluctuations of adjacent temperatures in the replica exchange scheme, few swap moves will be
accepted. Thus, as the system size increases, more and more intermediate temperature
replicas are needed to achieve a 50% acceptance ratio. This makes the method challenging to
apply to very large systems.

© M. S. Shell 2009 6/10 last modified 12/4/2019

Evaluating properties and distributions
Rigorously, all of the configurations at each temperature (regardless of whether they swapped
in from other temperatures) converge to a canonical distribution at that temperature. Thus, we
could compute the averages of properties at each temperature 𝑇𝑗 by simply averaging over the
configurations at that temperature. A trajectory is therefore considered the evolution of
configurations at a given temperature.

However, we can also compute an average at any arbitrary temperature between the minimum
and maximum in our simulation by using a reweighting approach. When a replica exchange
simulation is performed, we can collect histograms or histories of each potential energy visited
in each temperature, 𝑈𝑖𝑗 . We can then use the reweighting equations to compute the free
energies and configurational weights at each temperature. With the computed configurational
weights, we can express the average of any property as a weighted sum over all temperature
trajectories:

𝑒 −𝛽𝑈𝑖𝑗
𝑤𝑖𝑗 =
∑𝐽𝑙=1 𝑒 𝛽𝑙𝐴𝑙−𝛽𝑙𝑈𝑖𝑗

∑𝐽𝑗=1 ∑𝑛𝑖=1 𝑤𝑖𝑗 𝑋𝑖𝑗

⟨𝑋⟩ =
∑𝐽𝑗=1 ∑𝑛𝑖=1 𝑤𝑖𝑗
𝐽 𝑛

℘(𝑋) ∝ ∑ ∑ 𝑤𝑖𝑗 𝛿𝑋𝑖𝑗=𝑋

𝑗=1 𝑖=1

Implementation
To maintain the different replicas, one typically uses a parallel computing scheme, in which one
computer node (or processor) is assigned to each replica. When swaps are performed, all
nodes stop evolving the system in time and a master or head node sorts through and makes the
swap moves. This involves parallel communication between the different nodes and the head
node.

At each interval where swaps are performed, any number of swaps could be attempted. Typi-
cally, one chooses the number of swap attempts to be the number of temperature replicas.
The swaps can be attempted in serial order, from lowest/highest to highest/lowest tempera-
ture, or in random order.

In a MC simulation, the evolution of the system in between swaps can be governed by any
number of MC moves that accomplish changes in the potential energy, such as single particle
displacements.

© M. S. Shell 2009 7/10 last modified 12/4/2019

In a MD simulation, the evolution in between swaps can be performed using short MD trajecto-
ries. These moves can be considered hybrid MC/MD moves with good energy conservation
such that they are always accepted. After a round of swap moves, one picks random velocities
for each atom at every temperature 𝑗 from a Boltzmann distribution at the corresponding 𝑇𝑗 .
These velocities are used to start the short MD trajectory before another round of swap moves
is performed.

An alternative approach in MD would be to rescale the velocities in the configurations to the

new temperature after a swap move. This approach was first derived by Sugita and Okamoto
[1999]. Here, random new velocities are not picked at any time. Instead, after swapping
configurations, the momenta are scaled by a factor √𝑇new ⁄𝑇old , where 𝑇new is the new temper-
ature into which the configuration was swapped and 𝑇old is the temperature from which it
came. This approach requires a thermostat to be used during the MD trajectories. While the
velocity rescaling approach is more frequently used in the literature than the hybrid approach
(with random resampling), it is not necessarily more efficient and either approach is valid.

Variants
Here, we considered a replica exchange simulation where the replicas differed in temperature.
This is often also called parallel tempering for its similarity to simulated tempering. However,
we are not limited to temperature. One can have replicas differ in chemical potential or
pressure if the individual simulations are GCMC or 𝑁𝑃𝑇 simualtions, respectively. Two-
dimensional replica exchange methods allow replicas to differ in both temperature and chemi-
cal potential or pressure.

Another way to facilitate sampling is to modify the potential energy function itself. Here, one
wants to perturb the energy function so that the underlying energy landscape is smoother and
easier to sample. For example, one systematically scale the partial charges from zero to their
full values as one moves from replica to replica. This approach is often called Hamiltonian
exchange. The appropriate form of the acceptance criterion is derived as before, but with
different potential energy functions in each simulation 𝑈𝑗 :

𝑁 𝑁 𝑁 𝑁
acc
𝑃12,swap = min [1, 𝑒 𝛽1 𝑈1 (𝐫1 )+𝛽2 𝑈2 (𝐫2 )−𝛽1 𝑈1 (𝐫2 )−𝛽2 𝑈2 (𝐫1 ) ]

In evaluating this acceptance criterion, one must evaluate the cross-energies 𝑈2 (𝐫1𝑁 ) and
𝑈1 (𝐫2𝑁 ) between two replicas with each swap move, since the energy function is no longer
constant between them.

© M. S. Shell 2009 8/10 last modified 12/4/2019

In all of these cases, we can use multiple histogram reweighting techniques to express the
average or distribution of any property as a weighted average over all temperatures 𝑗 and
trajectory configurations 𝑖.

Extended ensemble molecular dynamics

One way to enhance the exploration of phase space in MD at lower temperatures is bias the
ensemble probabilities to artificially broad fluctuations in potential energy. This enables the
system to hop over higher-energy barriers more frequently. Ultimately, the behavior of the
system under normal conditions can be recovered by standard reweighting procedures. We
have already discussed this approach in the context of MC simulations. Now, we consider an
MD implementation.

We add to our MD simulation a weighting function that biases the configurational probabilities,
in the same way that we did during our discussion of biased sampling:

𝑈 w (𝐫 𝑁 ) = 𝑈(𝐫 𝑁 ) − 𝑘𝐵 𝑇𝜂(𝑈(𝐫 𝑁 ))

Here, 𝜂(𝑈) is a weighting function with dimensionless units, and it only depends on the poten-
tial energy. The equations of motion are derived as before. For an atom 𝑖:

𝑑𝑈 w
𝐟𝑖w = −
𝑑𝐫𝑖
𝑑𝑈 𝑑𝜂(𝑈) 𝑑𝑈
=− + 𝑘𝐵 𝑇
𝑑𝐫𝑖 𝑑𝑈 𝑑𝐫𝑖
𝑑𝜂(𝑈)
= 𝐟𝑖 (1 − 𝑘𝐵 𝑇 )
𝑑𝑈
In other words, the force on each atom in the weighted ensemble is scaled by a term involving
the derivative of the weighting function evaluated at the energy of the current configuration. If
the weighting function is a constant, no scaling occurs.

In order for this approach to work, the weighting function must be a continuous function. One
typically uses splines or other mathematical constructs so that an analytical derivative can be
computed.

How do we determine the weighting function? We can fit it to estimates of the probability
distribution function. In the canonical ensemble, we have the relationship

℘(𝑈)
∝ 𝑒 −𝜂(𝑈)
℘w (𝑈)

If we want the weighted distribution to be approximately flat, we can choose

© M. S. Shell 2009 9/10 last modified 12/4/2019

 𝜂(𝑈) = − ln ℘(𝑈) + const

Notice that we will have to fit the weighting function to discrete histogram calculation, since 𝜂
must be continuous. It is for this reason that it becomes difficult to evaluate 𝜂 to high statistical
accuracy using MD techniques. Still, broadening the energy distribution even a little bit may
speed equilibration in MD simulations.

© M. S. Shell 2009 10/10 last modified 12/4/2019