005: Jiggling Atoms
2026-10-12
Moving Figures
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
{
if (!navigator.gpu) throw new Error("WebGPU is not supported in this browser.");
// 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 (640x480 buffer, scaled to fill the column)
const canvas = DOM.canvas(640, 480);
canvas.style.width = "100%";
canvas.style.height = "auto";
const adapter = await navigator.gpu.requestAdapter();
const device = await adapter.requestDevice();
const context = canvas.getContext("webgpu");
const format = navigator.gpu.getPreferredCanvasFormat();
context.configure({ device, format });
// Reduced units: sigma = epsilon = mass = k_B = 1. The box is 4:3 so it fills
// the canvas, and a particle (diameter sigma) is about 16 pixels across. At
// this density (about 0.35) the gas and the liquid can coexist.
const Lx = 40;
const Ly = 30;
const count = 420;
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. Drop particles at random, rejecting any that land
// closer than 0.9 sigma (periodically) to one already placed, so no pair
// starts with a huge repulsive force.
const positions = new Float32Array(count * 2);
for (let k = 0; k < count; ) {
const x = Math.random() * Lx;
const y = Math.random() * Ly;
let ok = true;
for (let j = 0; j < k && ok; j++) {
let dx = x - positions[2 * j];
let dy = y - positions[2 * j + 1];
dx -= Lx * Math.round(dx / Lx);
dy -= Ly * Math.round(dy / Ly);
ok = dx * dx + dy * dy > 0.81;
}
if (ok) {
positions[2 * k] = x;
positions[2 * k + 1] = y;
k++;
}
}
// Gaussian random velocities, drift removed, rescaled to the starting temperature
// (in 2D the mean of v^2 per particle is 2T)
const gauss = () =>
Math.sqrt(-2 * Math.log(1 - Math.random())) * Math.cos(2 * Math.PI * Math.random());
const velocities = Float32Array.from({ length: count * 2 }, gauss);
for (let axis = 0; axis < 2; axis++) {
let mean = 0;
for (let i = 0; i < count; i++) mean += velocities[2 * i + axis] / count;
for (let i = 0; i < count; i++) velocities[2 * i + axis] -= mean;
}
const meanSquare = velocities.reduce((s, v) => s + v * v, 0) / count;
const scale = Math.sqrt((2 * temperatureInput.value) / meanSquare);
for (let k = 0; k < velocities.length; k++) velocities[k] *= scale;
// Storage buffers: positions, velocities and the accelerations between passes
const storage = (data) => {
const buffer = device.createBuffer({
size: data.byteLength,
usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST,
});
if (data.length) device.queue.writeBuffer(buffer, 0, data);
return buffer;
};
const positionBuffer = storage(positions);
const velocityBuffer = storage(velocities);
const accelerationBuffer = storage(new Float32Array(count * 2));
// One random-number-generator state (a u32) per particle, seeded from the CPU
const rngBuffer = storage(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,
});
// Vertex buffer: one unit square (two triangles), shared by every particle
const corners = new Float32Array([-1, -1, 1, -1, 1, 1, -1, -1, 1, 1, -1, 1]);
const cornerBuffer = device.createBuffer({
size: corners.byteLength,
usage: GPUBufferUsage.VERTEX | GPUBufferUsage.COPY_DST,
});
device.queue.writeBuffer(cornerBuffer, 0, corners);
// 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;
// PCG hash: advance the state, scramble it into a random u32
fn pcg(state: ptr<function, u32>) -> u32 {
*state = *state * 747796405u + 2891336453u;
let word = ((*state >> ((*state >> 28u) + 4u)) ^ *state) * 277803737u;
return (word >> 22u) ^ word;
}
// Two independent standard normal numbers (Box-Muller)
fn gaussian2(state: ptr<function, u32>) -> vec2f {
let u1 = (f32(pcg(state) >> 8u) + 0.5) / 16777216.0; // in (0, 1)
let u2 = (f32(pcg(state) >> 8u) + 0.5) / 16777216.0;
let radius = sqrt(-2.0 * log(u1));
let angle = 6.2831853 * u2;
return radius * vec2f(cos(angle), sin(angle));
}
@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 shaders: a disc per particle, coloured by its speed
const referenceSpeed = 1.5; // speed at which a particle is fully "hot"
const renderModule = device.createShaderModule({
code: `
const L = vec2f(${Lx}, ${Ly});
@group(0) @binding(0) var<storage, read> positions: array<vec2f>;
@group(0) @binding(1) var<storage, read> velocities: array<vec2f>;
struct Frag {
@builtin(position) position: vec4f,
@location(0) uv: vec2f,
@location(1) @interpolate(flat) colour: vec3f,
}
@vertex
fn vs(
@location(0) corner: vec2f,
@builtin(instance_index) i: u32,
) -> Frag {
// Box coordinates to clip space; the radius is sigma / 2
let centre = positions[i] / L * 2.0 - 1.0;
let radius = 1.0 / L;
let t = clamp(length(velocities[i]) / ${referenceSpeed}, 0.0, 1.0);
let cold = vec3f(0.35, 0.5, 0.8);
let hot = vec3f(1.0, 0.282, 0.204); // #FF4834
var out: Frag;
out.position = vec4f(centre + corner * radius, 0.0, 1.0);
out.uv = corner;
out.colour = mix(cold, hot, t);
return out;
}
@fragment
fn fs(in: Frag) -> @location(0) vec4f {
// Square to disc: drop the corners
if (length(in.uv) > 1.0) { discard; }
return vec4f(in.colour, 1.0);
}
`,
});
const renderPipeline = device.createRenderPipeline({
layout: "auto",
vertex: {
module: renderModule,
entryPoint: "vs",
buffers: [
{
arrayStride: 2 * 4,
attributes: [{ shaderLocation: 0, offset: 0, format: "float32x2" }],
},
],
},
fragment: {
module: renderModule,
entryPoint: "fs",
targets: [{ format }],
},
primitive: { topology: "triangle-list" },
});
const renderBindGroup = device.createBindGroup({
layout: renderPipeline.getBindGroupLayout(0),
entries: [
{ binding: 0, resource: { buffer: positionBuffer } },
{ binding: 1, resource: { buffer: velocityBuffer } },
],
});
// 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() {
if (!running) return;
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 = encoder.beginRenderPass({
colorAttachments: [
{
view: context.getCurrentTexture().createView(),
clearValue: { r: 0.102, g: 0.102, b: 0.102, a: 1.0 }, // #1a1a1a
loadOp: "clear",
storeOp: "store",
},
],
});
pass.setPipeline(renderPipeline);
pass.setBindGroup(0, renderBindGroup);
pass.setVertexBuffer(0, cornerBuffer);
pass.draw(6, count); // 6 vertices per square, `count` instances
pass.end();
device.queue.submit([encoder.finish()]);
requestAnimationFrame(frame);
}
// `invalidation` resolves when this cell re-runs or the page goes away: stop the
// loop and free the device.
let running = true;
invalidation.then(() => {
running = false;
device.destroy();
});
requestAnimationFrame(frame);
yield canvas;
}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