Gas, liquid, crystal

Every atom here obeys the same rule, the Lennard-Jones potential. Two atoms closer than about one diameter σ shove each other apart hard; a little further out they attract weakly, with the strongest pull at 2^1/6σ, where the energy bottoms out at −ε. Past 2.5σ the model ignores the pull altogether, which costs almost nothing and lets the GPU look only at nearby atoms.

Nothing in the code says “liquid” or “crystal”. Gas, droplets, surface tension and hexagonal crystals all come out of this one curve and a temperature. That is why the potential, written in this 12-6 form in 1931, is still the first model taught in molecular simulation.

Colour by speed, and each atom is blue when slow and amber when fast. Temperature is nothing more than the average kinetic energy: in two dimensions each atom carries kT of it on average, so the readout divides the summed kinetic energy by the number of atoms.

Speeds are never all equal. Collisions spread them into the Maxwell–Boltzmann distribution, which in two dimensions is f(v) = (mv/kT)·exp(−mv²/2kT), peaking at √(kT/m). The histogram beside the box counts the real speeds on the GPU ten times a second and lays them over that curve for the measured temperature. Watch them agree.

A thermostat holds the temperature. Every step it scales all the velocities by one shared factor, nudging the total kinetic energy toward the setting, with a pinch of randomness so the box samples the right spread of energies rather than pinning them. Because every atom is scaled alike, it changes how hot the box is without slowing its flow, so droplets still fall freely. Drop the temperature into the liquid band and the measured value sits a little abovethe setting while droplets form: condensing atoms fall into each other's wells and release that energy as heat, which is latent heat happening in front of you.

In reduced units (σ, ε and the atom's mass all set to one) this exact model, two-dimensional and cut and shifted at 2.5σ, has its critical point at kT/ε = 0.459. Above it there is no boiling: liquid and gas are one fluid that simply thins out. Below roughly 0.42, the triple point, the atoms skip the liquid and freeze straight out of the vapour.

So the liquid lives in a narrow band. Set the temperature to 0.45 and the uniform gas breaks up into a sponge of droplets that slowly merge. At 0.12, try colouring by neighbours: an atom with exactly six neighbours sits in a triangular lattice and turns green, the packing you get when you push coins together on a table. The lattice spacing is 1.113σ, a hair inside the potential's minimum, because each atom also feels the pull of atoms further out.

For scale, argon is well described with σ = 3.4 Å and ε/k = 120 K, which is how the temperature slider converts to kelvin, and one unit of time τ is 2.15 picoseconds. Real argon is three-dimensional, though, and its triple point is far higher, so read the kelvin as a scale, not a prediction.

Each atom is one GPU thread. The engine is written for WebGPU compute shaders, and every step runs four of them:

1. Move.The thermostat's shared scale factor, half a kick from the current forces and a full drift: velocity Verlet, the workhorse of molecular dynamics. Atoms that cross a wall bounce, and the momentum they hand the wall is added up: that sum, per unit length and time, is the pressure readout.
2. Bin.The box is cut into cells at least 2.5σ wide, and each atom atomically bumps its cell's count.
3. Sort.A prefix sum turns counts into offsets and a scatter writes every atom's index into its cell's slice, so each cell's atoms sit together in memory.
4. Force. Each atom reads only the nine cells around it, sums the Lennard-Jones forces from the ones inside 2.5σ, and finishes the step with the second half kick, adding its kinetic energy to the total the thermostat reads next step.

The work per atom stays the same however many atoms there are, so the cost grows in step with the count: twenty thousand atoms take a few hundred thousand steps a minute on a laptop, and a million still move. Only a few hundred bytes of sums come back to the CPU, ten times a second, for the readouts. This is the foundation the rest of this section builds on.

• It is two-dimensional.Phase lines move in 3D: the same model's critical point there is above 1.
• It does not settle how 2D crystals melt. Theory allows a strange in-between hexatic phase. Large simulations find Lennard-Jones crystals melt that way only at high temperature and density; near the conditions here melting is an ordinary jump. Seeing either needs far longer runs than a browser tab.
• It is classical. Atoms are points with a force law. Electrons, bonds and quantum effects are not in it yet; they come next in this section.

• Lennard-Jones 1931 · J. E. Lennard-Jones, "Cohesion", Proc. Phys. Soc. 43, 461 (1931): the 12-6 form of the potential
• Verlet 1967 · L. Verlet, Phys. Rev. 159, 98 (1967): the integration scheme and the neighbour list behind most molecular dynamics
• Rahman 1964 · A. Rahman, Phys. Rev. 136, A405 (1964): the first molecular dynamics of liquid argon
• Smit & Frenkel 1991 · B. Smit and D. Frenkel, J. Chem. Phys. 94, 5663 (1991): the vapour–liquid critical point of the 2D Lennard-Jones fluid
• Li & Pica Ciamarra 2020 · Y.-W. Li and M. Pica Ciamarra, Phys. Rev. E 102, 062101 (2020): how 2D Lennard-Jones crystals melt
• Swope et al. 1982 · W. C. Swope, H. C. Andersen, P. H. Berens and K. R. Wilson, J. Chem. Phys. 76, 637 (1982): the velocity form of the Verlet integrator
• Bussi, Donadio & Parrinello 2007 · G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys. 126, 014101 (2007): canonical sampling through velocity rescaling, the thermostat used here
• Allen & Tildesley · M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, 2nd ed., Oxford (2017)
• LAMMPS units · LAMMPS documentation: Lennard-Jones reduced units and the 0.005τ default step