GPU accelerated Monte Carlo simulations of lattice spin models

We consider Monte Carlo simulations of classical spin models of statistical mechanics using the massively parallel architecture provided by graphics processing units (GPUs). We discuss simulations of models with discrete and continuous variables, and using an array of algorithms ranging from single-spin flip Metropolis updates over cluster algorithms to multicanonical and Wang-Landau techniques to judge the scope and limitations of GPU accelerated computation in this field. For most simulations discussed, we find significant speed-ups by two to three orders of magnitude as compared to single-threaded CPU implementations.

💡 Research Summary

The paper investigates the use of modern graphics processing units (GPUs) to accelerate Monte Carlo simulations of classical lattice spin models, covering a range of algorithms from local single‑spin Metropolis updates to non‑local cluster methods and generalized‑ensemble techniques such as multicanonical (MUCA) and Wang‑Landau (WL) sampling. The authors first motivate the need for high‑performance computing in statistical‑physics simulations, noting historic reliance on special‑purpose machines (e.g., cluster processors, the Janus FPGA system) and the recent emergence of general‑purpose GPU computing as a more accessible alternative. They adopt NVIDIA’s CUDA toolkit for implementation, citing its maturity over OpenCL at the time of the study.

For the Metropolis algorithm, the authors treat O(n) spin models (Ising, XY, Heisenberg) on square (2D) and simple‑cubic (3D) lattices with nearest‑neighbour interactions. Parallelism is achieved through a double‑checkerboard decomposition combined with a two‑level hierarchical tiling scheme. Large tiles (e.g., 16 × 16 spins) are loaded into shared memory, and within each tile the two sub‑lattices are updated concurrently, separated by synchronization barriers. To amortize the cost of loading tiles, a multi‑hit parameter k (e.g., k = 100) is used, allowing each tile to be updated repeatedly before moving to the next sub‑lattice. Additional optimizations include pre‑tabulating Boltzmann factors as a texture, generating random numbers even when not strictly required to reduce thread divergence, and employing simple 32‑bit linear‑congruential generators per thread (acceptable for the precision required). Performance benchmarks on a Tesla C1060 and a newer GTX 480 (Fermi architecture) against a single‑threaded Intel Core 2 Quad Q9650 show that the GTX 480 achieves about 0.03 ns per spin flip for the 2D Ising model—a 235‑fold speed‑up. For continuous‑spin models, mixed‑precision calculations (single‑precision spins, double‑precision accumulators) and CUDA fast‑math intrinsics yield speed‑ups exceeding 1000× relative to the CPU. Multi‑spin coding further reduces spin‑flip times to the picosecond regime.

Cluster updates, exemplified by the Swendsen‑Wang algorithm, present a greater challenge because the cluster‑identification step is intrinsically non‑local, especially near criticality where percolation occurs. The authors adopt a two‑stage approach: first, clusters are identified independently within tiles, ignoring bonds that cross tile boundaries; second, a global consolidation step merges these partial clusters using a union‑find data structure. Various labeling schemes (Hoshen‑Kopelman, breadth‑first search, self‑labeling) were tested; self‑labeling proved most efficient for tile sizes up to T ≈ 16, despite its O(T³) operation count, because the high degree of parallelism outweighs the algorithmic cost. This implementation yields up to a 20‑fold speed‑up over the CPU reference, which, while modest compared with local updates, demonstrates that cluster algorithms can still benefit from GPU acceleration when carefully engineered.

Generalized‑ensemble methods (MUCA and WL) require knowledge of a global reaction coordinate (energy) before each update, which serializes the computation. To exploit GPU parallelism, the authors employ a “windowing” strategy: the total energy range is divided into overlapping windows, each processed by a separate thread block. Within a window, each thread simulates an independent replica, accumulating histograms and density‑of‑states estimates in shared memory. All calculations are performed in single precision, with results matching double‑precision CPU references. For the 2D Ising model, MUCA achieves a 128‑fold speed‑up, comparable to local algorithms, while WL attains a 46‑fold improvement. The lower WL speed‑up is attributed to the stochastic nature of the WL modification factor, which leads to thread divergence and idle cores; the authors suggest that load‑balancing schemes could mitigate this limitation.

Table 1 summarizes spin‑flip times for various models, lattice sizes, and algorithms across the three platforms (CPU, Tesla C1060, GTX 480). Notably, the GPU implementations maintain high efficiency only when the problem size is large enough to fully occupy the available cores (≈240 cores for C1060, ≈480 cores for GTX 480). The paper concludes that GPUs provide a powerful, cost‑effective platform for Monte Carlo simulations of lattice spin systems, delivering speed‑ups of two to three orders of magnitude for local updates and respectable gains for non‑local or ensemble methods. The authors emphasize that careful attention to memory hierarchy, thread synchronization, random‑number generation, and numerical precision is essential for achieving optimal performance, and they anticipate that the presented techniques can be extended to more complex systems such as spin glasses, quantum spin models, and long‑range interacting systems.