Optimization of Parallel Tempering Monte-Carlo Algorithm on 2D Ising Model Alexey Rybina,b, Dmitrii Kapitana,b, Konstantin Nefedeva,b, Petr Andriushchenkoa,b,c and Vitalii Kapitana,b a Far Eastern Federal University, Vladivostok, 690922 Russia b Institute of Applied Mathematics, Far Eastern Branch, Russian Academy of Science, Vladivostok, 690041 Russia c ITMO University, Kronverksky Pr. 49, bldg. A, St. Petersburg, 197101, Russia Abstract We consider efficient algorithms for thermodynamical characteristics calculation of 2D Ising model. We discuss optimization algorithms for temperature in order to improve effectiveness of replica exchange. We implement and test algorithm on a two-dimensional square Ising lattice. Keywords 1 Parallel tempering, Ising model, Monte-Carlo 1. Introduction Considering the growth of data volumes, there is a need for new research in the field of magnetic data carriers. The researchers use the Monte Carlo method to simulate various spin structures. However, this method has a drawback: in the phase transition region, the simulation process slows down. To combat this effect, replica exchange is used, which allows simulating several systems with different temperatures in parallel, and also exchanges configurations between neighboring systems. For the best efficiency parameters, it is proposed to distribute the temperatures of the systems non- linearly. However, it is not known which of the distributions will give the best result. Therefore, an urgent task is to study ways to optimize the replica-exchange Monte Carlo method. An optimized set of temperature values increases the efficiency of the algorithm by making more frequent replica visits to the temperature extremum. However, for efficient operation, careful tuning of parameters is required to ensure optimal execution time [1]. In simulations with parallel tempering, the replica exchange rate strongly depends on the simulated statistical ensemble, that is, on the selected temperature points {𝑇1 , 𝑇2 , . . . , 𝑇𝑀 }. 2. 2D Ising model Consider the Ising model on a flat square lattice. The probability of any configuration of the investigated models is described by the Gibbs distribution [2]. It is well known that knowing of the statistical sum for a system of interacting spins allows one to strictly calculate all possible mean physical quantities fully describing the state of the system at given parameters. Currently, the studies of the phenomenon of magnetic transformations (transitions) mainly use numerical methods VI International Conference Information Technologies and High-Performance Computing (ITHPC-2021), September 14–16, 2021, Khabarovsk, Russia EMAIL: rybin.ae@dvfu.ru (Aleksey Rybin) ©️ 2020 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings (CEUR-WS.org) 155 (canonical and multicanonical MC) [3]. The Hamiltonian of the Ising spin system in the external magnetic has the form (1). 𝐻 = βˆ’βˆ‘π½ 𝑆 𝑆 βˆ’ β„Žβˆ‘π‘  , (1) 𝑖𝑗 𝑖 𝑗 𝑖 βŒ©π‘–,𝑗βŒͺ 𝑖 where 𝑆𝑖 𝑆𝑗 are spins of the system, 𝑖, 𝑗 denotes summarizing over the lattice with size 𝑁, β„Ž - external magnetic field. 2.1. Replica exchange Monte-Carlo The Monte Carlo method with Markov chains is used as the main component of this research project. The probability of any possible configuration is determined by the Gibbs distribution: 1 (2) 𝑃= exp (βˆ’π›½πΈ β€² (π‘₯)) 𝑍(𝛽) In this model, each spin interacts only with its nearest four neighbors through a direct ferromagnetic exchange interaction randomly distributed in the lattice nodes, provided that βˆ‘π‘§=4 𝑖=1 𝐽𝑖 = 0 Let us recall that the acceptance probability 𝑃𝑓𝑙𝑖𝑝 of the Metropolis-Hastings algorithm is determined by the formula: 1 (3) 𝑃𝑓𝑙𝑖𝑝 = min{1, 𝑒 βˆ’π›½βˆ†π» } ; 𝛽 = π‘˜π΅ 𝑇 At low temperatures 𝛽 - a very large positive number. If we propose a spin flip with a positive energy difference, i.e. π›₯𝐻 > 0, we have: π›½βˆ†π» ≫ 0 β‡’ 𝑒 βˆ’π›½βˆ†π» β‰ˆ 0 β‡’ 𝑃𝑓𝑙𝑖𝑝 β‰ˆ 0 (4) Most likely, it will go through a sequence of spin flips with a negative energy difference, forcing the system back to the energy minimum. Therefore, it is not possible to create states according to the Boltzmann distribution, resulting in biased sampling [4]. Replica exchange serves to improve the convergence of the Metropolis-Hastings algorithm in the problem at hand. A number of systems initialized with different temperatures of the Metropolis- Hastings algorithm exchange configurations during a loop performing value sampling [5]. This is done to allow configurations at high temperatures to move to systems with low temperatures when the simulation process continues, and to save low temperatures from falling into unwanted stable states with a minimum of energy. The algorithm presented on Figure 1. 156 Figure 1: Replica exchange Monte-Carlo algorithm 3. Optimization Hukushima [6] proposed a method, which is simpler than the feedback method, for determining replica temperature values. The scheme starts with an initial arbitrary temperature distribution, captures the extreme temperatures, and iteratively corrects the intermediate temperatures so that the probability of replica exchange for all neighboring temperatures is the same. Since this scheme is based on estimating the energy in each replica as a function of its inverse temperature, we call this method the β€œenergy” method [7,8,9]. This method belongs to the category of approaches that strive for a uniform rate of exchange between replicas. Let 𝛽𝑖 and 𝐸𝑖 refer to the inverse temperature and average energy of replica 𝑖 respectively. The goal of the energy method is to adjust 𝛽𝑖 so that the replica exchange probabilities of neighboring temperatures are equal. β„™(πΈπ‘–βˆ’1 , π›½π‘–βˆ’1 ↔ 𝐸𝑖 , 𝛽𝑖 ) = β„™(𝐸𝑖 , 𝛽𝑖 ↔ 𝐸𝑖+1 , 𝛽𝑖+1 ). ( 7) More precisely, the replicas are divided into two groups: odd and even. Fixing the inverse temperatures of one group, the inverse temperatures of the other group are then corrected one by one [10]. The detailed procedure is outlined in Figure 2. 157 Figure 2: Optimized replica exchange Monte-Carlo algorithm 4. Results The speeds of the conventional and optimized algorithm were compared on systems with sizes L=100, L=1600, L=4900, L=10000 particles. The results are shown in Figure 3. Also, to evaluate the performance of the algorithm, heat capacity plots (Figure 4) were plotted for systems with different sizes. 18000 Optimized Normal 16000 14000 12000 t, sec 10000 8000 6000 4000 2000 0 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 L Figure 3: Conventional and optimized algorithm speed comparison 158 Figure 4: Heat capacity comparison of 2 methods 5. Conclusion In this paper, optimization algorithm was considered for replica exchange Monte-Carlo method. We considered β€œenergy” method of optimization. A program was also written to compare optimized and regular algorithms within the framework of the conditions we are interested in. On the basis of the obtained results, it can be concluded that choosing optimized temperature set leads to advantage in execution speed on larger lattices. It will help in further studies of ferromagnetic and antiferromagnetic spin models. 6. Acknowledgements This work was financially supported by the state task of the Ministry of Science and Higher Education of the Russian Federation β„–075-00400-19-01. The studies were carried out using the resources of the Center for Shared Use of Scientific Equipment "Center for Processing and Storage of Scientific Data of the Far Eastern Branch of the Russian Academy of Sciences", funded by the Russian Federation represented by the Ministry of Science and Higher Education of the Russian Federation under project No. 075-15-2021-663. Additional computing resources were provided by the FEFU supercomputer cluster (cluster.dvfu.ru). 7. References [1] Rozada, Ignacio et al. β€œEffects of Setting Temperatures in the Parallel Tempering Monte Carlo Algorithm.” Physical Review E 100.4 (2019): n. pag. Crossref. Web. [2] Landau L. D., Lifshitz E. M. and Pitaevskii L. P. Statistical physics / by L.D. Landau and E.M. Lifshitz; translated from the Russian by J.B. Sykes and M.J. Kearsley. β€” 3d rev. and enl. ed. β€” Oxford; New York: Pergamon Press, 1980. β€” Translation of Statisticheskaia fizika. [3] Newman M. E. J., Barkema G. T. Monte Carlo methods in statistical physics / M.E.J. Newman and G.T. Barkema. β€” Oxford : Clarendon Press, 1999. [4] Tawn, Nicholas & Roberts, Gareth & Rosenthal, Jeffrey. (2019). Weight-preserving simulated tempering. Statistics and Computing. 10.1007/s11222-019-09863-3. 159 [5] Katzgraber, Helmut G et al. β€œFeedback-Optimized Parallel Tempering Monte Carlo.” Journal of Statistical Mechanics: Theory and Experiment 2006.03 (2006): P03018–P03018. Crossref. Web. [6] K. Hukushima, Domain-wall free energy of spin-glass models: Numerical method and boundary conditions, Phys. Rev. E 60 (1999). [7] Hamze, Firas, Neil Dickson, and Kamran Karimi. β€œRobust Parameter Selection For Parallel Tempering.” International Journal of Modern Physics C 21.05 (2010): 603–615. Crossref. Web. [8] David J. Earl, Michael W. Deem Β«Optimal Allocation of Replicas to Processors in Parallel Tempering SimulationsΒ» The Journal of Chemical Physics, 2005 [9] P. Atkins and J. De Paula, Physical Chemistry for the Life Sciences (Oxford University Press, 2006). [10] H. G. Katzgraber, S. Trebst, D. A. Huse, and M. Troyer, Feedback-optimized parallel tempering Monte Carlo, J. Stat. Mech. P03018 (2006). 160