005: Jiggling Atoms
2026-10-12
Building on Compute Shader, the particles now push and pull on each other. Every pair interacts through the Lennard-Jones potential,
V(r) = 4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right],
which is strongly repulsive up close and weakly attractive further out. A first compute pass sums the force on each particle from every other one, and a second integrates the motion (leapfrog, which follows the same trajectory as velocity Verlet). The walls are replaced by periodic boundaries: leave on one side, re-enter on the other.
Condensing lowers the potential energy, and with nothing else to take it, that energy becomes motion, so on its own the box would just warm up. A Langevin thermostat takes it away: every step each particle is slightly slowed and given a small random kick, which holds the temperature at the slider’s value. The kicks need random numbers on the GPU, so each particle carries its own generator state. Particles are coloured by speed. Hot and they stay a gas; cool them and the attraction wins, and droplets of liquid condense out of it.
Show code
{
// Shared WebGPU boilerplate: setup, buffers and the animation loop (lib/gpu.js)
const gpu = await import(new URL("../lib/gpu.js", document.baseURI).href);
const md = await import(new URL("../lib/md.js", document.baseURI).href);
const ui = await import(new URL("../lib/ui.js", document.baseURI).href);
// The temperature is read live from its input element (referencing `viewof x`
// gives the element, not the value), so moving the slider does not restart the
// simulation.
const temperatureInput = viewof temperature;
// Setup, as before, but the canvas is 16:9 and the drawing buffer now follows its
// displayed size (times the screen's pixel ratio) instead of a fixed size
const { canvas, device, context, format } = await gpu.setup({
width: 640,
height: 360,
resolution: "auto",
});
// Reduced units: sigma = epsilon = mass = k_B = 1. The box is 16:9 so it fills
// the canvas, and a particle (diameter sigma) is about 15 pixels across. At
// this density (about 0.35) the gas and the liquid can coexist.
const Lx = 64;
const Ly = 36;
const count = 800;
const dt = 0.005;
const stepsPerFrame = 5;
const friction = 0.5; // thermostat coupling: how fast velocities are refreshed
const decay = Math.exp(-friction * dt); // per-step velocity damping
// Initial state: a gas of random, non-overlapping positions, with Gaussian
// velocities at the starting temperature (see lib/md.js)
const { positions, velocities } = md.randomGas({
count,
box: [Lx, Ly],
temperature: temperatureInput.value,
});
// Storage buffers: positions, velocities and the accelerations between passes
const positionBuffer = gpu.storageBuffer(device, positions);
const velocityBuffer = gpu.storageBuffer(device, velocities);
const accelerationBuffer = gpu.storageBuffer(device, new Float32Array(count * 2));
// One random-number-generator state (a u32) per particle, seeded from the CPU
const rngBuffer = gpu.storageBuffer(device, crypto.getRandomValues(new Uint32Array(count)));
// Uniform buffer: the strength of the thermostat's random kicks, which follows
// the temperature slider
const thermostatBuffer = device.createBuffer({
size: 16, // uniform buffers are happiest in multiples of 16 bytes
usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
});
// Compute shaders. Two entry points, so each pass only touches what it owns:
// `forces` reads every position but writes only its own acceleration, and
// `integrate` reads and writes only its own particle. No thread ever reads
// something another thread is writing in the same dispatch.
const computeModule = device.createShaderModule({
code: `
const L = vec2f(${Lx}, ${Ly}); // box size
const DT = ${dt};
const RC2 = 6.25; // cutoff radius 2.5 sigma, squared
const DECAY = ${decay}; // thermostat: velocity kept each step
struct Thermostat { noise: f32 } // kick size: sqrt(T (1 - DECAY^2))
@group(0) @binding(0) var<storage, read_write> positions: array<vec2f>;
@group(0) @binding(1) var<storage, read_write> velocities: array<vec2f>;
@group(0) @binding(2) var<storage, read_write> accelerations: array<vec2f>;
@group(0) @binding(3) var<storage, read_write> rng: array<u32>;
@group(0) @binding(4) var<uniform> thermostat: Thermostat;
${md.rngWGSL}
@compute @workgroup_size(64)
fn forces(@builtin(global_invocation_id) id: vec3u) {
let i = id.x;
if (i >= arrayLength(&positions)) { return; }
let here = positions[i];
var a = vec2f(0.0);
for (var j = 0u; j < arrayLength(&positions); j++) {
if (j == i) { continue; }
var d = here - positions[j];
// Minimum image: use the nearest periodic copy of particle j
d -= L * round(d / L);
let r2 = dot(d, d);
if (r2 < RC2) {
let s2 = 1.0 / r2; // (sigma / r)^2
let s6 = s2 * s2 * s2;
// force / r = 24 (2 s^12 - s^6) / r^2
a += d * (24.0 * s2 * s6 * (2.0 * s6 - 1.0));
}
}
accelerations[i] = a;
}
@compute @workgroup_size(64)
fn integrate(@builtin(global_invocation_id) id: vec3u) {
let i = id.x;
if (i >= arrayLength(&positions)) { return; }
// Leapfrog: kick the velocity...
var v = velocities[i] + accelerations[i] * DT;
// ...thermostat: slow it a little and add a random kick, which holds the
// mean kinetic energy at the target temperature...
var state = rng[i];
v = v * DECAY + thermostat.noise * gaussian2(&state);
rng[i] = state;
// ...then drift the position
var p = positions[i] + v * DT;
p -= L * floor(p / L); // wrap back into the box
velocities[i] = v;
positions[i] = p;
}
`,
});
// With layout "auto" each pipeline only gets the bindings its entry point uses,
// so each needs its own bind group.
const forcePipeline = device.createComputePipeline({
layout: "auto",
compute: { module: computeModule, entryPoint: "forces" },
});
const integratePipeline = device.createComputePipeline({
layout: "auto",
compute: { module: computeModule, entryPoint: "integrate" },
});
const forceBindGroup = device.createBindGroup({
layout: forcePipeline.getBindGroupLayout(0),
entries: [
{ binding: 0, resource: { buffer: positionBuffer } },
{ binding: 2, resource: { buffer: accelerationBuffer } },
],
});
const integrateBindGroup = device.createBindGroup({
layout: integratePipeline.getBindGroupLayout(0),
entries: [
{ binding: 0, resource: { buffer: positionBuffer } },
{ binding: 1, resource: { buffer: velocityBuffer } },
{ binding: 2, resource: { buffer: accelerationBuffer } },
{ binding: 3, resource: { buffer: rngBuffer } },
{ binding: 4, resource: { buffer: thermostatBuffer } },
],
});
// Render: a disc per particle, coloured by its speed (see lib/md.js)
const discs = md.discs(device, { format, positionBuffer, velocityBuffer, box: [Lx, Ly] });
// One frame: `stepsPerFrame` simulation steps in a single compute pass, then
// draw. Dispatches in a pass run in order, and each sees the previous one's writes.
// The temperature slider is read here every frame.
const workgroups = Math.ceil(count / 64);
function frame() {
const noise = Math.sqrt(temperatureInput.value * (1 - decay * decay));
device.queue.writeBuffer(thermostatBuffer, 0, new Float32Array([noise]));
const encoder = device.createCommandEncoder();
const compute = encoder.beginComputePass();
for (let s = 0; s < stepsPerFrame; s++) {
compute.setPipeline(forcePipeline);
compute.setBindGroup(0, forceBindGroup);
compute.dispatchWorkgroups(workgroups);
compute.setPipeline(integratePipeline);
compute.setBindGroup(0, integrateBindGroup);
compute.dispatchWorkgroups(workgroups);
}
compute.end();
const pass = gpu.beginRenderPass(encoder, context);
discs.draw(pass, count);
pass.end();
device.queue.submit([encoder.finish()]);
}
gpu.run(invalidation, device, frame);
// The canvas with the controls laid over its corner
yield ui.stage(canvas, [temperatureInput]);
}The first of Feynman’s lectures, Atoms in Motion, describes this picture: everything is made of atoms that jiggle, attract each other when a little way apart and repel when squeezed together.
Now the jiggling motion is what we represent as heat: when we increase the temperature, we increase the motion.
— Richard Feynman, The Feynman Lectures on Physics, Vol. I, ch. 1