Molecular Dynamics: The Lennard-Jones Potential Explained
Every molecular dynamics simulation, from a toy argon gas of a few hundred atoms to a protein folding in explicit water with millions of particles, ultimately reduces to one question: given the positions of all atoms, what force does each one feel? For simple, non-bonded atoms the answer is astonishingly compact — a single formula called the Lennard-Jones potential, discovered by John Lennard-Jones in 1924, still powers force fields a century later.
1. The Lennard-Jones Potential
Two neutral atoms separated by a distance r experience two competing effects: a weak, long-range attraction from fluctuating induced dipoles (the van der Waals / London dispersion force), and a strong, short-range repulsion once their electron clouds start to overlap (Pauli exclusion). The Lennard-Jones (LJ) 12-6 potential packages both effects into a single closed-form expression:
ε = depth of the potential well (bond strength, in energy units) σ = distance at which U(r) = 0 (roughly the atomic diameter) r = distance between the two particle centres
The attractive term −(σ/r)⁶ is not arbitrary — it is the actual leading-order result of quantum-mechanical perturbation theory for the induced-dipole/induced-dipole interaction between two polarisable atoms. The repulsive (σ/r)¹² term has no equally rigorous derivation; the twelfth power was chosen largely because it is cheap to compute (it is just the square of the sixth power) and reproduces experimental compressibility data reasonably well.
The potential has a minimum at rmin = 21/6·σ, where U(rmin) = −ε. This is the natural equilibrium spacing between two atoms with no other forces acting — the distance at which attraction and repulsion exactly balance.
2. Deriving the Force from the Potential
Molecular dynamics integrates Newton's second law, so we need the force, not the potential energy. The force between the pair is the negative gradient of U with respect to separation:
F(r) > 0 → repulsive (pushes atoms apart, at short range) F(r) < 0 → attractive (pulls atoms together, at long range)
In three dimensions, the force vector on particle i due to particle j points along the separation vector rij = ri − rj:
For N particles, the total force on particle i is the pairwise sum over every other particle:
This naive sum costs O(N²) pairwise evaluations per step — fine for a few hundred particles in a browser demo, but the reason production codes (LAMMPS, GROMACS) use cell lists and neighbour lists to reduce the cost to roughly O(N) once a cutoff radius is introduced (see Section 4).
3. Velocity Verlet Integration
With forces known, we must advance positions and velocities through time. The workhorse algorithm in molecular dynamics is velocity Verlet, prized for being time-reversible, symplectic (it conserves phase-space volume and, to a very good approximation, total energy over long runs), and cheap — one force evaluation per step:
a⃗(t+Δt) = F⃗(r⃗(t+Δt)) / m (recompute forces at new positions)
v⃗(t+Δt) = v⃗(t) + ½·[a⃗(t) + a⃗(t+Δt)]·Δt
Contrast this with plain Euler integration (v += a·Δt; r += v·Δt), which is not time-reversible and systematically drifts energy upward or downward over thousands of steps — it will visibly heat up or freeze an LJ fluid that should be at constant energy. Velocity Verlet's local error is O(Δt⁴) per step (O(Δt²) global), versus O(Δt²) local for Euler, at essentially the same computational cost.
4. Periodic Boundaries and Cutoffs
A simulation box of a few hundred atoms is dominated by surface effects — most atoms would sit at a wall rather than in a bulk-like environment. Periodic boundary conditions (PBC) solve this: the box is treated as one cell of an infinite tiling, and any particle leaving one face reappears at the opposite face.
Because U(r) decays as r⁻⁶, its contribution beyond a few σ is negligible. Production codes apply a cutoff radius rc (commonly 2.5σ) beyond which the interaction is simply set to zero, combined with the minimum image convention: for each pair, only the nearest periodic copy of particle j is considered when computing rij.
- Cutoff + neighbour lists: a Verlet list caches which particles are within rc + skin distance, rebuilt only every few steps — reducing the force loop to approximately O(N) for a fixed density.
- Cell lists: divide the box into cells of size ≥ rc; a particle only needs to check its own cell and 26 neighbouring cells (in 3D), giving the same asymptotic scaling with a simpler data structure.
- Tail corrections: because the cutoff discards a small attractive contribution, a long-range tail correction is usually added analytically to pressure and energy, assuming a uniform density beyond rc.
Molecular Dynamics Simulation
Watch an LJ gas condense into a liquid droplet or crystallise into a lattice as you tune density and temperature in real time.
5. Thermostats and Temperature Control
Velocity Verlet on its own conserves total energy (the microcanonical, NVE ensemble). Most physical questions — what happens at 300 K? at the melting point? — require holding temperature fixed instead (the canonical, NVT ensemble), which calls for a thermostat.
Instantaneous temperature is defined via the equipartition theorem from the kinetic energy of all N particles in 3D:
- Velocity rescaling: the crudest thermostat — multiply every velocity by √(Ttarget/Tcurrent) each step. Simple but does not sample the canonical ensemble correctly (no realistic fluctuations).
- Berendsen thermostat: rescales velocities gradually toward the target with a relaxation time τ, avoiding the harsh discontinuities of instant rescaling, but still not exactly canonical (it suppresses temperature fluctuations too strongly — the "flying ice cube" problem in extreme cases).
- Nosé-Hoover thermostat: introduces an extra dynamical variable representing a heat-bath degree of freedom, coupled to the physical particles through an extended Lagrangian. This produces a mathematically correct canonical (NVT) distribution and is the standard choice in production codes.
- Langevin thermostat: adds friction and a random (stochastic) force to each particle, mimicking collisions with an implicit solvent or heat bath. Simple to implement and physically intuitive, at the cost of altering the dynamics (it is not momentum-conserving in the way Nosé-Hoover is).
6. Why It Still Matters
Noble Gases and Simple Liquids
The LJ potential remains the textbook model for noble gases (Ar, Kr, Xe) and is the standard benchmark for testing new integrators, thermostats, and parallel MD algorithms, precisely because its phase diagram (triple point, critical point, melting curve) is known to high precision from decades of simulation studies.
Force Fields for Biomolecules
Modern all-atom force fields such as AMBER, CHARMM, and OPLS still use an LJ term for every non-bonded atom pair, combined with Coulomb electrostatics, bond/angle/dihedral terms. The 12-6 form persists mainly for historical and computational reasons — replacing it with a more physically accurate exponential repulsion (Buckingham potential) is more expensive per evaluation and rarely changes results enough to justify the cost at scale.
Materials Science
LJ and its variants model noble-gas solids, some simple metals (via embedded-atom extensions), and are frequently used as a coarse-grained proxy potential in method development for nucleation, fracture, and grain-boundary studies before moving to more expensive quantum or many-body potentials.
Frequently Asked Questions
Why do MD codes use a 12-6 power law instead of a more accurate repulsion?
The r⁻¹² repulsion is computationally convenient (it is simply the square of the r⁻⁶ term, so both can share the same distance calculation), and it reproduces bulk thermodynamic properties of simple fluids well enough for most purposes. More rigorous exponential (Buckingham-type) repulsions exist and are used when accuracy at very short range matters, but they cost extra transcendental function evaluations per pair, which adds up across billions of pair interactions in a large simulation.
What happens if the timestep is too large?
If Δt is too large relative to the fastest vibration in the system, the discrete update overshoots the true trajectory. Energy is no longer conserved — often the simulation visibly "explodes", with particles gaining unphysical velocity each step until the whole system diverges numerically. A rough stability criterion is that Δt should be well under a tenth of the shortest vibrational period present in the system.
How is pressure computed in a Lennard-Jones simulation?
Pressure is computed from the virial theorem, which relates the ensemble-averaged kinetic energy and the sum of pairwise force dotted with separation over all particle pairs: P·V = N·k_B·T + (1/3)·⟨Σᵢ<ⱼ F⃗ᵢⱼ·r⃗ᵢⱼ⟩. The first term is the ideal-gas contribution; the second, the virial term, captures the correction from interparticle forces and can be positive (net repulsion, at high density) or negative (net attraction, at low density near the critical point).
Who was John Lennard-Jones and when did he propose this potential?
Sir John Edward Lennard-Jones was a British mathematician and theoretical chemist who introduced the 12-6 potential in 1924 while studying the equation of state of gases at Bristol, later becoming the first professor of theoretical chemistry at Cambridge. His original work used the potential to fit the second virial coefficient of gases like argon, decades before the first digital MD simulations (which began in the 1950s with hard-sphere models and adopted LJ interactions in the 1960s, notably in the seminal Rahman 1964 simulation of liquid argon).
What is the difference between molecular dynamics and Monte Carlo simulation?
Molecular dynamics integrates Newton's equations of motion, producing a genuine time-ordered trajectory — useful for studying dynamic properties like diffusion coefficients, viscosity, and reaction rates. Monte Carlo methods (see our companion article on the Metropolis algorithm) instead generate a sequence of configurations sampled according to a target probability distribution, without any notion of physical time — well suited to computing equilibrium thermodynamic averages more efficiently when dynamics are not of interest.
Why does an LJ system crystallise at low temperature?
At low temperature the kinetic energy per particle (of order k_BT) is small compared to the depth of the potential well ε. Particles settle near the energy minimum r_min = 2^(1/6)σ of their neighbours, and the densest way to pack spheres at fixed low energy is a close-packed lattice (FCC or HCP). As temperature rises, thermal motion overcomes the well depth, the lattice melts into a liquid, and eventually (at low density) the system becomes a dilute gas — the same triple-point structure seen in real noble gases.
Can the Lennard-Jones potential model chemical bonds?
Not directly. LJ describes non-bonded, van der Waals-type interactions between atoms that are not chemically bonded to each other. Covalent bonds within a molecule are instead modelled with harmonic (or Morse) bond-stretching potentials, angle-bending terms, and dihedral torsion terms; the LJ term is reserved for interactions between atoms in different molecules, or between atoms in the same molecule that are separated by three or more bonds (the so-called 1-4 and beyond non-bonded interactions).
What units are typically used in LJ simulations?
To avoid working with tiny SI numbers, MD codes commonly use "reduced" Lennard-Jones units, where length is measured in units of σ, energy in units of ε, and mass in units of the particle mass m. Time is then naturally measured in units of τ = σ√(m/ε), and temperature in units of ε/k_B. In these units, ε = σ = m = 1, and all quantities become dimensionless numbers of order unity, which is both numerically convenient and makes results directly transferable between different substances by simple rescaling.