Optimal phase space sampling for Monte Carlo simulations of Heisenberg spin systems

We present an adaptive algorithm for the optimal phase space sampling in Monte Carlo simulations of 3D Heisenberg spin systems. Based on a golden rule of the Metropolis algorithm which states that an acceptance rate of is ideal to efficiently sample the phase space, the algorithm adaptively modifies a cone-based spin update method keeping the acceptance rate close to . We have assessed the efficiency of the adaptive algorithm through four different tests and contrasted its performance with that of other common spin update methods. In systems at low and high temperatures and anisotropies, the adaptive algorithm proved to be the most efficient for magnetization reversal and for the convergence to equilibrium of the thermal averages and the coercivity in hysteresis calculations. Thus, the adaptive algorithm can be used to significantly reduce the computational cost in Monte Carlo simulations of 3D Heisenberg spin systems.


Introduction
In recent years, atomistic models of magnetic materials have allowed the study of new magnetic phenomena, especially in nanomaterials. These models are ideal to study complex magn etic behaviors such as exchange bias and surface aniso tropy because they take into account finite size effects and changes in the magnetization at atomic scale [1]. The clas sical 3D Heisenberg model is the most common prototype of an atomistic model of magnetic materials with continuous degrees of freedom. It has been successfully employed to accurately describe the magnetic behavior of many physically interesting magnetic systems [2].
In systems represented by classical spin models, such as the classical 3D Heisenberg model, numerical simulations are usually used because of the difficulty to solve by analytical approaches their partition function. The thermal averages of a system can be estimated employing Markov chain Monte Carlo methods. For the generation of new states of the systems in the Markov chain, local and cluster spin update algorithms have been proposed. Cluster update algorithms, such as the Wolff [3] and the Swendsen-Wang [4] algorithms, are useful to reduce the critical slowing down generated close to the critical temperature (T c ). However, the applicability of these algorithms is limited according to the terms considered in the Hamiltonian. In its basic form, the Wolff and the Swendsen-Wang cluster algorithms are only applicable to Hamiltonians that just include exchange interaction contributions. If other contributions (e.g. magnetic anisotropy and/or externally applied magnetic field) are considered, the cluster algorithms must be modified correspondingly, increasing the complexity of the algorithms and its implementation. Furthermore, the efficiency of cluster algorithms can be even lower than that of local update algorithms for small systems [2]. On the other hand, local update algorithms have shown high flexibility which allows its application to a great variety of systems. Therefore, much effort have been made to enhance the perfor mance of these algorithms, which relies on their efficiency to sample the phase space and avoid long correlation times that slow the convergence of the thermal averages to equilibrium.
In the Monte Carlo Metropolis algorithm new states of the system are generated from a previous state by randomly selecting a single spin and updating it to a new trial position (trial move). Then, a transition probability, which depends on the energy difference between the initial and final state, is used to accept or reject the new state. Because the choice of the trial move has a direct influence on the efficiency of the algorithm and its physical interpretation, different types of trial moves have been proposed. Common choices include the spin flip, the random and the small step moves [5]. In the spin flip move, the direction of the spin is reversed. This trial move is the same employed for Ising spin systems. However, for classical spin systems it does not satisfy the condition of ergodicity because it can not sample the whole phase space. In the random move a new spin direction is chosen at random independently of the spin initial direction, sampling the whole phase space. In the small step move, the new spin direction is generated from a uni form probability distribution within a cone with a given opening angle around the initial spin direction, hence each spin can only move by small angular changes limited by the opening angle.
Hinzke and Nowak [5] evaluated the efficiency of these trial moves by running simulations of the thermally activated reversal of the magnetization of a ferromagnetic particle. Their results showed that for high anisotropy systems, where the reversal mechanism is nucleation, the least efficient trial move is the small step move while the most efficient is the spin flip move. For low anisotropy systems, where the reversal mechanism is coherent rotation, the least efficient trial move is the random move while the most efficient is the small step move. Independently of the anisotropy value of the system, the small step move is the least efficient move at high temper atures. They also showed that a trial move which combines the three types of trial moves has good efficiency independently of the anisotropy and the temperature values of the system.
In particular, the small step move is not efficient for high anisotropy systems because it can not produce large angular changes in the spins required to overcome the anisotropy energy barrier. Therefore, the acceptance rate of the new states is low. At high temperatures, although the acceptance rate is high, the sampling of the phase space is made by very small steps. A golden rule when using the Metropolis algorithm states that an acceptance rate of 50% is ideal to efficiently sample the phase space of the system [6,7]. High or low acceptance rates generate a sampling of the phase space by very small steps or a rejection of almost every new state, respectively. A varia tion of the small step move known as the Gaussian move [1] generates the new spin direction in the vicinity of the initial spin direction employing a Gaussian distribution and guaran teeing that all possible states are accessible. Unlike the other trial moves, the Gaussian move straightforwardly allows to control the acceptance rate by adjusting the value of the cone width, while allowing the new spin state to point to any direc tion on the unit sphere with nonzero probability.
In this work, we present a new adaptive algorithm for the Gaussian move (adaptive move) in 3D Heisenberg Monte Carlo simulations. The adaptive move keeps the acceptance rate close to 50%, enhancing the efficiency of the phase space sampling and generating low correlation times. We have made several tests to assess the performance of the algorithm, showing that it works efficiently at all temperatures and in systems with high and low anisotropy values. Our algorithm significantly improves the computational efficiency of Monte Carlo spin simulations allowing faster or higher quality statistical results for the same number of Monte Carlo steps. We have imple mented our algorithm in two open source software packages, VAMPIRE [1] and VEGAS [8], to provide sample code and easy access to the community for use in magnetic simulations.

Heisenberg Hamiltonian and the Monte Carlo method
We considered a system with ferromagnetic spin moments located on a cubic lattice with periodic boundary conditions. The size of the system is 10 × 10 × 10 for a total of N = 1000 spin moments. Each spin moment was modeled using the clas sical 3D Heisenberg model. The Hamiltonian, which includes nearest neighbor exchange interaction, uniaxial anisotropy and magnetic field interactions, is given by where i<j means the sum over the nearest neighbors pairs, J is the exchange interaction constant, S i and S j are the spins of the magnetic sites labeled i and j, respectively, k v is the aniso tropy constant, k is the canonical vector in the z direction and B is the magnetic field intensity. Both the uniaxial anisotropy and the magnetic field are directed along the z axis. The values for k v and B were normalized by J in order to perform generic simulations. Some of the simulations were made for k v /J = 0.001 and 1.0 to recreate low and high anisotropy systems, where the critical temperatures are k B T/J = 1.485 427 and 1.753 970, respectively.
The total magnetization (M) was computed according to where N is the total number of spins and M = N i S i .

Trial moves
For comparison purposes, we implemented other common trial moves besides the presented adaptive move. Figure 1 shows schematics of the different types of trial moves and the visualization of their Monte Carlo sampling on the unit sphere. The spin flip move reverses the spin direction according to where S i and S i are the initial and new spin directions, respec tively. In the random move the new spin direction is chosen at random, according to where Γ is a Gaussian distributed random vector. In the small step move, the new spin direction is also generated at random but within a cone with a given opening angle around the initial spin direction. And, in the Gaussian move the new spin direc tion is generated in the vicinity of the initial spin direction according to [1] where σ is proportional to the width of a cone around the initial spin direction. Unlike in the small step move, in the Gaussian move the sampling on the unit sphere is not limited by σ (see figure 1). When using the Gaussian move, the acceptance rate can be adjusted by varying the value of σ. Figure 2(a) shows the acceptance rate (R) as a function of σ at low temperature for kv J = 1.0. As indicated in the figure, there is an optimum cone width (σ opt ) at which the acceptance rate is 50%. Hence, if we want to keep an acceptance rate of 50% at all temperatures for a given system, σ should be a function of parameters such as the temperature, the exchange coupling, the anisotropy and the field intensity. Figure 2 of σ opt for different anisotropy values. As expected, at higher temperatures it is necessary to use a higher cone width in order to keep an acceptance rate of 50%. However, above T c , it is not possible to keep this acceptance rate because almost all the new states are accepted independently of the cone width value, then the acceptance rate is always higher than 50%. As σ increases, the Gaussian move tends to a random move as shown in figure 3, where the distribution of the polar angle (θ), the azimuthal angle (φ) and the cosine of the polar angle (cos(θ)) of spins generated using the Gaussian move (see equation (5)) are shown. Such as in the case of a random move, the mean value of θ tends to 90 • and the mean value of cos(θ) tends to 0, while the mean value of φ is always close to 180 • .

Adaptive algorithm for the Gaussian move
Because σ opt is characteristic of every system at a given temper ature and given Hamiltonian parameters, to find a gen eral equation for σ opt that works in any system would be a very complex task. Therefore, we have developed an adaptive algorithm for the Gaussian move that changes the cone width adaptively to keep an acceptance rate close to 50%.
The adaptive move is developed as follows: at each temper ature, the simulation starts using a high cone width (σ = 60) in the first Monte Carlo step (MCS). From then on, every MCS, the cone width is recalculated by multiplying the current cone width by a factor obtained according to the acceptance rate in the previous MCS. The selection of the factor is made such that the cone width approaches values close to the optimum cone width. From results as those shown in figure 2(a), it is possible to observe that a good approximation for the factor ( f ) as a function of the acceptance rate at all temperatures is of the form Therefore, when the acceptance rate R = 50%, the cone width is multiplied by 1, and when the acceptance rate is high (low) the cone width is multiplied by a large (small) factor approaching the optimum cone width. Figure 4 shows the time dependence of the acceptance rate and the cone width using the adaptive move at different temperatures when the system is initially ordered (all the spin moments pointing in the z direction) and disordered (all the spin moments pointing in a random direction). Independently of the spin moments' initial state, the acceptance rate converges to a specific value. At k B T/J = 0.1, because the simulation starts with a high cone width, the acceptance rate is initially very low when the system is initially ordered (see figure 4(a)). Then, the adaptive move keeps decreasing the cone width to a very low value, according to the equation (6), increasing the acceptance rate which stabilizes close to 50% within few MCS. On the other hand, when the system is initially disordered (see figure 4(d)), the acceptance rate is initially higher because large angular changes in the direction of the spin moments are taking place to order the system. Then, the acceptance rate starts to decrease because smaller angular changes are required as the system is being ordered. However, once the the cone width has decreased enough, the acceptance rate starts to increase approaching a value close to 50%. Figures 4(b) and (e) show the time dependence of the acceptance rate at a high temperature below T c (k B T/J = 1.48). In these cases, the adaptive move requires more MCS to reach equilibrium and the cone width For cone width values between 0 and 100, 10 000 spins were generated and the mean of their polar angles, azimuthal angles and cosine of the polar angles was calculated. For cone width values higher than 60 (shaded region), it is reasonable to assume that the distributions are already stabilized.
stabilizes at a higher value than the previous cases. Above T c (k B T/J = 2.0), when the system is initially ordered (see figure 4(c)), the acceptance rate is low because the exchange energy prevents the acceptance of large angular changes in the direction of the spin moments. For this reason, the cone width initially decreases in order to increase the acceptance rate. As mentioned before, at temperatures above T c , it is not possible to reach an acceptance rate of 50%. Then, the acceptance rate keeps increasing past 50% and, consequently, the adaptive move starts to increase the cone width, trying to decrease the acceptance rate. When the system is initially disordered (see figure 4(f)), the acceptance rate is initially high and stabilizes at a lower value. Independently of the spin moments initial state, the cone width keeps increasing indefinitely when the acceptance rate stabilizes above 50%. Therefore, we reset the cone width to 60 every time it reaches higher values because, at this value, the Gaussian move works as the random move (see figure 3) and employing a higher value would produce the same results. For this reason, it is expected that the adaptive move has the same efficiency as the random move at temper atures above T c , where the cone width stabilizes at 60.

Results and discussion
In order to assess the performance of the adaptive move, we have made four tests comparing its efficiency with that of other common trial moves: the spin flip, random, small step and Gaussian moves. Also, we considered a combinational move which includes three of the aforementioned trial moves. As implemented by Hinzke and Nowak [5], we considered the small step move with a fixed opening angle of 30 • and the combinational move as a combination of the spin flip, small step and random moves. In the combinational move, one of the three trial moves is selected from a set, composed by three random, one small step and one spin flip moves, at each MCS. All the tests were carried out in systems with low and high anisotropy values. Before making any of the tests, it is important to guarantee that the computed thermal averages are the same indepen dently of the trial move employed. Figure 5 shows the thermal dependence of the magnetization and the energy using the dif ferent trial moves for low and high anisotropy values. Results of Landau-Lifshitz-Gilbert (LLG) spin dynamics simula tions are also shown because in one of the tests the efficiency of the adaptive move is compared to that of spin dynamics. At high anisotropy (k v /J = 1.0), the spin moments remain more ordered as temperature increases than at low anisotropy (k v /J = 0.001), which increases the critical temperature (see figure 5(a)). This behavior also generates an important dif ference in the energy between both systems, specially at low temperatures (see figure 5(b)). For both high and low aniso tropy values, all the trial moves produced the same results, indicating that all of them can correctly sample the phase space and relax the system to equilibrium.

Convergence to equilibrium and integrated relaxation times
When computing averages of the thermal properties of a system, it is necessary to ensure that the system is already in  equilibrium. Depending on the characteristics of the system, such as its size, temperature and anisotropy constant, the time required to reach equilibrium can increase greatly. Usually, a large number of MCS needs to be discarded to ensure system equilibrium. Therefore, speeding up the convergence to equi librium and enhancing the statistical quality of the thermal averages is important.
In this test, we assessed the performance of the different trial moves when the system is both approaching and is already in equilibrium. First, we evaluated the convergence to equilib rium of the magnetization and the energy as a function of time at low and high temperatures below T c , as shown in figure 6. Because at low temperatures the thermal energy is low, the spin moments are very ordered and its direction tends to vary by small angular changes. Therefore, the small step move has good efficiency at low temperatures (see figures 6(a) and (c)). However, at high temperatures (see figures 6(b) and (d)), the small step move is the least efficient move because the thermal energy can overcome the exchange and anisotropy energies, generating disorder in the system and high fluctuations in the direction of the spin moments. Conversely, the random move is very efficient at high temperatures but is the least efficient at low temperatures. In the case of the combinational move, although it has good efficiency at all temperatures, it is not the best move in any of the cases. The adaptive move is the most efficient move at low temperatures and presents the same efficiency of the random move at high temperatures. Then, in general, the adaptive move is the most efficient move to relax both thermal averages to equilibrium independently of the anisotropy value of the system.
In the second part of the tests, we computed the integrated relaxation time (τ ) for the different trial moves. The statistical error of the thermal averages depends on the number of sta tistically independent configurations generated in the simula tion, and this is the total number of configurations divided by τ [2,9]. Thus, low τ values reduce the statistical error of the thermal averages. τ is obtained from the equilibrium relaxa tion function φ MM (t) according to the equation [2] φ MM (t) = 1 where N mcs is the total number of MCS and t is the number of MCS employed for relaxation. The time dependence of τ for the different trial moves is shown in figure 7. These results present similar behavior to that of the first part of the test. At very low temperatures, the small step move is the most effi cient move along with the adaptive move, while the random move is the least efficient. The adaptive move is the only move that produces low relaxation times at all temperatures. Specially at low temperatures, the relaxation times of the adaptive move are several orders of magnitude smaller than those of the random and the combinational moves, which indi cates that it requires less MCS to produce results with similar statistics. As expected, all the relaxation times diverge at the critical temperature.
Independently of the temperature and the anisotropy values, the adaptive move is the most efficient move relaxing the system to equilibrium and producing low correlation times when the system is already in equilibrium.

Susceptibility: a slow convergence problem
The magnetic susceptibility is a measure of the thermal fluc tuations of the magnetization. In a cubic system with uniaxial anisotropy along z at a given temperature, it is expected that the x and y spatial components of the susceptibility (χ x and χ y , respectively) converge to the same value. However, in Monte The convergence to equilibrium of these thermal averages is different for each trial move, specially at low temperatures. The adaptive move is efficient both at low and high temperatures. Carlo simulations, these spatial components of the suscepti bility typically require a large number of MCS to converge, particularly at high temperatures below T c .
In this test we evaluated the convergence of χ x and χ y in a system with uniaxial anisotropy along z at k B T/J = 1.48. The spatial components χ x and χ y are given by We computed the spatial components at each MCS after relax ation taking averages of the magnetization data before that MCS, according to where t is the number of MCS, The time dependence of the spatial components χ x and χ y and their difference |χ x − χ y | is shown in figure 8. As expected, χ x and χ y tend to the same value as time increases for all the trial moves, and consequently |χ x − χ y | converges to zero. Because the system is at a high temperature below T c and the thermal energy is considerable, the spin moments are more disordered in the low (see figure 8(a)) than in the high anisotropy system (see figure 8(b)). Therefore, χ x and χ y conv erge slower to the same value in the low anisotropy system, specially using the small step and the combinational move. Although the random move presents similar behavior to that of the adaptive move, the convergence of |χ x − χ y | is faster when the adaptive move is employed for both low and high anisotropy systems.

Convergence of the coercivity
The calculation of hysteresis loops allows the characteriza tion of a great variety of magnetic properties of materials. Particularly, a great number of studies deal with calculations of the coercivity, mainly because high coercivity materials have several applications as permanent magnets. Monte Carlo simulations are not considered efficient for the sim ulation of hysteresis loops because the collective behavior of the interacting moments is not taken into account, there fore requiring a huge number of MCS [10]. Thus, proper ties like the coercivity can have a very slow convergence to equilibrium. Other integration methods such as LLG spin dynamics are more commonly employed to make this kind of calculations.
In this test, we estimated the convergence of the coercivity using both LLG spin dynamics and Monte Carlo simulations in a system with parameters similar to those of cobalt. We considered a magnetic moment µ s = 1.72µ B , where µ B is the Bohr magneton. The hyster esis loops were obtained at k B T/J Co ≈ 0.023 (T ≈ 10 K) for k v /J = 0.001 (k v = 0.037 85 meV/atom) and 1.0 (k v = 37.85 meV/atom). For the LLG spin dynamics simu lations, we considered a damping parameter λ = 1, time step ∆t = 0.3 fs and magnetic field steps µ s ∆B = 0.015 and 0.15 meV for the low and high anisotropy systems, respectively. The time step is in the stable region for the atomistic model at low temperatures (T T c ) and only valid for ferromag nets (see figure 2 in [1]). To calculate the coercivity, we made independent simulations of the hysteresis loop employing in each simulation a specific value of time steps (MCS and LLG time steps per field step for the Monte Carlo and LLG spin dynamics simulations, respectively). Then, the coercivity was calculated in each independent simulation.  Figure 9 shows the coercivity ( B c ) as a function of the time steps per field step for the Monte Carlo and LLG spin dynamics simulations. Independently of the integration method and the trial move used, the coercivity is expected to converge to the same value as the number of time steps increases. At both low (see figure 9(a)) and high ( figure  9(b)) anisotropies, the coercivity converges faster with the adaptive move and slower with the random move. Then, the adaptive move requires less runtime to simulate a complete hysteresis loop than the random and combinational moves, and the LLG spin dynamics. Furthermore, according to the results shown in section 3.1, the adaptive move is expected to be significantly more efficient than the random move at lower temperatures and still more efficient at higher temper atures below T c .

Magnetization reversal
Different mechanisms of magnetization reversal present in magnetic materials make them ideal for technological applica tions such as magnetic recording. In Monte Carlo simulations, the reversal mechanism and the time required to reverse the spin moments depend on the trial move employed.
In this final test, we simulated thermally activated reversal processes, as simulated by Hinzke and Nowak [5]. To begin the test, we considered a system with an initial spin moments configuration where all of them were pointing up and an external magnetic field pointing down. After some time, the external magnetic field energy overcomes the anisotropy and exchange energy and eventually the spin moments, which are in a metastable state, reverse their magnetization. The time required to produce a magnetization in z equal to zero, i.e. M z (τ ) = 0, is known as the metastable lifetime (τ ml ). For low temperatures, τ ml is given by [5] τ ml = τ 0 e ∆E k B T (10) where ∆E is the energy barrier which depends on the reversal mechanism. The temperature dependence of τ ml for the dif ferent trial moves is shown in figure 10. It is observed that the results are very similar independently of the anisotropy value of the system. Such as in the results of the integrated relax ation times (figure 7), the small step and the random moves are the most and the least efficient moves at low temperatures, respectively. While at high temperatures, the adaptive and random move are the most efficient and the small step move is the least efficient. Overall, the adaptive move presents a good efficiency at all temperatures, indicating that it would require less time to reverse the magnetization of the system. Although we performed the simulations on a simple cubic lattice with nearest neighbor interactions, the adaptive move can be applied to any lattice and long range interactions. There is no difference in terms of sample for the more general case, as each spin will have an effective exchange energy, and the spins simply sample the relevant phase space with an optimal sampling method. Moreover, in systems composed of two or more magn etic phases with distinct magnetic properties, such as hard/soft core/shell nanoparticles, using the adaptive move with a unique adaptive cone width could lead to a poor phase space sampling for some or all the phases due to the different optimal cone width for a given temperature. In this case, it is more suitable to use independent adaptive cone widths for each magnetic phase.

Conclusions
In summary, we have developed an adaptive algorithm for the optimal phase space sampling in Monte Carlo simulations of 3D Heisenberg spin systems. The proposed adaptive algo rithm modifies adaptively a conebased spin update method keeping the acceptance rate close to 50%. We have shown that the adaptive algorithm is more efficient than other common spin update methods independently of the temperature and the anisotropy of the system in consideration. Low correlation times and a faster convergence to equilibrium of the thermal averages and the coercivity were obtained when the adaptive algorithm was used. Also, the adaptive algorithm showed good efficiency for magnetization reversal at all temperatures, even for high anisotropy systems. A generalization of the adaptive algorithm could be made to enhance the efficiency of Monte Carlo simulations of other kinds of systems.