008: Solid, Liquid, Gas

2026-10-13

Building on Jiggling Atoms and on the measuring tools of Is Energy Conserved?, the box starts as a triangular crystal that has been stretched out to a density of 0.5, far thinner than a crystal wants to be. It cannot hold together as one piece, and soon breaks into crystalline islands floating in vapour. We look at how the atoms are arranged, in two ways.

Each atom is coloured by how hexagonal its neighbourhood is. For every bond to a neighbour closer than 1.4 σ, take the bond’s angle \theta and average e^{6i\theta} over the bonds. In a triangular crystal all six bonds line up and the magnitude, |\psi_6|, is close to 1; in a liquid the angles are scrambled and it falls towards 0. Atoms bright orange are crystalline, blue are not. (The same idea is behind the hexatic phase.) Switch the colouring to speed to see the same atoms the way the earlier posts drew them.

The left chart is the radial distribution function g(r): how likely you are to find a neighbour at distance r, compared with a gas of the same density. A crystal has sharp, repeating peaks at the lattice spacings; a liquid has a first peak and a couple of fading ripples; a gas is flat at 1. The histogram is built on the GPU. Every thread adds one to a bin with atomicAdd while it is already looping over its neighbours, and only the 120 counts come back. The right chart is the fraction of atoms with |\psi_6| > 0.8, over time.

At the starting temperature the islands stay crystalline, with a thin vapour between them. Warm it and the islands melt into liquid droplets: the peaks wash out and the orange drains away. Cool it and the liquid freezes into crystallites, usually with defects and grain boundaries rather than a perfect lattice.

Temperature alone does not tell the whole story, because the phase also depends on density. The density slider expands or compresses the box, and the atoms with it. At the crystal’s own density, 0.87, there is no room for a gas: even a hot box is just a dense liquid. Spread the atoms thinly and keep it warm, and the colours go blue, the first peak shrinks and g(r) flattens towards 1.

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 controls are read live from their input elements (referencing `viewof x`
  // gives the element, not the value), so moving them does not restart anything.
  const temperatureInput = viewof temperature;
  const densityInput = viewof density;
  const orderInput = viewof byOrder;

  const { canvas, device, context, format } = await gpu.setup({
    width: 640,
    height: 360,
    resolution: "auto",
  });

  // Reduced units: sigma = epsilon = mass = k_B = 1. A triangular lattice of
  // 52 x 34 atoms at spacing 1.15 sigma fills a 16:9 box exactly (34 rows is even,
  // so it tiles across the periodic boundary).
  const nx = 52;
  const ny = 34;
  const spacing = 1.15;
  const { positions, velocities, box } = md.triangularLattice({
    nx,
    ny,
    spacing,
    temperature: temperatureInput.value,
  });
  const [Lx, Ly] = box; // the box at the crystal's density; the slider expands it
  const count = nx * ny;
  const density0 = count / (Lx * Ly);

  // Start at the slider's density: stretch the lattice to it now, so the first
  // frame does not have to jump. The box scale is the same in x and y.
  let scale = Math.sqrt(density0 / densityInput.value);
  for (let k = 0; k < positions.length; k++) positions[k] *= scale;
  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

  // g(r): 120 bins of width 0.05 sigma out to 6 sigma
  const bins = 120;
  const rMax = 6;
  const binWidth = rMax / bins;
  const sampleEvery = 30; // read back every this many frames
  const keep = 240; // samples kept in the time series

  // Storage buffers: positions, velocities and the accelerations between passes,
  // plus the per-particle order parameter and the g(r) histogram (u32 counts)
  const positionBuffer = gpu.storageBuffer(device, positions);
  const velocityBuffer = gpu.storageBuffer(device, velocities);
  const accelerationBuffer = gpu.storageBuffer(device, new Float32Array(count * 2));
  const orderBuffer = gpu.storageBuffer(device, new Float32Array(count));
  const histogramBuffer = gpu.storageBuffer(device, new Uint32Array(bins));
  const rngBuffer = gpu.storageBuffer(device, crypto.getRandomValues(new Uint32Array(count)));

  // Readback (lib/gpu.js): the histogram, then the orders, back to back
  const histogramBytes = bins * 4;
  const orderBytes = count * 4;
  const readback = gpu.readback(device, [
    { buffer: histogramBuffer, bytes: histogramBytes },
    { buffer: orderBuffer, bytes: orderBytes },
  ]);

  // Uniform buffers, rewritten every frame. `params`: the thermostat's kick size
  // (f32), whether this frame builds the histogram (u32), the box scale (f32) and
  // the change in scale since last frame (f32). `box`: the current box, for drawing.
  const paramBuffer = device.createBuffer({
    size: 16, // uniform buffers are happiest in multiples of 16 bytes
    usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
  });
  const boxBuffer = device.createBuffer({
    size: 16,
    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 order,
  // and `integrate` reads and writes only its own particle.
  const computeModule = device.createShaderModule({
    code: `
      const L0 = vec2f(${Lx}, ${Ly});  // box at the crystal's density
      const DT = ${dt};
      const RC2 = 6.25;                // force cutoff 2.5 sigma, squared
      const SHELL2 = 1.96;             // first neighbour shell, 1.4 sigma, squared
      const DECAY = ${decay};          // thermostat: velocity kept each step
      const BIN_WIDTH = ${binWidth};
      const RMAX2 = ${rMax * rMax};    // g(r) range, squared

      // scale: the box is L0 * scale; ratio: scale now over scale last frame
      struct Params { noise: f32, sample: u32, scale: f32, ratio: f32 }

      @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> params: Params;
      @group(0) @binding(5) var<storage, read_write> orders: array<f32>;
      @group(0) @binding(6) var<storage, read_write> histogram: array<atomic<u32>>;

      ${md.rngWGSL}
      @compute @workgroup_size(64)
      fn forces(@builtin(global_invocation_id) id: vec3u) {
        let i = id.x;
        if (i >= arrayLength(&positions)) { return; }

        let L = L0 * params.scale;
        let here = positions[i];
        var a = vec2f(0.0);
        var bond = vec2f(0.0); // sum of e^(6 i theta) over first-shell neighbours
        var shell = 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));
          }

          // Bond-orientational order: the angle of each bond to a close neighbour
          if (r2 < SHELL2) {
            let theta = atan2(d.y, d.x);
            bond += vec2f(cos(6.0 * theta), sin(6.0 * theta));
            shell += 1.0;
          }

          // g(r): on a sampling frame, every thread counts its pairs into bins.
          // Many threads may hit the same bin, so the increment must be atomic.
          if (params.sample == 1u && r2 < RMAX2) {
            atomicAdd(&histogram[u32(sqrt(r2) / BIN_WIDTH)], 1u);
          }
        }
        accelerations[i] = a;
        // An atom with one neighbour is trivially "ordered" (its single bond lines
        // up with itself), so only count atoms with a real neighbourhood
        orders[i] = select(0.0, length(bond) / shell, shell >= 4.0);
      }

      // The density slider changed the box: stretch every position with it
      @compute @workgroup_size(64)
      fn rescale(@builtin(global_invocation_id) id: vec3u) {
        let i = id.x;
        if (i >= arrayLength(&positions)) { return; }
        positions[i] = positions[i] * params.ratio;
      }

      @compute @workgroup_size(64)
      fn integrate(@builtin(global_invocation_id) id: vec3u) {
        let i = id.x;
        if (i >= arrayLength(&positions)) { return; }

        let L = L0 * params.scale;

        // Leapfrog: kick the velocity...
        var v = velocities[i] + accelerations[i] * DT;

        // ...thermostat: slow it a little and add a random kick...
        var state = rng[i];
        v = v * DECAY + params.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 rescalePipeline = device.createComputePipeline({
    layout: "auto",
    compute: { module: computeModule, entryPoint: "rescale" },
  });
  const rescaleBindGroup = device.createBindGroup({
    layout: rescalePipeline.getBindGroupLayout(0),
    entries: [
      { binding: 0, resource: { buffer: positionBuffer } },
      { binding: 4, resource: { buffer: paramBuffer } },
    ],
  });
  const forceBindGroup = device.createBindGroup({
    layout: forcePipeline.getBindGroupLayout(0),
    entries: [
      { binding: 0, resource: { buffer: positionBuffer } },
      { binding: 2, resource: { buffer: accelerationBuffer } },
      { binding: 4, resource: { buffer: paramBuffer } },
      { binding: 5, resource: { buffer: orderBuffer } },
      { binding: 6, resource: { buffer: histogramBuffer } },
    ],
  });
  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: paramBuffer } },
    ],
  });

  // Render: a disc per particle, coloured by its speed or by its order (lib/md.js)
  const speedDiscs = md.discs(device, {
    format,
    positionBuffer,
    velocityBuffer,
    boxBuffer,
  });
  const orderDiscs = md.discs(device, {
    format,
    positionBuffer,
    boxBuffer,
    colour: {
      buffer: orderBuffer,
      range: [0.3, 0.95],
      low: [0.28, 0.42, 0.8],
      high: [1.0, 0.8, 0.52],
    },
  });

  // Charts: one Observable Plot per host, redrawn after every readback
  const [rdfHost, orderHost] = charts.children;
  const draw = ui.plotter(Plot);
  const rdf = new Float32Array(bins); // smoothed over samples
  const series = []; // {t, fraction}: atoms with |psi6| > 0.8
  let simTime = 0;

  function drawCharts() {
    const rows = Array.from(rdf, (g, k) => ({ r: (k + 0.5) * binWidth, g }));

    draw(rdfHost, {
      x: { label: "distance r (σ)", domain: [0, rMax] },
      y: { label: "g(r)", grid: true, domain: [0, 6], clamp: true },
      marks: [
        Plot.ruleY([1], { stroke: "currentColor", strokeOpacity: 0.4, strokeDasharray: "4 3" }),
        Plot.line(rows, { x: "r", y: "g", stroke: "#FF4834", clip: true }),
      ],
    });

    draw(orderHost, {
      x: { label: "time (reduced units)" },
      y: { label: "fraction crystalline", domain: [0, 1], grid: true },
      marks: [Plot.lineY(series, { x: "t", y: "fraction", stroke: "#FF4834" })],
    });
  }

  // Turn one readback into chart data
  function record(buffer, time, steps, density) {
    const hist = new Uint32Array(buffer, 0, bins);
    const orders = new Float32Array(buffer, histogramBytes, count);

    // g(r): pair counts per particle per step, over the number an ideal gas of the
    // same density would put in the ring [r, r + dr]
    for (let k = 0; k < bins; k++) {
      const r = (k + 0.5) * binWidth;
      const ideal = count * steps * density * 2 * Math.PI * r * binWidth;
      const g = hist[k] / ideal;
      rdf[k] = rdf[k] === 0 ? g : 0.7 * rdf[k] + 0.3 * g;
    }

    let crystalline = 0;
    for (let i = 0; i < count; i++) if (orders[i] > 0.8) crystalline++;
    series.push({ t: time, fraction: crystalline / count });
    if (series.length > keep) series.shift();

    drawCharts();
  }

  // One frame: `stepsPerFrame` simulation steps in a single compute pass, then
  // draw. Every `sampleEvery` frames the histogram is also built and read back.
  const workgroups = Math.ceil(count / 64);
  const zeros = new Uint32Array(bins);
  const paramBytes = new ArrayBuffer(16);
  const paramView = new DataView(paramBytes);
  let frameCount = 0;

  function frame() {
    const sample = !readback.pending && ++frameCount % sampleEvery === 0;

    // Density rho means a box scale of sqrt(rho0 / rho) in each direction; `scale`
    // is last frame's value, so `ratio` is how much to stretch the positions now
    const nextScale = Math.sqrt(density0 / densityInput.value);
    const ratio = nextScale / scale;
    scale = nextScale;

    paramView.setFloat32(0, Math.sqrt(temperatureInput.value * (1 - decay * decay)), true);
    paramView.setUint32(4, sample ? 1 : 0, true);
    paramView.setFloat32(8, scale, true);
    paramView.setFloat32(12, ratio, true);
    device.queue.writeBuffer(paramBuffer, 0, paramBytes);
    device.queue.writeBuffer(boxBuffer, 0, new Float32Array([Lx * scale, Ly * scale]));
    if (sample) device.queue.writeBuffer(histogramBuffer, 0, zeros);

    const encoder = device.createCommandEncoder();

    const compute = encoder.beginComputePass();
    compute.setPipeline(rescalePipeline);
    compute.setBindGroup(0, rescaleBindGroup);
    compute.dispatchWorkgroups(workgroups);
    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();
    simTime += stepsPerFrame * dt;

    const pass = gpu.beginRenderPass(encoder, context);
    (orderInput.value ? orderDiscs : speedDiscs).draw(pass, count);
    pass.end();

    if (sample) readback.copy(encoder);
    device.queue.submit([encoder.finish()]);

    if (sample) {
      const time = simTime;
      const rho = densityInput.value;
      readback.read((buffer) => record(buffer, time, stepsPerFrame, rho));
    }
  }

  gpu.run(invalidation, device, frame);

  // The canvas with the controls laid over its corner
  yield ui.stage(canvas, [temperatureInput, densityInput, orderInput]);
}

Moving Figures