Monte Carlo in Statistical Mechanics: The Metropolis Algorithm

1. The Sampling Problem

Statistical mechanics says that a system in thermal equilibrium at temperature T visits configuration (microstate) i with probability given by the Boltzmann distribution:

For anything beyond a toy system, Z is an astronomically large sum — a lattice of just 20×20 Ising spins already has 2⁴⁰⁰ ≈ 10¹²⁰ microstates, more than there are atoms in the observable universe. Computing thermodynamic averages by literally summing over every microstate is therefore impossible. Ordinary Monte Carlo — sampling configurations uniformly at random and weighting by their Boltzmann factor — is also hopeless: at low temperature, overwhelmingly the random samples land on high-energy states with essentially zero weight, wasting almost all computational effort.

2. The Metropolis Algorithm

The Metropolis algorithm builds a Markov chain — a sequence of configurations where each one depends only on the previous one — whose long-run visiting frequency converges to the Boltzmann distribution. Each step:

The acceptance rule always accepts moves that lower the energy, and accepts energy-raising moves only occasionally, with probability falling exponentially with ΔE — exactly the balance needed to spend most computational effort near the thermodynamically important, low-energy configurations while still occasionally exploring higher-energy states (which is what allows the system to escape local minima and correctly sample thermal fluctuations).

3. Detailed Balance and Why It Works

The algorithm's correctness rests on detailed balance: at equilibrium, the flow of probability from state i to state j must exactly equal the reverse flow from j to i:

Splitting the transition probability into a proposal probability g(i→j) (chosen to be symmetric, g(i→j) = g(j→i), for a simple proposal such as flipping a random spin) and an acceptance probability A(i→j), detailed balance requires:

The Metropolis choice A(i→j) = min[1, exp(−ΔE/k_BT)] is one valid solution to this ratio (among several possible acceptance rules — the Glauber/heat-bath rule A = 1/(1+exp(ΔE/k_BT)) is another). Metropolis's choice is preferred in practice because it maximises the acceptance rate subject to satisfying detailed balance, making the Markov chain explore configuration space as efficiently as possible.

4. Practical Considerations

• Equilibration (burn-in): the chain starts from an arbitrary configuration and needs a number of steps to "forget" the initial state and start sampling the true equilibrium distribution. These early steps are discarded before computing averages.
• Autocorrelation time: successive configurations in the chain are correlated (each differs from the last by only one small move). Samples must be spaced by roughly the autocorrelation time τ before being treated as statistically independent for error estimation.
• Critical slowing down: near a phase transition (e.g. the Curie point of the Ising model), τ diverges, requiring vastly more steps to sample correctly. Cluster algorithms (Wolff, Swendsen-Wang) were developed specifically to overcome this by flipping entire correlated clusters of spins in a single move.
• Acceptance ratio tuning: for continuous variables (e.g. particle displacement), the proposal step size is tuned so roughly 30-50% of moves are accepted — too small a step wastes time with tiny moves; too large a step is almost always rejected.

5. Worked Example: The Ising Model

The 2D Ising model — a grid of spins sᵢ = ±1 with energy E = −J·Σ⟨ij⟩ sᵢsⱼ over nearest-neighbour pairs — is the canonical Metropolis test case. Flipping a single spin sₖ changes the energy by a purely local amount:

Because ΔE depends only on the four (or six, or eight, depending on lattice geometry) neighbours of the flipped spin, evaluating the acceptance criterion costs O(1) regardless of system size — this locality is what makes Metropolis Monte Carlo for lattice models like Ising extremely fast per step, even though many millions of steps may be needed to reach equilibrium near the critical temperature.

6. Beyond Physics

Computational Chemistry and Molecular Simulation

Monte Carlo methods complement molecular dynamics (see our companion article on the Lennard-Jones potential) whenever the quantity of interest is a thermodynamic average rather than a time-dependent property. Metropolis moves — displacing a molecule, rotating it, or attempting to insert/delete one (grand-canonical Monte Carlo) — are the workhorse of simulating phase equilibria, adsorption, and solvation free energies.

Bayesian Statistics and Machine Learning

The Metropolis-Hastings generalisation (allowing asymmetric proposals) underlies Markov chain Monte Carlo (MCMC) methods used throughout Bayesian statistics to sample posterior distributions that have no closed form — from cosmological parameter estimation to Bayesian neural network inference.

Combinatorial Optimisation

Simulated annealing — gradually lowering an artificial "temperature" parameter in a Metropolis loop applied to a cost function rather than a physical energy — is a general-purpose optimisation heuristic for problems like the travelling salesman problem, circuit layout, and protein structure prediction, where the Metropolis acceptance rule allows uphill moves early on to escape poor local optima.