Simulating spin models on GPU

Over the last couple of years it has been realized that the vast computational power of graphics processing units (GPUs) could be harvested for purposes other than the video game industry. This power, which at least nominally exceeds that of current CPUs by large factors, results from the relative simplicity of the GPU architectures as compared to CPUs, combined with a large number of parallel processing units on a single chip. To benefit from this setup for general computing purposes, the problems at hand need to be prepared in a way to profit from the inherent parallelism and hierarchical structure of memory accesses. In this contribution I discuss the performance potential for simulating spin models, such as the Ising model, on GPU as compared to conventional simulations on CPU.

💡 Research Summary

The paper presents a detailed study of how to exploit modern graphics processing units (GPUs) for high‑performance Monte Carlo simulations of spin systems, focusing on the ferromagnetic Ising model in two and three dimensions. After a brief motivation that highlights the growing computational demands of statistical‑physics problems such as spin glasses and protein folding, the author reviews the hierarchical memory architecture of NVIDIA GPUs: fast per‑core registers, on‑chip shared memory (16 KB on Tesla, 48 KB on Fermi), large but high‑latency global memory, and constant/texture caches. The key design principle is to minimize global‑memory traffic and to keep as many threads as possible active to hide memory latency.

To achieve this, the author introduces a double‑checkerboard decomposition. The lattice is first divided into coarse B × B tiles; each tile is then split again into fine T × T sub‑tiles arranged in a checkerboard pattern. This allows the assignment of one coarse tile to a CUDA thread block, with the even and odd sub‑lattices updated in separate phases. All spins of a tile, together with a one‑site boundary halo, are loaded cooperatively by the block’s threads into shared memory, ensuring coalesced global‑memory accesses. Once resident in shared memory, each thread updates its assigned spins using the Metropolis acceptance rule, synchronizes, and then updates the opposite sub‑lattice. The process is repeated k times (the “multi‑hit” technique) before the next tile is processed, thereby amortizing the cost of the initial data transfer.

Random numbers are generated per thread using simple 32‑bit linear congruential generators (LCGs) with distinct seeds. Although the period (2³² ≈ 10⁹) is far shorter than the total number of random draws in large simulations, the author reports that for the 2D Ising model the results (energy, specific heat) remain statistically correct. However, when using disjoint subsequences of a 64‑bit LCG, significant systematic errors appear, prompting the recommendation of higher‑quality generators such as lagged‑Fibonacci for production runs.

Performance measurements were carried out on an NVIDIA Tesla C1060 (4 GB memory) and compared with a 3.0 GHz Intel Core 2 Quad CPU (4 MB and 6 MB cache configurations). For lattice sizes L = 16–1024, optimal fine‑tile sizes were found to be T = 4 for L ≤ 64, T = 8 for L = 128, and T = 16 for larger L. With shared‑memory loading and multi‑hit updates (k = 5, 10, 100), the GPU achieved single‑spin‑flip times as low as 0.1 ns in 2D and 0.24 ns in 3D, corresponding to throughputs exceeding 100 GFLOP s⁻¹—roughly 10 % of the theoretical peak of the C1060. Compared with a carefully optimized CPU implementation, speed‑up factors of up to 100× were observed for the 2D case, and nearly 300× for the 3D case at L = 256. The author notes that these gains are fragile: cache effects, multi‑core utilization, and the choice of random‑number generator can significantly alter the observed speed‑up.

The paper also discusses extensions beyond the simple Ising model. Continuous‑spin systems such as the Heisenberg model would benefit from the GPU’s native single‑precision floating‑point performance. Disordered systems can be simulated in parallel by processing many disorder realizations simultaneously, and combining this with asynchronous multi‑spin coding yields performance of about 0.15 ps per spin flip for the Edwards‑Anderson spin glass. Future work will address more sophisticated pseudo‑random generators, hybrid schemes that combine Metropolis updates with cluster algorithms near criticality, and adaptations to newer CUDA‑capable architectures (Fermi, Kepler, etc.).

In summary, the study demonstrates that a careful mapping of the algorithmic structure onto the GPU’s memory hierarchy—via double checkerboard tiling, shared‑memory staging, and multi‑hit updates—allows Monte Carlo spin‑model simulations to achieve order‑of‑magnitude speed‑ups over optimized CPU codes. This provides a practical blueprint for researchers in statistical physics and related fields who wish to harness commodity GPUs for large‑scale lattice simulations.