What Is Molecular Dynamics?
Molecular dynamics (MD) numerically integrates Newton's equations of motion for a system of N interacting particles. Given initial positions and velocities, we repeatedly compute forces from a potential energy function, then advance positions and velocities in time via a symplectic integrator. Over millions of timesteps, equilibrium properties emerge from first-principles mechanics — no statistical assumptions needed.
The Lennard-Jones Potential
The Lennard-Jones (LJ) potential is the canonical model for noble-gas interactions (Ar, Ne, Kr) and the benchmark for MD algorithms:
Using reduced units ε = σ = m = k_B = 1 simplifies everything: temperature becomes T* = k_BT/ε, energy is in units of ε, length in σ, time in σ√(m/ε).
r⁻¹² Repulsion
Hard-core: models Pauli exclusion when electron clouds overlap. Steeply diverges as r→0.
−r⁻⁶ Attraction
Long-range van der Waals / London dispersion. The r⁻⁶ arises from dipole-dipole quantum fluctuations.
Well depth ε
Minimum V = −ε at r = 2^(1/6)σ ≈ 1.12σ — the equilibrium separation in a dimer.
Cutoff r_c = 2.5σ
Truncating at 2.5σ neglects only 1.6% of well depth. A long-range tail correction is applied.
Velocity-Verlet Integration
The simulation uses velocity-Verlet, a time-reversible symplectic integrator that conserves the shadow Hamiltonian to O(dt²):
x(t + Δt) = x(t) + Δt · v(t + ½Δt)
→ recompute F(t + Δt)
v(t + Δt) = v(t + ½Δt) + ½Δt · F(t + Δt)/m
Timestep Δt = 0.006 τ (reduced units). At this resolution total energy drift is <0.1% per million steps.
Berendsen Thermostat
To maintain a target temperature T₀ while thermodynamically equilibrating, velocities are rescaled every τ_c steps:
where T = ⟨mv²⟩/(2Nk_B) is the instantaneous kinetic temperature. In this simulation τ_c ≈ 40–60 steps, gentle enough to preserve short-time dynamics while steering toward equilibrium.
Phase Transitions in 2D
| Phase | Reduced T* | Reduced ρ* | Characteristics |
|---|---|---|---|
| Gas | > 1.2 | < 0.3 | Random motion, rare collisions, no structure |
| Liquid | 0.5 – 1.2 | 0.4 – 0.7 | Short-range order, diffusive dynamics, clustering |
| Solid (crystal) | < 0.4 | > 0.8 | Hexagonal close-packed lattice, phonon vibrations |
| Supercritical | > 1.3, high ρ | > 0.3 | No distinct liquid/gas boundary; single dense fluid phase |
In 2D the LJ triple point occurs near T* ≈ 0.4, ρ* ≈ 0.85. The hexatic phase — with quasi-long-range orientational but exponentially decaying positional order — may appear during melting (Kosterlitz-Thouless-KTHNY theory).
Periodic Boundary Conditions & Minimum Image
Particles interacting across the box boundary are handled by the minimum image convention: for each particle pair, choose the image of particle j (among its periodic copies) that is nearest to particle i:
This eliminates artificial surface effects and models bulk behaviour with only N ≈ 64–100 particles.
Maxwell-Boltzmann Distribution
Particle colours reflect speed: blue = slow, green = medium, red = fast. In equilibrium the speed distribution obeys the 2D Maxwell-Boltzmann distribution:
Most-probable speed v_p = √(k_BT/m). You can verify this by enabling the velocity overlay and noticing the distribution at different temperatures.
Summary Table
| Property | Symbol | Value / Method |
|---|---|---|
| Potential | V(r) | Lennard-Jones, cutoff 2.5σ |
| Integrator | — | Velocity-Verlet, Δt = 0.006 |
| Thermostat | τ_c | Berendsen, τ_c = 40–60 steps |
| Boundary | PBC | Periodic + minimum image |
| Force loop | O(N²) | N = 64–100, ~5–12 ms per frame |
| Units | σ,ε,m | Reduced LJ units |
| Temperature | T* | Instantaneous kinetic: T = ⟨mv²⟩/2Nk_B |
| Phase displayed | — | Inferred from KE/N threshold |
Frequently Asked Questions
Explore More Physics Simulations
Continue your exploration of emergent behaviour and statistical mechanics.