Equation of State Calculations by Fast Computing Machines

A gas of N interacting particles at temperature T has thermodynamic properties determined by the partition function — a sum over every possible arrangement of all N particles. Each arrangement has an energy E, and its contribution is weighted by the Boltzmann factor $\exp(-E / kT)$.

For a system of just 100 particles, the number of possible configurations is effectively infinite. Direct computation is hopeless. By 1953 physicists had been using Monte Carlo integration: pick random configurations, compute their energies, and average. But uniform random sampling wastes nearly all its effort on high-energy configurations whose Boltzmann weight is essentially zero. It's like searching for needles in a haystack by examining random straws.

The key idea is importance sampling with a random walk. Instead of choosing configurations uniformly and then weighting them, generate configurations that are already distributed according to the Boltzmann distribution. Then the average of any property A is just a simple unweighted mean of the values you see along the walk.

But how do you generate samples from the Boltzmann distribution when you can't even compute the partition function? This is the chicken-and-egg problem the Metropolis algorithm solves: you don't need to know the normalizing constant — only ratios of probabilities, which are ratios of Boltzmann factors, and those simplify to $\exp(-\Delta E / kT)$.

The Metropolis algorithm is strikingly simple. Start from any configuration. Then repeat:

1. Propose — pick a particle at random and move it by a small random displacement.
2. Evaluate — compute the energy change ΔE caused by the move.
3. Decide — if ΔE ≤ 0 (energy went down or stayed the same), accept the move. If ΔE > 0 (energy went up), accept with probability $\exp(-\Delta E / kT)$ — draw a uniform random number u ∈ [0, 1] and accept if $u < \exp(-\Delta E / kT)$.

That's the entire algorithm. The accepted configurations form a Markov chain whose stationary distribution is the Boltzmann distribution. After a burn-in period, every configuration along the chain is a sample from the target distribution.

Think of the acceptance rule as a doorman at a club. If you're bringing the energy down (the party gets better), you always get in. If you'd raise the energy (make it worse), you might get in — the doorman flips a weighted coin, and the heavier the energy cost the less likely you pass. At high temperature the coin is nearly fair (the doorman is lenient), so the walk explores freely. At low temperature the coin is heavily loaded against uphill moves, so the walk hugs the valleys.

The magic lies in a property called detailed balance. Imagine traffic flowing between two cities. If in the long run the same number of cars flow from A to B as from B to A, neither city gains or loses population — the system is in equilibrium.

Detailed balance says exactly this for the Markov chain: for any two states i and j, the flow of probability from i to j equals the flow from j to i. Formally:

$\pi(i) \, T(i \to j) = \pi(j) \, T(j \to i)$

where π is the target distribution and T is the transition probability. When this holds, π is the stationary distribution of the chain. The Metropolis acceptance rule is designed to satisfy this equation when π is the Boltzmann distribution. You can verify: the ratio $T(i \to j) / T(j \to i)$ equals $\pi(j)/\pi(i)$ by construction.

Temperature T plays a crucial role in balancing exploration and exploitation:

• At high temperature, $\exp(-\Delta E / kT) \approx 1$ for most moves, so nearly everything is accepted. The walk explores broadly but doesn't concentrate on any particular region — like a heated gas where particles fly everywhere.
• At low temperature, $\exp(-\Delta E / kT) \approx 0$ for uphill moves, so only downhill moves are accepted. The walk becomes trapped in the nearest energy minimum — like a frozen solid where particles barely vibrate.
• At just right temperature, the walk explores enough to find deep valleys but lingers in them long enough to sample them well.

This temperature knob inspired simulated annealing (Kirkpatrick et al., 1983): start at high temperature to explore, then slowly cool to concentrate on the global minimum — mimicking the physical annealing of metals.

The original Metropolis algorithm assumed the proposal distribution is symmetric: the probability of proposing a move from state i to state j equals the probability of proposing the reverse. In 1970, W. K. Hastings generalized the acceptance probability to handle asymmetric proposals:

$\alpha(i \to j) = \min\!\left(1,\; \frac{\pi(j)\, q(i \mid j)}{\pi(i)\, q(j \mid i)} \right)$

where $q(j \mid i)$ is the proposal probability. When the proposal is symmetric — $q(j \mid i) = q(i \mid j)$ — the q terms cancel and you recover the original Metropolis rule. This generalization opened MCMC to a vast range of proposal strategies and made the algorithm applicable far beyond physics.

The Markov chain will eventually converge to the target distribution — this is guaranteed by ergodicity (the chain can reach any state from any other state) and detailed balance. But "eventually" can be a long time.

Two practical concerns dominate:

• Burn-in — the initial samples reflect the starting point, not the target distribution. Discard the first N samples (the "burn-in" period) to remove this bias.
• Mixing — if the chain moves in tiny steps, consecutive samples are highly correlated and the chain takes a long time to "mix" (explore the full distribution). If steps are too large, most proposals are rejected and the chain barely moves. The proposal scale must be tuned: a common rule of thumb targets an acceptance rate of about 23–44% in high dimensions.

Metropolis and the Rosenbluths implemented the algorithm on the MANIAC computer at Los Alamos for a two-dimensional system of 224 rigid disks (hard spheres that can't overlap). Each "move" displaced one disk by a small random amount within a square of side 2α. If the new position caused disks to overlap, the move was rejected. Otherwise, it was always accepted (since all non-overlapping configurations have the same energy for hard spheres).

They ran the chain for hundreds of cycles (one cycle = 224 attempted moves, one per particle) and computed the equation of state — the relationship between pressure and density. Their results agreed with the theoretical four-term virial expansion at low density and revealed behavior at high density that no analytical method could predict.

The Metropolis algorithm is to sampling what gradient descent is to optimization: a simple, universal workhorse that launched an entire field. Every modern MCMC variant — Hamiltonian Monte Carlo, slice sampling, reversible-jump MCMC — descends from the same core idea: propose, evaluate the ratio, accept or reject.