⚗️ Tutorial · Canvas · Physics
📅 July 2026 ⏱ ≈ 3 hours 🎓 Beginner-Intermediate

Build a Molecular Gas Simulation in One Evening

You don't need a physics engine, a game framework, or WebGL to watch a gas condense into a liquid before your eyes. With plain canvas 2D and about 150 lines of JavaScript, you can build a Lennard-Jones molecular gas simulation from scratch — complete with realistic forces, stable integration, and a live temperature control — in a single evening.

1. Project Setup

A single HTML file with a canvas and a script tag is all we need — no build step, no dependencies: 

<!doctype html>
<html>
<body style="margin:0; background:#0a0a0f;">
<canvas id="c" width="800" height="600"></canvas>
<input id="temp" type="range" min="0" max="3" step="0.05" value="0.6"
       style="position:fixed; bottom:20px; left:20px; width:300px;">
<script src="gas.js"></script>
</body>
</html>

Everything else lives in gas.js, built up in the sections below.

2. Particle Data and Initial Conditions

Store particle state in flat typed arrays rather than an array of objects — this avoids per-particle allocation overhead and keeps memory access contiguous, which matters once you have a few hundred particles updating 60 times a second: 

const N = 300;                 // number of particles
const px = new Float64Array(N), py = new Float64Array(N);
const vx = new Float64Array(N), vy = new Float64Array(N);
const fx = new Float64Array(N), fy = new Float64Array(N);

const W = 800, H = 600;
const SIGMA = 8, EPSILON = 40;  // Lennard-Jones parameters (pixel units)

// Seed particles on a grid with small random jitter and random velocity
let idx = 0;
outer:
for (let row = 0; row < 30 && idx < N; row++) {
  for (let col = 0; col < 30 && idx < N; col++) {
    px[idx] = 40 + col * 22 + (Math.random() - 0.5) * 4;
    py[idx] = 40 + row * 22 + (Math.random() - 0.5) * 4;
    const speed = 40;
    const angle = Math.random() * Math.PI * 2;
    vx[idx] = Math.cos(angle) * speed;
    vy[idx] = Math.sin(angle) * speed;
    idx++;
    if (idx >= N) break outer;
  }
}
Why a grid, not random positions? Placing particles fully at random risks two particles landing almost on top of each other, producing a near-infinite LJ repulsive force that launches a particle off-screen at absurd speed on the very first frame. A jittered grid guarantees a safe minimum spacing.

3. Lennard-Jones Forces

For each pair of particles within range, compute the LJ force (see our companion article on the Lennard-Jones potential for the full derivation) and accumulate it onto both particles (Newton's third law: equal and opposite): 

const CUTOFF = SIGMA * 2.5;
const CUTOFF_SQ = CUTOFF * CUTOFF;

function computeForcePair(i, j) {
  let dx = px[i] - px[j];
  let dy = py[i] - py[j];
  const r2 = dx * dx + dy * dy;
  if (r2 > CUTOFF_SQ || r2 < 1e-6) return;

  const sr2 = (SIGMA * SIGMA) / r2;
  const sr6 = sr2 * sr2 * sr2;
  const sr12 = sr6 * sr6;
  // F(r)/r, so we can multiply directly by (dx, dy) without normalising
  const forceOverR = (24 * EPSILON / r2) * (2 * sr12 - sr6);

  fx[i] += forceOverR * dx;  fy[i] += forceOverR * dy;
  fx[j] -= forceOverR * dx;  fy[j] -= forceOverR * dy;
}

Note the shortcut: since the true force direction is (dx, dy)/r and the LJ magnitude formula already divides by r once, we divide by r² instead of r and skip the square root entirely — Math.sqrt is one of the more expensive operations in a tight inner loop, so avoiding it whenever possible is worth the small extra algebra.

4. Spatial Hashing for Speed

Calling computeForcePair for every pair costs O(N²) — fine for 300 particles (about 45,000 pairs), but it will not scale to a few thousand. Since the LJ force is essentially zero beyond the cutoff radius, we only need to check particles in nearby cells of a uniform grid: 

const CELL = CUTOFF;
const cols = Math.ceil(W / CELL), rows = Math.ceil(H / CELL);

function buildGrid() {
  const grid = new Map();  // key: "col,row" -> array of particle indices
  for (let i = 0; i < N; i++) {
    const col = Math.floor(px[i] / CELL);
    const row = Math.floor(py[i] / CELL);
    const key = col + ',' + row;
    if (!grid.has(key)) grid.set(key, []);
    grid.get(key).push(i);
  }
  return grid;
}

function computeAllForces() {
  fx.fill(0); fy.fill(0);
  const grid = buildGrid();
  for (const [key, cellParticles] of grid) {
    const [col, row] = key.split(',').map(Number);
    // Check this cell and the 8 neighbours (3x3 block)
    for (let dc = -1; dc <= 1; dc++) {
      for (let dr = -1; dr <= 1; dr++) {
        const neighbourKey = (col + dc) + ',' + (row + dr);
        const neighbourParticles = grid.get(neighbourKey);
        if (!neighbourParticles) continue;
        for (const i of cellParticles) {
          for (const j of neighbourParticles) {
            if (j > i) computeForcePair(i, j);  // avoid double-counting
          }
        }
      }
    }
  }
}
Why 3×3, not just the current cell? A particle near the edge of its cell can still be within the cutoff radius of a particle in the adjacent cell. Checking the current cell plus all 8 neighbours guarantees no interacting pair within the cutoff distance is ever missed, as long as CELL ≥ CUTOFF.

5. Velocity Verlet Integration

With forces known, advance every particle with velocity Verlet (see our companion article on molecular dynamics for why this is preferred over simple Euler integration): 

const DT = 0.01;
const MASS = 1;

function step() {
  computeAllForces();
  const ax = new Float64Array(N), ay = new Float64Array(N);
  for (let i = 0; i < N; i++) {
    ax[i] = fx[i] / MASS;
    ay[i] = fy[i] / MASS;
    px[i] += vx[i] * DT + 0.5 * ax[i] * DT * DT;
    py[i] += vy[i] * DT + 0.5 * ay[i] * DT * DT;
  }

  computeAllForces();  // forces at NEW positions
  for (let i = 0; i < N; i++) {
    const axNew = fx[i] / MASS, ayNew = fy[i] / MASS;
    vx[i] += 0.5 * (ax[i] + axNew) * DT;
    vy[i] += 0.5 * (ay[i] + ayNew) * DT;
  }

  applyWalls();
}
Why two force computations per step? Velocity Verlet needs the acceleration at both the old position (to update position) and the new position (to update velocity). This doubles the cost of the force loop compared to plain Euler, but the dramatically better energy conservation and stability is well worth it — an Euler-integrated LJ gas visibly heats up or explodes within seconds.

6. Walls, Rendering, and a Temperature Slider

Bounce particles elastically off the four walls of the canvas: 

function applyWalls() {
  const R = SIGMA * 0.5;  // visual particle radius
  for (let i = 0; i < N; i++) {
    if (px[i] < R)      { px[i] = R;     vx[i] = Math.abs(vx[i]); }
    if (px[i] > W - R)  { px[i] = W - R; vx[i] = -Math.abs(vx[i]); }
    if (py[i] < R)      { py[i] = R;     vy[i] = Math.abs(vy[i]); }
    if (py[i] > H - R)  { py[i] = H - R; vy[i] = -Math.abs(vy[i]); }
  }
}

Draw each particle, coloured by its speed (a cheap visual proxy for local kinetic temperature): 

const canvas = document.getElementById('c');
const ctx = canvas.getContext('2d');

function render() {
  ctx.fillStyle = '#0a0a0f';
  ctx.fillRect(0, 0, W, H);
  for (let i = 0; i < N; i++) {
    const speed = Math.hypot(vx[i], vy[i]);
    const hue = 220 - Math.min(speed * 2, 200);  // blue (cold) to red (hot)
    ctx.fillStyle = `hsl(${hue}, 90%, 60%)`;
    ctx.beginPath();
    ctx.arc(px[i], py[i], SIGMA * 0.5, 0, Math.PI * 2);
    ctx.fill();
  }
}

Finally, wire up the temperature slider using a simple velocity rescaling thermostat (see our companion Metropolis/statistical mechanics article for more rigorous alternatives): 

const tempSlider = document.getElementById('temp');
function applyThermostat(targetKineticPerParticle) {
  let totalKE = 0;
  for (let i = 0; i < N; i++) totalKE += 0.5 * MASS * (vx[i]**2 + vy[i]**2);
  const currentAvg = totalKE / N;
  if (currentAvg < 1e-6) return;
  const scale = Math.sqrt(targetKineticPerParticle / currentAvg);
  for (let i = 0; i < N; i++) { vx[i] *= scale; vy[i] *= scale; }
}

function loop() {
  step();
  applyThermostat(parseFloat(tempSlider.value) * 500);
  render();
  requestAnimationFrame(loop);
}
loop();

That's it — roughly 150 lines total. Slide the temperature down and watch the gas condense into liquid droplets; slide it up and watch the droplets evaporate back into a chaotic gas.

Next Steps

Frequently Asked Questions

Why does the gas sometimes explode instead of condensing?

Usually this means DT is too large for the given SIGMA/EPSILON combination, or two particles started too close together, producing a huge initial repulsive force. Try reducing DT, increasing the initial grid spacing, or clamping the maximum force magnitude per pair as a safety measure while you debug.

Why use a Map for the spatial grid instead of a 2D array?

A Map keyed by "col,row" strings only allocates memory for cells that actually contain particles, which matters when particles cluster in one region and leave large parts of the domain empty — a fixed 2D array would waste memory (and iteration time, if you loop over every cell) on those empty regions. For a small, bounded domain a typed 2D array is a viable and slightly faster alternative; the Map approach generalises better to larger or unbounded domains.

How much faster is the spatial hashing version than the naive O(N²) loop?

For N=300 with a cutoff covering roughly 5% of the domain area, the grid approach checks only a small constant number of neighbours per particle rather than all 299 others — typically a 5-10x speedup at this scale, growing to orders of magnitude as N increases into the thousands, since the naive approach scales as N² while the grid approach scales close to linearly with N (for roughly constant particle density).

Can this same approach be extended to 3D?

Yes — replace the 2D position/velocity/force arrays with 3D equivalents, extend the grid key to "col,row,layer", and check a 3×3×3 block of neighbouring cells instead of 3×3. Rendering becomes more involved (you'll want a simple orthographic or perspective projection, or a library like Three.js), but the physics and integration code are essentially unchanged.

Why does the thermostat rescale ALL velocities by the same factor?

This is the simplest possible thermostat (see our Metropolis algorithm article for more principled alternatives like Nosé-Hoover). Multiplying every velocity by the same scale factor preserves the direction of each particle's motion while adjusting the overall kinetic energy (and hence temperature) toward the target — it is crude and does not reproduce realistic thermal fluctuations, but it is more than adequate for an interactive visual demo.

Frequently Asked Questions

What will I learn in this tutorial?

Step-by-step tutorial: build a 2D Lennard-Jones molecular gas simulation from scratch in canvas JavaScript in one evening — velocity Verlet, spatial hashing, and a live temperature slider.

What topics are covered in this tutorial?

This tutorial covers: Project Setup, Particle Data and Initial Conditions, Lennard-Jones Forces, Spatial Hashing for Speed, Velocity Verlet Integration, Walls, Rendering, and a Temperature Slider.

What prerequisites do I need before starting?

This is a Beginner-Intermediate-level tutorial — no special preparation beyond basic JavaScript is assumed.