Physics and statistical mechanics

Why start with a ball on a hill

A molecular dynamics simulation of a protein or a crystal follows millions of atoms, but each atom obeys the same laws as a single ball rolling in a valley. If you understand one particle moving in one dimension, its position, velocity, momentum, force and energy, and the picture called phase space, you understand what every MD code computes, and why its results can be trusted. This chapter builds that picture.

Describing motion

The position $x(t)$ says where the particle is at time $t$. Its rate of change is the velocity, and the rate of change of velocity is the acceleration:

Newton recognised that the quantity a force changes is not velocity alone but velocity weighted by how much matter is moving, the momentum:

The second law is a differential equation: given the position and velocity now, and the force as a function of position, it determines the entire future and past. This determinism is what lets a computer predict motion step by step.

Forces from a landscape

For atoms, forces are not arbitrary: they come from a potential energy $V(x)$, a landscape that depends on where the particles are (the chapter on potential energy surfaces shows where atomic landscapes come from). The force is the downhill slope:

Three landscapes appear in the lab. A spring, $V = \tfrac12 kx^2$, gives a force $F = -kx$ that always pulls back towards the centre (Hooke’s law), the model of a chemical bond near equilibrium. A double well has two valleys separated by a barrier, the simplest model of a chemical reaction, a conformational change or a phase transition. A pendulum, $V = 1 - \cos x$, is a model of a rotating chemical group with a periodic barrier.

Energy, and why it is conserved

Watch the energy bar while the particle moves. As it rolls downhill, potential energy turns into kinetic energy $K = \tfrac12 mv^2 = p^2/2m$; as it climbs, kinetic turns back into potential; and without friction their sum never changes.

Conservation gives a lot for free. The particle can only be where $V(x) \le E$; the points where $V(x) = E$ are turning points, where it momentarily stops and reverses. With $E$ below the barrier of the double well, it is trapped in one valley forever; above it, it crosses back and forth. In molecular dynamics, the fact that a correct integrator keeps the total energy constant is the first and most important test of a simulation (the MD chapter).

Add friction in the lab: energy now leaks away, and the particle settles into a minimum. That is what an energy minimisation does to a structure, and what a thermostat has to counteract when it adds the right amount of random kicks back.

Phase space

Here is the central picture of this volume. The state of the particle at any instant is not its position alone, because two particles at the same place moving in opposite directions have different futures. The complete state is the pair $(x, p)$, a point in a plane called phase space. Newton’s law moves that point along a curve, the trajectory, and because the state determines the future uniquely, trajectories never cross.

Since energy is conserved, the trajectory stays on a curve of constant $E$: the energy contour $p = \pm\sqrt{2m(E - V(x))}$. In the lab:

• the spring gives ellipses: the particle oscillates forever around the centre;
• the double well shows two families of loops, one inside each valley, enclosed by a figure-of-eight called the separatrix; outside it, large loops circle both valleys;
• the pendulum shows closed loops (swinging) inside wavy open curves (rotating all the way round);
• friction turns closed loops into spirals that end at a minimum.

For $N$ atoms in three dimensions phase space has $6N$ dimensions, three positions and three momenta per atom, and a molecular dynamics run traces one trajectory through it. Statistical mechanics, the subject of the chapters after next, is the art of saying what fraction of that enormous space a system visits.

Hamilton’s formulation

Writing the energy as a function of position and momentum, $H(x, p) = p^2/2m + V(x)$, called the Hamiltonian, puts Newton’s law in a beautifully symmetric form:

Hamiltonian flow has a deep property, Liouville’s theorem: a small blob of initial conditions changes shape as it flows but never changes its area (volume, in many dimensions). This is why the integrators used in molecular dynamics are built to be symplectic, preserving phase-space volume exactly, and it is the foundation on which statistical mechanics builds its probabilities.

References and further reading

Why a simulation chapter needs quantum mechanics

Molecular dynamics moves nuclei with Newton’s laws, as in the previous chapter. Yet the forces between atoms come from electrons, which obey quantum mechanics, and even the nuclei are quantum particles. This chapter answers two questions every simulator should be able to answer. Where does quantisation come from? And when is it safe to forget about it?

Waves of matter

In 1924 de Broglie proposed that every particle with momentum $p$ behaves like a wave of wavelength $\lambda = h/p$, where $h = 6.626 \times 10^{-34}$ J s is Planck’s constant. For a cricket ball the wavelength is absurdly small; for an electron in an atom it is about the size of the atom, and wave behaviour dominates. The state of a quantum particle is described by a wavefunction $\psi(x)$, and the probability of finding the particle near $x$ is $|\psi(x)|^2\,\mathrm dx$ (Born’s rule).

States of definite energy satisfy the time-independent Schrödinger equation, in which the classical energy $p^2/2m + V$ becomes an operator:

The equation itself is not mysterious: it says that a wave curves more sharply where the kinetic energy $E - V$ is larger. What creates quantisation is the boundary condition. A confined wave has to fit inside its container, like the standing wave on a guitar string, and only a discrete set of wavelengths do.

The particle in a box

Trap a particle between two impenetrable walls, $V = 0$ for $0 < x < L$ and $V = \infty$ outside. The wavefunction must vanish at both walls, so it must be a sine wave with a whole number of half-wavelengths in the box:

Three consequences matter for materials. Energy is quantised. The lowest energy is not zero: squeezing a particle into a small space forces short wavelengths and therefore kinetic energy, the reason electrons do not collapse into the nucleus. And the spacing falls as $1/L^2$, so in a big box the levels are so close that the energy looks continuous. Quantum dots exploit the $1/L^2$ law: the smaller the nanocrystal, the larger its gap and the bluer its light.

The harmonic oscillator

Near the bottom of any valley the potential is a parabola, $V = \tfrac12 m\omega^2 x^2$, so the harmonic oscillator describes every chemical bond vibration and every phonon in a crystal. Its levels are equally spaced:

Zero-point energy is measurable and important. It is why liquid helium does not freeze at normal pressure, it shifts reaction barriers and bond lengths, and it makes deuterated molecules react more slowly than ordinary ones (the kinetic isotope effect), because the heavier isotope has a lower $\tfrac12\hbar\omega$.

From quantum to classical

Highlight a high level in the lab and switch to $|\psi|^2$. The quantum probability oscillates rapidly, but its average follows the teal dashed curve, the classical probability density: a classical oscillator moves fastest through the centre and slowest near the turning points, so it spends most of its time near the edges,

At high quantum number the two agree, as Bohr’s correspondence principle demands. Quantum effects matter when the level spacing is comparable to the energy available, and that brings temperature in.

Temperature populates the levels

At temperature $T$ a system in contact with its surroundings occupies level $n$ with a probability given by the Boltzmann distribution, derived in the chapter on Boltzmann and Maxwell:

Slide the temperature in the lab. When $k_BT \ll \hbar\omega$ almost everything sits in the ground state, and the oscillator is frozen: it cannot absorb heat, and its thermal energy is pinned at the zero-point value. When $k_BT \gg \hbar\omega$ many levels are populated, and the mean energy approaches the classical value $k_BT$ (for the box, $k_BT/2$), as the statistics row shows.

The numbers for real molecules are revealing. A C–H stretch vibrates at about 3000 cm⁻¹, so $\hbar\omega/k_BT \approx 14$ at room temperature and the ground state holds 99.9999 % of the molecules: that vibration is completely quantum. A soft torsion or lattice mode at 200 cm⁻¹ has $\hbar\omega/k_BT \approx 1$, and is close to classical. Einstein used exactly this freezing-out to explain in 1907 why the heat capacity of solids falls at low temperature.

Is classical molecular dynamics allowed?

Classical MD treats nuclei as points obeying Newton’s laws, so it gives each vibration $k_BT$ of energy whether or not it should. A useful yardstick is the thermal de Broglie wavelength, the typical quantum blur of a particle at temperature $T$:

The verdict is that classical dynamics of nuclei is accurate for heavy atoms and for structure, diffusion and most thermodynamics at room temperature. It fails for light atoms (hydrogen, lithium), for high-frequency vibrations, for low temperatures and for heat capacities. The remedies are to constrain stiff bonds (as most biomolecular force fields do with bonds to hydrogen), to apply quantum corrections afterwards, or to use path-integral MD, which represents each quantum nucleus as a ring of classical beads.

References and further reading

The question statistical mechanics answers

A glass of water contains about $10^{25}$ molecules. Following every one is hopeless, and also unnecessary: we want the temperature, the pressure, the heat capacity, whether ice melts. Statistical mechanics connects the microscopic world of quantum levels and phase-space trajectories to those few macroscopic numbers. The whole subject rests on a single idea, which this chapter makes concrete by counting.

Microstates and macrostates

A microstate is a complete specification of the system: which quantum level every particle occupies, or every position and momentum. A macrostate is what we measure: total energy, volume, number of particles. One macrostate corresponds to an enormous number of microstates, written $\Omega$ (the multiplicity).

The single assumption of statistical mechanics, the fundamental postulate, is that an isolated system in equilibrium is equally likely to be found in any of its accessible microstates. Everything else, including the second law of thermodynamics, follows by counting.

Counting energy in a solid

Einstein’s model of a solid treats each atom’s vibration as quantum harmonic oscillators with the same frequency, so energy comes in identical quanta $\varepsilon = \hbar\omega$. A microstate says how many quanta each oscillator holds. With $N$ oscillators and $q$ quanta, the number of microstates is the number of ways to drop $q$ identical balls into $N$ boxes:

The ten microstates of the example are easy to list: (3,0,0) in three arrangements, (2,1,0) in six, and (1,1,1) in one. Multiplicities grow at a ferocious rate: 100 oscillators sharing 100 quanta already have about $5 \times 10^{58}$ microstates.

Two solids in contact

Now put two solids, A and B, in thermal contact so that quanta can hop between them while the total $q$ stays fixed. A macrostate is now “how many quanta are in A”, and since every microstate of the whole is equally likely, the probability of a macrostate is proportional to its multiplicity,

Press the button and watch. Every hop moves one quantum from a random oscillator to another random oscillator, and nothing in the rule prefers either solid. Yet the energy flows out of A until A’s share fluctuates around the peak of the probability curve, and it never comes back. For the default settings (60 and 40 oscillators, 100 quanta) the peak is at $q_A = 60$ with a spread of about $\pm 7$, while the probability of finding all 100 quanta back in A is about $5 \times 10^{-15}$.

That is the second law of thermodynamics, and nothing more: energy flows to the macrostate with overwhelmingly more microstates, simply because there are overwhelmingly more ways to be there. Increase the sizes in the lab and the peak becomes narrower in relative terms, its width falling as $1/\sqrt N$. For $10^{23}$ oscillators the fluctuations are invisible and the macrostate behaves as if it were exact.

Entropy

Multiplicities multiply when systems are combined, $\Omega = \Omega_A\Omega_B$. It is more convenient to work with a quantity that adds, so Boltzmann defined the entropy as the logarithm of the multiplicity:

This equation is carved on Boltzmann’s tombstone in Vienna. Entropy is not a measure of “disorder” in any vague sense: it is the logarithm of a count.

Temperature

What decides where the peak lies? At the maximum of $\ln\Omega_A + \ln\Omega_B$, moving one quantum from B to A leaves the total unchanged, so

The two solids are in equilibrium when this quantity is the same for both. That shared quantity is what a thermometer measures, and it defines temperature:

The lab prints the temperature of each solid, computed from the exact counts. Start with all energy in A: A is hot and B is cold. As quanta flow, the two temperatures move towards each other and meet where the probability peaks. Energy flows from hot to cold because the cold body gains more entropy per quantum than the hot body loses; the total entropy rises until it can rise no further.

From counting to thermodynamics

Everything in classical thermodynamics follows from $S(U, V, N)$. Differentiating with respect to volume and particle number gives pressure and chemical potential,

which is the thermodynamic identity $\mathrm dU = T\,\mathrm dS - P\,\mathrm dV + \mu\,\mathrm dN$. This isolated-system picture, fixed $N$, $V$ and $E$, is the microcanonical ensemble, and it is exactly what an MD simulation without a thermostat samples (the NVE ensemble). For classical particles, $\Omega$ becomes the volume of the shell of phase space with energy $E$, measured in units of $h^{3N}$.

References and further reading

From an isolated system to one at a temperature

The previous chapter treated an isolated system: fixed energy, every microstate equally likely. Real samples, and most simulations, are instead in contact with surroundings at a fixed temperature, which can lend and borrow energy. What is the probability of finding a small system in a particular microstate then? The answer, the Boltzmann distribution, is used in every corner of chemistry and materials science.

Deriving the Boltzmann factor

Let a small system share energy with a huge reservoir, with total energy $E_{\text{tot}}$ fixed. When the system is in a particular microstate $s$ with energy $E_s$, the reservoir has energy $E_{\text{tot}} - E_s$, and the combined system is still isolated, so every joint microstate is equally likely. The probability of $s$ is therefore proportional to the number of reservoir microstates compatible with it:

Because the reservoir is huge, $E_s$ is a tiny change for it, and a first-order Taylor expansion is exact in the limit: $S_R(E_{\text{tot}} - E_s) \approx S_R(E_{\text{tot}}) - E_s\,\partial S_R/\partial E = S_R(E_{\text{tot}}) - E_s/T$. The constant first term cancels on normalising, which leaves one of the great results of physics:

This is the canonical ensemble (fixed $N$, $V$, $T$), the ensemble a thermostatted MD run samples (NVT). Its consequences are everywhere: vibrational populations, the Arrhenius law $k \propto e^{-E_a/k_BT}$ for reaction rates, vacancy concentrations in crystals, the exponential thinning of the atmosphere with height.

The partition function knows everything

$Z$ looks like a mere normalising constant, but all of thermodynamics can be extracted from it by differentiation. With $\beta = 1/k_BT$:

For a classical particle the sum over states becomes an integral over phase space, $Z = \frac{1}{h^{3N}N!}\int e^{-H(\mathbf r, \mathbf p)/k_BT}\,\mathrm d\mathbf r\,\mathrm d\mathbf p$. Because $H = K(\mathbf p) + V(\mathbf r)$, the integral splits into a momentum part, which can always be done exactly, and a configurational part, which is the hard part that MD and Monte Carlo exist to sample.

The speeds of molecules

The momentum part gives a striking prediction. Each velocity component has a Boltzmann weight $e^{-mv_x^2/2k_BT}$, a Gaussian. Combining the components and counting how many velocity vectors have the same speed (a circle of circumference $2\pi v$ in 2D, a sphere of area $4\pi v^2$ in 3D) gives the Maxwell–Boltzmann distribution of speeds:

Maxwell derived this in 1860 from symmetry arguments alone, before anyone could measure molecular speeds. The lab lets you watch it happen.

Restart and look at the histogram. At first every particle has the same speed, a single spike. Within a few hundred collisions the spike spreads into exactly the Maxwell–Boltzmann shape, though nothing in the program knows that curve: collisions conserve energy and momentum, and the distribution is simply the one with the most microstates. This is the second law happening in front of you. Boltzmann’s H-theorem (1872) proved that collisions drive any initial distribution towards this one.

Equipartition

Now tick “two species”. The orange particles are four times heavier but start with the same speed, so each carries four times the kinetic energy of a light particle. Watch the readout of mean kinetic energy: within a short time the two species have the same average kinetic energy, and the heavy particles move half as fast. This is the equipartition theorem:

Equipartition is how MD codes measure temperature. They compute the instantaneous kinetic temperature $T = \frac{2\langle K \rangle}{N_{\text{df}}\,k_B}$, with $N_{\text{df}}$ the number of degrees of freedom (about $3N$, less any constraints), and a thermostat nudges velocities to keep it near the target. It also explains the classical-MD error discussed earlier: equipartition gives every vibration $k_BT$, even those that quantum mechanics would freeze.

Pressure from collisions

Pressure is the momentum delivered to the walls per unit time per unit area. A particle bouncing off a wall reverses its normal velocity and delivers momentum $2m|v_x|$. Averaging over the Maxwell–Boltzmann distribution gives the ideal gas law:

The lab measures pressure by adding up the momentum transferred to the walls. Compare it with the ideal-gas value: the measured pressure is somewhat higher, because the disks have a finite size and the free area is smaller than the box. That is the excluded-volume correction in van der Waals’s equation, and the starting point of the theory of real gases and liquids. In MD, pressure is computed from the virial, $PV = Nk_BT + \frac{1}{d}\left\langle \sum_i \mathbf r_i \cdot \mathbf F_i \right\rangle$ in $d$ dimensions, which adds the interparticle forces to the kinetic term (the MD chapter).

References and further reading

Ensembles

An experiment never controls the microscopic state of a sample. It controls a few macroscopic quantities, such as the temperature, the pressure or the number of molecules, and the sample explores all the microstates compatible with those constraints. Statistical mechanics describes this exploration with an ensemble: an imagined collection of infinitely many copies of the system, each in a different microstate, with the copies distributed according to the probability of each microstate under the given conditions. The average of any quantity over the ensemble is what an experiment measures.

The chapter on entropy used the microcanonical ensemble, which describes an isolated system with fixed number of particles $N$, volume $V$ and energy $E$, and in which every accessible microstate is equally likely. The chapter on the Boltzmann distribution derived the canonical ensemble for a system at fixed $N$, $V$ and temperature $T$, in contact with a heat bath, and showed that a microstate of energy $E_s$ occurs with probability proportional to $e^{-E_s/k_BT}$. This chapter develops the canonical ensemble fully, because it is the ensemble of most experiments and of most simulations, and then introduces the ensembles at constant pressure and at constant chemical potential.

Different ensembles describe different experimental conditions, but for large systems they give the same thermodynamics, because the fluctuations that distinguish them are tiny compared with the averages; this is called the equivalence of ensembles. The choice is therefore a matter of convenience: one picks the ensemble in which the calculation is easiest or the simulation most natural.

The partition function

In the canonical ensemble the probability of microstate $s$ is

The partition function $Z$ looks like a mere normalising constant, but it is the central object of statistical mechanics. Its name is apt: it describes how the system’s probability is partitioned among its states. At low temperature only the lowest states have appreciable weight and $Z$ is close to the degeneracy of the ground state. As the temperature rises, more states acquire weight and $Z$ grows; roughly, $Z$ counts the number of states the system can effectively reach at temperature $T$.

The figure below makes this concrete. Each energy level contributes a weight $g_i e^{-\beta E_i}$, drawn as a gold bar; the weights add up to $Z$, and each level’s share of the total is its probability. Lower the temperature and the weights of the excited levels collapse exponentially; raise it and they approach their degeneracies.

Partition functions factorise in a way that makes them practical. If a system consists of independent parts, its energy is a sum of their energies, so each Boltzmann factor is a product and $Z = z_1 z_2 \cdots$. A molecule’s partition function factorises approximately into translational, rotational, vibrational and electronic parts, $z = z_{\text{trans}}\,z_{\text{rot}}\,z_{\text{vib}}\,z_{\text{elec}}$, each computed separately; this is how thermodynamic tables of gases are calculated from spectroscopic data. For $N$ identical, non-interacting molecules, $Z = z^N/N!$, where the $N!$ corrects for the fact that exchanging identical molecules does not produce a new microstate.

Thermodynamics from the partition function

Every thermodynamic quantity follows from $\ln Z$ by differentiation. The derivations are short, and they explain what each quantity means.

The mean energy formula follows because differentiating $Z = \sum e^{-\beta E_s}$ with respect to $\beta$ brings down a factor $-E_s$ in each term, so $-\partial \ln Z/\partial\beta = \sum E_s e^{-\beta E_s}/Z = \langle E\rangle$. Differentiating again gives the variance of the energy, which is why the heat capacity equals the energy fluctuations divided by $k_BT^2$: a system can absorb heat only if it has states at nearby energies to move into, and the same states allow its energy to fluctuate. This fluctuation formula is how heat capacities are computed from simulations, and the chapter on fluctuations develops it.

The entropy formula $S = -k_B\sum p_s \ln p_s$, the Gibbs entropy, generalises Boltzmann’s $S = k_B \ln \Omega$. If all $\Omega$ accessible states are equally likely, each $p_s = 1/\Omega$ and the sum reduces to $k_B \ln\Omega$. If the probabilities are uneven, the entropy is lower, because the system is effectively confined to fewer states. The same expression, without $k_B$, is Shannon’s measure of information, which is why entropy is often described as missing information about the microstate.

The figure’s bottom panel shows these quantities for each system. For the two-level system, the heat capacity rises from zero, peaks and falls again. This Schottky anomaly occurs because at low temperature the excited level is out of reach, and at high temperature both levels are already equally populated, so further heating changes nothing; only in between can the system absorb energy. The peak lies at $k_BT \approx 0.42\,\Delta$ for a level spacing $\Delta$, and such peaks are measured in magnetic materials and used to determine their energy levels. For the harmonic oscillator, the heat capacity rises from zero to the classical value $k_B$, the freezing-out of vibrations described in the chapter on quantum energy levels.

Free energy: energy against entropy

The Helmholtz free energy $F = U - TS$ is the quantity a system at constant temperature minimises. Its two terms pull in opposite directions. The energy term favours the states of lowest energy: atoms packed in a crystal, a protein in its folded form, a molecule in its deepest potential well. The entropy term, weighted by temperature, favours whatever has the most microstates: a gas filling its container, a disordered chain, a broad, shallow well rather than a narrow, deep one. Which wins depends on the temperature, and that competition is the origin of melting, boiling, unfolding and most other changes of state.

The cleanest way to see this is a system with two states, A and B. The probability of being in state A is the sum of the Boltzmann weights of all microstates belonging to A, so it is natural to define a partition function for each state, $Z_A$ and $Z_B$, and a free energy for each, $F_A = -k_BT\ln Z_A$. The ratio of probabilities then depends only on the free-energy difference:

The figure below samples a particle moving in a one-dimensional potential with two wells by the Metropolis Monte Carlo method (its own chapter). Well A is deep and narrow; well B is shallower but much broader. At low temperature the particle stays almost entirely in A, the state of lower energy. Raise the temperature and B gains, until above some temperature the particle spends more time in the shallow well than in the deep one. Nothing about the energies has changed. B wins because it offers more room: its larger entropy, multiplied by the higher temperature, outweighs its energy penalty. The readout decomposes the exact free-energy difference into its energy and entropy parts.

This is not a toy effect. Rubber contracts when heated because a stretched polymer chain has fewer conformations than a coiled one, so the entropic force pulling it back grows with temperature. Water and oil separate largely because water molecules next to an oily surface lose orientational freedom, an entropic penalty; this hydrophobic effect drives protein folding and the self-assembly of surfactants and nanoparticle ligands. And every first-order phase transition, from ice melting to a crystal changing its structure, is the point where the free energies of two phases cross.

Free-energy profiles and the potential of mean force

The same idea extends from two states to a continuous variable. Choose a collective variable $\xi$ that describes the progress of some process: the distance between two molecules, a dihedral angle, the number of atoms in a crystalline cluster. Its probability distribution $P(\xi)$, obtained by summing the Boltzmann weights of all microstates with a given value of $\xi$, defines a free-energy profile:

In the two-well figure, the orange curve is exactly this: a histogram of where the particle has been, turned into a free energy. The difference between it and the blue potential-energy curve is entropy. Free-energy profiles are central to molecular simulation. The barrier height $\Delta F^\ddagger$ of a profile along a reaction coordinate controls the rate of the process through a factor $e^{-\beta\Delta F^\ddagger}$, and the depth of a minimum controls binding constants and solubilities. When the barrier is many $k_BT$ high, as in the figure at low temperature, an unbiased simulation almost never crosses it and the histogram cannot be measured on the far side; the readout warns when this happens. The chapter on free-energy methods describes umbrella sampling and metadynamics, which bias the simulation to overcome this.

The classical partition function and the ideal gas

For a classical system of $N$ atoms, the sum over microstates becomes an integral over phase space, the $3N$ positions and $3N$ momenta:

Because the kinetic energy is a sum of squares of momenta, the momentum integrals are products of Gaussian integrals and can always be done exactly. All the difficulty, and all the interesting physics of interacting matter, lies in the remaining configurational integral over the potential energy. It is a $3N$-dimensional integral that cannot be computed on a grid, which is why it is sampled by Monte Carlo or molecular dynamics.

For an ideal gas, $U = 0$ and the configurational integral is simply $V^N$, so $Z = V^N/(N!\,\Lambda^{3N})$. Everything follows from this one line. The pressure is $P = k_BT\,\partial\ln Z/\partial V = Nk_BT/V$, the ideal gas law. The entropy, using Stirling’s approximation $\ln N! \approx N\ln N - N$, is the Sackur–Tetrode equation,

For argon at 298.15 K and 1 bar, the thermal wavelength is 0.16 Å, $V/(N\Lambda^3) \approx 10^7$, and the formula gives a molar entropy of 154.8 J mol⁻¹ K⁻¹, which matches the value measured calorimetrically. That an absolute entropy can be calculated from Planck’s constant and the mass of an atom, and agree with experiment, was one of the early triumphs of quantum statistical mechanics: the factor $h^{3N}$ and the $N!$ both matter.

Other ensembles

Experiments are more often done at constant pressure than at constant volume. In the isothermal–isobaric ensemble (NPT), the volume fluctuates as the system exchanges work with a barostat, a microstate with energy $E_s$ and volume $V$ has weight $e^{-\beta(E_s + PV)}$, and the corresponding free energy is the Gibbs free energy $G = F + PV$, minimised at constant temperature and pressure. Phase equilibria, reaction equilibria and solubilities at ordinary conditions are governed by $G$. Molecular dynamics simulations of liquids and solids are usually run in this ensemble, with a thermostat and a barostat (the chapter on thermostats).

In the grand canonical ensemble ($\mu VT$), the system exchanges particles as well as energy with a reservoir at fixed chemical potential $\mu$, the free-energy cost of adding one particle. A microstate with $N$ particles and energy $E_s$ has weight $e^{-\beta(E_s - \mu N)}$, and the grand potential $\Omega = -k_BT\ln\Xi = -PV$ follows from the grand partition function $\Xi$. This is the natural ensemble for adsorption: grand canonical Monte Carlo simulations of gases in porous materials such as zeolites and metal–organic frameworks predict adsorption isotherms by letting molecules be inserted and deleted at a chemical potential fixed by the gas pressure.

The pattern is always the same. Each ensemble has a weight for each microstate, a partition function that normalises the weights, and a thermodynamic potential equal to $-k_BT$ times its logarithm. The potentials are related by Legendre transforms, $F = U - TS$, $G = F + PV$, $\Omega = F - \mu N$, each swapping a controlled extensive variable for its conjugate intensive one, and each minimised at equilibrium under its own constraints.

References and further reading

Averages are not the whole story

Thermodynamics treats the energy of a beaker of water in a thermostat, or the number of molecules in a cubic millimetre of air, as definite numbers. Statistical mechanics says they are not. A system in contact with a heat bath constantly exchanges energy with it, so its energy fluctuates about its average; a region of a gas constantly exchanges molecules with its surroundings, so the number inside it fluctuates too. For a macroscopic sample these fluctuations are far too small to notice, which is why thermodynamics works. For the systems studied in molecular simulation, a few thousand atoms, and for the nanoparticles, single molecules and small droplets of modern experiments, they are large enough to see and to matter.

Fluctuations are not merely noise. Their size is fixed by the same statistical mechanics that fixes the averages, and it turns out to be directly related to how the system responds when it is disturbed. The heat capacity, the compressibility and the magnetic susceptibility can all be written as the size of a spontaneous fluctuation. This is one of the most useful results in statistical physics, because it means that response functions can be computed from an equilibrium simulation simply by recording the fluctuations, without ever disturbing the system.

Energy fluctuations and the heat capacity

In the canonical ensemble, the mean energy is $U = -\partial \ln Z/\partial\beta$, as the chapter on ensembles showed. Differentiating once more with respect to $\beta$ gives the variance of the energy, and since $\partial/\partial\beta = -k_BT^2\,\partial/\partial T$, the variance is proportional to the temperature derivative of the energy, which is the heat capacity:

The figure below samples the total energy of $N$ independent two-level units, each with an excited level at energy $\Delta$, in contact with a heat bath. Each unit is excited with probability $p = 1/(1 + e^{\Delta/k_BT})$, so the number of excited units fluctuates from moment to moment. The readout computes the heat capacity from nothing but the variance of the sampled energies and compares it with the exact result. Change the temperature: the fluctuations, and the heat capacity, are largest at intermediate temperatures where the units are most uncertain about their state, the Schottky peak seen in the ensembles topic.

The formula also explains why heat capacities are always positive: a variance cannot be negative. A system with a negative heat capacity would heat up when it lost energy, which cannot happen in a stable system in contact with a bath. (Isolated, self-gravitating systems such as stars can have negative heat capacities, and small isolated clusters near a melting transition can show them in the microcanonical ensemble, a reminder that the equivalence of ensembles holds only for large systems.)

The thermodynamic limit

How large are these fluctuations? For $N$ independent units, the energies add, and so do the variances (the chapter on probability), so the mean grows as $N$ and the standard deviation as $\sqrt N$. The relative fluctuation therefore shrinks:

This holds well beyond independent units: whenever correlations between parts of a system extend over a finite distance, the system behaves like many independent blocks, and the same $1/\sqrt N$ law applies. In the figure, increasing $N$ by a factor of 100 narrows the distribution tenfold. For a mole of matter, $N \approx 6 \times 10^{23}$, the relative fluctuation is about $10^{-12}$, utterly negligible: this is the thermodynamic limit, in which averages become sharp and thermodynamics becomes exact. The exception is near a critical point, where correlations extend over macroscopic distances, the parts of the system cease to be independent, and fluctuations grow enormously; the chapter on phase transitions shows this.

For a simulation of 1000 atoms, by contrast, the relative energy fluctuation is a few per cent, which is easily visible in a plot of the energy against time. That is not a malfunction of the thermostat; it is the canonical ensemble behaving as it should.

Density fluctuations and compressibility

The same reasoning applies to other quantities. Consider a small region inside a fluid, open to the rest. Molecules wander in and out, so the number $n$ inside fluctuates. In the grand canonical ensemble, where the region exchanges particles with a reservoir at fixed chemical potential, differentiating the grand partition function twice with respect to $\mu$ gives the variance of $n$, and thermodynamic identities relate it to the isothermal compressibility $\kappa_T = -\frac1V(\partial V/\partial P)_T$:

Switch the figure to the gas. Molecules are scattered at random through a box, and the shaded sub-volume holds a fraction of them. The ratio of the variance to the mean of the count comes out close to 1, as the formula predicts for an ideal gas, reduced by the sub-volume fraction because the total number in the whole box is fixed. In a liquid, molecules repel each other at short range, so they cannot crowd into one region, the count fluctuates much less than for a gas, and $\rho k_BT\kappa_T$ is small; for water at room temperature it is about 0.06. Liquids are nearly incompressible precisely because their density fluctuations are suppressed.

Density fluctuations are directly observable. They scatter light, X-rays and neutrons: the intensity scattered at small angles is proportional to $\langle(\delta n)^2\rangle/\langle n\rangle$, which is the structure factor $S(q)$ at $q \to 0$ discussed in the chapter on liquids. Near the liquid–gas critical point the compressibility diverges, density fluctuations grow to the size of the wavelength of light, and an otherwise transparent fluid turns milky, the phenomenon of critical opalescence explained by Einstein and Smoluchowski in 1908–1910.

Fluctuations at constant pressure

In the isothermal–isobaric ensemble, the volume of the system is not fixed but fluctuates, and the corresponding formulas are the ones used in most molecular dynamics simulations, which run at constant temperature and pressure:

where $H = E + PV$ is the enthalpy. Recording the box volume and the enthalpy during a simulation therefore gives the compressibility and the constant-pressure heat capacity without any additional calculation. A further formula of the same kind gives the thermal expansion coefficient from the correlation between volume and enthalpy fluctuations. These estimators are widely used to validate force fields against experiment, and their accuracy depends on long simulations, because the variance of a fluctuation converges more slowly than the mean.

The fluctuations of a quantity depend on the ensemble even when its average does not. In a microcanonical (NVE) simulation the total energy does not fluctuate at all, while the kinetic and potential energies fluctuate against each other; in a canonical simulation the total energy fluctuates as above. A thermostat or barostat that gets the average temperature or pressure right but the fluctuations wrong produces the wrong ensemble and wrong response functions. The Berendsen thermostat is the classic example, discussed in the chapter on thermostats.

The fluctuation–dissipation theorem

The formulas above are special cases of a general principle. Suppose a small external field $h$ is applied that couples to some quantity $A$ of the system, adding $-hA$ to the energy: a magnetic field coupling to the magnetisation, an electric field coupling to the dipole moment, a pressure coupling to the volume. The system responds by shifting the average of $A$. To first order in $h$, statistical mechanics gives

The proof is two lines: $\langle A\rangle = \sum_s A_s e^{-\beta(E_s - hA_s)}/Z$, and differentiating with respect to $h$ brings down $\beta A_s$ in each term, giving $\beta(\langle A^2\rangle - \langle A\rangle^2)$. The dielectric constant of water in simulations is computed exactly this way, from the fluctuations of the total dipole moment of the box.

The dynamic version of this result is the fluctuation–dissipation theorem, formulated by Callen and Welton in 1951 and by Kubo in 1957. It states that the way a system relaxes back to equilibrium after a small disturbance is the same as the way its spontaneous fluctuations decay, a principle Onsager had stated in 1931 as the regression hypothesis. Quantitatively, transport coefficients, which measure dissipation, are integrals of equilibrium time-correlation functions, which measure fluctuations. The diffusion coefficient is the time integral of the velocity autocorrelation function, the viscosity is the integral of the autocorrelation of the stress, and the electrical conductivity is the integral of the autocorrelation of the current. These Green–Kubo relations are how transport coefficients are computed from equilibrium simulations, and the chapter on transport coefficients develops them. Einstein’s 1905 relation between the diffusion coefficient of a Brownian particle and the friction it feels, $D = k_BT/\zeta$, was the first example of the theorem: the same molecular collisions that make the particle jiggle also slow it down when it is pushed.

Consequences for simulations

Several practical lessons follow. Fluctuations in a small simulation are real and expected: a 1000-atom canonical simulation at 300 K shows temperature fluctuations of about 3 %, and a pressure that fluctuates by hundreds of bar is normal for a small box of liquid. Response functions computed from fluctuations need long runs, since a variance converges more slowly than a mean and correlated samples count less than they seem (the chapter on probability explains effective sample sizes). Fluctuation formulas are only valid if the simulation samples the correct ensemble, so the choice of thermostat and barostat matters. And finite-size effects are themselves a kind of fluctuation effect: quantities that depend on long-wavelength fluctuations, such as the compressibility near a critical point, converge slowly with system size and should be checked by repeating the calculation in a larger box.

References and further reading

Phases and order parameters

Water freezes at 0 °C, iron loses its magnetism above 770 °C, and the copper and zinc atoms of β-brass, arranged in an ordered alternating pattern at room temperature, mix randomly on the same lattice above about 460 °C. These are phase transitions: changes in the collective organisation of matter that happen abruptly as a control parameter, usually temperature or pressure, passes a definite value. They are the most dramatic consequence of the competition between energy and entropy described in the chapter on free energy. Below the transition, the ordered arrangement has the lower free energy because its energy is lower; above it, the disordered arrangement wins because its entropy, multiplied by the temperature, outweighs the energy cost.

A phase transition is described by an order parameter, a quantity that is zero in the disordered phase and non-zero in the ordered one. For a magnet it is the magnetisation; for an alloy, the excess of one kind of atom on one sublattice; for the liquid–gas transition, the difference between the liquid and gas densities; for a nematic liquid crystal, the average alignment of the molecules. Transitions fall into two broad classes. In a first-order or discontinuous transition, such as melting or boiling, the order parameter jumps, the two phases coexist at the transition, a latent heat is absorbed, and the system can be superheated or supercooled. In a continuous or second-order transition, such as the loss of ferromagnetism at the Curie temperature, the order parameter falls smoothly to zero, there is no latent heat, but the heat capacity, the susceptibility and the size of fluctuations become very large near the critical point.

The Ising model

The simplest model that shows a phase transition was proposed by Wilhelm Lenz and solved in one dimension by his student Ernst Ising in 1925. Each site $i$ of a lattice carries a spin $s_i$ that can point up ($+1$) or down ($-1$). Neighbouring spins prefer to align, and an external field $h$ favours one direction:

Its simplicity is deceptive, because the same mathematics describes very different physical systems. Reading “up” as “site occupied by a molecule” and “down” as “empty site” turns the Ising model into a lattice gas, with $J$ an attraction between neighbouring molecules and $h$ related to the chemical potential; its transition is then the condensation of a gas into a liquid. Reading “up” and “down” as two kinds of atom turns it into a binary alloy, whose transition is ordering or phase separation. Lattice models of this kind, with $J$ fitted to first-principles calculations through a cluster expansion, are used today to predict the ordering of real alloys and of ions in battery electrodes.

Ising found that the one-dimensional chain has no transition at any non-zero temperature: a single break in an aligned chain costs a fixed energy $2J$ but can be placed anywhere, so its entropy always wins. He wrongly concluded that the model was useless in higher dimensions too. In 1944 Lars Onsager solved the two-dimensional square-lattice model exactly, a celebrated feat of mathematical physics, and showed that it orders spontaneously below a critical temperature

The spontaneous magnetisation below $T_c$, published by Chen Ning Yang in 1952, is $m = [1 - \sinh^{-4}(2J/k_BT)]^{1/8}$. These exact results make the 2D Ising model the standard test of approximate theories and of simulation methods.

Simulating the model

The model is easily simulated with the Metropolis Monte Carlo method (its own chapter explains why it samples the canonical ensemble). Repeatedly pick a spin, compute the energy change $\Delta E$ that flipping it would cause, which depends only on its four neighbours, and accept the flip with probability $\min(1, e^{-\Delta E/k_BT})$. One attempted flip per spin is called a sweep.

Start at a temperature well below $T_c$, say $k_BT = 1.5J$, with random spins. The lattice does not order all at once. Small aligned domains of each sign form, and their boundaries move and smooth out, so that domains coarsen over time; eventually one sign usually takes over, although stripes spanning the periodic box can survive for a long time. Above $T_c$, at $k_BT = 3J$, the spins remain a fine-grained mixture and the magnetisation fluctuates around zero. Close to $T_c$, the picture is remarkable: clusters of every size, from single spins to the size of the lattice, coexist and fluctuate, and the magnetisation wanders widely. This is the critical state, in which the system has no characteristic length scale.

To rebuild the phase diagram, set a temperature, wait until the magnetisation has settled, and press Record. Points well below and well above $T_c$ fall on Onsager’s curve. Near $T_c$ they do not: a finite lattice cannot have a sharp transition, so the recorded $\langle|m|\rangle$ is rounded off and stays above zero slightly beyond $T_c$, and the rounding shrinks as the lattice grows. Recording the same temperatures on 32, 64 and 128 spin lattices shows this finite-size effect directly. Close to $T_c$ the simulation also slows down dramatically, because large clusters take many sweeps to change, a phenomenon called critical slowing down. Cluster algorithms such as those of Swendsen and Wang or Wolff, which flip whole correlated clusters at once, were invented to overcome it.

The readouts show the susceptibility and heat capacity estimated from the fluctuations of the magnetisation and energy, using the formulas of the chapter on fluctuations. Both peak near $T_c$, where fluctuations are largest.

Mean-field theory

Onsager’s solution is a mathematical tour de force that does not generalise to three dimensions or to other models. A simple approximate theory that works in any dimension is mean-field theory, introduced by Pierre Weiss in 1907 for ferromagnets. Its idea is to replace the fluctuating neighbours of each spin by their average. Each spin then feels an effective field, the external field plus $J$ times the sum of its $z$ neighbours’ average values, $h_{\text{eff}} = h + zJm$. A single spin in a field $h_{\text{eff}}$ has average value $\tanh(\beta h_{\text{eff}})$, and since every spin is equivalent, this average must equal $m$ itself:

The equation cannot be solved in closed form, but it can be solved by drawing. Plot the straight line $y = m$ and the curve $y = \tanh((zJm + h)/k_BT)$; the solutions are where they cross.

At high temperature the tanh curve is shallow and crosses the line only at $m = 0$: the disordered phase. As the temperature falls, the curve steepens, and when its slope at the origin, $zJ/k_BT$, exceeds 1, it crosses the line at two more points, $\pm m_0$; the solution $m = 0$ becomes unstable and the system orders spontaneously. Mean-field theory thus predicts a transition at $k_BT_c = zJ$, which for the square lattice is $4J$, well above the exact $2.269J$. It overestimates $T_c$ because it ignores fluctuations, which disrupt order; it even predicts a transition in the one-dimensional chain, where there is none. It improves as the number of neighbours grows, and becomes exact when every spin interacts with every other.

Landau theory

Lev Landau generalised the mean-field idea in 1937 into a theory that applies to any transition, by writing the free energy as a function of the order parameter and expanding it in powers of $m$ near the transition, keeping only terms allowed by the symmetry of the problem. For the Ising magnet, flipping every spin leaves the energy unchanged when $h = 0$, so only even powers can appear:

The figure plots this free energy. For the continuous transition, as the temperature passes $T_c$ the single minimum at $m = 0$ flattens and splits into two symmetric minima that move apart continuously; a small field tilts the curve and selects one. The predicted magnetisation grows as $(T_c - T)^{1/2}$ below the transition. Switch to the first-order form: now, over a range of temperatures, the free energy has three minima, one at $m = 0$ and two at large $\pm m$. The equilibrium state jumps from one to the other when their depths cross, which is the discontinuous transition, while the shallower minimum remains as a metastable state that the system can be trapped in. That is the origin of supercooling, superheating and hysteresis, and escaping a metastable minimum requires crossing a free-energy barrier, which is what nucleation of a new phase means.

Critical phenomena and universality

Near a continuous transition, many quantities follow power laws in the distance from the critical temperature, $t = (T - T_c)/T_c$, characterised by critical exponents: the order parameter vanishes as $m \propto (-t)^\beta$, the susceptibility diverges as $\chi \propto |t|^{-\gamma}$, the heat capacity as $C \propto |t|^{-\alpha}$, and the correlation length $\xi$, the typical size of correlated regions, diverges as $\xi \propto |t|^{-\nu}$. The divergence of $\xi$ is the heart of the matter: at the critical point fluctuations are correlated over all length scales, which is what the simulation shows as clusters of every size.

Mean-field theory gets the exponents wrong in two and three dimensions because it neglects exactly these long-range fluctuations. The most striking discovery about critical exponents is universality: they do not depend on the details of a system, only on its dimension and the symmetry of its order parameter. The liquid–gas critical point of carbon dioxide, the critical point of a binary liquid mixture and the Curie point of a uniaxial magnet all have the exponents of the 3D Ising model, even though the forces involved are completely different. The explanation came from the renormalisation group developed by Kenneth Wilson in the early 1970s, for which he received the 1982 Nobel Prize in Physics: near a critical point, as the system is viewed on longer and longer length scales, microscopic details become irrelevant and only symmetry and dimension survive.

Phase transitions in simulations

Simulations of phase transitions must contend with the finite size of the simulated system. A box of $N$ particles cannot show true singularities, so transitions are rounded and shifted, and the critical properties of the infinite system are obtained by finite-size scaling: repeating the simulation for several sizes and extrapolating using the known way in which quantities scale with $L/\xi$. A practical tool is the Binder cumulant of the order parameter, whose curves for different system sizes cross near $T_c$.

First-order transitions pose the opposite problem. The two phases are separated by a free-energy barrier that grows with system size, so a simulation started in one phase can remain there far beyond the true transition temperature: a simulated crystal heated at constant rate often melts well above its melting point, because there is no surface from which melting can start. Locating first-order transitions reliably therefore requires free-energy calculations that compare the two phases directly, coexistence simulations in which both phases share a box with an interface, or biased sampling methods that drive the system over the barrier, described in the subject on enhanced sampling.

References and further reading

Order without a lattice

A crystal is easy to describe: give the lattice and the contents of one unit cell, and every atomic position follows. A gas is easy too: the atoms are almost independent. Liquids, glasses, amorphous solids and dense fluids lie in between. They have no long-range order, yet they are far from random: each atom is surrounded by a well-defined first shell of neighbours at nearly fixed distance, a less sharply defined second shell, and so on, with order fading over a few molecular diameters. The theory of liquids, developed from the 1930s by John Kirkwood, Joseph Mayer, Jean-Pierre Hansen, John Barker, Douglas Henderson and many others, and tested at every stage by computer simulation, describes this statistical order with a few correlation functions. The most important is the radial distribution function.

The radial distribution function

The radial distribution function $g(r)$ measures how the density of particles around a given particle varies with distance, relative to the average density:

The figure runs a two-dimensional Lennard-Jones fluid by molecular dynamics and accumulates $g(r)$ as it goes. At low density, $g(r)$ shows a single weak peak just beyond the repulsive core, reflecting the attraction of neighbouring pairs, and is otherwise close to 1. In the liquid, it develops a strong first peak at about one particle diameter, a minimum, a second peak and weaker oscillations that decay within a few diameters. Compress further and the peaks sharpen and split as the system crystallises into a triangular lattice. The tagged atom (rose) and its first-shell neighbours (teal) show the shell that the first peak counts.

Coordination numbers and thermodynamics from g(r)

Integrating $\rho g(r)$ over the first peak gives the coordination number, the average number of nearest neighbours,

with the integral stopped at the first minimum $r_{\min}$. Simple liquids near their triple point have about 12 neighbours in three dimensions, close to the 12 of a close-packed crystal, which is why melting changes the density of simple substances by only around 10 %; water has only about 4.5 neighbours, the signature of its open, hydrogen-bonded network. In two dimensions, as in the figure, the dense liquid has close to 6 neighbours.

For particles interacting through a pair potential $u(r)$, $g(r)$ determines the thermodynamics. The average potential energy per particle and the pressure are

and the isothermal compressibility follows from the compressibility route, $\rho k_BT\kappa_T = 1 + \rho\int_0^\infty [g(r) - 1]\,4\pi r^2\,\mathrm dr$, which links $g(r)$ to the density fluctuations of the chapter on fluctuations. Because these routes are exact only for the exact $g(r)$, comparing them is a sensitive test of approximate theories such as the Percus–Yevick and hypernetted-chain integral equations.

The potential of mean force

The probability of finding two particles at distance $r$, relative to an ideal gas, defines a free energy as a function of their separation:

Switch the figure to the potential of mean force. At low density the orange curve follows the blue pair potential. At liquid density it develops extra minima at the second and third neighbour distances and a barrier between the first two: there is no direct interaction at those distances, but the surrounding particles make some separations more probable than others. This is the same idea as the free-energy profiles of the chapter on ensembles, applied to a distance. It is important in practice: the potential of mean force between two ions in water shows a contact minimum and a solvent-separated minimum, with a barrier between them corresponding to a water molecule entering or leaving the gap, and its depth determines ion pairing and the stability of colloids and nanoparticles in solution. Coarse-grained models are often built by matching potentials of mean force, by iterative Boltzmann inversion for example.

Scattering and the structure factor

$g(r)$ is not only a simulation concept; it is measured. X-rays and neutrons scattered by a liquid produce an intensity pattern whose variation with the scattering vector $q$ is the static structure factor, related to $g(r)$ by a Fourier transform:

Comparing simulated and measured structure factors is one of the main tests of force fields and MLIPs for liquids, glasses and amorphous materials. For multicomponent systems there are partial structure factors and partial $g_{\alpha\beta}(r)$ for each pair of species, and X-ray and neutron experiments, which weight the species differently, together constrain them; for water, the oxygen–oxygen $g(r)$ from X-ray diffraction is the standard benchmark.

Hard spheres and the van der Waals picture

The structure of simple dense liquids is dominated by repulsion. Simulations of hard spheres by Berni Alder and Thomas Wainwright (1957) and by William Wood and J. D. Jacobson showed that even particles with no attraction at all crystallise at high density, an entropy-driven transition, and that their liquid $g(r)$ looks remarkably like that of real simple liquids. This led to the modern van der Waals picture of liquids, developed by John Weeks, David Chandler and Hans Andersen (1971) and by Barker and Henderson: the structure is fixed by packing of the repulsive cores, and the attraction mainly provides a uniform cohesive background that holds the liquid together. Perturbation theories built on this idea predict liquid thermodynamics accurately and underlie much of the theory of fluids.

Computing g(r) in a simulation

The recipe is simple, and the figure follows it: for every pair of particles compute the distance, applying the minimum-image convention in a periodic box; histogram the distances into bins of width $\Delta r$; average over many frames; and normalise each bin by the number of pairs an ideal gas would have in that shell, $\tfrac12 N\rho\cdot 4\pi r^2\Delta r$ (counting each pair once) or $N\rho\cdot 4\pi r^2\Delta r$ (counting both orders). Distances beyond half the box length are unreliable, because shells there are cut off by the periodic boundary. Bins must be fine enough to resolve the first peak but coarse enough to collect statistics, and the result should be averaged over frames separated by more than the correlation time. Analysis libraries such as MDAnalysis, MDTraj and OVITO compute $g(r)$, coordination numbers and structure factors directly from trajectories, as the chapter on trajectory analysis describes.

References and further reading