Molecular Dynamics Simulation

Molecular dynamics is the in-silico realization of classical statistical mechanics: given a potential-energy surface $U(\mathbf r_1,\dots,\mathbf r_N)$ and initial conditions, one integrates Newton's equations forward in time. The resulting trajectory is a single realization of the classical $N$ -body problem; because that problem is chaotic for $N\gtrsim3$ , individual trajectories are exponentially sensitive to their initial data (Lyapunov instability), yet the statistical distribution they visit is robust. This is exactly the regime where statistical mechanics is meant to apply — we never ask where one trajectory goes, only how often it visits each region of phase space. The conceptual engine is ergodicity. If the trajectory explores the constant-energy surface uniformly with respect to the microcanonical measure, then long time averages of an observable equal its ensemble average, and thermodynamics becomes a sampling problem: run long enough, average, and read off pressure, energy, radial distribution functions and heat capacities. Beyond statics, molecular dynamics has a unique strength — genuine dynamics: the Green–Kubo and Einstein relations express diffusion constants, viscosities and thermal conductivities as time integrals of equilibrium correlation functions, which molecular dynamics computes directly. That side of the story — the fluctuation–dissipation content — is owned by the langevin-equation and fluctuation-dissipation-theorem articles, recalled here by name. The algorithmic heart is the integrator, and its single most important property is not raw local accuracy but geometric fidelity. The velocity-Verlet scheme succeeds because it is symplectic — it exactly conserves a nearby "shadow" Hamiltonian — and time-reversible, so its energy error oscillates within a bounded band forever rather than accumulating. A naively "more accurate" but non-symplectic integrator like explicit Euler suffers secular energy drift and destroys the sampled ensemble over long runs. This contrast is the article's owned thesis, and it is why a single conceptual choice separates a simulation you can run for $10^8$ steps from one that explodes. Bare Newtonian dynamics samples the microcanonical ( $NVE$ ) ensemble; extending it to fixed temperature or pressure is the job of thermostats and barostats, owned by the thermostats-md-ensembles article.

Key takeaways

Molecular dynamics (MD) integrates Newton's equations $m_i\ddot{\mathbf r}_i=-\nabla_i U$ for $N$ interacting particles, producing a deterministic phase-space trajectory.
Time averages along the trajectory estimate equilibrium ensemble averages by the ergodic hypothesis; time-correlation functions yield transport coefficients via Green–Kubo relations.
The standard integrator is velocity Verlet: symplectic, time-reversible, conserving a shadow Hamiltonian, with global error $O(\Delta t^2)$ and bounded (non-secular) energy error — the decisive property for long runs.
A non-symplectic scheme like explicit Euler is locally comparable but suffers secular energy drift, silently biasing the sampled ensemble; geometric fidelity beats per-step accuracy.
Bare Newtonian MD samples the microcanonical ( $NVE$ ) ensemble; temperature is a diagnostic via equipartition $\langle\tfrac12 m v^2\rangle=\tfrac32 k_BT$ , and pressure comes from the virial. Thermostats/barostats extend it to $NVT$ / $NPT$ .
Key idea: integrate deterministic, time-reversible, symplectic dynamics on a fixed potential; trust ergodicity to convert time averages into ensemble averages and time-correlation functions into transport coefficients — and pick a symplectic integrator so the energy error stays bounded forever.

The core insight

Integrate deterministic, time-reversible, symplectic dynamics on a fixed potential; trust ergodicity to convert time averages into ensemble averages, and time-correlation functions into transport coefficients. The symplectic velocity-Verlet integrator conserves a shadow Hamiltonian so its energy error is bounded, whereas non-symplectic schemes drift secularly and corrupt the sampled ensemble.

Where this lives today

Molecular dynamics is the microscope of computational science. From the trajectory one builds the radial distribution function $g(r)$ (whose Fourier transform is the measured structure factor $S(q)$ ), computes transport coefficients from Green–Kubo integrals of the velocity, stress and heat-current autocorrelations (water's self-diffusion $\sim2.3\times10^{-9}\,\mathrm{m^2/s}$ at 298 K), and simulates protein folding, ligand binding and membrane dynamics on microsecond–millisecond scales with empirical force fields (AMBER, CHARMM, OPLS) on GPUs and special-purpose hardware. Constraint algorithms (SHAKE/RATTLE) freeze the stiffest bond vibrations to permit $\Delta t\sim2$ fs, and particle-mesh Ewald handles long-range electrostatics at $O(N\log N)$ . In materials science, MD with embedded-atom or machine-learned potentials resolves dislocation motion, crack propagation, melting and nucleation — atomistic mechanisms invisible to continuum models. Coupled to thermostats and barostats it spans ensembles; coupled to enhanced sampling it reaches rare events.

The physical picture

Think of molecular dynamics as releasing $N$ billiard-like particles on a smooth potential landscape and filming them. Any single frame-to-frame path is unpredictable past a fraction of a nanosecond (chaos), but the histogram of where the particles spend their time converges to the Boltzmann distribution. Statistical mechanics is the claim that this histogram, not the path, is what physics cares about — and that is what makes deterministic chaos a virtue rather than a bug. The Verlet integrator works because it respects geometry, not arithmetic. It does not minimize the local truncation error; it preserves phase-space volume (symplecticity) and time-reversal symmetry. As a result the discrete trajectory exactly conserves a "shadow Hamiltonian" close to the true one, so the energy wanders within a bounded band forever instead of drifting away. Explicit Euler, though of the same nominal order per step, is neither symplectic nor time-reversible; its energy climbs secularly, and a plot of $E(t)$ versus $t$ for the two integrators is the single clearest picture of why symplecticity matters. Two observables are read directly off the trajectory. Temperature is not an input but a diagnostic: equipartition $\langle\tfrac12 m_i v_i^2\rangle=\tfrac32 k_BT$ lets you infer $T$ from the instantaneous kinetic energy (fixing its mean and fluctuations is what a thermostat does). Pressure comes from the virial, $PV=Nk_BT+\tfrac13\langle\sum_{i<j}\mathbf r_{ij}\cdot\mathbf F_{ij}\rangle$ , a mechanical route that needs only forces. And transport is fluctuation: the same spontaneous velocity correlations that decay in equilibrium, integrated over time, give the diffusion constant — a direct microscopic Green–Kubo measurement, which the langevin-equation and fluctuation-dissipation articles ground.

Figure

Energy E(t): symplectic velocity-Verlet stays bounded, explicit Euler drifts secularly — the punchline

Figure

Velocity-Verlet as symmetric kick–drift–kick: time-reversible and symplectic

The core equations

The dynamics are Newton's second law on the potential-energy surface; everything else follows from integrating it.

m_i\,\ddot{\mathbf r}_i = \mathbf F_i = -\nabla_{\mathbf r_i} U(\mathbf r_1,\dots,\mathbf r_N) \qquad (1)

where: $m_i$ — mass of particle $i$ ; $\mathbf r_i$ — position; $U$ — total potential (force field); $\mathbf F_i$ — net force on $i$ . The workhorse integrator is velocity Verlet, a kick–drift–kick scheme that is symplectic and time-reversible. It is the discrete map whose long-time behaviour this article is about.

\begin{aligned}\mathbf r_i(t+\Delta t)&=\mathbf r_i(t)+\mathbf v_i(t)\Delta t+\tfrac12\mathbf a_i(t)\Delta t^2\\ \mathbf v_i(t+\Delta t)&=\mathbf v_i(t)+\tfrac12\big[\mathbf a_i(t)+\mathbf a_i(t+\Delta t)\big]\Delta t\end{aligned} \qquad (2)

where: $\mathbf v_i$ — velocity; $\mathbf a_i=\mathbf F_i/m_i$ — acceleration; $\Delta t$ — timestep. Temperature is inferred from the kinetic energy by equipartition — a diagnostic, not an input. The canonical ensemble that fixes $T$ is owned by the SM2 canonical-ensemble article and by thermostats-md-ensembles.

T = \frac{2}{(3N-N_c)k_B}\left\langle \sum_{i=1}^{N}\tfrac12 m_i v_i^2\right\rangle \qquad (3)

where: $k_B$ — Boltzmann constant; $N_c$ — number of constraints (conserved momenta/rigid bonds subtracted). Pressure comes from the virial — a mechanical route needing only forces, which is why it is natural in molecular dynamics.

P = \frac{N k_B T}{V} + \frac{1}{3V}\left\langle \sum_{i<j}\mathbf r_{ij}\cdot\mathbf F_{ij}\right\rangle \qquad (4)

where: $V$ — volume; $\mathbf r_{ij}\cdot\mathbf F_{ij}$ — pair virial. The ergodic hypothesis identifies the time average with the ensemble average — the bridge that makes the trajectory a thermodynamic sampler.

\bar A = \frac{1}{M}\sum_{n=1}^{M} A\big(\Gamma(n\Delta t)\big)\;\xrightarrow[M\to\infty]{}\;\langle A\rangle_{\rm ens} \qquad (5)

where: $\Gamma$ — phase-space point $\{\mathbf r_i,\mathbf p_i\}$ ; $M$ — number of sampled steps. The discrete map does not conserve $H$ exactly, but it conserves a nearby shadow Hamiltonian to all orders in $\Delta t$ — the mathematical reason the energy error stays bounded.

H_{\rm shadow} = H + \Delta t^2\,H_2 + \Delta t^4\,H_4 + \cdots \qquad (6)

where: $H_{\rm shadow}$ — conserved quantity of the symplectic map; $H_2,H_4$ — Baker–Campbell–Hausdorff corrections.

Derivation AVelocity Verlet from symmetric Trotter splitting

The purpose of this derivation is to expose why velocity Verlet is symplectic and time-reversible — properties that are invisible if one simply writes down the update rules. Starting from the exact Liouville propagator and applying a symmetric (Strang/Trotter) operator splitting, each factor turns out to be an exactly solvable symplectic map, and the symmetric ordering makes the whole step invariant under time reversal. These two structural facts are what guarantee a conserved shadow Hamiltonian and hence bounded energy error, the central claim of the article.

Step 1. Write the Liouville operator generating the phase-space flow:

i\mathcal L=\sum_i\left(\mathbf v_i\cdot\frac{\partial}{\partial\mathbf r_i}+\frac{\mathbf F_i}{m_i}\cdot\frac{\partial}{\partial\mathbf v_i}\right),

whose exact propagator is $e^{i\mathcal L t}$ . Why: $e^{i\mathcal L\Delta t}$ advances any phase-space function by one exact step; the integrator approximates this exponential.

Step 2. Split into a position-streaming part $i\mathcal L_r=\sum_i\mathbf v_i\cdot\partial_{\mathbf r_i}$ (drift) and a velocity-kick part $i\mathcal L_v=\sum_i(\mathbf F_i/m_i)\cdot\partial_{\mathbf v_i}$ (kick). Each sub-propagator is exactly solvable because it moves only one set of variables.

Step 3. Apply the symmetric factorization

e^{i\mathcal L\Delta t}=e^{i\mathcal L_v\Delta t/2}\,e^{i\mathcal L_r\Delta t}\,e^{i\mathcal L_v\Delta t/2}+O(\Delta t^3).

Reading the factors right to left gives: a half-velocity kick $\mathbf v\to\mathbf v+\tfrac12\tfrac{\mathbf F(t)}{m}\Delta t$ ; a full position drift $\mathbf r\to\mathbf r+\mathbf v\Delta t$ ; recompute forces; a second half-kick $\mathbf v\to\mathbf v+\tfrac12\tfrac{\mathbf F(t+\Delta t)}{m}\Delta t$ — exactly the velocity-Verlet equations (2).

Step 4. Read off the structural properties. Each factor is a shear map with unit Jacobian, hence symplectic (phase-space-volume-preserving); the symmetric ordering is invariant under $\Delta t\to-\Delta t$ with momentum reversal, hence time-reversible. Why this guarantees stability: a symplectic map conserves a shadow Hamiltonian (6) to all orders, so the true energy oscillates within a bounded band — no secular drift. Verified: the update reproduces velocity Verlet, and the symplectic/time-reversible energy stability vs Euler's secular drift is demonstrated numerically on the SHO in verify_aq-sm-molecular-dynamics.py.

Go deeper · Derivation B · 1 minVirial expression for the pressure

Derivation BVirial expression for the pressure

This derivation gives molecular dynamics its mechanical route to pressure — computing $P$ from interparticle forces rather than from a thermodynamic volume derivative. It starts from the Clausius virial theorem, separates the external (wall) and internal (pairwise) forces, and uses Newton's third law to collapse the double sum to distinct pairs. The result is the estimator (4) used at every step of a production run, and it reduces correctly to the ideal-gas law when interactions vanish.

Step 1. Virial theorem. For a bounded equilibrium system the average of $\sum_i\mathbf r_i\cdot\mathbf p_i$ is constant, so its time derivative vanishes, giving

\left\langle\sum_i \mathbf r_i\cdot\mathbf F_i^{\rm tot}\right\rangle = -3Nk_BT.

Why: $\tfrac{d}{dt}\sum_i\mathbf r_i\cdot\mathbf p_i=\sum_i(\dot{\mathbf r}_i\cdot\mathbf p_i+\mathbf r_i\cdot\mathbf F_i)$ ; the first term averages to $2\langle K\rangle=3Nk_BT$ by equipartition, and the whole time-derivative averages to zero.

Step 2. Split the forces into wall (external) and internal parts, $\mathbf F_i^{\rm tot}=\mathbf F_i^{\rm wall}+\mathbf F_i^{\rm int}$ . The wall part integrates to $-3PV$ (the confining pressure acting inward over the boundary).

Step 3. Collapse the internal sum with Newton's third law and $\mathbf r_{ij}=\mathbf r_i-\mathbf r_j$ :

\sum_i\mathbf r_i\cdot\mathbf F_i^{\rm int}=\sum_{i<j}\mathbf r_{ij}\cdot\mathbf F_{ij},

the pairwise virial. Why pairs: each action–reaction pair contributes once through the relative coordinate.

Step 4. Assemble: $-3PV-3Nk_BT+\langle\sum_{i<j}\mathbf r_{ij}\cdot\mathbf F_{ij}\rangle=0$ , i.e.

PV = Nk_BT + \tfrac13\left\langle\sum_{i<j}\mathbf r_{ij}\cdot\mathbf F_{ij}\right\rangle.

The first term is the ideal-gas pressure; the second is the interaction correction, which vanishes for the ideal gas ( $PV=Nk_BT$ ). Verified: the virial pressure reduces to $PV=Nk_BT$ for the ideal gas and the pairwise-sum collapse is checked in verify_aq-sm-molecular-dynamics.py.

Limiting cases

Case	Result	Meaning
Timestep stability	$\Delta t\lesssim\tfrac1{20}\tau_{\min},\ \tau_{\min}=2\pi/\omega_{\max}$	The step must resolve the fastest vibrational period (C–H stretch $\sim10$ fs → $\Delta t\sim0.5$ – $1$ fs; constraints allow larger).
Classical limit	$\lambda_{\rm dB}=h/\sqrt{2\pi m k_BT}\ll a$	Classical MD needs the thermal de Broglie wavelength small versus spacing $a$ ; fails for light atoms (H, He) at low $T$ — use path-integral MD.
Finite size	$L\gg\xi,\ L\gg$ (force range)	The box must exceed the correlation length and force cutoff; periodic boundaries suppress but do not remove finite-size artifacts.
Ergodicity / sampling	$t_{\rm sim}\gg\tau_{\rm relax}$	Time average → ensemble average only if the run far exceeds the slowest relaxation time; glassy/folding systems break this.
Non-symplectic integrator	$\Delta E\sim t$ (secular)	Explicit Euler or RK4 drift in energy over long runs, biasing the ensemble — the failure velocity Verlet avoids.

Worked examples

Example 1 — One velocity-Verlet step for the harmonic oscillator. Take the SHO $\ddot x=-\omega^2 x$ (unit mass) with $\omega=1$ , $x_0=1$ , $v_0=0$ , $\Delta t=0.1$ . Compute $x_1,v_1$ explicitly and the energy after one step.

Step 1. Acceleration at $t=0$ : $a_0=-\omega^2 x_0=-1$ . Position update: $x_1=x_0+v_0\Delta t+\tfrac12 a_0\Delta t^2=1+0-\tfrac12(1)(0.01)=\mathbf{0.9950}$ .

Step 2. Acceleration at $t=\Delta t$ : $a_1=-\omega^2 x_1=-0.9950$ . Velocity update: $v_1=v_0+\tfrac12(a_0+a_1)\Delta t=0+\tfrac12(-1-0.9950)(0.1)=\mathbf{-0.099750}$ .

Step 3. Energy: $E_1=\tfrac12 v_1^2+\tfrac12\omega^2 x_1^2=\tfrac12(0.099750)^2+\tfrac12(0.9950)^2=0.004975+0.495013=\mathbf{0.49999}$ , versus $E_0=0.5$ . The energy changed by only $\sim10^{-5}$ in one step and, crucially, will oscillate around $0.5$ rather than drift — the symplectic guarantee.

Example 2 — Symplectic versus non-symplectic energy drift over a long run. Integrate the same SHO for $10^4$ steps with velocity Verlet and with explicit Euler, and compare the energy behaviour.

Step 1. Velocity Verlet. The energy oscillates within a band of relative width $O(\Delta t^2)=O(0.01)$ and shows no net drift over $10^4$ steps — the mean energy stays at $E_0$ . This is the shadow-Hamiltonian conservation in action.

Step 2. Explicit Euler ( $x_{n+1}=x_n+v_n\Delta t$ , $v_{n+1}=v_n-\omega^2 x_n\Delta t$ ). Each step multiplies the phase-space area by $(1+\omega^2\Delta t^2)>1$ , so the energy grows geometrically: after $10^4$ steps of $\Delta t=0.1$ the energy has blown up by a factor $\approx(1.01)^{10^4/2}$ , astronomically large. Same nominal order, opposite fate: the non-symplectic scheme's secular drift makes the ensemble meaningless, while the symplectic one is stable indefinitely.

Origins · 1 minloup verlet and how this was discovered

loup verlet

1931–2019 · French physicist whose 1967 study of the Lennard-Jones fluid introduced the position "Verlet" integration scheme and the Verlet neighbour list — two of the most-used tools in all of molecular simulation. His paper computed the radial distribution function and thermodynamics of argon, establishing molecular dynamics as a quantitative method. Later in life he turned to the history and philosophy of science. The velocity-Verlet and leapfrog variants that dominate modern codes are direct descendants of his algorithm.

History

1957

— Alder and Wainwright run the first molecular dynamics of a hard-sphere fluid

1964

— Rahman performs the first MD of a realistic (Lennard-Jones argon) liquid

1967

— Verlet introduces the Verlet integrator and neighbour list for the Lennard-Jones fluid

1970

— Alder and Wainwright observe the long-time tail of the velocity autocorrelation function

1977

— Ryckaert, Ciccotti and Berendsen introduce the SHAKE bond-constraint algorithm

1982

— Swope, Andersen, Berens and Wilson formulate the velocity-Verlet scheme (Andersen's constant-pressure MD is 1980)

Going further

References. Pathria & Beale, Statistical Mechanics, 3rd ed., ch. on computer simulation · Frenkel & Smit, Understanding Molecular Simulation, 2nd ed. · Allen & Tildesley, Computer Simulation of Liquids, 2nd ed. · Sethna, Statistical Mechanics: Entropy, Order Parameters, and Complexity · Tuckerman, Statistical Mechanics: Theory and Molecular Simulation · Verlet, Phys. Rev. 159, 98 (1967).

Provenance. A verifier should recompute: (i) the LJ force by differentiating $u(r)$ ; (ii) the velocity-Verlet update from the symmetric Trotter factorization of $e^{i\mathcal L\Delta t}$ and its symplectic/time-reversible energy stability versus explicit Euler's secular drift (this article's verify script); (iii) the virial pressure from the Clausius virial theorem, reducing to $PV=Nk_BT$ for the ideal gas. The water self-diffusion value ( $\sim2.3\times10^{-9}\,\mathrm{m^2/s}$ at 298 K) is a standard experimental number. The Green–Kubo/fluctuation–dissipation content is owned by the langevin-equation and fluctuation-dissipation-theorem articles; the canonical ensemble by SM2 and thermostats-md-ensembles.

Builds on: Hamiltonian and Newtonian mechanics · Microcanonical and canonical ensembles · Kinetic Theory and the Boltzmann Equation · Fluctuation–Dissipation Theorem. Leads to: Thermostats and Ensemble Control in Molecular Dynamics · Monte Carlo Method · Langevin Equation and Stochastic Differential Equations · Green–Kubo and transport coefficients.

Why does chaos not ruin molecular dynamics?

Because MD does not need individual trajectories to be accurate — only their statistics. The classical $N$ -body problem is chaotic (Lyapunov instability), so a computed trajectory diverges exponentially from the "true" one within a fraction of a nanosecond. But statistical mechanics asks only how often phase space is visited, and that distribution is robust to the chaotic scrambling. The shadow-Hamiltonian argument makes this rigorous for symplectic integrators: the computed trajectory is the exact trajectory of a slightly perturbed Hamiltonian, so it samples a slightly perturbed but well-defined ensemble.

What is the shadow Hamiltonian, concretely?

A symplectic integrator does not conserve the true Hamiltonian $H$ , but by the Baker–Campbell–Hausdorff expansion of the split propagator it conserves a nearby $H_{\rm shadow}=H+\Delta t^2 H_2+\Delta t^4 H_4+\cdots$ to all orders. Because $H_{\rm shadow}$ is exactly conserved by the discrete map and differs from $H$ by only $O(\Delta t^2)$ , the true energy $H$ stays within an $O(\Delta t^2)$ band of its initial value forever — it oscillates rather than drifts. This is the precise statement behind "bounded energy error."

How does MD extract a diffusion constant?

Two equivalent routes. The Green–Kubo route integrates the velocity autocorrelation function, $D=\tfrac13\int_0^\infty\langle\mathbf v(0)\cdot\mathbf v(t)\rangle\,dt$ ; the Einstein route fits the long-time slope of the mean-squared displacement, $\langle\Delta r^2\rangle=6Dt$ . They are equivalent because the MSD is the double time-integral of the velocity correlation. The langevin-equation and fluctuation-dissipation-theorem articles own the fluctuation–dissipation content that makes these equilibrium correlations equal to a transport coefficient; MD simply computes the correlation directly from the trajectory.

Why can't I just rescale velocities to control temperature?

Because naive velocity rescaling produces a "perfectly conserved" kinetic energy that is NOT canonical — it suppresses the physical kinetic-energy fluctuations that a real heat bath would impose, and it drains energy into center-of-mass motion (the "flying ice cube" artifact). The correct way to fix temperature is a thermostat that reproduces the right distribution, not just the right mean. That machinery — Nosé–Hoover, Andersen, Langevin thermostats — is owned by the thermostats-md-ensembles article, which this one recalls by name.

When does classical MD break down?

When quantum effects on the nuclei matter. Classical MD treats nuclei as point particles obeying Newton's laws, valid only when the thermal de Broglie wavelength $\lambda_{\rm dB}=h/\sqrt{2\pi m k_BT}$ is small compared to the interparticle spacing. For light atoms (hydrogen, helium) at low temperature, zero-point motion and delocalization are significant, and one must switch to path-integral MD (which maps each quantum particle onto a classical ring polymer) or apply quantum corrections. Electronic degrees of freedom are always treated separately, via the force field or an ab-initio potential.

What is the role of neighbour lists and Ewald summation?

They make MD scale to large systems. Evaluating all pairwise forces naively costs $O(N^2)$ ; for short-range potentials, neighbour/cell lists restrict each particle to nearby partners, giving $O(N)$ scaling. Long-range Coulomb interactions cannot be cut off, so they are handled by Ewald summation — splitting the sum into a rapidly converging real-space part and a reciprocal-space part — with particle-mesh Ewald (PME) achieving $O(N\log N)$ via the fast Fourier transform. These are engineering choices that do not change the physics but determine whether a simulation is feasible.