Skip to content

The Lennard-Jones pair force, live

The force between two particles along the line between them, positive apart, minus the derivative of the 12-6 potential V = 4ε ((σ/r)¹² − (σ/r)⁶): F = 24ε (2 (σ/r)¹² − (σ/r)⁶) / r, repulsive inside r = 2^(1/6) σ, where the potential is deepest, attractive beyond it, and fading to nothing far away. Below, two argon atoms (ε = 1.65 × 10⁻²¹ J, σ = 3.4 Å), as their separation and its rate with the reduced mass, are sent toward each other from 1.5 σ at 200 m/s and stepped by the library's velocity Verlet (lusine.integrator.kick_drift, the force lusine.laws.pair_force over the mass, lusine.integrator.kick), two femtoseconds a step and two hundred steps a second, on your GPU as you read. They meet, exchange kinetic and potential energy, and part; the same program adds the two energies up at every step and checks that their sum stays where it started.

This page runs the program on your GPU through WebGPU, so it needs a browser with WebGPU (a recent Chrome, Edge, Safari or Firefox). Without it the program does not run and nothing is drawn below.

What it shows

  • The strip above the controls is green while every check holds, and turns red for good once one fails.
  • The pair: the first chart draws the kinetic energy of the separation's rate, the potential and their sum, in units of the well's depth ε, over the last 2000 steps: as the atoms meet the two trade places and the sum stays flat. The second draws their separation in units of σ. Every 2000 steps (four picoseconds) the pair is sent in again. Below them, the separation, the closest the atoms have come, and the force between them.
  • The checks: the chart draws the energy's drift from its start at each step as a share of what keeps its energy allows, 0.1% of the start, on a log scale, so the check holds while it stays below 1. Below it, the force at the potential's minimum, 2^(1/6) σ; zero at the minimum, true while that force has stayed within binary32's cancellation there of nothing; repels, attracts, fades, while the force has been apart just inside the minimum, together just beyond it, and less than a ten-millionth at 20 σ of what it is at 1.2 σ.

Every number here is computed by the program on the GPU; the page only reads the values back and draws them. Argon's parameters, the reduced mass, the start, the step and the run's length are the program's settings, which the host may write; this page offers no control to change them yet.

The program

The pair, its energy and the checks, the panels the page draws, and the one procedure, main, run each step: one velocity Verlet step, the pair sent in again where a run starts, and the checks after it.

// The Lennard-Jones pair force's live page: two argon atoms, as their separation and its rate with the reduced mass,
// sent toward each other and stepped by the library's velocity Verlet under the pair force, their energy kept as they
// meet and part, with the page's checks computed here on the GPU and shown by the panels the program declares.

namespace pair_force;

lusine.version version = { { number = 1 } };

// The reader's settings: argon's well depth and size, the pair's reduced mass (half an argon atom's), where they start
// and how fast they close, and the step (two femtoseconds), and the steps a run lasts.
host struct {
  float<J> epsilon = 0.00000000000000000000165 J; float<m> sigma = 0.00000000034 m;
  float<kg> mass = 0.00000000000000000000000003315 kg; float<m> start = 0.00000000051 m;
  float<m/s> closing = -200.0 m/s;
  float<s> dt = 0.000000000000002 s; integer steps = 2000;
} settings;

// The pair: their separation, its rate, its acceleration and the force between them; the steps since they were sent.
struct { float<m> r; float<m/s> velocity; float<m/s2> acceleration; float<N> force; integer step; } pair;

// The pair's energy, the kinetic energy of the separation's rate and the potential V = 4 eps ((s/r)^12 - (s/r)^6),
// each in units of the well's depth, against what it started with; the separation in units of sigma, and the closest
// it has come; and the force at the potential's minimum, 2^(1/6) sigma, inside it and beyond it, near and far.
struct {
  float kinetic; float potential; float total; float start; float drift_share; float separation; float closest = 100.0;
  float<N> at_minimum; float<N> inside; float<N> beyond; float<N> near; float<N> far;
} energy;

// The checks: the energy within 0.1% of its start at every step; the force zero at the potential's minimum to
// binary32's cancellation there; and repulsive inside it, attractive beyond it, and nothing far away.
struct { bool keeps_energy = true; bool minimum = true; bool sides = true; } checks;

// A strip of lamps the page draws, lit green while every check holds and red once one fails.
struct lamp { vector<float, 4> colour; };
lattice<lamp, 2> lamps(spacing = 1 m, boundary = wall);

vector<float, 4> lamp_colour(bool holds)
{
  return holds ? vector<float, 4>(0.15, 0.6, 0.3, 1.0) : vector<float, 4>(0.75, 0.15, 0.15, 1.0);
}

// Newton's second law: the force over the mass.
float<m/s2> newton(float<N> force, float<kg> mass) { return force / mass; }

// The Lennard-Jones potential, in units of its well's depth: 4 ((s/r)^12 - (s/r)^6).
float potential(float<m> r, float<m> sigma)
{
  let s = sigma / r;
  let s6 = s * s * s * s * s * s;
  return 4.0 * (s6 * s6 - s6);
}

float kinetic(float<kg> mass, float<m/s> velocity, float<J> epsilon)
{
  return 0.5 * mass * velocity * velocity / epsilon;
}

output table<lusine.widget> panels = {
  lusine.series { { panel = "the pair", chart = "energy (in eps)", quantity = energy.kinetic, steps = 2000,
                    label = "kinetic" } },
  lusine.series { { panel = "the pair", chart = "energy (in eps)", quantity = energy.potential, steps = 2000,
                    label = "potential" } },
  lusine.series { { panel = "the pair", chart = "energy (in eps)", quantity = energy.total, steps = 2000,
                    label = "total" } },
  lusine.series { { panel = "the pair", chart = "separation (in sigma)", quantity = energy.separation, steps = 2000,
                    label = "r / sigma" } },
  lusine.value { { panel = "the pair", quantity = energy.separation, label = "separation (in sigma)" } },
  lusine.value { { panel = "the pair", quantity = energy.closest, label = "closest (in sigma)" } },
  lusine.value { { panel = "the pair", quantity = pair.force, label = "force" } },
  lusine.series { { panel = "the checks", chart = "share of the tolerance", quantity = energy.drift_share, steps = 2000,
                    log = true, label = "energy drift" } },
  lusine.value { { panel = "the checks", quantity = energy.at_minimum, label = "force at 2^(1/6) sigma" } },
  lusine.value { { panel = "the checks", quantity = checks.keeps_energy, label = "keeps its energy" } },
  lusine.value { { panel = "the checks", quantity = checks.minimum, label = "zero at the minimum" } },
  lusine.value { { panel = "the checks", quantity = checks.sides, label = "repels, attracts, fades" } },
};

void main()
{
  // Every run's steps the pair is sent toward each other again from the start.
  pair.r = pair.step == 0 ? settings.start : pair.r;
  pair.velocity = pair.step == 0 ? settings.closing : pair.velocity;
  pair.force = lusine.laws.pair_force(pair.r, settings.epsilon, settings.sigma);
  pair.acceleration = newton(pair.force, settings.mass);
  energy.start = pair.step == 0 ? kinetic(settings.mass, velocity = settings.closing, settings.epsilon)
                                  + potential(r = settings.start, settings.sigma)
                                : energy.start;
  (pair.r, pair.velocity) = lusine.integrator.kick_drift(position = pair.r, velocity = pair.velocity,
                                                         acceleration = pair.acceleration, dt = settings.dt);
  pair.force = lusine.laws.pair_force(pair.r, settings.epsilon, settings.sigma);
  pair.acceleration = newton(pair.force, settings.mass);
  pair.velocity = lusine.integrator.kick(pair.velocity, pair.acceleration, settings.dt);
  pair.step = pair.step + 1 >= settings.steps ? 0 : pair.step + 1;
  energy.kinetic = kinetic(settings.mass, pair.velocity, settings.epsilon);
  energy.potential = potential(pair.r, settings.sigma);
  energy.total = energy.kinetic + energy.potential;
  energy.drift_share = abs(energy.total - energy.start) / (0.001 * abs(energy.start));
  energy.separation = pair.r / settings.sigma;
  energy.closest = energy.separation < energy.closest ? energy.separation : energy.closest;
  checks.keeps_energy = checks.keeps_energy && energy.drift_share <= 1.0;
  energy.at_minimum = lusine.laws.pair_force(r = 1.122462 * settings.sigma, settings.epsilon, settings.sigma);
  energy.inside = lusine.laws.pair_force(r = 0.95 * 1.122462 * settings.sigma, settings.epsilon, settings.sigma);
  energy.beyond = lusine.laws.pair_force(r = 1.1 * 1.122462 * settings.sigma, settings.epsilon, settings.sigma);
  energy.near = lusine.laws.pair_force(r = 1.2 * settings.sigma, settings.epsilon, settings.sigma);
  energy.far = lusine.laws.pair_force(r = 20.0 * settings.sigma, settings.epsilon, settings.sigma);
  checks.minimum = checks.minimum
                   && abs(energy.at_minimum) <= 8.0 * 0.00000011920929 * 24.0 * settings.epsilon / settings.sigma;
  checks.sides = checks.sides && energy.inside > 0.0 N && energy.beyond < 0.0 N
                 && abs(energy.far) < 0.0000001 * abs(energy.near);
  lamps.colour = lamp_colour(holds = checks.keeps_energy && checks.minimum && checks.sides);
}

The law's other checks (the force against a central difference of the potential, in double precision) are on its proof page.