GPU accelerated population annealing algorithm

Abstract

Population annealing is a promising recent approach for Monte Carlo simulations in statistical physics, in particular for the simulation of systems with complex free-energy landscapes. It is a hybrid method, combining importance sampling through Markov chains with elements of sequential Monte Carlo in the form of population control. While it appears to provide algorithmic capabilities for the simulation of such systems that are roughly comparable to those of more established approaches such as parallel tempering, it is intrinsically much more suitable for massively parallel computing. Here, we tap into this structural advantage and present a highly optimized implementation of the population annealing algorithm on GPUs that promises speed-ups of several orders of magnitude as compared to a serial implementation on CPUs. While the sample code is for simulations of the 2D ferromagnetic Ising model, it should be easily adapted for simulations of other spin models, including disordered systems. Our code includes implementations of some advanced algorithmic features that have only recently been suggested, namely the automatic adaptation of temperature steps and a multi-histogram analysis of the data at different temperatures. The program calculates the internal energy, specific heat, several magnetization moments, entropy and free energy of the 2D Ising model on square lattices of edge length L with periodic boundary conditions as a function of inverse temperature β. Solution method: The code uses population annealing, a hybrid method combining Markov chain updates with population control. The code is implemented for NVIDIA GPUs using the CUDA language and employs advanced techniques such as multi-spin coding, adaptive temperature steps and multi-histogram reweighting. Restrictions: The system size and size of the population of replicas are limited depending on the memory of the GPU device used. Unusual features: Additional comments: Code repository at https://github.com/LevBarash/PAising. Running time: For the default parameter values used in the sample programs, L = 64, θ = 100, β0 = 0, β f = 1, ∆β = 0.005, R = 20 000, a typical run time on an NVIDIA Tesla K80 GPU is 151 seconds for the single spin coded (SSC) and 17 seconds for the multi-spin coded (MSC) program (see Sec. 2 for a description of these parameters). PROGRAM SUMMARY

Estimating the density of states of frustrated spin systems

“…The new algorithm introduced here, which we call EPA, is not based on Markov chains but on the sequential Monte Carlo method. PA was first studied in [53,54] and more recently developed further in [55][56][57][58][59][60][61]. It is based on the initialization of a population of replicas drawn from the equilibrium distribution at high temperatures, which is then subsequently cooled to lower and lower temperatures.…”

“…During this process, a combination of population control and spin flips is used to ensure that the ensemble remains in equilibrium. The simulation entails the following steps [55,58].…”

“…temperature steps were used on increasing β i >β i−1 [55], it turns out that a better, parameter-free method consists of choosing β i to ensure a certain overlap of the energy distributions between the two temperatures [58]. This overlap can be computed from the resampling factors,…”

“…The method is not very sensitive to the precise value of α * , and we choose α * =0.86 in the runs below. While the algorithm described above is just PA [55] improved by adaptive temperature steps [58,60], the possibility of sampling the entropy arises from a combination of the method with multi-histogram techniques [62]. An estimator of the free energy follows directly from the resampling factors [55],…”

“…This approach is naturally suited for (moderately or massively) parallel calculations as the R replicas are simulated independently of each other and the only interaction occurs during resampling. An efficient GPU implementation was discussed in [58]. Importantly for our application, EPA does not require any prior knowledge of the range of realized energies.…”