Molecular simulation

Atoms move on a landscape

Almost every atomistic simulation, classical or machine-learned, rests on one approximation. Nuclei are thousands of times heavier than electrons, so on the timescale of nuclear motion the electrons relax instantly to their ground state. This Born–Oppenheimer approximation separates the problem: for each fixed arrangement of nuclei $\mathbf R = (\mathbf r_1, \dots, \mathbf r_N)$, solving the electronic Schrödinger equation gives one number, the energy $E(\mathbf R)$. As the nuclei move, that number traces out a $3N$-dimensional landscape, the potential energy surface (PES).

The force on atom $i$ is the downhill slope of the landscape:

Everything else follows from $E(\mathbf R)$. Minima are stable structures: molecules, crystals, adsorption geometries. First-order saddle points connect minima and are transition states; their height sets reaction rates. The curvature at a minimum, the Hessian $\partial^2 E/\partial\mathbf r_i\partial\mathbf r_j$, gives vibrational frequencies and phonons. Molecular dynamics moves the atoms across the surface under these forces.

Computing $E(\mathbf R)$ from quantum mechanics, typically with density functional theory (DFT), is accurate but expensive: the cost grows roughly as $N^3$ with the number of electrons, which limits routine DFT molecular dynamics to a few hundred atoms for tens of picoseconds. A force field or interatomic potential replaces the quantum calculation with a cheap analytic function of the positions. This lesson is about the classical forms; the lesson on interatomic potentials replaces them with machine-learned ones.

Pair potentials

The simplest model assumes the energy is a sum over pairs that depends only on distance, $E = \sum_{i<j} u(r_{ij})$. The force between a pair is then directed along the line joining them, with magnitude $F(r) = -\mathrm du/\mathrm dr$.

Lennard-Jones. For closed-shell atoms such as argon,

The $r^{-6}$ attraction has a physical origin: it is the London dispersion interaction between fluctuating dipoles. The $r^{-12}$ repulsion, standing in for Pauli repulsion between overlapping electron clouds, was chosen mainly because it is the square of $r^{-6}$ and cheap to compute; the true repulsion is closer to exponential. The minimum lies at $r_{\min} = 2^{1/6}\sigma \approx 1.12\,\sigma$ with depth $-\varepsilon$, and the potential crosses zero at $r = \sigma$.

Morse. For a chemical bond that can dissociate,

with dissociation energy $D_e$, equilibrium length $r_0$ and width parameter $a$. Unlike a spring, it levels off at large $r$: bonds break.

In the lab, notice three things. The force is zero exactly at the energy minimum and largest in magnitude where the energy curve is steepest. Repulsion at short range is far steeper than attraction at long range. And the harmonic approximation (dashed) is good only near the minimum: real bonds are anharmonic, softer when stretched and stiffer when compressed, which is why thermal expansion exists.

The harmonic approximation and vibrations

Expanding any potential to second order about its minimum gives a spring,

and a diatomic of reduced mass $\mu = m_1m_2/(m_1+m_2)$ vibrates with angular frequency $\omega = \sqrt{k/\mu}$, usually quoted as a wavenumber $\tilde\nu = \omega/(2\pi c)$. For the Morse potential, $k = 2D_e a^2$. The H₂ parameters in the lab give about 4400 cm⁻¹, close to the measured harmonic frequency; the C–H bond gives about 2900 cm⁻¹, the familiar C–H stretch region of an infrared spectrum. These numbers matter for simulation: the fastest vibration in a system sets the largest stable time step (the molecular dynamics lesson).

Anatomy of a classical force field

Biomolecular and soft-matter force fields such as AMBER, CHARMM, OPLS-AA and GROMOS write the energy as a sum of simple terms:

The first line, the bonded terms, holds the covalent skeleton together. The second, the non-bonded terms, acts between all atoms not closely bonded: van der Waals interactions and electrostatics between fixed partial charges. Parameters for unlike pairs usually come from combining rules, for example Lorentz–Berthelot: $\sigma_{ij} = (\sigma_i + \sigma_j)/2$, $\varepsilon_{ij} = \sqrt{\varepsilon_i\varepsilon_j}$.

Every parameter is fitted, to quantum-chemical calculations (geometries, torsion scans, electrostatic potentials) and to experiments (densities, heats of vaporisation, hydration free energies). A force field is therefore a model with a domain of validity, and it can fail outside it. Two failures recur.

• Fixed charges cannot polarise. A molecule next to a metal surface or an ion has its electron density distorted; fixed point charges cannot respond. Polarisable force fields add induced dipoles or Drude oscillators, at a cost. For gold surfaces, for example, the GolP-CHARMM force field adds a mobile image charge to every gold atom.
• Combining rules are approximations. In concentrated solutions, the default cross-interactions between unlike groups can be wrong enough to cause artificial aggregation. Pair-specific corrections (NBFIX in CHARMM) are fitted to experimental osmotic coefficients to fix them.

Other material classes need other functional forms. Metals are poorly described by pair potentials because metallic bonding depends on the local electron density; the embedded-atom method (EAM) adds a many-body embedding term $F_i(\bar\rho_i)$ for each atom. Covalent solids such as silicon need angle-dependent terms (Stillinger–Weber, Tersoff). Reactive force fields such as ReaxFF let bond orders change continuously so that bonds can form and break.

In practice

Check any potential you implement against a finite-difference derivative. For Lennard-Jones argon:

import numpy as np

eps, sig = 0.997, 3.405                    # argon: kJ/mol, Å
def U(r):                                  # Lennard-Jones energy
    s6 = (sig / r) ** 6
    return 4 * eps * (s6 * s6 - s6)
def F(r):                                  # force = -dU/dr (positive = repulsive)
    s6 = (sig / r) ** 6
    return 24 * eps * (2 * s6 * s6 - s6) / r

r = np.linspace(3.2, 8.0, 2001)
rmin = r[np.argmin(U(r))]
print(f"minimum at r = {rmin:.3f} Å   (2^(1/6) σ = {2**(1/6)*sig:.3f} Å)")
print(f"U(rmin) = {U(rmin):.3f} kJ/mol (= -ε)")

h = 1e-5                                   # check the analytic force by finite differences
print(f"F(4.0) analytic {F(4.0):.6f}, numerical {-(U(4.0+h)-U(4.0-h))/(2*h):.6f}")

minimum at r = 3.822 Å   (2^(1/6) σ = 3.822 Å)
U(rmin) = -0.997 kJ/mol (= -ε)
F(4.0) analytic -0.544033, numerical -0.544033

References and further reading

Origins

The averages of statistical mechanics are integrals over all configurations, weighted by the Boltzmann factor, in a space of thousands of dimensions. The chapter on integrals showed why such integrals cannot be done on a grid and why random sampling, with an error falling as $1/\sqrt N$ whatever the dimension, is the only practical approach. The question is how to sample. Drawing configurations uniformly at random is useless: almost all of them have overlapping atoms and a Boltzmann weight of essentially zero. The samples must be drawn from the Boltzmann distribution itself.

The solution was found at Los Alamos at the dawn of electronic computing. Stanislaw Ulam and John von Neumann had used random sampling for neutron-diffusion calculations in the late 1940s, and Nicholas Metropolis named the approach after the Monte Carlo casino; Metropolis and Ulam described it in 1949. In 1953 Metropolis, Arianna and Marshall Rosenbluth, and Augusta and Edward Teller published “Equation of state calculations by fast computing machines”, in which they computed the equation of state of a two-dimensional system of hard disks on the MANIAC computer. Their algorithm, now called the Metropolis algorithm, generates a sequence of configurations that visits each one with a frequency proportional to its Boltzmann weight, without ever computing the partition function. It is among the most influential algorithms of the twentieth century.

Markov chains and detailed balance

The algorithm generates a Markov chain: a sequence of configurations in which each new one is chosen at random with a probability that depends only on the current one, $\pi(o \to n)$ from old to new. If the chain is run long enough, it settles into a stationary distribution $p(\mathbf x)$ that no longer changes from step to step. The design problem is to choose $\pi$ so that the stationary distribution is the Boltzmann distribution. A sufficient condition is detailed balance: in the stationary state, the flow of probability from any configuration to any other is exactly balanced by the reverse flow,

The transition is split into two stages: propose a trial configuration $n$ with probability $\alpha(o \to n)$, then accept it with probability $\mathrm{acc}(o \to n)$, otherwise keep $o$. If the proposal is symmetric, $\alpha(o \to n) = \alpha(n \to o)$, detailed balance only constrains the ratio of acceptance probabilities, and the Metropolis choice satisfies it in the simplest way:

The last point is often missed and causes real errors: when a move is rejected, the current configuration must be counted again in all averages. Discarding rejected steps biases the averages towards configurations that are easy to leave.

The algorithm in practice

For a system of particles, a Monte Carlo step consists of choosing a particle, displacing it by a random vector drawn uniformly from a small cube of side $2\Delta$, computing the change in energy, and accepting or rejecting by the Metropolis rule. Only the interactions of the moved particle change, so the energy difference costs $O(N)$ (or $O(1)$ with neighbour lists), not a full energy evaluation. One attempted move per particle is a cycle or sweep.

The figure runs the algorithm for a single coordinate in a double-well potential. Each arrow is a trial move: green if accepted, rose if rejected. The histogram of visited positions converges to the Boltzmann distribution.

The step size $\Delta$ controls efficiency. Very small steps are almost always accepted but move the system hardly at all, so successive configurations are nearly identical. Very large steps jump far but are almost always rejected, so the system stays put. Between these extremes lies an optimum, which the figure measures through the integrated autocorrelation time $\tau_{\text{int}}$: the number of steps needed to produce one effectively independent sample, as defined in the chapter on probability. A common rule of thumb targets an acceptance rate of 30–50 %, but the real criterion is the smallest $\tau_{\text{int}}$ for the quantities of interest, and in a double well the optimal step may be the one that can occasionally jump the barrier. Codes adjust $\Delta$ during equilibration to reach a target acceptance, then keep it fixed during production, because changing the move during sampling breaks detailed balance.

Monte Carlo in other ensembles

Molecular dynamics naturally samples the microcanonical ensemble, and other ensembles require thermostats and barostats. Monte Carlo is the opposite: ensembles are chosen simply by choosing which moves to make and how to accept them, each derived from detailed balance in that ensemble.

In the isothermal–isobaric (NPT) ensemble, the volume becomes a variable. A volume move rescales the box and all positions from $V$ to $V'$ and is accepted with

where the last term comes from the change in the volume available to $N$ scaled coordinates (for moves uniform in $V$).

In the grand canonical ($\mu VT$) ensemble, particles are inserted at random positions and deleted at random, with

where $\Lambda$ is the thermal wavelength. Grand canonical Monte Carlo is the standard method for computing adsorption isotherms of gases in porous materials such as zeolites and metal–organic frameworks, with the chemical potential set from the gas pressure through an equation of state; codes such as RASPA are built around it.

The Gibbs ensemble method of Athanassios Panagiotopoulos (1987) simulates two coexisting phases in two separate boxes without an interface, exchanging volume and particles between them so that pressure and chemical potential equalise; it gives liquid–vapour coexistence curves directly.

Smarter moves

The freedom to invent moves is Monte Carlo’s great strength: moves need not be physical, only correct. Inserting a long chain molecule at a random position almost always produces an overlap and is rejected. Configurational-bias Monte Carlo (J. Ilja Siepmann and Daan Frenkel; Frenkel, G. C. A. M. Mooij and Berend Smit, 1992) grows the chain segment by segment, choosing each segment among several trial positions with probability proportional to its Boltzmann weight, and corrects for the bias in the acceptance rule with a Rosenbluth weight; it made simulations of alkanes in zeolites and of polymer melts practical. Cluster moves, such as the Swendsen–Wang (1987) and Wolff (1989) algorithms for spin models, flip whole correlated clusters at once and largely eliminate critical slowing down near phase transitions (the chapter on phase transitions). Hybrid Monte Carlo (Simon Duane, Anthony Kennedy, Brian Pendleton and Duncan Roweth, 1987) proposes a new configuration by running a short molecular dynamics trajectory with random initial velocities, and accepts it with a Metropolis test on the change in total energy, combining the large collective moves of dynamics with the exactness of Monte Carlo. Replica-exchange moves between simulations at different temperatures are also Monte Carlo moves (the chapter on replica exchange).

Monte Carlo or molecular dynamics?

Both sample equilibrium ensembles, and for equilibrium properties the choice is practical. Molecular dynamics provides genuine time evolution and therefore transport properties and kinetics; it moves all particles collectively and parallelises naturally; and it needs forces, which every modern potential supplies. Monte Carlo provides no dynamics (its “time” is the number of trials), but it can sample ensembles with fluctuating particle numbers, use unphysical moves that jump over barriers, handle discontinuous potentials such as hard spheres, and needs only energy differences. Systems with discrete degrees of freedom, such as lattice models, alloy orderings and adsorption, are natural for Monte Carlo; liquids, biomolecules and anything requiring dynamics are natural for MD. Kinetic Monte Carlo, discussed in the chapter on rates, is a different method that uses rates to generate true stochastic dynamics.

References and further reading

The idea

Molecular dynamics (MD) solves Newton’s equations of motion for every atom,

one small time step at a time, starting from some positions and velocities. The result is a trajectory, a movie of the atoms. Averages over the trajectory give thermodynamic and structural properties; correlations in time give transport properties. The method dates from Alder and Wainwright’s hard-sphere simulations (1957) and Rahman’s simulation of liquid argon (1964); today it runs routinely on millions of atoms.

Two facts make MD work in practice. First, a good integrator keeps the total energy constant for millions of steps even though each step has an error. Second, statistical mechanics tells us that a long trajectory samples the same states as an ensemble average (the ergodic hypothesis), so time averages from one simulation give equilibrium properties.

Integrating the equations: Verlet

Expand the position forward and backward in time by a step $\Delta t$:

Adding the two expansions cancels the odd terms and gives the Verlet algorithm (Verlet, 1967),

The equivalent form used in most codes, velocity Verlet, carries velocities explicitly:

One force evaluation per step, which is where nearly all the cost lies.

Why Verlet and not something more accurate? The obvious scheme, explicit Euler ($\mathbf r \mathrel{+}= \mathbf v\Delta t$, $\mathbf v \mathrel{+}= \mathbf a\Delta t$), has a smaller per-step error than you might think, yet its energy grows without limit. Verlet is time-reversible and symplectic: it preserves the volume of phase space exactly, as the true dynamics does. A deep consequence is that a symplectic integrator exactly conserves a slightly perturbed “shadow” Hamiltonian $\tilde H = H + \mathcal O(\Delta t^2)$. The energy therefore fluctuates around a constant value with an amplitude proportional to $\Delta t^2$, with no systematic drift for a well-chosen step. Higher-order, non-symplectic methods such as Runge–Kutta are more accurate over a few steps but drift over millions of steps.

Things to try in the lab:

1. NVE with Verlet at Δt = 0.005. Total energy is flat to four decimal places; kinetic and potential energy trade places.
2. Euler. The total energy rises steadily: the integrator is pumping energy in. The system heats up until it boils.
3. Verlet with Δt = 0.02–0.03. Fluctuations grow roughly as $\Delta t^2$ until close collisions are no longer resolved and the energy jumps. The NumPy code below shows the same trend quantitatively.
4. Langevin thermostat at T* = 0.4 and density 0.85, then Restart. The atoms lock into a close-packed (triangular) arrangement, with defects because 100 atoms do not tile a square box perfectly. In the g(r) plot the second peak splits into two, at about $\sqrt3\,\sigma$ and $2\sigma$, the signature of hexagonal order. At T* = 2 and density 0.3 the system is a dense gas with a single weak peak.

Choosing the time step

The step must resolve the fastest motion in the system, with about 10 to 20 steps per period. The fastest motions in molecules are bond stretches involving hydrogen: a C–H stretch at about 3000 cm⁻¹ has a period of about 11 fs, so unconstrained all-atom MD uses $\Delta t \approx 0.5$–1 fs. Constraining bonds to hydrogen with SHAKE or LINCS removes these motions and allows 2 fs, the standard in biomolecular simulation; hydrogen mass repartitioning (moving mass from heavy atoms to their hydrogens) allows about 4 fs. Coarse-grained models, which have no fast vibrations, use 10–40 fs. Monoatomic Lennard-Jones argon is comfortable at about 5 fs, $0.002$–$0.005\tau$ in reduced units.

Simulating bulk matter with a few thousand atoms

Periodic boundary conditions. A box of $10^4$ atoms would be mostly surface. Instead the box is surrounded by infinite copies of itself: an atom leaving through one face re-enters through the opposite one. Each atom interacts with the nearest image of every other atom (the minimum-image convention), which is why the box must be more than twice the interaction cut-off. Periodicity suppresses fluctuations longer than the box, and some properties, notably diffusion coefficients, need finite-size corrections.

Cut-offs and neighbour lists. Short-range interactions are truncated at a cut-off $r_c$, typically 1.0–1.2 nm in biomolecular force fields or $2.5\sigma$ for Lennard-Jones. Shifting the potential so that it reaches zero at $r_c$ avoids a jump in energy, and long-range dispersion corrections can be added analytically. Rebuilding a list of neighbours within $r_c$ plus a skin every few steps makes the cost proportional to $N$ instead of $N^2$.

Electrostatics cannot be cut off. The Coulomb interaction decays as $1/r$ and its sum over a periodic system converges only conditionally. Ewald summation splits it into a short-range part summed in real space and a smooth long-range part summed in reciprocal space; the particle-mesh Ewald (PME) method does the reciprocal sum with fast Fourier transforms in $\mathcal O(N\log N)$ time. Every serious biomolecular or electrolyte simulation uses PME.

Temperature and pressure control

Plain MD conserves energy and samples the microcanonical (NVE) ensemble. Experiments are usually at constant temperature and pressure. The instantaneous temperature comes from equipartition,

where $N_{\text{df}}$ is the number of degrees of freedom ($3N - 3$ in 3D after removing centre-of-mass motion, $2N - 2$ in the 2D lab). A thermostat modifies the dynamics so that the temperature fluctuates correctly around a target.

• Berendsen rescales velocities so that $T$ relaxes exponentially to the target with time constant $\tau_T$. It is robust for equilibration but suppresses the natural fluctuations of the kinetic energy, so it does not generate the canonical ensemble and should not be used for production runs.
• Nosé–Hoover adds a fictitious heat-bath variable with its own equation of motion. It is deterministic and samples the canonical ensemble, though a single thermostat can be non-ergodic for small or stiff systems (Nosé–Hoover chains fix this).
• Stochastic velocity rescaling (Bussi, Donadio and Parrinello) is a Berendsen-like scheme with a correctly designed random term. It samples the canonical ensemble and is the default in GROMACS (v-rescale).
• Langevin dynamics adds friction and a random force to every atom, $m\dot{\mathbf v} = \mathbf F - \gamma m\mathbf v + \sqrt{2\gamma m k_BT}\,\boldsymbol\xi(t)$. It is canonical and very stable, but the friction slows diffusion, so dynamical properties depend on $\gamma$.

Barostats do the same for pressure by letting the box volume (and shape) fluctuate: Parrinello–Rahman and the Martyna–Tobias–Klein (MTK) equations are canonical for production; Berendsen and the stochastic cell-rescaling method are used for equilibration or, in the latter case, also production.

A typical workflow

1. Build the system: structure, topology (force-field parameters), solvent and ions.
2. Minimise the energy to remove bad contacts.
3. Equilibrate at constant volume and temperature (NVT) for tens to hundreds of picoseconds, often with restraints on the solute.
4. Equilibrate at constant pressure (NPT) until the density is stable.
5. Production: a long run with canonical thermostat and barostat, saving coordinates every few picoseconds.
6. Analyse, checking equilibration and statistical errors (the lesson on trajectory analysis).

In GROMACS each stage is a preprocess-then-run pair (each .mdp file holds the stage’s settings):

gmx grompp -f em.mdp  -c solvated.gro -p topol.top -o em.tpr
gmx mdrun  -deffnm em                        # energy minimisation
gmx grompp -f nvt.mdp -c em.gro  -r em.gro  -p topol.top -o nvt.tpr
gmx mdrun  -deffnm nvt                       # NVT equilibration
gmx grompp -f npt.mdp -c nvt.gro -t nvt.cpt -r nvt.gro -p topol.top -o npt.tpr
gmx mdrun  -deffnm npt                       # NPT equilibration
gmx grompp -f md.mdp  -c npt.gro -t npt.cpt -p topol.top -o md.tpr
gmx mdrun  -deffnm md -nb gpu                # production, non-bonded work on the GPU
gmx energy -f md.edr -o density.xvg          # pick "Density" from the menu

The same physics in LAMMPS, the classic Lennard-Jones melt in reduced units:

units           lj
atom_style      atomic
lattice         fcc 0.8442
region          box block 0 10 0 10 0 10
create_box      1 box
create_atoms    1 box
mass            1 1.0
velocity        all create 1.44 87287 loop geom
pair_style      lj/cut 2.5
pair_coeff      1 1 1.0 1.0 2.5
neighbor        0.3 bin
fix             1 all nve
thermo          50
run             250

MD in NumPy, and the energy test

A complete velocity-Verlet code fits on one screen. Run it to see how the energy drift depends on the time step (64 atoms in 2D, reduced units):

import numpy as np

def lj_forces(x, L, rc=2.5):
    """Forces and potential energy for a periodic Lennard-Jones system (reduced units)."""
    d = x[:, None, :] - x[None, :, :]            # all pair vectors, shape (N, N, dim)
    d -= L * np.round(d / L)                      # minimum-image convention
    r2 = (d ** 2).sum(-1)
    np.fill_diagonal(r2, np.inf)
    mask = r2 < rc * rc
    inv6 = np.where(mask, r2 ** -3, 0.0)
    ec = 4 * (rc ** -12 - rc ** -6)               # shift so U(rc) = 0
    pe = 0.5 * np.sum(np.where(mask, 4 * (inv6 ** 2 - inv6) - ec, 0.0))
    fmag = np.where(mask, 24 * (2 * inv6 ** 2 - inv6) / r2, 0.0)
    return (fmag[:, :, None] * d).sum(1), pe

def run(dt, nsteps=2000, seed=1):
    rng = np.random.default_rng(seed)
    n, rho = 8, 0.6                                # 8 x 8 = 64 atoms in 2D
    L = np.sqrt(n * n / rho)
    g = (np.arange(n) + 0.5) * L / n
    x = np.array([[a, b] for a in g for b in g])
    v = rng.normal(size=x.shape); v -= v.mean(0)
    v *= np.sqrt(1.0 / (v ** 2).sum(1).mean() * 2)   # T* = 1 (2 degrees of freedom per atom)
    f, pe = lj_forces(x, L)
    energies = []
    for step in range(nsteps):                     # velocity Verlet
        v += 0.5 * dt * f
        x += dt * v
        x %= L
        f, pe = lj_forces(x, L)
        v += 0.5 * dt * f
        energies.append(0.5 * (v ** 2).sum() + pe)
    e = np.array(energies) / len(x)
    return e[-1] - e[0], e.std()

for dt in (0.002, 0.005, 0.01, 0.02):
    drift, fluct = run(dt)
    print(f"dt = {dt:.3f}:  drift of E/N = {drift:+.5f},  fluctuation = {fluct:.5f}")

dt = 0.002:  drift of E/N = -0.00007,  fluctuation = 0.00003
dt = 0.005:  drift of E/N = -0.00078,  fluctuation = 0.00022
dt = 0.010:  drift of E/N = -0.00163,  fluctuation = 0.00071
dt = 0.020:  drift of E/N = +0.03473,  fluctuation = 0.01650

Halving the step from 0.010 to 0.005 cuts the fluctuation by about a factor of three, close to the $\Delta t^2$ scaling expected of Verlet; at 0.020 close encounters are no longer resolved and the energy jumps. Checking energy conservation in a short NVE run is the first test of any new force field, MLIP or integrator. It catches wrong forces (forces that are not the gradient of the energy) immediately.

With ASE, the same physics for 108 three-dimensional argon atoms:

import numpy as np
from ase.lattice.cubic import FaceCenteredCubic
from ase.calculators.lj import LennardJones
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution, Stationary
from ase.md.verlet import VelocityVerlet
from ase import units

# 108 argon atoms on an fcc lattice; LJ parameters in eV and Å
atoms = FaceCenteredCubic("Ar", size=(3, 3, 3), latticeconstant=5.26)
atoms.calc = LennardJones(epsilon=0.0104, sigma=3.40, rc=8.5, smooth=True)
MaxwellBoltzmannDistribution(atoms, temperature_K=60, rng=np.random.default_rng(0))
Stationary(atoms)                                  # remove centre-of-mass motion

dyn = VelocityVerlet(atoms, timestep=5 * units.fs)
def report():
    n = len(atoms)
    epot, ekin = atoms.get_potential_energy() / n, atoms.get_kinetic_energy() / n
    print(f"t = {dyn.get_time() / units.fs:6.0f} fs   Epot {epot:.4f}   Ekin {ekin:.4f}   "
          f"Etot {epot + ekin:.5f} eV/atom   T = {ekin / (1.5 * units.kB):5.1f} K")
dyn.attach(report, interval=200)
dyn.run(1000)

t =      0 fs   Epot -0.0789   Ekin 0.0080   Etot -0.07089 eV/atom   T =  62.2 K
t =   1000 fs   Epot -0.0748   Ekin 0.0040   Etot -0.07088 eV/atom   T =  30.6 K
t =   2000 fs   Epot -0.0748   Ekin 0.0039   Etot -0.07088 eV/atom   T =  30.5 K
...

Notice the temperature falling from 60 K to 30 K while the total energy stays constant to $10^{-5}$ eV/atom. The run started from a perfect lattice, at the bottom of the potential well; half of the initial kinetic energy flows into potential energy as the atoms start to vibrate, exactly as the equipartition theorem predicts for a harmonic solid. This is why simulations are equilibrated before anything is measured.

References and further reading

Why control temperature

Newton’s equations conserve energy, so plain molecular dynamics samples the microcanonical (NVE) ensemble: constant number of particles, volume and energy. Experiments are almost always done at constant temperature, and usually at constant pressure. Simulations also need temperature control for practical reasons: to equilibrate a system prepared far from equilibrium, to remove the heat produced by a driving force (as in shear or current flow), and to compensate for the slow energy drift caused by numerical integration errors over long runs.

A thermostat modifies the equations of motion so that the system samples the canonical (NVT) ensemble at a chosen temperature $T_0$. The temperature it controls is the kinetic temperature, defined through equipartition (the chapter on Boltzmann and Maxwell) as

where $N_{\text{df}}$ is the number of degrees of freedom ($3N$ minus constraints and conserved momenta). In the canonical ensemble the kinetic energy is not constant: it fluctuates, with a gamma distribution whose relative width is $\sqrt{2/N_{\text{df}}}$, about 3 % for a thousand atoms. A thermostat that holds $K$ exactly fixed, or suppresses its fluctuations, does not produce the canonical ensemble, and every quantity derived from fluctuations, such as the heat capacity (the chapter on fluctuations), will be wrong. A good thermostat must get the whole distribution right, not just the average, and ideally disturb the dynamics as little as possible.

Rescaling velocities, and why it is not enough

The simplest idea is to multiply all velocities at every step by the factor that makes the kinetic temperature exactly $T_0$, $\lambda = \sqrt{T_0/T}$. This velocity rescaling (Leslie Woodcock, 1971) holds the temperature perfectly constant and therefore removes all its fluctuations. The Berendsen thermostat (Herman Berendsen and colleagues, 1984) softens it by relaxing the temperature exponentially towards $T_0$ with a coupling time $\tau$:

Berendsen coupling was for decades the most popular thermostat because it is stable and gentle, but its incorrect ensemble has consequences. Beyond wrong fluctuations, it can transfer energy systematically between degrees of freedom: in simulations of solvated biomolecules, kinetic energy drained from fast internal motions into overall translation and rotation, so that the molecule “froze” while flying through the box, an artefact named the flying ice cube by Stephen Harvey, Robert Tan and Thomas Cheatham (1998). Berendsen coupling is now recommended only for equilibration.

Stochastic thermostats: Andersen and Langevin

A thermostat that samples the canonical ensemble exactly is obtained by mimicking a heat bath with random collisions. In the Andersen thermostat (Hans Andersen, 1980), each particle is occasionally, with probability $\nu\Delta t$ per step, given a new velocity drawn from the Maxwell–Boltzmann distribution at $T_0$, as if it had collided with a bath particle. Between collisions the dynamics is Newtonian. The method samples the canonical ensemble exactly, but the collisions randomise momentum and disturb dynamical properties such as diffusion.

The Langevin thermostat adds a friction force and a random force to every particle:

Langevin dynamics is the physical model of a particle in a solvent, the equation Paul Langevin wrote down in 1908 for Brownian motion, and as a thermostat it is robust, simple and exactly canonical. Its friction slows diffusion and other dynamics, increasingly so as $\gamma$ grows, so weak coupling (a friction time $1/\gamma$ of about a picosecond or longer) is used when dynamics matters.

The Nosé–Hoover extended system

A deterministic route to the canonical ensemble was found by Shuichi Nosé in 1984 and reformulated by William Hoover in 1985. The system is coupled to an extra dynamical variable, a friction coefficient $\xi$ with a fictitious mass $Q$, that speeds up or slows down all particles depending on whether the system is too cold or too hot:

Nosé–Hoover dynamics is deterministic, time-reversible and samples the canonical ensemble if it is ergodic. For a single harmonic oscillator, or for stiff, weakly coupled systems, it is not: the trajectory stays on a subset of phase space and the distribution is wrong. Glenn Martyna, Michael Klein and Mark Tuckerman (1992) fixed this with Nosé–Hoover chains, in which the thermostat variable is itself thermostatted by another, and so on; chains of a few thermostats are ergodic in practice and are a standard choice in many codes.

Stochastic velocity rescaling

In 2007 Giovanni Bussi, Davide Donadio and Michele Parrinello combined the simplicity of rescaling with exact canonical sampling in the canonical sampling through velocity rescaling (CSVR) thermostat. Instead of rescaling the kinetic energy to its target value, they rescale it to a value drawn stochastically, so that the kinetic energy follows exactly the canonical distribution:

CSVR samples the canonical ensemble for any coupling time, is as stable as Berendsen, and perturbs dynamical properties very little because it acts only on the total kinetic energy. It is the default in GROMACS (tcoupl = v-rescale) and widely used elsewhere.

The figure runs a small two-dimensional Lennard-Jones fluid under each thermostat and histograms the kinetic energy against the canonical gamma distribution.

Velocity rescaling collapses the distribution to a spike. Berendsen produces a distribution that is too narrow, increasingly so for short coupling times. NVE dynamics also gives a narrower distribution, because kinetic energy can only be exchanged with the potential energy of a small isolated system. Andersen, Langevin, Nosé–Hoover (this fluid is chaotic enough to be ergodic) and CSVR all reproduce the canonical fluctuations.

Barostats

Controlling pressure requires letting the volume, and possibly the shape, of the simulation box change. The instantaneous pressure is computed from the virial, $P = \frac{1}{dV}\left(2K + \sum_i \mathbf r_i\cdot\mathbf F_i\right)$ in $d$ dimensions (with the pairwise form used for periodic systems), and in the isothermal–isobaric ensemble the volume fluctuates with a variance set by the compressibility, as the chapter on fluctuations explained.

The Berendsen barostat scales coordinates and box to relax the pressure exponentially towards its target; like its thermostat, it equilibrates well but does not sample the correct ensemble, suppressing volume fluctuations. Andersen’s barostat (1980) treats the volume as a dynamical variable with a fictitious mass, driven by the difference between internal and external pressure, and samples the NPH ensemble correctly. Parrinello–Rahman (Michele Parrinello and Aneesur Rahman, 1981) extends this to the full cell matrix, allowing the box shape to change, which is essential for studying solid–solid phase transitions and anisotropic materials; their simulations of structural transformations in crystals under stress were the first of their kind. Martyna–Tobias–Klein (MTK, 1994) equations combine extended-system thermostats and barostats in a form that samples the NPT ensemble exactly. Stochastic cell rescaling (Mattia Bernetti and Giovanni Bussi, 2020) applies the stochastic idea of CSVR to the volume, giving a first-order barostat that is robust like Berendsen but samples correct volume fluctuations; it is available in GROMACS as pcoupl = c-rescale. Monte Carlo volume moves, as in the chapter on Monte Carlo, are another exact alternative.

Practical recommendations

For equilibration from a poor starting structure, Berendsen thermostats and barostats (or stochastic cell rescaling) are robust and quick. For production runs in which ensemble properties matter, use a thermostat that samples the canonical ensemble: CSVR, Langevin, or Nosé–Hoover chains; and a barostat that samples the correct volume distribution: Parrinello–Rahman, MTK or stochastic cell rescaling. Coupling times of 0.1–1 ps for thermostats and 1–5 ps for barostats are typical. When transport properties such as diffusion or viscosity are computed, use weak coupling or a global thermostat (CSVR, Nosé–Hoover), or run the measurement in NVE after equilibrating, because local stochastic thermostats such as Langevin or Andersen with strong coupling alter the dynamics being measured. And for systems with large, stiff, weakly coupled degrees of freedom, such as a solid at low temperature, check ergodicity: a single Nosé–Hoover thermostat may fail where chains or stochastic thermostats succeed.

References and further reading

Averages, and why the error bar is the hard part

Statistical mechanics says an equilibrium property is an ensemble average, $\langle A\rangle = \int A(\mathbf x)\,p(\mathbf x)\,\mathrm d\mathbf x$, over all microstates $\mathbf x$ weighted by their Boltzmann probability. A simulation replaces this with a time average over $N$ saved frames,

which converges to $\langle A\rangle$ if the trajectory is long enough to visit the relevant states (ergodicity). Computing $\bar A$ is trivial. Knowing how wrong it is is not, because successive frames are correlated: a frame saved 10 fs after the last one is nearly the same configuration. Treating $N$ correlated frames as independent gives an error bar that is too small, sometimes by an order of magnitude.

Correlation time and effective sample size

The autocorrelation function of the fluctuations $\delta A_t = A_t - \bar A$ measures how long the signal remembers itself:

For a correlated series of variance $\sigma^2$, the variance of the mean is

with the integrated autocorrelation time $\tau_{\text{int}}$ in units of the sampling interval ($\tau_{\text{int}} = \tfrac12$ for independent samples). Equivalently, the data contain only

independent samples. Saving frames more often does not help once the interval is shorter than the correlation time; it only adds disk space.

Block averaging

Summing $C(k)$ directly is noisy at long lags. The block-averaging method of Flyvbjerg and Petersen (1989) avoids it. Divide the series into $M$ blocks of length $b$, average within each block, and treat the block means as samples:

For $b = 1$ this is the naive standard error. As $b$ grows beyond the correlation time, the block means become independent and $\sigma_{\bar A}(b)$ rises to a plateau, the honest error. If no plateau appears before the number of blocks becomes small, the simulation is too short to estimate its own error.

Use the lab to see three lessons.

1. The naive error is wrong by $\sqrt{2\tau_{\text{int}}}$. At $\tau = 40$ steps it is about 9 times too small; the block-averaged plateau matches the exact value.
2. The plateau needs many correlation times. At $\tau = 400$ steps there are only about 10 independent samples in 8192 steps, the plateau is noisy and the error itself is uncertain.
3. Equilibration bias is not a statistical error. With the far-from-equilibrium start, the mean is biased however carefully you compute its error bar. Discarding the first stretch removes the bias. Chodera’s (2016) method chooses the cut-off automatically by maximising the number of effectively independent samples that remain.

The same calculation in NumPy, on a series with $\tau = 50$:

import numpy as np

def block_average(x, min_blocks=8):
    """Standard error of the mean of a correlated series, for block sizes 1, 2, 4, ..."""
    x = np.asarray(x, float)
    out = []
    b = 1
    while len(x) // b >= min_blocks:
        m = len(x) // b
        means = x[: m * b].reshape(m, b).mean(axis=1)
        out.append((b, means.std(ddof=1) / np.sqrt(m)))
        b *= 2
    return out

# a correlated test series: AR(1) with correlation time tau = 50 steps
rng = np.random.default_rng(0)
tau, n = 50, 2 ** 15
phi = np.exp(-1 / tau)
x = np.empty(n); x[0] = rng.normal()
for i in range(1, n):
    x[i] = phi * x[i - 1] + np.sqrt(1 - phi ** 2) * rng.normal()

for b, se in block_average(x):
    print(f"block {b:5d}: SE = {se:.4f}")
tau_int = (1 + phi) / (2 * (1 - phi))
print(f"naive SE {x.std(ddof=1) / np.sqrt(n):.4f}, exact SE {np.sqrt(2 * tau_int / n):.4f}")

The block error climbs from 0.0056 (naive) to a plateau of about 0.053 at block sizes of 512–1024 steps, against an exact value of 0.055, then becomes noisy as the number of blocks falls. My browser tool Converge does this analysis, together with equilibration detection, for any GROMACS .xvg file.

Structure: the radial distribution function

The radial distribution function $g(r)$ is the density of particles at distance $r$ from a reference particle, relative to the bulk density $\rho$:

In practice: histogram all pair distances (minimum image) in shells of width $\Delta r$, divide each bin by the number expected for an ideal gas at the same density, $\rho\,N\,4\pi r^2\Delta r$, and average over frames. $g(r) \to 1$ at large $r$ in a fluid. Its first peak marks the nearest-neighbour shell, and integrating to the first minimum $r_1$ gives the coordination number

$g(r)$ is measurable: its Fourier transform is the structure factor $S(q)$ from X-ray or neutron scattering, which makes it a direct test of a force field against experiment. Integrals of $g(r) - 1$ over all space, the Kirkwood–Buff integrals

connect structure to thermodynamics: preferential solvation, activity coefficients and partial molar volumes. They converge slowly and need finite-size corrections in a closed simulation box; my tool KBkit implements them. The g(r) view of the MD lab builds the 2D version live.

Dynamics: diffusion

The mean-squared displacement of a particle grows linearly at long times in a fluid, and its slope gives the self-diffusion coefficient through the Einstein relation,

with $d$ the dimension (6Dt in 3D). Two practical rules: use unwrapped coordinates (not folded back into the box), and fit only the linear regime, after the short-time ballistic region. The equivalent Green–Kubo formula integrates the velocity autocorrelation function, $D = \frac1d\int_0^\infty\langle\mathbf v(0)\cdot\mathbf v(t)\rangle\,\mathrm dt$; similar formulas give viscosity and conductivity from stress and current correlations. Diffusion coefficients in periodic boxes are systematically too small because of hydrodynamic self-interaction with the periodic images; Yeh and Hummer’s correction adds $k_BT\xi/(6\pi\eta L)$ with $\xi \approx 2.837$ for a cubic box of side $L$.

import numpy as np

def msd(pos, max_lag):
    """Mean-squared displacement from unwrapped positions of shape (frames, atoms, 3)."""
    return np.array([np.mean(np.sum((pos[lag:] - pos[:-lag]) ** 2, axis=-1)) if lag else 0.0
                     for lag in range(max_lag)])

# test on Brownian particles with a known diffusion coefficient
rng = np.random.default_rng(1)
D, dt, frames, atoms = 0.5, 0.01, 4000, 200           # nm^2/ns, ns
steps = rng.normal(scale=np.sqrt(2 * D * dt), size=(frames, atoms, 3))
pos = np.cumsum(steps, axis=0)
lags = np.arange(400)
m = msd(pos, 400)
slope = np.polyfit(lags[50:] * dt, m[50:], 1)[0]       # fit the linear regime
print(f"D from MSD = {slope / 6:.3f} nm^2/ns (true {D})")

This recovers $D = 0.509$ against a true value of 0.5. Testing analysis code on synthetic data with a known answer, as here, is good practice before trusting it on a real trajectory.

Free energies are not simple averages

Many quantities that matter most, such as binding free energies, solubilities, the free-energy profile of a molecule crossing a membrane or adsorbing on a surface, are not averages of any instantaneous observable. A free-energy difference involves the ratio of partition functions, and states separated by high barriers are rarely visited in a plain simulation. Umbrella sampling adds harmonic bias potentials (windows) along a reaction coordinate so that every region is sampled, and the weighted histogram analysis method (WHAM) or MBAR recombines the biased windows into an unbiased potential of mean force. Window placement and overlap decide whether the result can be trusted; my tool WindowWise plans and checks them. These methods produced the roughly 470,000 free-energy profiles used to train the adhesion model described on the research page.

References and further reading

From fluctuations to transport

Transport coefficients describe how a system responds to gradients: the diffusion coefficient relates a flux of particles to a concentration gradient (Fick’s law), the shear viscosity relates a momentum flux to a velocity gradient, the thermal conductivity relates a heat flux to a temperature gradient, and the electrical conductivity relates a current to an electric field. They can be measured in simulations by imposing a gradient and measuring the flux, but the fluctuation–dissipation theorem offers a more elegant route: in equilibrium, the spontaneous fluctuations of the fluxes already contain the transport coefficients. This chapter shows how, starting with the oldest example, diffusion.

Brownian motion and the mean-squared displacement

In 1827 the botanist Robert Brown observed pollen particles in water moving erratically, a motion he could not explain. In 1905 Albert Einstein showed that such motion is caused by collisions with the molecules of the liquid and that, although each trajectory is random, its statistics are precise: the mean-squared displacement grows linearly in time, with a slope set by the diffusion coefficient, and the diffusion coefficient is related to the friction the particle feels, $D = k_BT/\zeta$. Jean Perrin’s measurements of this growth in 1908–1909 confirmed Einstein’s theory and gave one of the first accurate values of Avogadro’s number, settling the long debate about the reality of atoms.

For a particle in a simulation the same law holds once its motion has become diffusive:

At short times the motion is ballistic, $\text{MSD} \approx \langle v^2\rangle t^2$, because a particle moves freely between collisions. In a dense liquid there is an intermediate regime where each particle rattles in a cage formed by its neighbours before escaping. Only at times long compared with these processes does the MSD become linear, and the diffusion coefficient is read from the slope in that regime, typically on a log–log plot to check that the exponent has reached 1.

Time-correlation functions and Green–Kubo

The displacement is the time integral of the velocity, so the MSD can be rewritten in terms of the velocity autocorrelation function (VACF), $C_v(t) = \langle \mathbf v_i(0)\cdot\mathbf v_i(t)\rangle$, the average product of a particle’s velocity now with its velocity a time $t$ later. Differentiating the Einstein relation twice gives the Green–Kubo formula for diffusion:

Melville Green (1952, 1954) and Ryogo Kubo (1957) showed that every transport coefficient has such a representation: the integral over time of the equilibrium autocorrelation function of the corresponding microscopic flux. The general pattern is

with $J$ the particle velocity for diffusion, an off-diagonal component of the pressure tensor for shear viscosity, $\eta = \frac{V}{k_BT}\int_0^\infty\langle P_{xy}(0)P_{xy}(t)\rangle\mathrm dt$, the heat current for thermal conductivity, and the electric current for electrical conductivity.

The figure runs a two-dimensional Lennard-Jones fluid and computes both estimates of the diffusion coefficient from the same trajectory, averaging over particles and over many time origins.

At low density the VACF decays smoothly and monotonically: a particle keeps moving in roughly the same direction until a collision randomises it. At high liquid densities it develops a negative dip: a particle moving into its cage of neighbours rebounds, so its velocity a short time later is anticorrelated with its initial velocity. This cage effect reduces the diffusion coefficient, because the negative part of the VACF subtracts from its integral. The two estimates agree within statistical error, as they must, but they converge differently: the Einstein estimate needs the MSD to be clearly linear, and the Green–Kubo integral needs the VACF to have decayed to zero, with noise in its tail accumulating in the integral.

Long-time tails and finite-size effects

In 1967–1970 Berni Alder and Thomas Wainwright, simulating hard spheres, discovered that the VACF does not decay exponentially at long times, as kinetic theory had assumed, but as a power law, $t^{-d/2}$. A moving particle sets up a vortex flow in the surrounding fluid that pushes it from behind, a hydrodynamic memory effect. This long-time tail was one of the first major discoveries made by computer simulation. In two dimensions the $t^{-1}$ tail makes the Green–Kubo integral diverge logarithmically, so a strictly defined diffusion coefficient does not exist in two-dimensional fluids, which is worth remembering when interpreting the figure; in three dimensions the $t^{-3/2}$ tail is integrable.

The same hydrodynamic interactions cause a systematic finite-size effect in periodic simulations: a particle interacts with the flows set up by its own periodic images, which lowers its diffusion coefficient. In three dimensions the correction derived by In-Chul Yeh and Gerhard Hummer (2004) is $D_\infty = D_{\text{PBC}} + \frac{\xi\,k_BT}{6\pi\eta L}$, with $\xi \approx 2.837$ for a cubic box of side $L$ and $\eta$ the shear viscosity. For water in a box of a few nanometres the correction is about 10–20 %, comparable to many of the differences between force fields that simulations try to resolve, so it should always be applied or checked by varying the box size.

Practical transport calculations

Several practical points determine whether a transport calculation is reliable. Use thermostats that do not disturb the dynamics, or measure in NVE after equilibration, as the chapter on thermostats explains; strong Langevin friction lowers diffusion coefficients directly. Unwrap trajectories and remove any drift of the centre of mass. Average over many time origins, but remember that overlapping windows are correlated, so estimate errors by splitting the trajectory into independent blocks or running independent replicas. For collective properties such as viscosity and conductivity, which involve a single flux for the whole system rather than one per particle, statistics are much poorer than for self-diffusion, and long runs, several independent simulations or the Einstein–Helfand (integrated) forms of the Green–Kubo relations are needed. Non-equilibrium methods, such as imposing shear flow (Lees–Edwards boundaries, or the reverse non-equilibrium method of Florian Müller-Plathe, 1999), are often more efficient for viscosity and thermal conductivity.

References and further reading