Umbrella Sampling: Monte Carlo Free-Energy Estimation

A Methodological Shift in Monte Carlo Simulations

This is a Method paper that introduces a novel computational technique for Monte Carlo simulations. It presents Umbrella Sampling, an importance sampling approach that uses non-physical distributions to calculate free energy differences in molecular systems.

The Sampling Gap in Phase Transitions

The paper addresses the failure of conventional Boltzmann-weighted Monte Carlo to estimate free energy differences.

The Problem: Free energy depends on the integral of configurations that are rare in the reference system. In a standard simulation, the relevant probability density $f_0(\Delta U^*)$ is too small to be sampled accurately by conventional Boltzmann-weighted Monte Carlo.
Phase Transitions: Conventional “thermodynamic integration” fails near phase transitions because it requires a path of integration where ensemble averages can be reliably measured, which is difficult in unstable regions.

Bridging States with Non-Physical Distributions

The authors introduce a non-physical distribution $\pi(q^N)$ to bridge the gap between a reference system (0) and a system of interest (1).

Arbitrary Weights: They generate a Markov chain with a limiting distribution $\pi(q^N)$ that differs from the Boltzmann distribution of either system. This distribution is written as $\pi(q’^N) = w(q’^N) \exp(-U_0(q’^N)/kT_0) / Z$, where $w(q^N) = W(\Delta U^*)$ is a weighting function chosen to favor configurations with values of $\Delta U^*$ important to the free-energy integral.
Reweighting Formula: The unbiased average of any property $\theta$ is recovered via the ratio of biased averages:

$$\langle\theta\rangle_{0}=\frac{\langle\theta/w\rangle_{w}}{\langle1/w\rangle_{w}}$$

Overlap: The method allows sampling a range of $\Delta U^*$ up to three times that of a conventional Monte Carlo experiment, enabling accurate determination of values of $f_0(\Delta U^*)$ as small as $10^{-8}$. If a single weight function cannot span the entire gap, additional overlapping umbrella-sampling experiments are carried out with different weighting functions exploring successively overlapping ranges of $\Delta U^*$.

Validation on Lennard-Jones Fluids

The authors validated Umbrella Sampling using Monte Carlo simulations of model fluids.

Experimental Setup

System Specifications: The study used a Lennard-Jones (LJ) fluid and an inverse-12 “soft-sphere” fluid.
System Size: Simulations were primarily performed with $N=32$ particles, with some validation runs at $N=108$ particles to check for size dependence.
State Points: Calculations covered a wide range of densities ($N\sigma^3/V = 0.50$ to $0.85$) and temperatures ($kT/\epsilon = 0.7$ to $2.8$), including the gas-liquid coexistence region.

Baselines

Baselines: Results were compared to thermodynamic integration data from Hansen, Levesque, and Verlet.
Quantitative Success:
- Agreement: The free energy estimates agreed with pressure integration results to within statistical uncertainties (e.g., at $kT/\epsilon=1.35$, Umbrella Sampling gave -3.236 vs. Conventional -3.25).
- Precision: Free energy differences were obtained with high precision ($\pm 0.005 NkT$ for $N=108$).
- Efficiency: A single umbrella run could replace the “numerous runs” required for conventional $1/T$ integrations.

Temperature Scaling via Reweighting

When the reference system has the same internal energy function as the system of interest (i.e., the same fluid at a different temperature), the free-energy expression simplifies to:

$$\frac{A(T)}{kT} = \frac{A(T_0)}{kT_0} - \ln \int f_0(U) \exp\left[-U\left(\frac{1}{kT} - \frac{1}{kT_0}\right)\right] dU$$

This is especially useful because a single determination of $f_0(U)$ over a wide energy range gives the free energy over a whole range of temperatures simultaneously. For 32 Lennard-Jones particles, only two umbrella-sampling experiments are needed to span the temperature range from the triple point ($kT/\epsilon = 0.7$) to twice the critical temperature ($kT/\epsilon = 2.8$). For 108 particles, four experiments suffice.

Mapping the Liquid-Gas Free Energy Surface

Methodological Utility: The method successfully mapped the free energy of the LJ fluid across the liquid-gas transition, a region where conventional methods face convergence problems.
N-Dependence: Comparison between $N=32$ and $N=108$ showed no statistically significant size dependence for free energy differences, suggesting small systems are sufficient for these estimates.
Comparison with Gosling-Singer Method: The paper contrasts its results with free energies derived from Gosling and Singer’s entropy estimation technique, finding discrepancies as large as $0.4N\epsilon$ (a 20% error in the nonideal entropy), equivalent to overestimating the configurational integral of a 108-particle system by a factor of $10^{16}$.
Generality: While demonstrated on energy ($U$), the authors note the weighting function $w$ can be any function of the coordinates, generalizing the technique beyond simple free energy differences.

Reproducibility

This 1977 paper predates modern code-sharing practices, and no source code or data files are publicly available. However, the paper provides sufficient algorithmic detail for reimplementation:

Constructing $W$: The paper does not derive $W$ analytically. It uses a trial-and-error procedure: start with a short Boltzmann-weighted experiment, then broaden the distribution in stages through short test runs, adjusting weights to flatten the probability density $f_w(\Delta U^*)$. The paper acknowledges this requires “interaction between the trial computer results and human judgment.”
Specific Weights: Table I provides the exact numerical weights used for the 32-particle soft-sphere experiment at $N\sigma^3/V = 0.85$, $kT/\epsilon = 2.74$, with values spanning from $W=1{,}500{,}000$ at the lowest energies down to $W=1.0$ at the center and back up to $W=16.0$ at the highest energies.
Potentials: The Lennard-Jones and inverse-twelve potentials are fully specified (Eqs. 8 and 9).
State Points: Densities and temperatures are enumerated in Tables II and III.
Block Averaging: Errors were estimated by treating sequences of $m$ steps as independent samples, where $m$ is determined by increasing block size until no systematic trends can be detected in either the average or the standard deviation of the mean.

Paper Information

Citation: Torrie, G. M., & Valleau, J. P. (1977). Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. Journal of Computational Physics, 23(2), 187-199. https://doi.org/10.1016/0021-9991(77)90121-8

Publication: Journal of Computational Physics, 1977

@article{torrie1977nonphysical,
  title={Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling},
  author={Torrie, Glenn M and Valleau, John P},
  journal={Journal of Computational Physics},
  volume={23},
  number={2},
  pages={187--199},
  year={1977},
  publisher={Elsevier},
  doi={10.1016/0021-9991(77)90121-8}
}

A Methodological Shift in Monte Carlo Simulations#

The Sampling Gap in Phase Transitions#

Bridging States with Non-Physical Distributions#

Validation on Lennard-Jones Fluids#

Experimental Setup#

Baselines#

Temperature Scaling via Reweighting#

Mapping the Liquid-Gas Free Energy Surface#

Reproducibility#

Paper Information#