Lennard-Jones Potential — The Foundation of Molecular Dynamics

1. The 12-6 Formula

The Lennard-Jones (LJ) potential between two atoms separated by distance r:

The (σ/r)¹² term is repulsive (always positive); the −(σ/r)⁶ term is attractive. Together they produce a potential well with a minimum at r = 2^(1/6)·σ.

2. Physical Basis

Attractive term — r⁻⁶

The r⁻⁶ dependence comes from London dispersion forces (induced dipole-induced dipole), derived by Fritz London in 1930 via quantum perturbation theory. Even non-polar atoms like argon experience these forces because instantaneous quantum fluctuations in electron density create temporary dipoles that induce dipoles in neighbours. The interaction scales as α² (polarisability squared) and decays as r⁻⁶.

Repulsive term — r⁻¹²

The true repulsion from Pauli exclusion (electron cloud overlap) decays approximately exponentially (Buckingham potential: A·e^(−Br)). The r⁻¹² is a convenient mathematical approximation — it has no deep physical justification at this precise power, but it is computationally advantageous because (σ/r)¹² = [(σ/r)⁶]², requiring only one extra multiplication over the r⁻⁶ term already computed.

3. Equilibrium Distance and Well Depth

ε/k_B is reported in Kelvin because ε is frequently compared to the thermal energy k_B·T. For argon, ε/k_B = 119.8 K means the LJ bond energy equals k_B × 119.8 K — argon condenses into a liquid only below this temperature (boiling point ≈ 87 K).

4. Force from the LJ Potential

Force is the negative gradient of potential. In 3D, the force on particle i due to particle j along the vector r_ij:

The force is zero at r = r_eq (equilibrium), repulsive at r < r_eq, and attractive at r > r_eq, becoming negligible for r > ~3σ.

5. Reduced Units

Molecular dynamics simulations conventionally use reduced (dimensionless) LJ units to eliminate ε, σ, and m from the equations. All quantities are expressed as multiples of the LJ parameters:

For argon, the time unit τ = σ√(m/ε) ≈ 2.1 ps. A typical MD timestep is 0.002τ ≈ 4 fs. Reduced units also make results transferable: two LJ fluids at the same T* and ρ* have the same structure and dynamics regardless of their σ and ε values (the principle of corresponding states).

6. Cutoff Radius and Tail Corrections

Computing all N(N-1)/2 pair interactions is O(N²). In practice the LJ potential is truncated at a cutoff radius r_c (typically 2.5σ), where V(r_c) ≈ −0.016ε << room temperature thermal energy:

// LJ with cutoff — fast inner loop
const rc = 2.5 * sigma;  // cutoff
const rc2inv = 1 / (rc * rc);
const Vc = 4 * eps * (Math.pow(rc2inv, 6) - Math.pow(rc2inv, 3));

function ljEnergy(r2) {  // r2 = r², avoids sqrt
  if (r2 > rc * rc) return 0;
  const r2inv = 1 / r2;
  const r6inv = r2inv * r2inv * r2inv;
  return 4 * eps * r6inv * (r6inv - 1) - Vc;  // shifted cutoff
}

function ljForce(r2) {  // returns |F|/r (multiply by r_ij vector)
  if (r2 > rc * rc) return 0;
  const r2inv = 1 / r2;
  const r6inv = r2inv * r2inv * r2inv;
  return 48 * eps * r2inv * r6inv * (r6inv - 0.5);
}

The shifted cutoff (subtracting V(r_c)) ensures the potential goes smoothly to zero at the cutoff without a discontinuous jump in energy — essential for conserving total energy in NVE simulations.

7. Implementation Notes

• Avoid sqrt: Compare r² against r_c² (not r against r_c) — saves a sqrt per pair.
• r6inv trick: Compute r6inv = (1/r²)³, then r12inv = r6inv². LJ force and energy use both — one multiply away from the other.
• Neighbour lists: Verlet neighbour list with skin radius r_s = r_c + 0.5σ reduces force calculations from O(N²) to O(N) for most timesteps.
• Periodic boundary conditions: Minimum image convention — use the nearest image of each particle for pair calculations. Critical for bulk properties.
• Linked cell lists: For large N, divide space into cells of size r_c. Only check pairs in neighbouring cells — typically 27 cells in 3D. O(N) complexity.

8. Beyond LJ — Other Potentials

• Buckingham potential: A·e^(−Br) − C/r⁶ — physically motivated repulsion. More accurate but slower and can blow up at r → 0.
• Morse potential: ε[e^(−2α(r−r_eq)) − 2e^(−α(r−r_eq))] — commonly used for covalent bonds and diatomic molecules.
• Embedded atom method (EAM): Metal potential accounting for embedding energy in electron density. Required for metals: aluminium, copper, gold.
• CHARMM / AMBER / OPLS: Force fields for proteins and biomolecules — combine LJ, Coulomb, bond, angle, dihedral terms for each atom type.

For simulating noble gases (argon, neon, krypton), LJ remains the gold standard — simple and accurate. For polar liquids (water), at minimum a partial-charge Coulomb term must be added (TIP3P, SPC/E). LJ is the foundation everything else builds on.