Accelerating spin-space sampling by auxiliary spin-dynamics and temperature-dependent spin-cluster expansion

TL;DR

This paper introduces a Hamiltonian Monte Carlo method utilizing auxiliary spin-dynamics and a temperature-dependent spin-cluster expansion to efficiently sample high-dimensional spin spaces in magnetic material simulations, compatible with electronic-structure methods.

Contribution

The authors develop a novel sampling framework that combines auxiliary spin-dynamics with a flexible, temperature-dependent spin-cluster expansion, improving efficiency in magnetic simulations.

Findings

Fast warm-up observed in classical Heisenberg model

Small dynamical critical exponent indicating efficiency

Successful application to magnetic phase transition in bcc iron

Abstract

Atomistic simulations of thermodynamic properties of magnetic materials rely on an accurate modelling of magnetic interactions and an efficient sampling of the high-dimensional spin space. Recent years have seen significant progress with a clear trend from model systems to material specific simulations that are usually based on electronic-structure methods. Here we develop a Hamiltonian Monte Carlo framework that makes use of auxiliary spin-dynamics and an auxiliary effective model, the temperature-dependent spin-cluster expansion, in order to efficiently sample the spin space. Our method does not require a specific form of the model and is suitable for simulations based on electronic-structure methods. We demonstrate fast warm-up and a reasonably small dynamical critical exponent of our sampler for the classical Heisenberg model. We further present an application of our method to the magnetic phase transition in bcc iron using magnetic bond-order potentials.

Equations

$$⟨O⟩=\frac{1}{Z}\int \prod _{i}d{\mathbf{s}}_i\pi (\{{\mathbf{s}}_i\})O(\{{\mathbf{s}}_i\}),$$

$$\pi (\{{\mathbf{s}}_i\})=\exp \left[-\frac{1}{k_{\text{B}}T}E(\{{\mathbf{s}}_i\})\right],$$

$$\pi (q,p)=\frac{1}{Q}\exp \left[-\frac{1}{k_{\text{B}}T}\mathcal{H}(q,p)\right],$$

$$\mathcal{H}=\frac{I}{2}\sum _i{\mathbf{\omega }}_i^T\cdot {\mathbf{\omega }}_i+E(\{{\mathbf{s}}_i\}),$$

$$\pi (\{{\mathbf{s}}_i\})=\int \prod _{i}d{\mathbf{\omega }}_i\pi (\{{\mathbf{\omega }}_i,{\mathbf{s}}_i\}),$$

$$\pi (\{{\mathbf{\omega }}_i,{\mathbf{s}}_i\})=\frac{1}{Q}\exp \left[-\frac{1}{k_{\text{B}}T}\mathcal{H}(\{{\mathbf{\omega }}_i,{\mathbf{s}}_i\})\right].$$

$$\frac{d{\mathbf{s}}_i}{dt}=\frac{{\gamma }_e}{m_i}\frac{\partial E}{\partial {\mathbf{s}}_i}\times {\mathbf{s}}_i,$$

$$I\frac{d{\mathbf{\omega }}_i}{dt}=\frac{\partial E}{\partial {\mathbf{s}}_i}\times {\mathbf{s}}_i,$$

$$\frac{d{\mathbf{s}}_i}{dt}={\mathbf{\omega }}_i\times {\mathbf{s}}_i.$$

$${\tilde{\mathbf{\omega }}}_i^{t+1}={\mathbf{\omega }}_i^t+\frac{1}{2}\frac{ϵ}{I}\cdot \frac{\partial E}{\partial {\mathbf{s}}_i^t}\times {\mathbf{s}}_i^t,$$

$${\mathbf{s}}_i^{t+1}=\mathbf{D}({\tilde{\mathbf{\omega }}}_i^{t+1},ϵ)\cdot {\mathbf{s}}_i^t,$$

$${\mathbf{\omega }}_i^{t+1}={\tilde{\mathbf{\omega }}}_i^{t+1}+\frac{1}{2}\frac{ϵ}{I}\cdot \frac{\partial E}{\partial {\mathbf{s}}_i^{t+1}}\times {\mathbf{s}}_i^{t+1},$$

$$\mathbf{D}({\tilde{\mathbf{\omega }}}_i^{t+1},ϵ)=\mathbf{I}+{\mathbf{W}}_i^{t+1}\sin ({\omega }_i^{t+1}ϵ)+({\mathbf{W}}_i^{t+1}{)}^2\left[1-\cos ({\omega }_i^{t+1}ϵ)\right].$$

$$p^{\mathrm{acc}}({\mathbf{x}}_{\mathrm{new}}|{\mathbf{x}}_{\mathrm{old}})=\min \left\{1,\exp \left(-\frac{\Delta {\mathcal{H}}_{\mathrm{old}\to \mathrm{new}}}{k_{\text{B}}T}\right)\right\},$$

$$p^{\mathrm{acc}}({\mathbf{x}}_{\mathrm{new}}|{\mathbf{x}}_{\mathrm{old}})$$

$$=\min \left[1,\exp \left(-\frac{\Delta {\mathcal{H}}_{\mathrm{old}\to \mathrm{new}}^{\mathrm{target}}}{k_{\text{B}}T}\right)\right]$$

$$=\min [1,\exp \left(\frac{\Delta {\mathcal{H}}_{\mathrm{old}}^{\mathrm{auxiliary}\to \mathrm{target}}}{k_{\text{B}}T}\right)\times$$

$$\exp \left(-\frac{\Delta {\mathcal{H}}_{\mathrm{old}\to \mathrm{new}}^{\mathrm{auxiliary}}}{k_{\text{B}}T}\right)\times$$

$$\exp \left(-\frac{\Delta {\mathcal{H}}_{\mathrm{new}}^{\mathrm{auxiliary}\to \mathrm{target}}}{k_{\text{B}}T}\right)].$$

$$E=-J\sum _{<ij>}{\mathbf{s}}_i\cdot {\mathbf{s}}_j,$$

$$\tau \approx L^z,$$

$$\varphi (t)=A{e}^{-t/\tau },$$

$$E(\{{\mathbf{r}}_i,{\mathbf{s}}_i\})=U_{\text{bond}}+U_{\text{trans}}+U_{\text{rep}}+U_{\text{C}}+U_{\text{X}}+U_{\text{ext}},$$

$${\bar{p}}^{\mathrm{acc}}\approx \frac{1}{2}\left[1+\exp \left(-\frac{\Delta E_{\mathrm{RMS}}}{k_{\mathrm{B}}T}\right)\right].$$

$${\Delta }_O=\sqrt{\frac{\mathrm{Var}O}{N}(1+2{\tau }_O)},$$

$$\Delta (t)=\frac{1}{2}\sum _i{\left[{\mathbf{s}}_i(t)-{\mathbf{s}}_i(0)\right]}^T\cdot \left[{\mathbf{s}}_i(t)-{\mathbf{s}}_i(0)\right].$$

$$\frac{d\Delta }{dt}<0,$$

$$\sum _i{\left[{\mathbf{s}}_i(t)-{\mathbf{s}}_i(0)\right]}^T\cdot [{\mathbf{\omega }}_i(t)\times {\mathbf{s}}_i(t)]<0,$$

Full text

I Introduction

Understanding and predicting thermodynamic properties of magnetic materials is not only of general scientific interest, but also crucial for many applications, e.g., magnetic refrigeration Pecharsky and Jr (1999), iron and steel design Hickel et al. (2012), and optimization of magnetic nanoparticles in biomedicineLu et al. (2007). In a spin system, the potential energy $E$ is defined as a function of the spin configuration $\{\mathbf{s}_{1},\mathbf{s}_{2}...\mathbf{s}_{N}\}$. Finite-temperature properties are then investigated using Boltzmann statistics,

[TABLE]

with the configurational probability density in spin space,

[TABLE]

where $k_{\text{B}}$ is the Boltzmann constant and $T$ the temperature. $O$ may be any spin-dependent observable, and $Z$ is the partition function in the classical limit and contains no longitudinal spin-fluctuationsDietermann et al. (2012); Ruban et al. (2007); Ma and Dudarev (2012a).

Frequently, model Hamiltonians are used for evaluating spin-dependent observables, such as the rigid-lattice Heisenberg model or its coordinate-dependent variants that account for spin-lattice coupling Ma et al. (2008). For materials design it would be desirable to work with material specific Hamiltonians that explicitly take into account the electronic structure. This would also facilitate an adequate treatment of spin fluctuations in itinerant-electron magnetsMoriya (1985), complex exchange-interactions in FeKvashnin et al. (2016) and magnon-phonon coupling in magnetic transition metalsKörmann et al. (2014). Recent years have seen significant progress in the development of electronic-structure based models, e.g., density functional theory for non-collinear magnetism Kubler et al. (1988); Stocks et al. (1998); Hobbs et al. (2000); Kurz et al. (2004); Peralta et al. (2007); Ma and Dudarev (2015), non-collinear magnetic tight-bindingMukherjee and Cohen (2001); Barreteau et al. (2016) or bond-order potentialsDrautz and Pettifor (2011); Ford et al. (2014, 2015); Drautz et al. (2015).

For model Hamiltonians there are many Monte Carlo (MC) sampling algorithms that efficiently sample the spin space, but unfortunately none of them is suitable for electronic-structure based models. We notice that there are three main differences between model Hamiltonians and electronic-structure based models. First, model Hamiltonians, at least most of them, contain only pair-wise interactions, while electronic-structure based models require many-spin interactions. Second, the range of the interaction is different. In model Hamiltonians typically only first and second nearest-neighbour interactions are taken into account, while in principle all spins are coupled in electronic-structure based models. Third, electronic-structure based models are orders of magnitude slower in the evaluation of the Hamiltonian. These differences prohibit application of many efficient MC algorithms. For example, the checkerboard MC algorithm Peczak and Landau (1990) is not applicable to electronic-structure based models as the system cannot be decomposed into non-interacting sublattices and the checkerboard-decomposition method is not applicable. The Swendsen-Wang Swendsen and Wang (1987) and WolffWolff (1989) cluster algorithms significantly reduce the correlation of samples, but they only work for models which may be mapped onto percolation models, which is difficult for electronic-structure based models. The over-relaxationBrown and Woch (1987); Creutz (1987) algorithm works well for the classical Heisenberg model, but its efficient implementation depends upon the checkerboard-decomposition method Landau and Binder (2005). The heat-bath spin dynamicsMa and Dudarev (2011); Tranchida et al. (2018) might be seen as a variant of the over-relaxation algorithm and suffers from the same problem. The Wang-Landau samplingWang and Landau (2001) can in principle overcome the critical slowing-down, but the convergence of the density of states usually requires millions of energy evaluations, which is not affordable for electronic-structure based methods. In a recent application of the Wang-Landau sampling to first-principle non-collinear magnetism 590,000 energy evaluations were performed yielding a reasonable estimate of the density of statesEisenbach et al. (2011), but despite the huge computational effort the results could not be fully converged.

Our work builds on recent progress in the Hamiltonian Monte Carlo (HMC)Neal (2011); Mark and Ben (2011); Beskos et al. (2013); Betancourt (2017); Betancourt et al. (2017); Wang et al. (2013); Hoffman and Gelman (2014) and efficient methods to accelerate first-principles thermodynamic calculationsHellman et al. (2013); Grabowski et al. (2009); Duff et al. (2015). For the former, a rigorous theoretical proof has been given that underpins the empirical success of HMC Betancourt et al. (2017) and the theoretical considerations Beskos et al. (2013) and techniquesNeal (2011); Betancourt (2017); Wang et al. (2013); Hoffman and Gelman (2014) developed for the automatic tuning of its hyper-parameters. In the latter, effective potentials are employed to speed up first-principles thermodynamic calculations, as the direct calculations with first-principles methods are too expensive. However, the conventional HMC algorithm is not applicable to spin systems as the spin length is not preserved in standard molecular dynamics that is used in HMC. The effective-potential methods in literature Hellman et al. (2013); Grabowski et al. (2009); Duff et al. (2015) are designed to accelerate calculations of free energies and cannot be applied to evaluate other thermodynamic quantities straightforwardly. These considerations form the basis for the methods developed in this work. First, we propose an auxiliary spin-dynamics as a basis for a HMC algorithm for spin systems that rigorously preserve spin lengths. Second, we propose a HMC framework in which the temperature-dependent spin-cluster expansion (SCE)Drautz and Fähnle (2004, 2005) is used as an auxiliary model to further accelerate the sampling of the spin space.

The paper is structured as follows: we first introduce the HMC algorithm for spin systems using auxiliary spin-dynamics, and discuss the automatic tuning of its hyper-parameters. Then we introduce the temperature-dependent SCE as an auxiliary model to accelerate the sampling of the spin space. In Sec. III, we employ the classical Heisenberg model to demonstrate the efficiency of our method, and apply our algorithm to sample the magnetic phase transition in bcc iron with magnetic bond-order potentialsDrautz and Pettifor (2011); Mrovec et al. (2011); Ford et al. (2015).

II Methodology

Our target is to draw efficiently independent samples according to the configurational probability density defined in spin space. To this end we extend HMC Neal (2011); Betancourt et al. (2017) for the sampling of spin space variables. HMC does not sample the configurational distribution directly but a joint distribution of positions $q$ and momenta $p$,

[TABLE]

where $Q$ is the partition function in the phase space. The marginal distribution of $q$ then restores the target distribution.

For sampling spin space using HMC we introduce auxiliary spin angular velocities $\bm{\omega}_{i}$ as canonical variables of the spin directions $\{\bm{\omega},\bm{s}\}$ in formal analogy to the classical canonical variables $\{\bm{p},\bm{q}\}$ and define the Hamiltonian as

[TABLE]

where $I$ is a fictitious mass that later is used as a parameter to optimize the efficiency of the HMC sampling.

The configurational probability density in spin space, Eq. (2), is restored by the marginal distribution of the configurational variables in phase space,

[TABLE]

with the joint probability density defined as

[TABLE]

We may thus sample the joint probability density Eq. (6) and then obtain the configurational probability density from Eq. (5), which leads to the correct thermal average Eq. (1). For sampling Eq. (6), we employ Hamiltonian dynamics together with Monte Carlo sampling, as outlined in the following.

Different from usual spin dynamics, which is based on a first-order differential equation in timeAntropov et al. (1996); Gilbert (2004); Ma et al. (2008)

[TABLE]

where $m_{i}$ is the magnitude of magnetic moments, and $\gamma_{e}$ is the gyromagnetic ratio for an electron spin, the auxiliary Hamiltonian, Eq. (4), dictates that the spins follow conventional Hamiltonian dynamics for rigid bodies

[TABLE]

It is evident that the spin dynamics described by Eq. (7) and Eq. (8) is different. Eq. (8) describes a completely fictitious dynamics of the spins that, however, by construction may be used to sample spin space according to the probability density $\pi(\{\mathbf{s}_{i}\})$. We denote Eq. (7) as the semi-classical spin dynamics and Eq. (8) as the auxiliary spin-dynamics in this paper.

While it is very difficult or impossible to find an effective symplectic time-reversible integrator for the semi-classical spin dynamics McLachlan et al. (2014); Krech et al. (1998), this is not a problem for the auxiliary spin-dynamics. We combine the velocity Verlet method Swope et al. (1982) and a spin rotation scheme to obtain an efficient numerical integration algorithm that is time-reversible, area-preserving and preserves spin length. There are three steps per update,

[TABLE]

where the index $i$ denotes spin, $t$ the current state and $t+1$ the next state. $\epsilon$ is the time step, and $\mathbf{D}(\bm{\tilde{\omega}}^{t+1}_{i},\epsilon)$ is a $3\times 3$ rotation matrix Omelyan et al. (2001),

[TABLE]

$\omega^{t+1}_{i}$ is the magnitude of $\bm{\tilde{\omega}}^{t+1}_{i}$, and $\mathbf{W}^{t+1}_{i}$ is a skew - symmetric matrix with $\mathbf{W}^{t+1}_{i,XY}=-\hat{\omega}^{t+1}_{i,Z}$, $\mathbf{W}^{t+1}_{i,XZ}=\hat{\omega}^{t+1}_{i,Y}$, $\mathbf{W}^{t+1}_{i,YZ}=-\hat{\omega}^{t+1}_{i,X}$, where ${\hat{\bm{\omega}}}^{t+1}_{i}$ is the directional vector of $\bm{\tilde{\omega}}^{t+1}_{i}$. $X$, $Y$ and $Z$ denote Cartesian components.

We may now employ standard HMC for sampling the spin space probability density $\pi(\{\mathbf{s}_{i}\})$ on the basis of the auxiliary spin-dynamics. There are three steps per MC update. The first step performs a Gibbs sampling of angular velocities, in which we fix the configurational variables $\{\mathbf{s}_{i}\}$ and sample the angular velocities according to their conditional, Gaussian distribution. In the second step the auxiliary spin-dynamics is run for a specific trajectory of length $L$. In the third step the Metropolis-Hastings acceptance-rejection is performed for the proposal state generated by auxiliary spin-dynamics in order to guarantee detailed balance. The acceptance probability is given by

[TABLE]

where $\mathbf{x}_{\mathrm{old}}$ and $\mathbf{x}_{\mathrm{new}}$ are the state variables containing both the spins and the angular velocities. Since in Hamiltonian dynamics the energy is conserved, the Metropolis-Hastings acceptance criterion only corrects for changes in energy due to numerical integration errors, resulting in high acceptance ratios. Repeating the three steps above leads to a Markov chain which obeys the joint distribution and is used to evaluate thermal averages.

There are three hyper-parameters in our algorithm, the mass $I$, the time-step $\epsilon$ for numerical integration and the trajectory length $L$. They have no effect on accuracy of the MC sampling but strongly influence the efficiency and in practical applications should be set automatically without user intervention. We adapt the time step such that the exponential moving average of the acceptance probability is in the range from 0.6 to 0.7, which is centered around the optimal value 0.651 suggested by Beskos et al.Beskos et al. (2013) More specifically, if the acceptance probability is smaller (larger) than 0.6 (0.7), the time step is decreased (increased) by a specific factor. For the tuning of the trajectory length $L$, we employ an empirical termination criterion, the U-turn terminationHoffman and Gelman (2014), as given in Appendix A. The hyper-parameters are only tuned in the warm-up phase and then fixed to leave the distribution function invariant in the sampling phase. We fix the mass $I$ to $1~{}\text{eV}\text{fs}^{2}$ in our work, which is an empirical optimal value according to our tests.

We denote the HMC based on auxiliary spin-dynamics as Algorithm I in this paper. We next incorporate the temperature-dependent SCE into Algorithm I to further accelerate the sampling for expensive spin models and denote the new method Algorithm II. This is based on the observation that the auxiliary spin-dynamics is only used to generate the proposal state and may be run with a cheaper auxiliary model instead of the expensive target one. This leaves the sampling correct as long as the detailed balance is guaranteed for the target model through the Metropolis-Hastings acceptance-rejection. The criterion for the auxiliary model is that it should be as close to the target one as possible since the acceptance probability is now determined by the difference between auxiliary and target models,

[TABLE]

In Eq. (12), the first and last exponentials contain the energy difference between the auxiliary and the target model for the old and new states, respectively. The second exponential arises from the numerical error of the integration of the auxiliary spin-dynamics with the auxiliary model, which is typically a small contribution.

For spin systems, an ideal auxiliary model is the spin-cluster expansionDrautz and Fähnle (2004, 2005) which may be fitted to accurately reproduce the target model. Here we propose to generate temperature-dependent SCEs due to two considerations. First, only a specific area in the spin space is explored with high probability at a given temperature, as shown in Fig. 1,

where we plot the magnetization and the potential energy of spin configurations at different temperatures for a classical ferromagnetic Heisenberg model. Clearly, the configurations at different temperatures are clustered into specific areas. This locality makes it easier to fit a temperature-dependent effective potential. Second, the temperature dependence of the exchange parameters is inherited in some models, e.g., in spin-density-functional-theory for itinerant-electron magnetsRosengaard and Johansson (1997). In practice, we collect the spin configurations in the warm-up phase and fit the temperature-dependent SCE at different temperatures separately. The extra cost required for fitting the temperature-dependent SCEs is normally negligible compared to the total computational time for sampling. The auxiliary spin-dynamics is then run with the temperature-dependent SCE to generate proposal states while the Metropolis-Hastings acceptance-rejection is performed for the target models in order to guarantee correct sampling. Implementation details are given in Appendix B. We note that no gradient calculations for the electronic-structure based models are required in this algorithm as the auxiliary spin-dynamics is run with the temperature-dependent SCE, which is another advantage of our method since the evaluation of gradients usually requires considerable extra computational cost.

III Results and Discussion

We demonstrate two applications in this section. In the first application, we perform efficiency tests for Algorithm I using the classical Heisenberg model. In the second application, we test and discuss Algorithm II using magnetic bond-order potentialsDrautz and Pettifor (2011); Ford et al. (2014, 2015); Drautz et al. (2015).

III.1 Application to the classical Heisenberg model

We first employ the classical ferromagnetic Heisenberg model,

[TABLE]

on a simple cubic lattice to perform an efficiency test for HMC, Algorithm I. $J$ is the exchange parameter, and $<ij>$ denotes the first-nearest-neighbour pairs without double counting. $J$ is chosen to be $1~{}\text{eV}$ in this test.

In Fig. 2, we measure the warm-up efficiency for a $6\times 6\times 6$ simple cubic lattice from the high (low)-energy to low (high)-energy states.

As expected, the efficiency from low-energy to high-energy states is higher than the other way round, and the former (latter) takes around 400 (1200) gradient calls.

We next fix the temperature to the critical temperature of the 3D classical Heisenberg model ($T_{\mathrm{C}}=1.4459~{}J/k_{\mathrm{B}}$) Peczak and Landau (1990), and estimate the dynamical critical exponent. The estimation is based on the dynamical finite-size scaling ansatzSuzuki (1977)

[TABLE]

where $L$ is the side length of the cubic simulation cell. $\tau$ is the relaxation time of magnetization and evaluated according to

[TABLE]

where $\phi(t)$ is the time auto-correlation function of magnetization. The relaxation time $\tau$ is measured in units of gradient calls instead of MC step as multiple gradient calls may be required for one MC step. In Fig. 3,

we show the log-log plot for the relaxation time versus the side length of the cubic simulation cell. The dynamical critical exponent is estimated to be around 2.23. This value is comparable to that of the checkerboard MC algorithm (1.96, cf. Ref. Peczak and Landau, 1990), but cannot compete with the cluster algorithms whose relaxation times are almost independent of the size Holm and Janke (1993). However, as we discussed in the introduction, these algorithms are only applicable for a small group of spin models, whereas HMC is generally applicable.

III.2 Application to electronic-structure based models

Next, we employ the magnetic bond-order potential (BOP)Drautz and Pettifor (2011); Mrovec et al. (2011); Ford et al. (2015) to demonstrate the application of our algorithm with an auxiliary model, Algorithm II. In algorithm I, we need to run both MC and auxiliary spin-dynamics with the target model, which requires too many energy and gradients calls to converge thermal averages for electronic-structure based models. A trajectory length of more than ten is usually needed in auxiliary spin-dynamics in order to decorrelate the current and the proposal state. The computational cost of this part can be dramatically reduced with algorithm II in which the auxiliary spin-dynamics is run with a temperature-dependent SCE. As the temperature-dependent SCE is orders of magnitude faster than the electronic-structure based models, this gives a significant speed-up.

The magnetic BOP is one of the simplest electronic-structure based models for magnetic transition metals. In this model, the potential energy is based on the electronic density of states. It is given as a function of atomic positions and spin orientations,

[TABLE]

where $U_{\text{bond}}$ is the bond energy, $U_{\text{trans}}$ the electron transfer energy, $U_{\text{rep}}$ the repulsion energy, $U_{\text{C}}$ the Coulomb energy, $U_{\text{X}}$ the exchange energy, and $U_{\text{ext}}$ the external energy. Readers are referred to Ref. Drautz and Pettifor, 2011 for a detailed discussion of this model.

The magnetic BOP gives a robust description of ferromagnetism, and more importantly, the real-material properties such as phase stability, elastic constants, and dislocations are described properly Mrovec et al. (2011); Ford et al. (2015, 2014). More specifically, we use the 9-moments magnetic BOP of iron by Mrovec et al. Mrovec et al. (2011) and the implementation in the BOPfox code Hammerschmidt et al. (2019) to calculate potential energies for spin configurations in the bcc lattice. The $5\times 5\times 5$ simulation cell contains 250 spins, and the atomic positions are fixed.

We employ a temperature-dependent SCE fitted for the BOP of iron as an auxiliary model to run the auxiliary spin-dynamics. As shown in Eq. (12), the energy difference between the electronic-structure based model and the auxiliary model should be as small as possible to guarantee a good acceptance probability. A rough estimation can be given for the relationship between the root-mean-squared (RMS) error $\Delta E_{\mathrm{RMS}}$ and the average acceptance probability $\bar{p}^{acc}$,

[TABLE]

For example, a RMS error of one $k_{\mathrm{B}}T$ corresponds to an average acceptance probability of around 0.68, which is a good value for MC acceptance-rejection. In Fig. 4,

we show a typical plot for the RMS error versus neighbour shells of pair-wise and three-spin interactions. The pair-wise spin clusters are not sufficient to reduce the root-mean-squared (RMS) error to less than one $k_{\mathrm{B}}T_{\mathrm{C}}$ and three-spin interactions are taken into account in order to further reduce the RMS error to 0.80 $k_{\mathrm{B}}T_{\mathrm{C}}$. In practice, we include pair-wise interactions up to the sixth nearest-neighbour shell and the first-nearest-neighbour three-spin interactions in our temperature-dependent SCE. The two-spin and three-spin interactions are sufficient to converge the SCE for the ideal lattice with fixed atomic positions that is used in this work. Breaking the geometric degeneracy by, e.g., vibrations or defects, will introduce distance- and environment-dependent interaction terms which rule the convergence of the SCE much more difficult or even intractable.

As discussed in Appendix B, we employ the temperature-independent and temperature-dependent SCE in the warm-up phase and sampling phase, respectively. In Fig. 5,

we show the average acceptance probability at different temperatures in both the warm-up and sampling phase. It is higher in the sampling phase than in the warm-up phase as the temperature-dependent SCE reproduces the energy landscape of the BOP better, as discussed in Sec. II.

We use 10,000 MC steps to evaluate thermal averages, and perform a binning analysis Ambegaokar and Troyer (2010) for the potential energy to check the convergence. As shown in Fig. 6, the autocorrelation times at all temperatures reach a plateau, which indicates that the calculation is fully converged. The autocorrelation times range from two to six depending on the temperature. The running error can be estimated according toAmbegaokar and Troyer (2010)

[TABLE]

where $O$ is the observable, $\mathrm{Var}O$ is its variance, $\tau_{O}$ is its autocorrelation time, and $N$ is the number of MC steps. Based on this, we plot in Fig. 7 the running errors in the estimation of thermal averages of potential energies at different temperatures. Clearly, more MC steps are required to reach convergence at temperatures closer to the critical point (1200 K in this case, as shown later) as the variance is larger. For our algorithm, an error of 0.5 meV/atom can be reached within 5000 MC steps. This excellent efficiency is due to the short autocorrelation time, as already discussed.

The magnetization and the magnetic contribution to the specific heat of bcc iron for magnetic BOP are plotted as a function of temperature in Fig. 8 and Fig. 9, respectively.

The magnetization curve we obtain here has similar features to those of the classical Heisenberg model Ma and Dudarev (2012b); Körmann et al. (2011a), which indicates that the exchange parameters are not influenced considerably by spin fluctuations in magnetic iron. We notice that there is a recent paper by Ruban and Peil Ruban and Peil (2018) who found that the exchange parameters are significantly influenced by atomic vibrations in magnetic iron. This effect is not considered here as we fix atomic positions in our simulation. There is residual magnetization up to 2000 K due to the finite-size effect. Based on the magnetic contribution to the specific heat (see Fig. 9), we estimate the Curie temperature to be around 1200 K, which is close to the experimental value (1043 K). The difference in $T_{\mathrm{C}}$ is attributed to the neglect of magnon-phonon coupling Ruban and Peil (2018) and limitations of the parametrization of the magnetic BOP used here. To our knowledge, this is the first result for the direct and high-fidelity simulations of the magnetic phase transition of bcc iron with an electronic-structure based model, in contrast to related works in literature that are based on parametrized model Hamiltonians Körmann et al. (2011b); Ma and Dudarev (2012b); Ruban and Peil (2018).

IV Conclusions

We have developed a Hamiltonian Monte Carlo (HMC) framework in order to efficiently sample the spin space for electronic-structure based models of magnetic materials. From auxiliary spin-dynamics we derive a HMC algorithm for spin systems, which is denoted as Algorithm I in this paper. The hyper-parameters (the time step and the trajectory length) are tuned with automatic tuning procedures without user intervention. Our tests with the classical Heisenberg model show that with the automatic tunning procedures of hyper-parameters this algorithm has a fast warm-up efficiency. The dynamical critical exponent is estimated to be $2.23\pm 0.02$, close to that of the checkerboard MC algorithm. While the latter is only applicable for a limited group of spin models, our algorithm is of general applicability. The utilization of the temperature-dependent spin-cluster expansion as an auxiliary model to run auxiliary spin-dynamics further accelerates the exploration of the spin space, and based on this we develop Algorithm II in this paper. Our application employing the magnetic bond-order potentials demonstrates the efficiency of our sampler. For a 250-atom supercell of bcc iron, the autocorrelation time of the potential energy is less than six MC steps, and the thermal average can be converged with an error of $\pm 0.5$ meV/atom within 6000 steps in the whole temperature range. We compute the magnetization curve, the magnetic contribution to the specific heat, and the Curie temperature of bcc iron as predicted by the BOP model with high fidelity. We conclude that this work paves the way towards atomistic simulations of magnetic materials with complex spin interactions, and look forward to seeing applications of our method for more complex models, such as non-collinear magnetic density-functional theory.

ACKNOWLEDGEMENT

NW acknowledges a PhD fellowship from the International Max Planck Research School for Surface and Interface Engineering (IMPRS-SurMat). The authors thank Martin Staadt and Tilmann Hickel for fruitful discussions.

APPENDIX

Appendix A U-turn termination criterion

The U-turn termination criterionHoffman and Gelman (2014) is an empirical estimate for the optimal length of Hamiltonian dynamics per MC step. The basic idea is to maximize the squared distance between the initial and final states. Its implementation in this work is sightly different from the one in the original paper as we do not use standard molecular dynamics.

We first define the half squared distance for spin systems,

[TABLE]

The U-turn termination criterion is then derived according to,

[TABLE]

and given as,

[TABLE]

where the equations of motion, Eq. (8), are employed.

Appendix B Implementation details

We first fit a temperature-independent spin-cluster expansion as a start-up auxiliary model. This involves 1000 spin configurations generated with the classical Heisenberg model at the critical point and the corresponding potential energies for electronic-structure based models (magnetic bond-order potential in our application). The warm-up phase is split into two stages. In the first stage, the system is thermalized purely with the start-up auxiliary model and the length is set to 1000 MC steps. The time-step and the trajectory length are automatically tuned with the methods in Section II. In the second stage, the start-up model is used to run auxiliary spin-dynamics while the Metropolis-Hastings acceptance-rejection is performed for the target electronic structure based models. The time-step is further tuned while the trajectory length is fixed to the mean value in the first stage. In our implementation for the magnetic bond-order potential, a typical time step is around 0.3 fs, and a typical trajectory length is around 18. The length of the second stage is set to be 500 MC steps. Once the warm-up phase is finished, we collect all the spin configurations and the potential energies in the second stage to fit a temperature-dependent SCE with least-squares fitting, which is then used to run auxiliary spin-dynamics in the sampling phase. Both the time-step and the trajectory length are fixed in this phase. We use 10,000 MC steps to evaluate thermal averages, which are sufficient to guarantee convergence due to the small autocorrelation time in our algorithm.