7 Molecular dynamics: principles

7.1 From local curvature to global exploration

Geometry optimisation finds minima on the potential energy surface. Phonon calculations characterise the curvature at those minima — how the energy changes for small displacements. Both assume atoms remain close to their equilibrium positions: optimisation iterates toward a minimum, and the harmonic approximation expands the energy around it.

At finite temperature, this picture becomes incomplete. Atoms have kinetic energy and explore the potential energy surface. For mildly anharmonic systems, we can extend lattice dynamics with higher-order terms — but this becomes expensive. For systems where atoms don’t remain near any single minimum — diffusing ions, phase transitions, liquids — no expansion around one point will suffice.

Molecular dynamics takes a different approach. Rather than expanding \(E(\mathbf{r})\) and truncating, we follow the actual motion of atoms over the potential energy surface. Given positions and velocities at some instant, we calculate the forces on all atoms, use these to update the velocities and positions a short time later, and repeat. The result is a trajectory: a sequence of configurations showing how the system evolves in time.

7.2 Equations of motion

The force on atom \(i\) is the negative gradient of the potential energy:

\[\mathbf{F}_i = -\nabla_i E(\mathbf{r})\]

Given the force, Newton’s second law gives the acceleration:

\[m_i \frac{\mathrm{d}^2 \mathbf{r}_i}{\mathrm{d}t^2} = \mathbf{F}_i\]

We can split this into two first-order equations:

\[\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t} = \mathbf{a} = \frac{\mathbf{F}}{m}\]

\[\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t} = \mathbf{v}\]

Written in integral form, these become:

\[\mathbf{v}(t) = \mathbf{v}(0) + \int_0^t \mathbf{a}(t') \, \mathrm{d}t'\]

\[\mathbf{r}(t) = \mathbf{r}(0) + \int_0^t \mathbf{v}(t') \, \mathrm{d}t'\]

If we knew how the acceleration varied with time, we could evaluate these integrals and predict the trajectory exactly. But the acceleration depends on the positions, which depend on the velocities, which depend on the accelerations — the problem is self-referential. For any realistic system, we cannot solve these equations analytically.

7.3 Numerical integration

Instead, we integrate numerically by advancing the system in small time steps. The simplest approach is to assume the acceleration is roughly constant over a short interval \(\Delta t\). The velocity integral then becomes straightforward:

\[\mathbf{v}(\Delta t) = \mathbf{v}(0) + \int_0^{\Delta t} \mathbf{a}(0) \, \mathrm{d}t' = \mathbf{v}(0) + \mathbf{a}(0) \Delta t\]

For the position, we need to integrate the velocity — but now the velocity is changing during the interval. Substituting \(\mathbf{v}(t') = \mathbf{v}(0) + \mathbf{a}(0) t'\) into the position integral:

\[\mathbf{r}(\Delta t) = \mathbf{r}(0) + \int_0^{\Delta t} \left[\mathbf{v}(0) + \mathbf{a}(0) t'\right] \mathrm{d}t' = \mathbf{r}(0) + \mathbf{v}(0) \Delta t + \frac{1}{2}\mathbf{a}(0) \Delta t^2\]

These are the familiar kinematic equations for motion under constant acceleration. At each timestep, we calculate the forces, update positions and velocities, then repeat. The approximation becomes exact as \(\Delta t \to 0\), but practical calculations require finite timesteps.

This simple approach has problems. The assumption of constant acceleration introduces errors at each step — small discrepancies between the predicted trajectory and the true one. These errors accumulate over many steps, causing the total energy to drift rather than remaining constant as it should for an isolated system.

There is also a more subtle problem. Newtonian mechanics is symmetric under time reversal: if we reverse all velocities, the system should retrace its path. But using only the acceleration at the start of each timestep breaks this symmetry — a forward step followed by reversing velocities and stepping again does not return exactly to the starting point.

7.4 Velocity Verlet

The velocity Verlet algorithm partially addresses these problems by treating the two halves of each timestep symmetrically.

The idea is to split the velocity update into two half-steps. In the first half, we advance the velocity from \(t\) to \(t + \frac{\Delta t}{2}\) using the acceleration at time \(t\):

\[\mathbf{v}(t + \frac{\Delta t}{2}) = \mathbf{v}(t) + \frac{1}{2}\mathbf{a}(t)\Delta t\]

Then we update positions using this half-step velocity:

\[\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \mathbf{v}(t + \frac{\Delta t}{2}) \Delta t\]

With the new positions, we calculate the forces to get \(\mathbf{a}(t + \Delta t)\). Finally, we complete the velocity update using this new acceleration:

\[\mathbf{v}(t + \Delta t) = \mathbf{v}(t + \frac{\Delta t}{2}) + \frac{1}{2}\mathbf{a}(t + \Delta t)\Delta t\]

The two half-steps are mirror images of each other: the first uses the acceleration at the start, the second uses the acceleration at the end. This symmetry is what makes the algorithm time-reversible.

In practice, these steps can be combined. Substituting the half-step velocity into the position update gives:

\[\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \mathbf{v}(t) \Delta t + \frac{1}{2}\mathbf{a}(t) \Delta t^2\]

And combining the two velocity half-steps:

\[\mathbf{v}(t + \Delta t) = \mathbf{v}(t) + \frac{1}{2}[\mathbf{a}(t) + \mathbf{a}(t + \Delta t)] \Delta t\]

The position equation is unchanged from the simple scheme. The velocity update now uses the average of the old and new accelerations — the trapezoidal rule rather than a rectangle.

Discretisation errors remain — the trajectory is still approximate. But the symmetry makes velocity Verlet time-reversible: run the system forward one step, reverse all velocities, run another step, and you return exactly to where you started. This eliminates systematic energy drift. Energy still fluctuates due to discretisation, but the fluctuations remain bounded rather than accumulating over time. In practice, this means we can run longer simulations without them becoming unphysical.

The timestep must still be small enough to capture the fastest motion in the system — typically the highest-frequency vibrations — which means \(\Delta t\) is usually around 1 femtosecond. A nanosecond trajectory therefore requires \(10^6\) timesteps, each requiring a force evaluation.

7.5 Phase space and trajectories

The state of a classical system is specified by the positions and momenta of all atoms — a point in phase space. As we integrate the equations of motion, this point traces out a trajectory through phase space.

For an isolated system (no external heat bath), the total energy E = K + V is conserved — kinetic energy converts to potential and back, but the sum remains constant. The trajectory is confined to a surface of constant energy in phase space. This is the microcanonical ensemble: all accessible states have the same total energy.

The temperature in such a system is related to the average kinetic energy:

\[\langle K \rangle = \frac{3}{2} N k_B T\]

where N is the number of atoms. This follows from the equipartition theorem: each quadratic degree of freedom contributes ½kT to the average energy, and kinetic energy is quadratic in the velocities.

7.6 Ensembles and what they mean

Real experiments are usually not performed at constant energy. A sample in a furnace exchanges heat with its surroundings; it has a well-defined temperature, not a fixed total energy. Different experimental conditions correspond to different statistical ensembles.

The microcanonical ensemble (NVE) has fixed particle number N, volume V, and energy E. This is what unmodified Newtonian dynamics samples. It describes an isolated system.

The canonical ensemble (NVT) has fixed N, V, and temperature T. The system exchanges energy with a heat bath at temperature T. States are weighted by the Boltzmann factor exp(−E/kT): lower-energy configurations are more probable, but higher-energy configurations are accessible with probability that decreases exponentially.

The isothermal-isobaric ensemble (NPT) has fixed N, pressure P, and temperature T. The system exchanges both energy and volume with its surroundings. This corresponds to most laboratory conditions — a sample at ambient pressure and controlled temperature.

Which ensemble we simulate determines what statistical distribution we sample. This matters for calculating thermodynamic properties and for comparing with experiment.

7.7 Thermostats and barostats

To simulate NVT or NPT, we need to modify the equations of motion so that the trajectory samples the correct statistical distribution. A thermostat couples the system to a heat bath; a barostat couples it to a pressure reservoir.

The details of how thermostats work are technical, but the key point is conceptual: a good thermostat is designed so that, over a long trajectory, configurations are visited with the correct Boltzmann weights. The time-average of any property over the trajectory then equals the ensemble average — the thermodynamic expectation value.

This is the ergodic hypothesis in practice: for a sufficiently long trajectory, the time spent in each region of phase space is proportional to its statistical weight. We can therefore compute ensemble averages — the quantities that connect to thermodynamics and to time-averaged experimental measurements — from a single long trajectory.

7.8 What you get

A molecular dynamics simulation produces a trajectory: atomic positions (and velocities) at each timestep. From this trajectory, we can calculate:

Equilibrium properties — averages over the trajectory. The mean potential energy, the average structure, the distribution of atomic positions. These are ensemble averages, directly comparable to time-averaged experimental measurements.

Dynamic properties — how the system evolves in time. How do atomic positions correlate between different times? How quickly do correlations decay? These time correlation functions contain information about dynamics that equilibrium measurements cannot access.

Transport properties — diffusion coefficients, viscosities, thermal conductivities. These emerge from how quickly the system loses memory of its initial state.

7.9 The \(E(r)\) choice

Every MD timestep requires calculating forces — one evaluation of \(E(r)\) and its gradient. A typical simulation might run for millions of timesteps. The cost of \(E(r)\) therefore determines what simulations are feasible.

Classical potentials are cheap to evaluate. Simulations of millions of atoms for nanoseconds to microseconds are routine. The limitation is accuracy: the functional form constrains what physics can be captured.

Ab initio MD (AIMD) calculates forces from DFT at each timestep. This is accurate but expensive — practical for hundreds of atoms and picoseconds. Many dynamical processes happen on longer timescales or require larger system sizes.

Machine-learned interatomic potentials approach DFT accuracy at near-classical cost. They are increasingly the method of choice when classical potentials aren’t accurate enough but AIMD is too expensive. The caveats discussed in Chapter @ref{mlips} apply: MLIPs inherit limitations from their training data, and can fail outside their training domain.

The choice is problem-dependent. For a quick estimate of dynamics in a well-characterised material, classical potentials may suffice. For quantitative comparison with experiment in a system where accuracy matters, MLIPs or AIMD may be necessary.