Hard spheres are particles that do not attract each other and only collide. Packed densely enough, they still form crystals, driven by entropy alone. In a team project of three students, we simulated them with event-chain Monte Carlo: how they crystallize, and how they melt in two dimensions. My parts were a faster search for the next collision and sedimentation under gravity.

Event-chain Monte Carlo

Standard Monte Carlo proposes a small random move and rejects it if two spheres would overlap. In a dense packing, most moves are rejected. Event-chain Monte Carlo rejects no moves: a sphere moves in a straight line until it touches another sphere, and that sphere continues the remaining distance. One chain moves many spheres, and the method still samples the correct equilibrium.

Finding the next collision

Each step of a chain needs the first sphere in the way. Checking all pairs of spheres costs time that grows with the square of their number. Instead, the box is divided into cells about one sphere wide, and each sphere is stored in the cell of its centre. A collision partner can only be in cells near the path.

A grid of square cells with spheres on it; the cells that hold a sphere’s centre are blue, although several spheres reach into neighboring cells

Searching along the path

The simple search checks shells of cells around the moving sphere, so its cost grows with the searched volume. I adapted Bresenham’s line algorithm, a standard method from computer graphics, to the three-dimensional grid. This ray tracing visits only the cells the path crosses and checks their direct neighbors.

Two grids side by side. Left, the naïve method: a circle around the start cell, as wide as the chain vector, fills a large square of orange cells. Right, ray tracing: only the cells along the arrow are blue, with a narrow band of orange cells around them
In two dimensions. Left: the shell search checks every cell within reach of the start (blue). Right: ray tracing walks the cells the path crosses (blue) and checks only their neighbors (orange).

Cells checked per step

Counted from geometry, ray tracing adds nine cells per step in three dimensions. After ten steps, the shell search has checked nearly 1,400 cells, ray tracing just over 100.

This helps when moves are long, in dilute systems. In a dense packing, moves end within the first cells, and both searches cost about the same.

Plot of the number of checked cells against the search step, one to ten: the naïve method curves up to nearly 1,400 cells, ray tracing rises in a straight line to just over 100
Cells checked after each step of the search, in three dimensions, counted from geometry.

Gravity

For sedimentation, the box got a fixed floor and ceiling. Gravity enters as an energy budget: for each upward move, a budget is drawn at random, and moving up uses it. When it runs out before a collision, a lifting event turns the move downward, as the ceiling does. No move is rejected, so the method stays rejection-free.

A sphere moving upward along a line meets, in turn, the point where its energy budget ΔU runs out (gravitation event), another sphere (collision event), the ceiling (wall event), or nothing (free movement)
Four ways an upward move can end: the energy budget runs out, another sphere, the ceiling, or no event at all.

Sedimentation

With 91,125 spheres, the bottom layers became a solid while the top stayed fluid. The transition between them spanned about three layers. A structure analysis classified each sphere by its neighbors: the solid was FCC and HCP crystal, with defects.

A view from below, in a second run with 64,000 spheres, shows small crystal patches that grow and merge into larger domains.

Eight square views from below of the lowest layers, frames 9 to 79: scattered green and red spheres at first, then patches of crystal that grow and merge into large green domains with red borders
The lowest layers of the 64,000-sphere run seen from below (0 ≤ z ≤ 6.5), frames 9 to 79; only spheres in crystal order are drawn. Green: FCC, red: HCP.

Validation against the equation of state

Before gravity was added, the code had to reproduce the known equation of state of hard spheres, the Carnahan–Starling equation. It matched in the fluid, showed a clear phase transition, and slightly overestimated the pressure in the solid.

With gravity, the weight of the spheres above sets the pressure at each height, so the pressure follows from the density profile. With one reference pressure chosen to match, it follows the equation reasonably well, slightly below it, and the phase transition is still visible.

Plot of pressure βP/ρ against packing fraction η: the simulation points with sedimentation lie slightly below the Carnahan–Starling curve through the gas and liquid region, flatten in the coexistence region, and rise along the Alder–Hoover–Young curve in the solid
Pressure under gravity, from the density profile (dots), against the Carnahan–Starling equation for the fluid and the Alder–Hoover–Young equation for the solid.

Crystal growth over time

Spheres in HCP order appeared first and levelled off after about 20 frames. FCC grew almost linearly and passed HCP after about 90 frames. FCC was still growing when the run ended. Next steps would be longer runs, smoother density profiles, and other strengths of gravity.

Plot of particle counts over 100 frames: particles without crystal order fall steadily; HCP rises quickly and levels off after about 20 frames; FCC grows almost linearly and passes HCP near frame 90
Spheres by structure over the run; one frame is a million chains (15.62 chains per particle).

Code

The simulation was written in C++ and later ported to Rust. The code was lost in a disk failure; the repository keeps the report and its figures. We exchanged code as files; one shared repository with feature branches would have kept it.

Context

Primary references

  1. Bernard, Krauth & Wilson (2009): Event-chain Monte Carlo algorithms for hard-sphere systems
  2. Michel, Kapfer & Krauth (2014): Generalized event-chain Monte Carlo: constructing rejection-free global-balance algorithms from infinitesimal steps
  3. Kampmann et al. (2021): Event-chain Monte-Carlo simulations of dense soft matter systems
  4. Verlet (1967): Computer “experiments” on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules
  5. Bresenham (1965): Algorithm for computer control of a digital plotter
  6. Carnahan & Starling (1969): Equation of state for nonattracting rigid spheres
  7. Alder, Hoover & Young (1968): Studies in molecular dynamics. V. High-density equation of state and entropy for hard disks and spheres
  8. Piazza, Bellini & Degiorgio (1993): Equilibrium sedimentation profiles of screened charged colloids: a test of the hard-sphere equation of state
  9. Ackland & Jones (2006): Applications of local crystal structure measures in experiment and simulation
  10. Marechal, Hermes & Dijkstra (2011): Stacking in sediments of colloidal hard spheres
  11. Annuss, Nienhaus & Lehner (2024): Hard sphere simulations, final report

Novelty scope: Hard-sphere sedimentation, event-chain Monte Carlo, and grid-based neighbor searches are established methods. The contributions are implementations in a team project: a Bresenham-style walk through the cell grid for the collision search, compared by counting the cells it checks; and gravity and fixed walls added as lifting events, checked against the Carnahan–Starling equation of state. The crystallization and two-dimensional parts of the report are my colleagues’ work.