Biased sampling and related free energy techniques ChE210D

Today's lecture: creating artificial ensembles in simulations that facilitate sam-

pling and the computation of free energetic quantities.

Biased or extended ensembles

So far, we have discussed simulations in standard statistical-mechanical ensembles. However,
we are free to perform a simulation in any ensemble of interest, including artificial ones of our
own design. These biased ensembles or extended ensembles have some benefits:

• We can achieve broader sampling of particular states that are rarely visited, in order to
compute properties of them. The usual ensembles give rise to sharply peaked distribution
functions whose relative width scales as 𝑁 −1⁄2 ; extended ensembles allow us to broaden
these fluctuations.

• We can achieve enhanced sampling of states in a way that promotes faster equilibration
times.

• By designing artificial ensembles that visit frequently states that are normally rarely-sam-
pled, we can compute free energies of them to much higher accuracy than would nor-
mally be the case.

• Artificial ensembles enable us to couple sampling with the calculation of free energies
through iterative feedback loops. We will discuss this approach in greater depth later.

By biased ensemble, we mean that we specify a configurational /microstate probability distribu-

tion ℘(𝐫 𝑁 ) that does not extend from a standard ensemble (𝑁𝑉𝐸, 𝑁𝑉𝑇, 𝑁𝑃𝑇, 𝜇𝑉𝑇). We use this
custom probability distribution to design a Monte Carlo simulation through the usual detailed
balance equation for determining acceptance probabilities. Alternatively, we could also design a
molecular dynamics simulation in a biased ensemble, but this is not as straightforward or simple.
For the remainder of this lecture, we will consider primarily Monte Carlo simulations, but we will
come back to biased MD runs in a later lecture.

Ultimately, reality exists in standard ensembles. Therefore, any simulation of a system that we
perform in a biased ensemble must somehow be related back to results for that system in usual
ensembles. To do this, we can use reweighting techniques. In fact, biased ensemble simulations
are closely tied to reweighting methods and the two are almost always performed hand-in-hand
in modern simulation studies.

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

Keep in mind that the specification of the ensemble is entirely separate from the specification of
the Monte Carlo move sets. We can perform a given MC move in any ensemble that we like, so
long as we properly derive the acceptance criterion using the ensemble distribution ℘(𝐫 𝑁 ).

General formalism and connection to statistical mechanics

For the sake of simplicity, we will use a discrete notation to indicate the probability distribution,
as we did in our initial discussion of MC methods:

℘𝑚

where 𝑚 denotes one configurational microstate of the system, e.g., a set of coordinates 𝐫 𝑁 .

Consider a MC simulation in the canonical ensemble, where we specify

℘u𝑚 ∝ 𝑒 −𝛽𝑈𝑚

where “u” signifies that we are in the unweighted ensemble. We now want to perform a biased
simulation in this ensemble where we introduce a weighting factor for the microstates in the
exponential:

−𝛽𝑈𝑚 +𝜂𝑚
℘w
𝑚 ∝𝑒

The weighting function 𝜂𝑚 modifies the canonical distribution so that some configurations have
higher or lower probabilities than would normally be expected. It gives a value that is dependent
on the particular configuration of interest. As we will see shortly, this function can take a variety
of forms.

We can think of the weighting function as giving rise to an effective "weighted" potential energy
function:
w
𝑒 −𝛽𝑈𝑚 +𝜂𝑚 ≡ 𝑒 −𝛽𝑈𝑚

or,

𝑈 w ≡ 𝑈 − 𝑘𝐵 𝑇𝜂

Notice that if 𝜂𝑚 = const, then we recover the usual canonical distribution because its effect is
removed via the normalization condition:

𝑒 −𝛽𝑈𝑚 +𝜂𝑚
℘w
𝑚 =
∑𝑚′ 𝑒 −𝛽𝑈𝑚′ +𝜂𝑚′

Comparing the weighted and unweighted ensembles,

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

 ℘u𝑚
∝ 𝑒 −𝜂𝑚
℘w𝑚

Imagine that we perform a MC simulation in the weighted ensemble. For symmetric moves, our
acceptance criterion would follow
acc
𝑃12 ℘w2
acc = w
𝑃21 ℘1
= 𝑒 −𝛽(𝑈2 −𝑈1 )+𝜂2 −𝜂1
= 𝑒 −𝛽Δ𝑈+Δ𝜂

With the Metropolis form,

acc
𝑃12 = min[1, 𝑒 −𝛽𝛥𝑈+𝛥𝜂 ]

Ultimately, the MC simulation would generate a trajectory of configurations according to the

distribution ℘w
𝑚 . We could compute the average of any property 𝐴 in the weighted ensemble,

1
⟨𝐴⟩w = ∑ 𝐴𝑖
𝑛
where the index 𝑖 = 1, … , 𝑛 is an index over the trajectory configurations.

What if we wanted to compute the average in the unweighted ensemble (the usual canonical
average)? We could still perform a simulation in the weighted ensemble, but our expression for
the unweighted average would need to take into account differences in ℘w u
𝑚 and ℘𝑚 :

℘u𝑖
∑ 𝐴𝑖
℘w 𝑖
⟨𝐴⟩u =
℘u𝑖
∑ w
℘𝑖

Substituting from above,

∑ 𝐴𝑖 𝑒 −𝜂𝑖
⟨𝐴⟩u =
∑ 𝑒 −𝜂𝑖

Here, 𝜂𝑖 gives the value of the weighting function for each configuration included in the trajectory
𝑖. The denominator in this expression serves as a normalization for the unweighting expression.

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

Practical and statistical issues
To perform such an unweighting of weighted simulation results, we need keep lists of the values
𝐴𝑖 and 𝜂𝑖 for each of the trajectory configurations. We also must rearrange the exponential so
that we don’t run out of precision:

∑ 𝐴𝑖 𝑒 −𝜂𝑖 +𝜂min
⟨𝐴⟩u = 𝜂min = min 𝜂𝑖
∑ 𝑒 −𝜂𝑖 +𝜂min 𝑖

Keep in mind that the error in this expression is determined by the correlation times of the ob-
servable in the weighted ensemble.

The error is also highly dependent on whether the distribution of 𝐴 in the unweighted ensemble
is similar to the distribution in the weighted ensemble. If ⟨𝐴⟩u is rarely sampled in the weighted
simulation, then the statistics of this unweighting procedure will produce a poor estimate of the
unweighted average.

Typically, however, one explicitly selects a weighting function that allows the simulation to ex-
plore a very wide range of 𝐴 to accumulate good statistics, much wider than runs in conventional
ensembles. In fact, one of the primary purposes of the artificial ensemble is to generate a broad
sampling of 𝐴.

Umbrella sampling
The umbrella sampling method was developed by Torrie and Valleau in 1977 and since has been
one of the major approaches for performing simulations along predetermined reaction coordi-
nates. It is readily applied to both molecular dynamics and Monte Carlo simulations. Here, we
use a simple example for demonstrating the application of the approach; however, keep in mind
that umbrella sampling can be applied to many, arbitrary reaction coordinates.

Example: reaction coordinate

In this example, we consider the interaction of a spherical solute in solution with a surface:

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

 𝑧

Though we have not shown it explicitly, the solute and the surface atoms both interact with the
solution atoms, which can translate, rotate, and move throughout the system volume. Here, we
assume the canonical ensemble.

Consider the probability distribution of the molecule from the surface, which relates to the po-
tential of mean force (PMF):

𝑁
℘(𝑧) ∝ ∫ 𝑒 −𝛽𝑈(𝐫 ) 𝛿[𝑧 − 𝑧(𝐫 𝑁 )]𝑑𝐫 𝑁

∝ 𝑒 −𝛽𝐹(𝑧)

If the solute interacts strongly and favorably with the surface, we might expect a probability dis-
tribution of the coordinate 𝑧 to look something like:

ln ℘(𝑧) 𝛽Δ𝐹 ‡

𝑧
𝑧‡

If the free energy barrier for a surface-bound solute for escaping is fairly large, we will expect an
extremely small probability for the system to adopt states other than the bound state. That is,

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

we would almost never see a solute molecule far from the surface. This would cause problems
if we wanted to compute the free energy 𝐹(𝑧) for large 𝑧 values since the statistics of our histo-
gram would be so few and poor there.

To attain better statistics, we can construct biased ensembles that allow us to sample the full
range of 𝑧 space to accrue many counts in our histogram bins. The umbrella sampling procedure
is as follows:

• Perform 𝐽 simulations 𝑗 of the same system.

• In each simulation, restrain the system to sample a small range of 𝑧 values centered
around 𝑧𝑗 . This restraint can be added as an energetic penalty for configurations that
venture far from the value 𝑧𝑗 . That is, we create an extended ensemble that biases us
towards configurations near 𝑧𝑗 :

−𝛽𝑈jw (𝐫 𝑁 ) = −𝛽𝑈(𝐫 𝑁 ) + 𝜂𝑗 (𝑧)

• Use a different target value 𝑧𝑗 for each simulation. The 𝑧𝑗 should span the entire range of
interest.

• Measure the weighted ensemble distribution ℘w

𝑗 (𝑧) for each simulation using histo-
grams. Let the variable 𝑐𝑗 (𝑧) denote the counts in simulation 𝑗 for different 𝑧 observa-
tions. Notice that one has to discretize the coordinate 𝑧 into bins since it is a continuous
variable.

• Unweight and stitch together all of the simulation ℘w

𝑗 (𝑧) to produce the unweighted un-
derlying free energy function 𝐹(𝑧) = −𝑘𝐵 𝑇 ln ℘(𝑧).

Form of the biasing potential

In each simulation, we need to bias the run to sample values of 𝑧 near 𝑧𝑗 . We want to energeti-
cally penalize configurations outside of this range so as to lower their probability weight. Many
forms of this penalty can be chosen. The most common is a harmonic potential:

𝑘 2
𝜂𝑗 (𝑧) = −𝛽 (𝑧 − 𝑧𝑗 )
2
such that the effective weighted potential is given by

𝑘 2
𝑈jw (𝐫 𝑁 ) = 𝑈(𝐫 𝑁 ) + (𝑧 − 𝑧𝑗 )
2
Here a value of the force constant must be specified. Too small a value will not sufficiently bias
the simulation. Too large a value will result in a very narrow distribution in the sampled 𝑧 in each

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

simulation; this will result in poor overlap between the 𝑧 distributions in each simulation and will
make hard to patch the unweighted results together with good statistical accuracy. In practice,
it can be a trial-and-error process to determine good values for 𝑘.

The advantage of the harmonic potential is that it enables us to apply either MD or MC algorithms
to compute the distribution ℘w𝑗 (𝑧). In MD, the harmonic term will add an extra force to the
solute in its z-coordinate. In MC, this term will modify the energy used in the acceptance crite-
rion.

Reweighting procedure
Each simulation is connected to the unweighted distribution via the following considerations:
𝑁 −𝜂𝑗 (𝑧)
℘(𝐫 𝑁 ) ∝ ℘w
𝑗 (𝐫 )𝑒

Here, we are omitting the superscript "u" from the unweighted distribution, and will continue to
do so throughout. Following this relationship, if we integrate this distribution over all of the
coordinates except for 𝑧:
−𝜂𝑗 (𝑧)
℘(𝑧) ∝ ℘w
𝑗 (𝑧)𝑒

Taking the logarithm,

𝐹(𝑧) = −𝑘𝐵 𝑇 ln ℘w
𝑗 (𝑧) + 𝑘𝐵 𝑇𝜂𝑗 (𝑧) + const

Notice that each simulation 𝑗 should return, in principle, the same distribution ℘(𝑧). We have
forced the different distributions to span a broad range of 𝑧 solely to attain accurate statistics
over the full range of 𝑧.

We could stitch together the 𝐽 estimates of 𝐹(𝑧) by shifting the unknown constants so as to
obtain overlap in the common regions:

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

 ℘(𝑧)

𝑗=1 𝑗=2 𝑗=3

ln ℘w
𝑗 (𝑧)
− 𝜂𝑗 (𝑧)
+ const

We need to have sufficient overlap between the distributions in each simulation in order to do
this effectively.

This approach may be very reminiscent of the histogram reweighting issues that we discussed in
the previous lecture. Indeed, this scenario extends from a generalization of that procedure to
arbitrary parameters (here, the reaction coordinate 𝑧).

Rather than visually shift these curves to overlap, a more robust approach to finding the under-
lying free energy function would be to use maximum likelihood arguments as we did earlier. In-
deed, the multiple histogram reweighting method (WHAM) can be derived for such cases. One
needs to take into account the different energy functions used in the 𝑗 simulations in this case.
The final equations for 𝐹(𝑧), which must be solved iteratively, are:

−𝛽𝐹(𝑧) = ln 𝑐𝑡𝑜𝑡 (𝑧) − ln 𝑛 − ln ∑ 𝑒 𝜂𝑗(𝑧)+𝛽𝐴𝑗 for all 𝑧

𝑗=1

−𝛽𝐴𝑗 = ln ∑ 𝑒 −𝛽𝐹(𝑧)+𝜂𝑗(𝑧) for all 𝑗

𝑧

Here, the 𝐴𝑗 give the weighted ensemble free energy of each simulation 𝑗.

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

Flat-histogram sampling
In the past two decades, so-called flat histogram methods have become major tools for compu-
ting free energies and phase equilibria. These methods are designed to construct an extended
ensemble that generates a uniform or flat distribution in one or more parameters or reaction
coordinates. By doing so, two goals are achieved:

• broad sampling and good statistics for the flat histogram parameters, which expands the
range of conditions at which reweighting can be effectively performed

• computation of free energies or entropies along the flat histogram parameters, through
a connection between them and the presence of a flat distribution

Such methods are beginning to supersede umbrella sampling approaches in Monte Carlo simula-
tions because they automatically determine the weight functions, instead of requiring a specific
form (e.g., a harmonic potential). Below we discuss a specific example in the grand-canonical
ensemble; however, these approaches are general to any simulation ensemble and parameter of
interest.

Example: Grand-Canonical Monte Carlo

Consider a grand canonical MC simulation with 𝜇 = 𝜇1 , where both the energy and particle num-
ber fluctuate. The probability of seeing a particular configuration is given by

𝑁
𝑒 −𝛽𝑈+𝛽𝜇1 𝑁
℘(𝐫 , 𝑁) ∝
Λ(𝑇)3𝑁 𝑁!

If we integrate this expression over all configurations at a particular particle number, we can find
the distribution in 𝑁:

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

𝑒 −𝛽𝑈+𝛽𝜇1 𝑁 𝑁
=∫ 𝑑𝐫
Λ(𝑇)3𝑁 𝑁!
𝑒 −𝛽𝑈
= 𝑒 𝛽𝜇1 𝑁 ∫ 𝑑𝐫 𝑁
Λ(𝑇)3𝑁 𝑁!
= 𝑒 𝛽𝜇1 𝑁−𝛽𝐴(𝑇,𝑉,𝑁)

That is, the particle number distribution involves the Helmholtz free energy 𝐴(𝑇, 𝑉, 𝑁). We will
actually drop the variables 𝑇, 𝑉 since these do not fluctuate during the simulation:

℘(𝑁) ∝ 𝑒 𝛽𝜇1 𝑁−𝛽𝐴(𝑁)

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

Now consider that we have specified a temperature and chemical potential that places us at co-
existence between a liquid and vapor phase. At coexistence, this distribution might look some-
thing like:

ln ℘(𝑁) 𝛽Δ𝐴‡

𝑁
𝑁‡

The two probability peaks correspond to the liquid and gas phase. Notice that there is an inter-
mediate value of 𝑁 with much lower probability. The probability that the system will visit this
value of 𝑁 is
‡
℘(𝑁 ‡ ) ∝ 𝑒 −𝛽Δ𝐴

Small free energy differences can make the probability of intermediate values of 𝑁 very, very tiny
due to the exponential. Thus, even though the system is at coexistence per the bimodal distri-
bution above, the implication is that fluctuations that traverse intermediate densities are very
rare.

In a GCMC simulation, this means that we would rarely see the system interconvert between low
and high particle numbers due to the low probability of visiting intermediate particle number
states between them. We would have an extremely difficult time equilibrating a GCMC simula-
tion at these conditions since it almost certainly would visit only one of the two phases during
the simulation run, and not both.

We can use an extended ensemble to enhance the probability of intermediate states. One way
to do this would be to modify the ensemble probabilities so that we would obtain a uniform
distribution or a flat histogram of particle numbers, between two predetermined limits 𝑁min and
𝑁max . Typically we pick 𝑁min = 1 and 𝑁max to be well beyond the particle number associated
with the liquid density.

To do this, we can add a weighting function to the ensemble probabilities that depends on 𝑁:

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

 w (𝐫 𝑁
𝑒 −𝛽𝑈+𝛽𝜇1 𝑁+𝜂(𝑁)
℘ , 𝑁) ∝
Λ(𝑇)3𝑁 𝑁!

This is the function that we use in the determination of our acceptance criteria. Following
through the detailed balance equation, we find that:

acc
𝑉 ′
𝑃12 = min [1, 𝑒 −𝛽Δ𝑈+𝛽𝜇1 +Δ𝜂 ] for insertions
𝑁+1

acc
𝑁 −𝛽Δ𝑈−𝛽𝜇′ +Δ𝜂
𝑃12 = min [1, 𝑒 1 ] for deletions
𝑉

In each case, Δ𝜂 = 𝜂(𝑁2 ) − 𝜂(𝑁1 ).

Choice of weighting function and reweighting procedure

How do we pick the function 𝜂(𝑁)? This function would be tabulated as an array in our simula-
tion and would be defined in the range [𝑁min , 𝑁max ]. We want to pick this function so that our
ultimate probability distribution in 𝑁 looks flat:

ln ℘w (𝑁)

𝑁min 𝑁max 𝑁

We can figure out what we need to pick for 𝜂(𝑁) by examining the expected distribution in the
weighted ensemble:

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

= 𝑒 𝛽𝜇1 𝑁−𝛽𝐴(𝑁)+𝜂(𝑁)

We want this distribution to be flat, ℘w (𝑁) = const. This gives us a way to solve for 𝜂(𝑁):

𝜂(𝑁) = 𝛽𝐴(𝑁) − 𝛽𝜇1 𝑁 + const

Notice two important aspects of this result:

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

 • We do not need to know the additive constant in 𝜂(𝑁), since the acceptance criterion
depends only upon differences in this function at different values of 𝑁. Ultimately this is
because additive shifts in 𝜂(𝑁) do not affect microstate probabilities, via the normaliza-
tion condition. Typically we simply set 𝜂(𝑁min ) = 0 as a working assumption.

• To find 𝜂(𝑁) that gives a flat histogram, we need to know the Helmholtz free energy as a
function of 𝑁. This may seem like a difficult task. In practice, however, this connection
provides us with a way to determine 𝐴(𝑁): given a simulation with a trial function 𝜂(𝑁),
if our biased simulation produces a uniform distribution in 𝑁, then we have computed
the true 𝐴(𝑁). Various flat histogram techniques discussed below enable us to compute
𝐴(𝑁) based on this connection.

Let’s say that we find some approximate 𝜂(𝑁) that gives rise to the distribution

ln ℘w (𝑁)

𝑁min 𝑁max 𝑁

Even though this distribution is not totally flat, we still sample the intermediate values of 𝑁 with
much higher probability and our simulation will likely alternate between liquid and vapor densi-
ties with greater frequency than would be present in the unweighted ensemble. As a result, our
simulation reaches equilibrium faster and we get good statistics at all particle numbers.

We measure the function ℘w (𝑁) from a histogram. Now, we want to convert the measured
℘w (𝑁) back to the unweighted ℘(𝑁), to compute the expected distribution (and averages) in
the normal grand-canonical ensemble:

℘(𝐫 𝑁 , 𝑁) ∝ 𝑒 −𝜂(𝑁) ℘w (𝐫 𝑁 , 𝑁)

Integrating over the particle positions,

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

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

Using this approach, we could reweight to different chemical potentials than the original chemi-
cal potential. In the unweighted ensemble,

℘(𝑁; 𝜇2 ) ∝ ℘(𝑁; 𝜇1 )𝑒 𝛽(𝜇2 −𝜇1)

Making the above substitutions:

℘(𝑁; 𝜇2 ) ∝ ℘w (𝑁; 𝜇1 )𝑒 −𝜂(𝑁)+𝛽(𝜇2 −𝜇1 )

The constant of proportionality is given by the normalization condition. This equation enables
the following general procedure:

• Specify a chemical potential 𝜇1 and approximate weighting function 𝜂(𝑁) ≈ 𝛽𝐴(𝑁) −

𝛽𝜇1 𝑁.

• Perform a simulation in the weighted ensemble and measure ℘w (𝑁; 𝜇1 ) using a histo-
gram.

• By using a weighted ensemble with a roughly flat distribution in 𝑁, we are accumulating

good statistics in our histogram for a wide range of 𝑁 values. Such would not be the case
if we performed a traditional GCMC simulation, which has a narrowly peaked distribution
of 𝑁.

• Use the reweighting equation to find the true ℘(𝑁; 𝜇2 ) at any arbitrary chemical poten-
tial 𝜇2 . The quality of the reweighted distribution is limited only by regions where the
measured ℘w (𝑁; 𝜇1 ) has few histogram entries.

Ultimately this approach enables us to find conditions of phase equilibrium to high accuracy: we
tune the reweighting chemical potential 𝜇2 until the weight under the two peaks in the probabil-
ity distribution is equal:

𝜇2 < 𝜇coexist (𝑇)

𝜇2 = 𝜇coexist (𝑇)
ln ℘(𝑁; 𝜇2 ) 𝜇2 > 𝜇coexist (𝑇)

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

The multicanonical method
How do we find the optimal weighting function 𝜂(𝑁) such that a histogram of 𝑁 in the weighted
ensemble is flat? Berg and Neuhaus in 1992 devised an iterative solution to this problem called
the multicanonical method. The basic idea is to perform a serial series of 𝐽 simulations in which
the weighting function is updated after each:

• Consider simulations number 𝑗 and 𝑗 + 1.

• We perform simulation 𝑗 using weighting function 𝜂𝑗 (𝑁) and measure ℘w

𝑗 (𝑁).

• If ℘w
𝑗 (𝑁) is flat, then 𝜂𝑗 (𝑁) has converged and we are done.

• If it is not flat, we use deviations of ℘w

𝑗 (𝑁) to update 𝜂𝑗 (𝑁) to get 𝜂𝑗+1 (𝑁).

• The process is repeated until convergence.

We can derive a simple update rule for extracting 𝜂𝑗+1 (𝑁) from 𝜂𝑗 (𝑁) by comparing the ex-
pected weighted ensemble distributions:

℘w
𝑗 (𝑁) ∝ ℘(𝑁)𝑒
𝜂𝑗 (𝑁)

℘w
𝑗+1 (𝑁) ∝ ℘(𝑁)𝑒
𝜂𝑗+1 (𝑁)

Dividing these two equations and taking the logarithm gives,

ln ℘w w
𝑗+1 (𝑁) − ln ℘𝑗 (𝑁) = 𝜂𝑗+1 (𝑁) − 𝜂𝑗 (𝑁) + const

We demand that ℘w
𝑗+1 (𝑁) = const so that we can find the optimal 𝜂𝑗+1 (𝑁) at the next iteration
of the weighting function. Applying this constraint and rearranging,

𝜂𝑗+1 (𝑁) = 𝜂𝑗 (𝑁) − ln ℘w

𝑗 (𝑁) + const

This equation provides us with an update rule for determining the next weighting function upon
each iteration. Note that,

• We can only know the weights to within an additive constant. This is not a problem,
however, because only weight differences appear in the acceptance criterion. Typically
we demand 𝜂(𝑁min ) = 0 with each of these updates.

• If ℘w
𝑗 (𝑁) = const, then every value in the weights is shifted by the same amount. Since
this effectively only serves to change the arbitrary additive constant, it therefore leaves
the weights unchanged. This is a crucial feedback property of the method: if we measure
a flat histogram, then we have converged to the desired 𝜂(𝑁).

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

Zero-entry bins
When we measure ℘w
𝑗 (𝑁), we use a histogram of counts 𝑐𝑗 (𝑁). However, this histogram can
lead to a zero probability in bins that have no counts. This would make the update equation
above ill-defined due to the logarithm term. A statistical analysis by Smith and Bruce [J. Phys. A,
1995] shows that the following update equation is a more reliable statistical estimator and avoids
this problem:

𝜂𝑗+1 (𝑁) = 𝜂𝑗 (𝑁) − ln[𝑐𝑗 (𝑁) + 1] + const

Convergence
Several iterations can be required for convergence. The following shows the results for four it-
erations using the Lennard-Jones system at 𝑉 = 125 and 𝜇 ′ = −3.7.

25
20 1 2 3 4
(N)

15
10
5
0
6
10
5
10
cj(N)

4
10
3
10
2
10
50 100
N
Note that, even if the 𝜂(𝑁) is not fully converged, a “flat enough” histogram can be used to
compute the unweighted average via

℘(𝑁; 𝜇2 ) ∝ ℘w (𝑁; 𝜇1 )𝑒 −𝜂(𝑁)+𝛽(𝜇2 −𝜇1 )

Error analysis
The multicanonical approach provides a way to estimate 𝜂(𝑁) that gives a fairly flat histogram.
In principle this method should give 𝜂(𝑁) = 𝛽𝐴(𝑁) − 𝛽𝜇1 𝑁 + const for a perfectly flat histo-
gram. Thus, we might be able to estimate from the final iteration

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

 𝐴(𝑁) = 𝜇1 𝑁 + 𝑘𝐵 𝑇𝜂𝐽 (𝑁) + const

In practice, the multicanonical method is not able to resolve 𝜂(𝑁) to sufficient statistical accuracy
so that we can determine 𝐴(𝑁) to high quality. Each iteration introduces statistical error into
this estimate since the histograms employed have a finite number of counts.

Wang-Landau sampling
The Wang-Landau method [Wang and Landau, PRL 86, 2050 (2001); PRE 64, 056161] is a very
recent approach to estimating 𝜂(𝑁) to high statistical accuracy. It overcomes some of the sta-
tistical problems with the multicanonical method and has been rapidly adopted as one of the
major tools in flat histogram algorithms. Our discussion here entails a slightly different presen-
tation than the original formulation, for clarity in connecting the approach to the current exam-
ple.

In this approach, we do not iterate over different simulations to find the optimal 𝜂(𝑁) that gives
a flat histogram. Rather, we modify this function at every MC step in a way that enforces a flat
distribution and allows it to converge to its optimal value. After every MC move, we update our
running estimate for 𝜂(𝑁) via

𝜂(𝑁) ← 𝜂(𝑁) − 𝑔

Here, 𝑔 is termed the modification factor. It is a number that we choose so as to gauge the rate
at which the weighting function is updated. Typically, simulations begin with 𝑔 = 1.

Think about what the behavior of this update would be:

• Values of 𝑁 that are oversampled will have their weights decreased, on average, more
than values of 𝑁 that are undersampled, since there will be more updates to them

• When all particle numbers are sampled with equal probability, on average, we will update
𝜂(𝑁) uniformly across all 𝑁. This results in a net additive shift in the weights, which ulti-
mately doesn’t affect the microstate probabilities. Thus, if we are sampling a uniform
distribution in 𝑁, the weights are not effectively modified. In any case, we always demand
that 𝜂(𝑁min ) = 0 by shifting the curve.

Thus, the Wang-Landau method enforces a kind of feedback loop between the sampled distribu-
tion of 𝑁 and the determination of the weighting function.

Modification factor schedule

Still, one must address two important issues:

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

 • Modifying the weighting factor at every step breaks detailed balance, because it changes
the state probabilities with time.

• We cannot resolve 𝜂(𝑁) to differences less than the modification factor 𝑔.

To resolve the issues, Wang and Landau proposed that 𝑔 → 0 over the course of the simulation
run. They suggested the following procedure:

1. Set 𝜂(𝑁) = 0 and 𝑔 = 1 initially.

2. Perform a simulation, updating 𝜂(𝑁) at every simulation step using the modification fac-
tor. Start a new histogram of 𝑁 and collect observations throughout the run.

3. When the histogram is “flat enough”, scale down the value of the modification factor ac-
cording to:

1
𝑔← 𝑔
2
4. Re-zero the histogram counts and return to step 2. Continue until 𝑔 is very small. Typi-
cally, we stop when 𝑔 < 10−6 − 10−8.

The histogram can be determined to be flat enough using the 80% rule: when the number of
counts in the least-visited histogram bin is no less than 80% of the average number of counts
over all bins, the histogram can be considered “flat”.

In this way, the simulation proceeds in stages of decreasing values of the modification factor.
Initial stages help to rapidly build a good estimate of 𝜂(𝑁), while later stages refine this calcula-
tion to increasing precision and satisfy detailed balance asymptotically.

Reweighting of results
At the end of the simulation, we have computed 𝜂(𝑁) to such high accuracy that we can use it
directly to perform reweighting. That is, we assume that ℘w (𝑁) = const:

℘(𝑁; 𝜇2 ) ∝ 𝑒 −𝜂(𝑁)+𝛽(𝜇2 −𝜇1 )

In fact, 𝜂(𝑁) provides a high-quality estimate of the free energy:

𝐴(𝑁) = 𝜇1 𝑁 + 𝑘𝐵 𝑇𝜂(𝑁) + const

Choice of initial state

Since the Wang-Landau method automatically and dynamically determines the weighting func-
tion, it is typical to choose the initial state of the simulation so that the only factor appearing in

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

the probability distribution is the weighting function itself. In this case, we would choose 𝜇1 = 0
so that:

𝐴(𝑁) = 𝑘𝐵 𝑇𝜂(𝑁) + const

and the reweighting equation becomes

℘(𝑁; 𝜇2 ) ∝ 𝑒 −𝜂(𝑁)+𝛽𝜇2

Transition matrix methods

The most recent methods to have emerged in biased simulations are those based on transition
matrix estimators. These have been shown to be very easy to implement and to provide very
high quality estimates of free energies, perhaps better than the Wang-Landau approach. Erring-
ton and coworkers have pioneered the application of these methods to fluid phase equilibria
[Erringron, JCP 118, 9915 (2003)]. Transition matrix approaches are applied to Monte Carlo sim-
ulations because they rely on the detailed balance equation.

The idea of transition matrix estimators is that we measure macrostate transition probabilities
and use these to compute underlying free energies. A macrostate transition probability is simply
the probability associated with transitions of the system between different values of some mac-
roscopic parameter. In this example, we measure the conditional probability that a system will
make a transition between one value 𝑁1 to another value 𝑁2 , given that it initially has 𝑁1 parti-
cles:

Π(𝑁1 → 𝑁2 )

Notice that, in a grand canonical simulation, we only make incremental steps in 𝑁. Thus,
Π(𝑁1 → 𝑁2 ) = 0 if 𝑁2 is not one of (𝑁1 − 1, 𝑁1 , 𝑁1 + 1).

Definition of the macroscopic transition probability

We can relate the macroscopic transition probability to the microscopic ones that we presented
in the initial discussion of MC simulations:

∑𝑚∈{𝑁1 } ∑𝑛∈{𝑁2 } ℘𝑚 𝜋𝑚𝑛

Π(𝑁1 → 𝑁2 ) =
∑𝑚∈{𝑁1 } ℘𝑚

Here, the sums over microstates 𝑚 and 𝑛 are performed for all states with 𝑁1 and 𝑁2 particles,
respectively.

If the microstate transition probabilities obey detailed balance:

℘𝑚 𝜋𝑚𝑛 = ℘𝑛 𝜋𝑛𝑚

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

then a simple summation of this equation over all 𝑚 and 𝑛 shows that the macrostate transition
probabilities also obey detailed balance:

℘(𝑁1 )Π(𝑁1 → 𝑁2 ) = ℘(𝑁2 )Π(𝑁2 → 𝑁1 )

Estimating free energies from transition matrices

The equation above allows us to estimate free energies from macroscopic transition probabili-
ties. We rearrange it according to

Π(𝑁1 → 𝑁2 ) ℘(𝑁2 )
ln = ln
Π(𝑁2 → 𝑁1 ) ℘(𝑁1 )
= 𝛽𝜇(𝑁2 − 𝑁1 ) − 𝛽[𝐴(𝑁2 ) − 𝐴(𝑁1 )]

Solving for the free energy difference,

Π(𝑁2 → 𝑁1 )
𝐴(𝑁2 ) − 𝐴(𝑁1 ) = 𝜇(𝑁2 − 𝑁1 ) + ln
Π(𝑁1 → 𝑁2 )

Looking at neighboring particle numbers,

Π(𝑁 + 1 → 𝑁)
𝐴(𝑁 + 1) − 𝐴(𝑁) = 𝜇 + ln
Π(𝑁 → 𝑁 + 1)

Using this equation, we could map out an entire 𝐴(𝑁) curve by computing free energy differ-
ences at each particle number. To do that, we would need to measure the relative probabilities
of seeing transitions between 𝑁1 and 𝑁2 . We can do this using a histogram:

𝑐𝑁1 ,𝑁2 = counts of observations of transitions from 𝑁1 to 𝑁2

Then,
𝑐𝑁1,𝑁2
Π(𝑁1 → 𝑁2 ) =
∑𝑖 𝑐𝑁1 ,𝑁𝑖

We can actually do better than this. Instead of tallying counts, we can tally the actual acceptance
probabilities computed for use in the Metropolis criterion:

𝑐𝑁1 ,𝑁2 = sums of observations of 𝑃𝑁acc

1 ,𝑁2

Computing the weighting function from transition probabilities

The transition probabilities provide us with a way to periodically update an estimate for 𝜂(𝑁):

𝜂(𝑁) = 𝛽𝐴(𝑁) − 𝛽𝜇𝑁 + const

where 𝐴(𝑁) is determined using the above procedure.

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

Note that, to determine 𝐴(𝑁) using the transition probabilities, we need to be in the unweighted
ensemble (since that is what we used to derive the relationship). If we add a weighting function,
then it would seem that we need to take this into account. However, the advantage of this ap-
proach is that we do not need to take this into account if we sum the acceptance probabilities in
the unweighted ensemble in our transition matrix, and not in the weighted ensemble. That is,
we compute two acceptance probabilities.

For particle additions:

acc
𝑉 ′
𝑃12 = min [1, 𝑒 −𝛽Δ𝑈+𝛽𝜇 +Δ𝜂 ] → used in acceptance of moves
𝑁+1

acc
𝑉 ′
𝑃12 = min [1, 𝑒 −𝛽Δ𝑈+𝛽𝜇 ] → used in updating sums in transition probabilities
𝑁+1

Notice that we don’t include the weight factor in the update of the transition probability matrix.
This enables us to compute transition probabilities as if we were still in the unweighted ensem-
ble, even though we are actually performing a simulation using a weighting function. A similar
case exists for particle deletions.

Other parameters
In all of the above examples, we presented flat histogram calculations in the context of a grand
canonical simulation. There, we computed a weighting function 𝜂(𝑁) that had a relation with
the 𝑁-dependence of the underlying Helmholtz free energy 𝐴(𝑁).

All of these methods presented can be used to compute any arbitrary free energy or entropy
function, including potentials of mean force like the one presented for umbrella sampling. In
general,

To compute the free energy along a given reaction coordinate, we need to bias
the simulation to perform flat-histogram sampling along that coordinate.

The table at the end of this section summarizes the microstate probabilities and reweighting ex-
pression that one might use in a number of different flat-histogram sampling. Keep in mind that
the microstate probabilities are used to determine acceptance criteria in Monte Carlo simula-
tions.

Keep in mind that continuous coordinates, like the energy or a distance-based reaction coordi-
nate, require us to discretize our weighting function, histograms, and reweighting procedure.

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

Stratification
When a flat histogram simulation is performed, the system experiences large fluctuations in the
flat histogram variable. As such, the correlation time for that variable can be quite large. In other
words, it can take the system a very long period of time to explore the complete range of interest
in the flat histogram. In our example, this means that the system has a long time scale for trav-
ersing values of 𝑁 between 𝑁min and 𝑁max . The time it takes for the system to perform a com-
plete walk of particle numbers between the two limits is called the tunneling time.

The tunneling time in flat histogram simulations can grow to be very long for a wide range 𝑁min
to 𝑁max . If the system performs a random walk in 𝑁, we expect

𝜏tunnel ∝ (𝑁max − 𝑁min )2

In reality, the walk is not completely random due to correlations, even if the system has a per-
fectly flat distribution in 𝑁. Typically the tunneling time grows with an exponent greater than 2.

One way around this problem is to use stratification to determine the underlying free energy or
entropy function that governs a flat histogram variable. In stratification, we break the entire
range of 𝑁 into smaller, overlapping sub-ranges or windows. Then we perform a separate sim-
ulation in each. At the end of the run, we know that the underlying free energy estimates should
be the same to within an unknown additive constant. We can shift these constants so as to obtain
overlap.

Consider that we perform multiple grand canonical simulations for a number of windows each
with different 𝑁min , 𝑁max . The simulations would reject moves that took their particle number
outside of the specified sub-range. From each, we would then obtain a high-quality estimate of
𝐴(𝑁):

𝐴(𝑁) = 𝜇𝑁 + 𝑘𝐵 𝑇𝜂(𝑁) + const

Then, we patch the 𝐴(𝑁) together by shifting to form a master curve:

𝜇𝑁 + 𝑘𝐵 𝑇𝜂(𝑁) 𝐴(𝑁)
+ const

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

However, the flat histogram approach makes this shifting much easier than the earlier un-
weighted case in which we used Ferrenberg-Swendsen reweighting. Here, because each value
of 𝑁 is sampled with the same frequency, we expect the same statistical error in our estimates
for 𝐴(𝑁) at each value of 𝑁. This means that we can use a simple least-squares procedure to
optimally shift each curve so as to obtain overlap.

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

Common flat-histogram ensembles
initial ensemble speci- flat his- microstate probabilities ideal weighting function
fied togram
variables
1 canonical 𝑇1 𝑈 ℘w (𝐫 𝑁 ) ∝ 𝑒 −𝛽1𝑈+𝜂(𝑈) 𝜂(𝑈) = 𝛽1 𝑈 − 𝑆(𝑈)

2 grand canonical 𝑇1 , 𝜇1 𝑈, 𝑁 ℘w (𝐫 𝑁 , 𝑁) ∝ 𝑒 −𝛽1𝑈+𝛽1 𝜇1𝑁+𝜂(𝑈,𝑁) 𝜂(𝑈, 𝑁) = 𝛽1 𝑈 − 𝛽1 𝜇1 𝑁 − 𝑆(𝑈, 𝑁)

3 isothermal-isobaric 𝑇1 , 𝑃1 𝑈, 𝑉 ℘w (𝐫 𝑁 , 𝑉) ∝ 𝑒 −𝛽1𝑈−𝛽1 𝑃1 𝑉+𝜂(𝑈,𝑉) 𝜂(𝑈, 𝑉) = 𝛽1 𝑈 + 𝛽1 𝑃1 𝑉 − 𝑆(𝑈, 𝑉)

4 grand canonical 𝑇1 , 𝜇1 𝑁 ℘w (𝐫 𝑁 , 𝑁) ∝ 𝑒 −𝛽1 𝑈+𝛽1 𝜇1 𝑁+𝜂(𝑁) 𝜂(𝑁) = −𝛽1𝜇1 𝑁 + 𝛽1 𝐴(𝑁; 𝑇1 )

5 isothermal-isobaric 𝑇1 , 𝑃1 𝑉 ℘w (𝐫 𝑁 , 𝑉) ∝ 𝑒 −𝛽1 𝑈−𝛽1𝑃1 𝑉+𝜂(𝑉) 𝜂(𝑉) = 𝛽1 𝑃1 𝑉 + 𝛽𝐴(𝑉; 𝑇1 )

canonical
6 (reaction coordi- 𝑇1 𝜉 ℘w (𝐫 𝑁 ) ∝ 𝑒 −𝛽1 𝑈+𝜂(𝜉) 𝜂(𝜉) = 𝛽𝐹(𝜉; 𝑇1 )
nate)

initial ensemble reweightable con- reweighting expression

ditions
1 canonical 𝑇2 ℘(𝑈; 𝑇2 ) ∝ ℘w (𝑈)𝑒 −(𝛽2−𝛽1 )𝑈−𝜂(𝑈)

2 grand canonical 𝑇2 , 𝜇2 ℘(𝑈, 𝑁; 𝑇2 , 𝜇2 ) ∝ ℘w (𝑈, 𝑁)𝑒 −(𝛽2 −𝛽1 )𝑈+(𝛽2 𝜇2−𝛽1 𝜇1 )𝑁−𝜂(𝑈,𝑁)

3 isothermal-isobaric 𝑇2 , 𝑃2 ℘(𝑈, 𝑉; 𝑇2 , 𝑃2 ) ∝ ℘w (𝑈, 𝑉)𝑒 −(𝛽2 −𝛽1)𝑈−(𝛽2 𝑃2−𝛽1 𝑃1)𝑉−𝜂(𝑈,𝑉)

4 grand canonical 𝜇2 (𝑇2 = 𝑇1 ) ℘(𝑁; 𝜇2 ) ∝ ℘w (𝑁)𝑒 (𝛽2 𝜇2−𝛽1 𝜇1)𝑁−𝜂(𝑁)

5 isothermal-isobaric 𝑃2 (𝑇2 = 𝑇1 ) ℘(𝑉; 𝑃2 ) ∝ ℘w (𝑉)𝑒 −(𝛽2 𝑃2 −𝛽1 𝑃1 )𝑉−𝜂(𝑉)

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

 canonical
6 (reaction coordi- none N/A
nate)
* 𝑆 denotes the dimensionless configurational entropy; 𝐴 denotes the configurational free energy
** All chemical potentials are relative by a factor of −3𝑘𝐵 𝑇 ln Λ(𝑇)

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