Simulate a Lennard-Jones Gas From Scratch

1. Why Lennard-Jones

Real atoms and small molecules attract each other weakly at moderate distance (van der Waals / London dispersion forces) and repel ferociously once their electron clouds start to overlap. The Lennard-Jones (LJ) potential, proposed by John Lennard-Jones in 1924, captures both effects with a single simple formula that is cheap to evaluate millions of times per frame — which is exactly why it remains the default test system for teaching, benchmarking, and studying generic phase behaviour (gas, liquid, solid) without committing to any one real chemical species.

2. The LJ Potential and Force

The 12-6 form balances a steep repulsive term (electron overlap, exponent 12) against a gentler attractive term (dispersion, exponent 6):

The minimum of U(r) sits at r = 2^(1/6) σ ≈ 1.122σ, the natural equilibrium spacing between two isolated particles. Below that distance the force is strongly repulsive; above it, weakly attractive, decaying to zero as r → ∞.

function ljForceMagnitude(r, epsilon = 1, sigma = 1) {
  const sr6  = Math.pow(sigma / r, 6);
  const sr12 = sr6 * sr6;
  return 24 * epsilon / r * (2 * sr12 - sr6); // positive = repulsive, negative = attractive
}

3. Reduced (LJ) Units

MD codes almost never simulate in SI units. Instead they set ε = σ = m = 1 (particle mass) and express everything else — time, temperature, pressure — as a dimensionless combination of these. This keeps numbers near unity (numerically well behaved) and makes one simulation applicable to any real substance simply by rescaling.

4. Velocity Verlet Integration

Velocity Verlet is the standard MD integrator: it is symplectic (conserves phase-space volume), time-reversible, and only needs one force evaluation per step despite being second-order accurate — critical when the force computation is the expensive part of every frame.

function velocityVerletStep(particles, forces, dt, computeForces) {
  // Step 1: update positions using old accelerations
  for (const p of particles) {
    p.x += p.vx * dt + 0.5 * p.ax * dt * dt;
    p.y += p.vy * dt + 0.5 * p.ay * dt * dt;
  }
  // Step 2: recompute forces at the new positions
  const oldAx = particles.map(p => p.ax);
  const oldAy = particles.map(p => p.ay);
  computeForces(particles); // sets p.ax, p.ay for the new configuration

  // Step 3: update velocities using the average of old and new accelerations
  particles.forEach((p, i) => {
    p.vx += 0.5 * (oldAx[i] + p.ax) * dt;
    p.vy += 0.5 * (oldAy[i] + p.ay) * dt;
  });
}

5. Periodic Boundaries and Minimum Image

A box of a few hundred particles is dominated by surface effects unless you remove the surface entirely: periodic boundary conditions tile the simulation box infinitely in every direction, so a particle exiting the right edge re-enters on the left. Force calculations then use the minimum image convention — for each pair, use whichever periodic copy of the second particle is closest.

function wrapPosition(p, boxSize) {
  p.x = ((p.x % boxSize) + boxSize) % boxSize;
  p.y = ((p.y % boxSize) + boxSize) % boxSize;
}

function minimumImageDelta(dx, dy, boxSize) {
  // shift delta into (-boxSize/2, boxSize/2] so pairs "see" the nearest periodic copy
  if (dx >  boxSize / 2) dx -= boxSize;
  if (dx <= -boxSize / 2) dx += boxSize;
  if (dy >  boxSize / 2) dy -= boxSize;
  if (dy <= -boxSize / 2) dy += boxSize;
  return { dx, dy };
}

6. Cutoff Radius and Neighbour Lists

Beyond a few multiples of σ, LJ attraction is negligible. Truncating the potential at a cutoff radius — conventionally r_c = 2.5σ — turns an O(N²) all-pairs force loop into one that only needs to check pairs within r_c, and a spatial grid (cell list) makes that check O(N) on average instead of O(N²).

function buildCellList(particles, boxSize, cutoff) {
  const nCells = Math.max(1, Math.floor(boxSize / cutoff));
  const cellSize = boxSize / nCells;
  const cells = Array.from({ length: nCells * nCells }, () => []);
  particles.forEach((p, i) => {
    const cx = Math.floor(p.x / cellSize) % nCells;
    const cy = Math.floor(p.y / cellSize) % nCells;
    cells[cy * nCells + cx].push(i);
  });
  return { cells, nCells, cellSize }; // only check the 3×3 block of neighbouring cells per particle
}

7. Thermostatting by Velocity Rescaling

Left alone, an MD system conserves total energy (microcanonical, NVE) but temperature drifts as kinetic and potential energy trade back and forth. To hold a target temperature (canonical, NVT), the simplest thermostat rescales every velocity by a common factor each step so the instantaneous kinetic temperature matches the target exactly:

function rescaleVelocities(particles, targetT) {
  const n = particles.length;
  const kineticSum = particles.reduce((s, p) => s + p.vx**2 + p.vy**2, 0);
  const currentT = kineticSum / n;
  if (currentT < 1e-9) return;
  const lambda = Math.sqrt(targetT / currentT);
  for (const p of particles) { p.vx *= lambda; p.vy *= lambda; }
}

8. The Full Simulation Loop

Assembling force computation, integration, boundaries, and thermostatting into a single per-frame function:

function computeForces(particles, boxSize, cutoff = 2.5, epsilon = 1, sigma = 1) {
  for (const p of particles) { p.ax = 0; p.ay = 0; }
  for (let i = 0; i < particles.length; i++) {
    for (let j = i + 1; j < particles.length; j++) {
      let dx = particles[j].x - particles[i].x;
      let dy = particles[j].y - particles[i].y;
      ({ dx, dy } = minimumImageDelta(dx, dy, boxSize));
      const r2 = dx * dx + dy * dy;
      if (r2 > cutoff * cutoff || r2 < 1e-6) continue;
      const r = Math.sqrt(r2);
      const fMag = ljForceMagnitude(r, epsilon, sigma) / r; // pre-divide to project onto dx, dy
      particles[i].ax -= fMag * dx; particles[i].ay -= fMag * dy;
      particles[j].ax += fMag * dx; particles[j].ay += fMag * dy;
    }
  }
}

function simulate(particles, boxSize, dt, targetT) {
  velocityVerletStep(particles, null, dt, p => computeForces(p, boxSize));
  for (const p of particles) wrapPosition(p, boxSize);
  if (targetT !== null) rescaleVelocities(particles, targetT);
}

Run this loop at low density and high T* and you get a disordered gas; lower the temperature or raise the density and particles condense into liquid droplets, then a crystalline solid — the same phase behaviour explored qualitatively in the site's van der Waals equation article and quantitatively in this LJ system's phase diagram.

9. Common Pitfalls

• Force divergence at r→0: the r⁻¹³ term in the LJ force blows up if two particles ever overlap exactly — always start from a non-overlapping lattice and use a small enough time step that particles never approach closer than about 0.8σ.
• Forgetting minimum image: computing raw dx = x_j − x_i without wrapping makes particles near opposite box edges feel a huge, wrong long-range force instead of the correct short-range one.
• Time step too large: LJ's steep repulsive core needs a small dt (typically 0.001–0.005 in reduced time units) — energy conservation drifting upward over a run is the classic symptom of too large a step.
• Skipping the potential shift: an unshifted cutoff introduces a small discontinuity in the force at exactly r = r_c, which slowly injects energy over thousands of steps.
• Rescaling every single step: constant hard rescaling suppresses natural temperature fluctuations; rescale only every N steps, or during an initial equilibration phase, then switch it off (or to a gentler thermostat) to sample the correct ensemble.