The Lennard-Jones potential is a mathematical model for the interaction energy between two neutral atoms or molecules as a function of their separation distance r. It has the form V(r) = 4ε[(σ/r)¹² - (σ/r)⁶], where ε is the depth of the potential well (binding energy) and σ is the finite distance at which the inter-particle potential is zero. The r⁻¹² repulsion term models the Pauli exclusion principle (electron clouds cannot overlap), while the r⁻⁶ attraction term models van der Waals (London dispersion) forces.
The Lennard-Jones potential has a minimum at r_min = 2^(1/6)σ, where the equilibrium separation occurs and the force is zero. At smaller distances the repulsion dominates and pushes particles apart; at larger distances the weak attraction pulls them together. By simulating many particles interacting via this pairwise potential using molecular dynamics (Newton's equations integrated numerically by the Verlet or leapfrog algorithm), one can observe the spontaneous formation of solid, liquid, and gaseous phases and transitions between them.
Molecular dynamics simulations using the Lennard-Jones potential have reproduced the equation of state of noble gases (argon, krypton) with high accuracy, validated the Lennard-Jones parameter set for rare gas atoms. The model correctly captures the liquid-gas critical point, triple point (solid-liquid-gas coexistence), diffusion coefficients, viscosity, and structural features like pair correlation functions. Despite its simplicity, it serves as the benchmark potential for testing new simulation algorithms, thermostat methods, and sampling techniques in computational chemistry and materials science.
The attractive r⁻⁶ term models London dispersion forces (van der Waals attraction): temporary fluctuating dipoles in one atom induce a dipole in a neighbouring atom, creating a weak but universal attractive force. The repulsive r⁻¹² term is a numerical approximation to the much stronger Pauli exclusion repulsion that prevents electron clouds from overlapping; the r⁻¹² form is computationally convenient but physically less rigorous than exponential repulsion.
The equilibrium separation (where force is zero and potential is minimum) is r_min = 2^(1/6)×sigma ≈ 1.122 sigma. At this distance, the energy is -epsilon (the well depth). For argon, sigma ≈ 3.4 Angstroms and epsilon/k_B ≈ 120 K, so the equilibrium separation is about 3.82 Angstroms and the binding energy per pair corresponds to about 120 K (very weakly bound, consistent with argon's low boiling point of 87 K).
Molecular dynamics (MD) integrates Newton's equations of motion numerically: given positions and velocities, forces are computed, and positions are updated step by step. It follows the true dynamical trajectory and gives time-dependent information (diffusion, velocity autocorrelations, spectra). Monte Carlo (MC) generates equilibrium configurations by proposing and accepting/rejecting random moves based on the Boltzmann factor; it samples equilibrium properties efficiently but contains no physical dynamics or time information.
Depending on temperature (relative to epsilon/k_B) and density, a Lennard-Jones system exhibits gas (low density, particles freely moving), liquid (intermediate density, short-range order, diffusive motion), and solid (close-packed FCC crystal, long-range order). The gas-liquid critical point, the triple point (where solid, liquid, and gas coexist), and melting/solidification can all be observed by varying temperature and pressure in the simulation.
The pair correlation function g(r) measures how particle density varies as a function of distance from a reference particle, relative to the average density. In a gas, g(r) ≈ 1 everywhere (no structure). In a liquid, g(r) shows peaks at preferred separations (first, second coordination shells) that decay to 1 at large r. In a crystal, sharp peaks persist at large r reflecting the regular lattice. The positions and heights of peaks in g(r) reveal bond lengths, coordination numbers, and structural order in the material.