Massively parallel Monte Carlo for many-particle simulations on GPUs

Current trends in parallel processors call for the design of efficient massively parallel algorithms for scientific computing. Parallel algorithms for Monte Carlo simulations of thermodynamic ensembles of particles have received little attention because of the inherent serial nature of the statistical sampling. In this paper, we present a massively parallel method that obeys detailed balance and implement it for a system of hard disks on the GPU. We reproduce results of serial high-precision Monte Carlo runs to verify the method. This is a good test case because the hard disk equation of state over the range where the liquid transforms into the solid is particularly sensitive to small deviations away from the balance conditions. On a Tesla K20, our GPU implementation executes over one billion trial moves per second, which is 148 times faster than on a single Intel Xeon E5540 CPU core, enables 27 times better performance per dollar, and cuts energy usage by a factor of 13. With this improved performance we are able to calculate the equation of state for systems of up to one million hard disks. These large system sizes are required in order to probe the nature of the melting transition, which has been debated for the last forty years. In this paper we present the details of our computational method, and discuss the thermodynamics of hard disks separately in a companion paper.

💡 Research Summary

The paper addresses a long‑standing challenge in computational statistical physics: how to parallelize Monte Carlo (MC) sampling, which is intrinsically sequential, on modern massively parallel hardware such as graphics processing units (GPUs). The authors develop a novel algorithm that respects detailed balance while allowing thousands of trial moves to be processed concurrently. Their approach hinges on a checkerboard domain decomposition. The simulation box is tiled into square cells that are colored alternately; cells of the same color are never neighbors and can therefore be updated in parallel without risking particle overlap across cell boundaries. Each CUDA thread block is assigned to a cell, loads the cell’s particle list into shared memory, and generates candidate displacements for its particles. Candidate moves are drawn from a uniform distribution of direction and distance, and a Metropolis acceptance test is performed after checking for overlaps with particles in the same cell and the eight neighboring cells.

To guarantee detailed balance, the algorithm alternates the color of the active cells after a fixed number of MC steps, effectively swapping the roles of black and white cells. In addition, a global re‑partitioning step periodically redistributes particles that have crossed cell borders, ensuring that the Markov chain remains reversible and that the stationary distribution is unchanged. Random numbers are produced by per‑thread XORWOW generators with distinct seeds, eliminating correlations that could bias the sampling.

Implementation details are carefully tuned for the NVIDIA Tesla K20 architecture. Particle data are stored in a structure‑of‑arrays layout to enable coalesced memory accesses. Shared memory is used to cache local particle coordinates, reducing global memory traffic, and atomic operations are avoided wherever possible. The authors benchmark their code against a high‑precision serial MC implementation running on an Intel Xeon E5540 core. The GPU version achieves more than one billion trial moves per second, a speed‑up factor of 148× over the single‑core baseline. Energy consumption is reduced by a factor of 13, and performance per dollar improves by 27×. Scaling tests show that efficiency remains high from 10⁴ up to 10⁶ particles; the latter size would be infeasible on a CPU in a reasonable time frame.

Validation is performed by reproducing the hard‑disk equation of state, especially in the narrow density window where the liquid‑to‑solid transition occurs. This regime is notoriously sensitive to violations of detailed balance, and the agreement with published serial results confirms that the parallel algorithm does not introduce systematic bias.

The scientific payoff of the accelerated simulation is demonstrated by calculating the equation of state for systems containing up to one million hard disks. Such large‑scale data are essential for resolving the decades‑old debate over the nature of two‑dimensional melting (whether it follows the Kosterlitz‑Thouless‑Halperin‑Nelson‑Young scenario or a first‑order transition). The authors defer the thermodynamic analysis to a companion paper, focusing here on the computational methodology.

In the discussion, the authors acknowledge limitations: the current implementation is optimized for short‑range, hard‑core interactions, and extending it to long‑range potentials (e.g., Lennard‑Jones, Coulomb) will require larger interaction neighborhoods and possibly a different domain‑decomposition strategy. Multi‑GPU scaling and efficient inter‑node communication are identified as future work, as is the adaptation of the framework to other MC ensembles such as grand‑canonical or Gibbs‑ensemble simulations.

Overall, the work provides a robust, detailed‑balance‑preserving blueprint for massively parallel MC simulations on GPUs, opening the door to high‑throughput studies of large particle systems in condensed‑matter physics, chemistry, and materials science.