Macroscopically constrained Wang-Landau method for systems with multiple order parameters and its application to drawing complex phase diagrams

TL;DR

This paper introduces a generalized, macroscopically constrained Wang-Landau method for efficiently computing the density of states in systems with multiple order parameters, enabling rapid phase diagram analysis.

Contribution

It develops a novel approach that decomposes multidimensional random walks into constrained one-dimensional walks, facilitating quick thermodynamic calculations across phase diagrams.

Findings

Successfully applied to a spin-crossover model with complex interactions.

Enabled rapid computation of phase diagrams and order-parameter distributions.

Demonstrated efficiency in systems with multiple macroscopic parameters.

Abstract

A generalized approach to Wang-Landau simulations, macroscopically constrained Wang-Landau, is proposed to simulate the density of states of a system with multiple macroscopic order parameters. The method breaks a multidimensional random-walk process in phase space into many separate, one-dimensional random-walk processes in well-defined subspaces. Each of these random walks is constrained to a different set of values of the macroscopic order parameters. When the multi-variable density of states is obtained for one set of values of field-like model parameters, the density of states for any other values of these parameters can be obtained by a simple transformation of the total system energy. All thermodynamic quantities of the system can then be rapidly calculated at any point in the phase diagram. We demonstrate how to use the multi-variable density of states to draw the phase diagram, as well as order-parameter probability distributions at specific phase points, for a model spin-crossover material: an antiferromagnetic Ising model with ferromagnetic long-range interactions. The field-like parameters in this model are an effective magnetic field and the strength of the long-range interaction.

Tables

Table 1. Table 1: Crude tests for the computational time to obtain g ​ ( E , M , M s ) 𝑔 𝐸 𝑀 subscript 𝑀 𝑠 g(E,M,M_{s}) , using WL with multi-dimensional random walk in phase space, with energy spectrum partitioned into windows plus replica-exchange Valentim  et al. ( 2015 ) . Parallel programming was performed with each window assigned 1 core and 5 walkers. The histogram is considered to be ‘flat’ if all the walkers satisfy the root-mean-square ‘flatness’ criterion Brown ( 2014 ) (Eq. ( 2 )). The time recorded includes the time for initializing the systems. The asterisks mark the cases in which too many windows were used for a small system. In general, changing the system size from L 𝐿 L to L + 2 𝐿 2 L+2 causes the simulation time to increase by a factor of 6 to 10.

$L\times L$	3 windows	5 windows	7 windows
$6\times 6$	$38$ s	$33$ s	*$56$ s
$8\times 8$	$636$ s	$430$ s	*$1128$ s
$10\times 10$	$6572$ s	$4541$ s	$2720$ s
$12\times 12$	$52728$ s	$29312$ s	$22052$ s

Equations

$$\mathcal{H}=J\sum _{⟨i,j⟩}{s}_i{s}_j-\frac{A}{2N}{M}^2-HM.$$

$$\sqrt{\frac{{\sum }_{E,M,M_s}|H(E,M,M_s)-H_{\mathrm{average}}{|}^2}{N_{E,M,M_s}}}<20\%,$$

$$g(E,M,M_s)=\frac{H(E,M,M_s)}{{\sum }_{M,M_s}H(E,M,M_s)}g(E).$$

$$g(V_2,V_1)$$

$$g(V_3,V_2,V_1)$$

$$⋮$$

$$g(V_K,...,V_1)$$

$$g(E,M,M_s)=\frac{g(E|M,M_s)}{{\sum }_{E}g(E|M,M_s)}g(M,M_s).$$

$$M$$

$$M_s$$

$$g(E,M,M_s)\to g(E-HM-\frac{A{M}^2}{2N},M,M_s).$$

$$M$$

$$M_s$$

$$g(E,M,M_s)=g(-E,-M_A+M_B,-M_A-M_B).$$

$$\sqrt{\frac{{\sum }_E|H(E)-H_{\mathrm{average}}{|}^2}{N_E}}<20\%,$$

$$Z_{m,m_s}=\sum _{E}g(E,m,m_s)e^{-E/T}.$$

$$Z_{\mathrm{all}}=\sum _{m,m_s}{Z}_{m,m_s}.$$

$$P(m,m_s)\Delta m\Delta m_s=\frac{Z_{m,m_s}}{Z_{\mathrm{all}}},$$

$$F(m,m_s)=-T\ln Z_{m,m_s}.$$

$$P(m)\Delta m$$

$$P(m_s)\Delta m_s$$

$$u_p=1-\frac{⟨(p-⟨p⟩{)}^4⟩}{3⟨(p-⟨p⟩{)}^2{⟩}^2},$$

$$var(p)={\chi }_{p}T=L^2(⟨p^2⟩-⟨p{⟩}^2).$$

$$F(m)$$

$$F(m_s)$$

$$M$$

$$M_s$$

$$g(M,M_s)$$

$$=$$

$$C_k^n=C_{k-1}^{n-1}+C_k^{n-1},1\le k\le n-1$$

$$C_0^n=C_n^n=1$$

$$\ln n!\approx n\ln n-n.$$

$$M$$

$$M_s$$

$$g(E,M,M_s)=g(-E,-M_A+M_B,-M_A-M_B).$$

Full text

I introduction

Classical spin models have found wide application to a vast array of problems in many branches of physics and other sciences. This is due to the relative simplicity of such models and the fact that the different spin states can be given many interpretations besides that of magnetic spins, including different kinds of atoms or molecules (“lattice-gas models”), opinions, biological species in an ecosystem, etc. A few examples of these diverse applications are magnetic materials Richards et al. (1995), high-energy physics Caselle et al. (2003), astrophysics Hasnaoui and Piekarewicz (2013), electrochemistry Brown et al. (1999), polymer science Luo and Huang (2003), network reliability problems Ren et al. (2016), and economics Sornette (2014). The archetypal member of this class of models is the $S=1/2$ ferromagnetic Ising model of binary spin variables placed at the sites of a lattice or a more general network and interacting via a simple Hamiltonian. Since its introduction for a one-dimensional system almost a century ago Ising (1925), this model has been joined by many antiferromagnetic or ferromagnetic generalizations to higher spatial dimensions, multiple local states and/or multidimensional order-parameter spaces, such as the $S=1$ or three-state Blume-Capel Blume (1966); Capel (1966) and Blume-Emery-Griffiths models Blume et al. (1971). Although the order parameter for the $S=1/2$ square-lattice Ising model in zero field has been obtained exactly Onsager (1944); Yang (1952), solution of this class of models under general conditions, including in nonzero field, is known to be NP hard Barahona (1982). As a consequence, much effort has been applied to developing accurate approximate and numerical solutions, including mean-field approximations Bragg and Williams (1934, 1935); Agra et al. (2006), series expansions Domb (1974), numerical transfer-matrix calculations Nightingale (1990), and a variety of Monte Carlo (MC) methods Landau and Binder (2015). Although many ingenious algorithms have been introduced, development of improved numerical methods to study equilibrium and nonequilibrium aspects of classical spin systems remains an active research area.

In this paper we present a generalization of the Wang-Landau (WL) MC method for calculating densities of states (DOS) Wang and Landau (2001a, b) to systems with multiple order parameters. Development of the method was inspired by applications to a class of molecular crystals known as spin-crossover materials Halcrow (editor), some of which can have competing antiferromagnetic-like short-range and ferromagnetic-like long-range interactions. A discrete-spin model of such a system was recently studied for a few values of two field-like model parameters (an effective external field and the strength of the long-range interaction) by a computationally intensive Metropolis importance-sampling MC method Rikvold et al. (2016). In order to obtain results for a wide range of model parameters with a manageable computational effort, a simulation method is needed that can produce three-dimensional DOS, $g(E,M,M_{s})$, where $E$ is the total system energy, $M$ can be interpreted as a total system magnetization, and $M_{s}$ as a staggered magnetization Endnote1 . The method can also in principle deal with higher-dimensional order-parameter spaces. The original time-consuming WL random-walk simulation in multiple-dimensional phase space is broken down into many stages. In each stage, many independent WL simulations perform one-dimensional walks, each with different constrained macroscopic parameters. For the lattice model we consider in this paper, exact combinatorial calculations can be applied to simplify the process, so that only one stage of simulation in $E$ is required. As the simulations are run independently, no special skills in parallel programming are required. From the one-dimensional random-walk WL simulations in $E$, performed separately over a grid in the order-parameter space at one single set of model parameters, the method can produce DOS for any value of model parameters and temperatures, using a simple transformation of the total system energy. This contrasts with both importance-sampling and the original WL MC methods, in which separate simulations must be performed for each set of model parameters of interest. From the resulting multidimensional DOS, properties such as phase diagrams, free-energy landscapes, and joint and marginal order-parameter distributions can be simply obtained.

Another advantage of this method is the ability to use symmetries in the order-parameter space to reduce the number of simulations needed. (For instance, in the case of the spin-crossover model, such symmetry considerations lead to an additional eight-fold reduction in the computational work.) For concreteness, the details of the method will be demonstrated here in the context of the model spin-crossover material of Ref. Rikvold et al. (2016). Further results for several sets of model parameters of experimental interest will be described in forthcoming papers Chan et al. (a, a).

The remainder of this paper is organized as follows. In Sec. II we introduce the model spin-crossover material Hamiltonian that inspired the method, and which we will use to illustrate the application of our algorithm. In Sec. III we first summarize the relevant basics of the WL algorithm, and then we discuss several ways of finding the joint DOS, $g(E,M)$, previously introduced in the literature, pointing out their weaknesses when applied to the situation studied here. In Sec. IV we discuss the macroscopically constrained WL algorithm in detail. Some calculations and symmetry considerations are discussed in Appendices A and B, and the detailed implementation of the method is described in Appendix C. In Sec. V we give numerical results for the model spin-crossover material, and we demonstrate how to use $g(E,M,M_{s})$ to draw and investigate its phase diagram. A discussion of methods to extend the system size is given in Sec. VI, and conclusions and a brief discussion of future work are given in Sec. VII.

II Two-dimensional Ising model with Antiferromagnetic Short-range and Ferromagnetic Long-range interactions (2D Ising-ASFL model)

To demonstrate the details and performance of the proposed method, and for its comparison with other methods, we will use a pseudo-spin model of a spin-crossover material with short-range antiferromagnetic-like and long-range ferromagnetic-like interactions, which was previously introduced and studied by Metropolis importance-sampling MC in Refs. Brown et al. (2014); Rikvold et al. (2016). It is defined by the Hamiltonian,

[TABLE]

The local variables are $s_{i}=\pm 1$, and $M=\sum_{i}s_{i}$ is the corresponding global “magnetization.” The explicit sum runs over all nearest-neighbor pairs on an $N=L\times L$ square lattice with periodic boundary conditions, and $J>0$ makes the local interactions antiferromagnetic. The second term models long-range, mean-field like ferromagnetic interactions of strength $A\geq 0$. In the third term, $H$ is the applied field (actually an effective field in the spin-crossover model Wajnflasz and Pick (1971); Miyashita et al. (2008); Nakada et al. (2011); Rikvold et al. (2016)), which breaks the symmetry between positive and negative $M$. $A$ and $H$ are the model’s two field-like parameters. Throughout this paper, we will use the notation $E$ to represent the total energy obtained through this Hamiltonian, including the contributions from the terms proportional to $A$ and $H$. $E$, $A$, and $H$ will be given in units of $J$, and temperatures in units of $J/k_{\rm B}$ with $k_{\rm B}$ being Boltzmann’s constant.

While this Hamiltonian includes only the ferromagnetic order parameter $M$, we note that both linear and nonlinear terms (in this case $M^{2}$) are included. Linear and/or nonlinear terms in the staggered magnetization $M_{s}$ (see below) could also be added as needed to model other particular systems. However, the form given here is sufficient to demonstrate and validate the algorithm.

III Basic Wang-Landau and current methods to find joint densities of states

In this section, we review some relevant basics of the WL method, followed by some current methods to obtain joint DOS, and we explain why these methods are not appropriate to obtain $g(E,M,M_{s})$ for the system defined by the Hamiltonian (1). These are mostly WL based methods.

III.1 Basic Wang-Landau Monte Carlo method

The WL method is a restricted random-walk method for finding the DOS of a system. Its idea is based on the observation that if one imposes an acceptance probability for proposed energy transitions in the random walk which is proportional to the reciprocal of the DOS, then the system will spend roughly equal times in all different energy states. If a histogram $H(E)$ is used to record the number of times that the walker visits each energy state, a ‘flat’ histogram will eventually be generated. The whole random walk process is divided into many sweeps. In each sweep, the estimated DOS is multiplied by a modification factor, $f>1$. The sweep is finished when a ‘flat’ histogram is obtained. Then, the next sweep starts with a smaller value of $f$. The whole process ends when $\ln f<10^{-8}$. To speed up the simulation, a wide energy spectrum may be divided into separate energy windows, with separate simulations performed in each window. The partitioned windows may be uniform Wang and Landau (2001a, b) or nonuniform Cunha-Netto et al. (2008). As errors are generated near the window boundaries Schulz et al. (2003); Wang and Landau (2001b) whenever proposed moves outside the window are rejected, neighboring energy windows should overlap by a certain fraction, and the estimated $g(E)$ obtained in neighboring windows should be joined at the point where their slopes with respect to $E$ are closest, so that a smooth $g(E)$ over the whole energy spectrum can be obtained. In complex systems with rugged energy landscapes, states that lie in the same energy window may sometimes be connected only by paths that go via a different window. Therefore, the Replica Exchange Wang-Landau (REWL) scheme Vogel et al. (2013); Li et al. (2014); Vogel et al. (2014) was proposed to ensure that all microstates are visited. The scheme allows two walkers that both have energies within the overlap region of adjacent windows to exchange their microstates with a certain probability, so that ergodicity is preserved. REWL is performed in parallel Yin and Landau (2012). Errors and convergence of WL have also been studied Zhou and Bhatt (2005); Lee et al. (2006); Belardinelli and Pereyra (2007a, b); Brown et al. (2011); Komura and Okabe (2012). It was found that the statistical errors in $\ln g(E)$ are proportional to $\sqrt{\ln f}$ Zhou and Bhatt (2005), and the fluctuations in the histogram are proportional to $1/\sqrt{\ln f}$ Lee et al. (2006). The accuracy of $g(E)$ may be increased by using the $1/t$ algorithm Belardinelli and Pereyra (2007a); Belardinelli et al. (2008); Belardinelli and Pereyra (2007b). A mathematical generalization of the WL algorithm is the Stochastic Approximation Monte Carlo (SAMC) algorithm Liang et al. (2007); Liang (2009), which has also involved the concept used in the $1/t$ algorithm. Recently, Junghans et al. have demonstrated that WL, statistical temperature molecular dynamics, and metadynamics are equivalent under consistent initial conditions and update rules Junghans et al. (2014).

III.2 Wang-Landau with multi-dimensional random walk in phase space

The basic WL method for finding the joint DOS of a system, $g(E,V_{1})$, is to perform a two-dimensional random walk in the $(E,V_{1})$ space Wang and Landau (2001b); Kwak (2012); Silva et al. (2006); Tsai et al. (2007); Valentim et al. (2015). However, this approach is quite slow. To speed up the simulation, the system could be divided into multiple energy windows, with each window containing all the compatible $V_{1}$, using replica-exchange to ensure that all the microstates are accessible in each energy window Valentim et al. (2015). However, if the joint DOS contains one more variable, like the $g(E,M,M_{s})$ we want to obtain here, the simulation will again become slow. Here we performed several crude tests on the antiferromagnetic Ising model with different system sizes using this method.

We first adopt a strict ‘flatness’ criterion similar to the original WL papers Wang and Landau (2001a, b), which considers a histogram ‘flat’ if for every state, the deviation in the histogram $H(E,M,M_{s})$ is less than $20\%$ from the average histogram. Using parallel programming with the energy spectrum divided into 5 energy windows, with each window assigned 5 random walkers and 1 core, a $6\times 6$ Ising system takes 23 min to finish the simulation, while an $8\times 8$ Ising system takes 972 min ($\sim$ 16 hours). This approach obviously does not scale well with system size.

A relaxed ‘flatness’ criterion Brown (2014) considers a histogram ‘flat’ if the root-mean-square of the deviation from the average histogram is less then $20\%$, i.e.,

[TABLE]

where $N_{E,M,M_{s}}$ is the number of accessible $(E,M,M_{s})$ in the window. Using this criterion significantly improves the convergence times, but it does not improve the scaling behavior of this approach (see Table 1). A very rough estimate for the time it would take for an $32\times 32$ system to finish is $20000s\times 5^{10}\sim 6000$ years.

III.3 Two-stage method

Another method to obtain a joint DOS breaks the simulation into two stages of random walk Gervais et al. (2009); Li et al. (2013); Kalyan et al. (2016). In the first stage, a normal WL process is carried out and $g(E)$ is obtained. Then, a second stage of random walk is performed with an acceptance probability of $1/g(E)$, but only $H(E,M,M_{s})$ is updated, i.e. $g(E)$ remains unchanged. The joint DOS is then obtained from

[TABLE]

We partition the first-stage WL process into 5 energy windows with replica exchange. Each window has 5 walkers, and a separate core is assigned to each window. The strict ‘flatness’ criterion requiring that every state $(E,M,M_{s})$ does not deviate more than $20\%$ from the average histogram is adopted Wang and Landau (2001a, b). In the second stage of random walk, each core is assigned a random walker in a random initial state, which can walk through the whole phase space. We set the number of time steps spent on stage 2 to be 10 times of that spent on stage 1. This method is much faster than the previous method: an $8\times 8$ system can be finished in a few seconds, and even $12\times 12$ can be finished in 170 seconds if the walkers do not get ‘stuck’ Yin and Landau (2012); Koh et al. (2013, 2015).

However, this apparently reasonable approach does not yield the correct DOS for the two-dimensional Ising-ASFL model. The joint DOS $g(M,M_{s})$ should be independent of $H$ and $A$. Figures 1 (b) to (d) show that the results obtained by this method change significantly when $H$ and $A$ change, and all of them are quite different from the result obtained by exact combinatorial calculation in (a).

The reason for the differences can be explained as follows. Consider applying the two-stage method to get the DOS in terms of only two macroscopic variables, $(E,V_{1})$. In the second stage of the process, suppose a walker is trying to move from a phase point $(E,V_{a})$ to a point $(E,V_{b})$, and then to a point $(E,V_{c})$. As these three points have the same energy, the moves will have the same acceptance probabilities, $1/g(E)$. In this sense, these moves will be similar to an unbiased random walk confined to an energy $E$, and with the histogram corresponding to $V_{1}$ being recorded. Thus, if the change in $g(E,V_{1})$ along $V_{1}$ for a given value of $E$ is not too large, the method should give good results. However, if $g(E,V_{1})$ changes significantly when $V_{1}$ changes, the method may not give correct statistics. Unbiased sampling works for extremely small systems where the differences in $g(E)$ are small, but the WL method is required if the differences are significant.

The method can in principle be improved by doing separate WL simulation for each energy $E$ to find the statistics corresponding to $V_{1}$, if ergodicity is not broken. This observation leads to the basic principle of our macroscopically constrained WL method. To preserve ergodicity and obtain $g(E,M,M_{s})$ for our system, we can simplify the process to just perform simulations in the energy space for fixed $M$ and $M_{s}$. We discuss the method in detail in the following sections.

III.4 Other methods

A different WL method to obtain $g(E,M)$ was proposed by Zhou et al. Zhou et al. (2006). A kernel function is applied when one tries to update the histogram. The method appears to save time, but it is quite complicated and tuning of kernel functions seems to be required for different systems.

Very recently, Zablotskiy et al. Zablotskiy et al. (2016) used the stochastic approximation MC method to obtain the joint DOS, $g(V_{2},V_{1})$, of a polymer model, and then used it to deduce the $g(E)$ of the system. The way they obtain $g(V_{2},V_{1})$ is similar to the method in Sec. III.2 but including the $1/t$ algorithm, and only small $\ln g(V_{1},V_{2})$ were considered.

Two papers have recently been published that use methods similar to, but less general to the one presented here. Lourenço and Dickman Lourenço and Dickman (2016) obtained the two-dimensional joint DOS for the square-lattice Ising antiferromagnet, $g(E,M)$ at a single phase point, using the tomographic entropic sampling method Dickman and Cunha-Netto (2011); Belardinelli et al. (2014). The method employs many random walkers, each starting in a different energy state, and their results are combined to get $g(E,M)$. From this, they obtained the critical line and canonical averages of the order parameter, $M_{s}$. However, as pointed out in Dickman and Cunha-Netto (2011), the tomographic entropic sampling method cannot give correct results when the system size is large. In a study motivated by a network reliability problem, Ren, Eubank, and Nath presented a method to obtain the joint DOS for the square-lattice Ising ferromagnet, $g(E,M)$, using parallel WL simulations at fixed $M$ Ren et al. (2016). While conceptually similar to the method we present here, we note that they do not discuss generalizations to higher-dimensional order-parameter spaces. Moreover, we will here discuss several methods to simplify the computational process.

IV Macroscopically constrained Wang-Landau

IV.1 Basic idea

Suppose one wants to obtain the joint DOS for a system with $K$ macroscopic variables, $g(V_{K},...,V_{2},V_{1})$. Instead of letting the random walker travel in a $K-$dimensional phase space ($V_{K}$,…,$V_{2}$,$V_{1}$), which would require a very long time to obtain a ‘flat’ histogram, the simulation can be broken into many simulations performed in smaller phase spaces as follows. First, obtain $g(V_{1})$ through normal WL simulation or direct calculation. Then, break the large phase space into smaller phase spaces, each with a different fixed value of $V_{1}$. For each value of $V_{1}$, a separate WL simulation is performed to obtain the DOS with respect to only one macroscopic variable $V_{2}$, denoted as $g(V_{2}|V_{1})$. Next, each phase space is broken into smaller phase spaces, each with a different fixed value of $(V_{2},V_{1})$. Again, separate simulations for different fixed values of $(V_{2},V_{1})$ are performed to obtain DOS with respect to only one macroscopic variable $V_{3}$, denoted as $g(V_{3}|V_{2},V_{1})$. Iterating the process, the joint DOS with $K$ variables, $g(V_{K},...,V_{2},V_{1})$, can be obtained as

[TABLE]

In general, the macroscopic variables should be arranged such that the more fundamental building blocks of $g(V_{K},...,V_{1})$, like $g(V_{1})$, $g(V_{2},V_{1})$ and $g(V_{3},V_{2},V_{1})$, can be obtained in the most accurate manner. Therefore, if the joint DOS for two macroscopic variables can be obtained directly by exact calculation, they should be chosen as $V_{1}$ and $V_{2}$, so that the joint DOS $g(V_{K},...,V_{1})$ involves an exact factor, $g(V_{2},V_{1})$. However, when partitioning the simulations into different stages, one must be careful that in each stage, a simple method can be found to let the walker walk through the whole confined phase space, such that ergodicity is not broken.

Each time the walker performs a WL process with only one free macroscopic variable, the constrained DOS, e.g., $g(V_{4}|V_{3},V_{2},V_{1})$, can be partitioned into different windows of $V_{4}$, and then joined together smoothly through choosing the contact point with the most similar slopes with respect to $V_{4}$ as in REWL Vogel et al. (2013); Li et al. (2014).

Breaking down the single WL processes into many independent processes like this works fast and is more accurate compared to the methods discussed in Sec. III. Furthermore, the algorithm itself is very suitable for parallelization on many independent processors.

IV.2 $g(E,M,M_{s})$ for the Ising-ASFL model

To obtain the joint DOS, $g(E,M,M_{s})$, for the Ising-ASFL model, we choose $V_{1}=M_{s}$, $V_{2}=M$, and $V_{3}=E$. This is convenient because we can calculate $g(M,M_{s})$ exactly through direct combinatorial calculation as shown in Appendix A. Therefore, we can directly arrive at Eq. (5) and write

[TABLE]

Separate independent WL simulations will be performed for different fixed values of $(M,M_{s})$, each obtaining a DOS in terms of one macroscopic variable ($E$), $g(E|M,M_{s})$.

Any square lattice can be simply broken down into two sublattices, $A$ and $B$. Every alternate site belongs to the same sublattice. The magnetization ($M$) and the staggered magnetization ($M_{s}$) can be written in terms of the magnetizations of these two sublattice, $M_{A}$ and $M_{B}$, as

[TABLE]

Exchanging spins (Kawasaki dynamics) independently on each sublattice will preserve $M_{A}$ and $M_{B}$, and thus also preserves the values of $M$ and $M_{s}$. Moreover, it will allow the walker to walk through all the possible configurations and energies corresponding to each ($M,M_{s}$), and thus preserve ergodicity. This is the method used here to perform the random walk in microstates.

IV.3 Advantages of obtaining $g(E,M,M_{s})$

If we want to know the DOS $g(E)$ under different conditions, say for different external magnetic fields $H$, and different long-range interaction strengths $A$, using the simple WL method or importance sampling MC, we would have to perform separate runs every time these conditions are changed. However, if we can obtain $g(E,M,M_{s})$, we only have to do WL for a single set of $H$ and $A$. For simplicity we choose zero field and zero long-range interaction, i.e., $H=A=0$. The results for other parameter values can be obtained by simply shifting the result obtained for $H=A=0$. This happens because all the microstates are equally shifted in energy when a field-like model parameter changes, as the field is coupled to a global property, such as $M$, according to Eq. (1). Therefore, we can shift the DOS result from $H=A=0$ to the DOS for arbitrary $H$ and $A$ through the transformation,

[TABLE]

This shifting approach saves a very large amount of work.

IV.4 Simplification through symmetry considerations

Through the use of the shifting approach described in Sec. IV.3, we only have to consider $H=A=0$. This enables further simplification through symmetry considerations. Consider a spin configuration (microstate) that belongs to the macrostate $(M,M_{s})$ lying in region 0 of Fig. 2. For $H=A=0$, if we flip all the spins on sublattice $A$, $E$ and $M_{A}$ of the system will be reversed. From Eqs. (8) and (9), $M$ and $M_{s}$ of this new microstate are related to the original $M_{A}$ and $M_{B}$ through

[TABLE]

Therefore, we have

[TABLE]

This means that if we have obtained $g(E,M,M_{s})$ at one sampling point $(M,M_{s})$ in region 0, we can directly obtain $g(E,M,M_{s})$ at another point in region 1. There are seven similar symmetry considerations, which correspond to regions 1-7 in Fig. 2. It is important to note that $g(E,M,M_{s})$ must not be double-counted along the four symmetry axes in Fig. 2 when combining the results. These symmetry considerations reduce the computational work by roughly a factor of eight. In Appendix B, we show explicitly how to map from region 0 to the other seven regions.

IV.5 Simplification through uniform sampling

In region 0 of Fig. 2, after we have chosen a data point at ($M,M_{s}$)=($0,0$), the next possible pairs are ($0,2$) and ($2,0$), i.e. the smallest increment is $M_{\rm const}=2$. For a big system, the number of possible $(M,M_{s})$ pairs is very large. To get the DOS for every pair of $(M,M_{s})$ would require a huge amount of computational resources. Therefore, data points in the $(M,M_{s})$ space are chosen with a constant increment, $M_{\rm const}$. Here, choosing $M_{\rm const}=32$ gives good results for a system size of $L=32$. Proper values for $M_{\rm const}$ for different system sizes will be further discussed in Sec. VI.1.

IV.6 Simulations in practice

With the simplifications introduced in the previous sections, WL processes can be carried out for separate ($M,M_{s}$). Here, we just list a few points adopted in our simulations. The detailed implementation is described in Appendix C.

First, every ($M,M_{s}$) pair has a different accessible range of energies, and it is known that the extreme energy states have spins aligned close to either a strip or a droplet shape Leung and Zia (1990) (Fig. 3). We first estimate the extreme energies and then decide how many energy windows shall be used for each ($M,M_{s}$) pair.

Second, starting from the extreme energy states, artificial spin exchange processes are carried out in the initialization stage, so that we can find more accessible energy states at the beginning of the simulation.

Third, the root-mean-square ‘flatness’ criterion Brown (2014) is used, such that a histogram is regarded as ‘flat’ according to the root-mean-square deviation criterion,

[TABLE]

where $N_{E}$ is the number of accessible energy levels in that energy window, and the summation runs over these energy levels. This relaxation of the ‘flatness’ criterion can make simulations finish much earlier, as already demonstrated in Sec. III.2.

Fourth, the small statistical fluctuations in the DOS found at $H=A=0$ may be magnified near a critical point, causing difficulties in locating it accurately. Here we reduce the statistical fluctuations by obtaining 10 different $g(E,M,M_{s})$ through independent simulations, and taking the ensemble average.

IV.7 Simulation time

For $L=32$, which is the largest system size we have considered, we have kept around 125 energy levels in edge windows and around 200 energy levels in non-edge windows. Most non-edge windows can finish simulations in a few minutes, but the edge windows, which contain energy levels with low density of states, may take 20 min to more than one hour to finish. Therefore, most pairs of ($M,M_{s}$) considered can finish the simulation within a few minutes to a few hours. However, some edge windows may get ‘stuck’ at energy levels with low DOS Yin and Landau (2012); Koh et al. (2013, 2015) and do not converge after several days, especially when the sampling points include $M_{s}=0$. Therefore, we reject a simulation that does not finish in two days and re-start the run. Some sampling points may have to be rejected and re-started several times. With around two hundred computing cores, all the simulations for the $L=32$ system (including data for 10 different ensembles) could be finished within one week, with most of the time devoted to these ‘stuck’ sampling points. If one intends to take the ensemble average of 10 different $g(E)$ as we do here, one may submit more than 10 identical jobs for the sampling points with $M_{s}=0$ at the beginning, and reject all the runs for the remaining jobs after 10 of them have finished. This can make the simulation finish earlier. The ‘stuck’ problem is further discussed in Sec. VI.2.

V Application of $g(E,M,M_{s})$

V.1 Density of states for arbitrary $H$ and $A$

Figure 4(a) shows the joint DOS $g(E,M,M_{s})$ for $H=A=0$ obtained from our simulation for the Ising-ASFL model (Sec. II). By shifting the energy as stated in Sec. IV.3, we obtain $g(E,M,M_{s})$ for $(H,A)=(3,7)$ as shown in Fig. 4(b). Figure 5 shows that by summing over one component of $g(E,M,M_{s})$, we can obtain $g(E,M)$ and $g(E,M_{s})$, which give very smooth results. Indeed, $g(E,M,M_{s})$, $g(E,M)$, $g(E,M_{s})$ and $g(E)$ for arbitrary $H$ and $A$ can be obtained in this way.

V.2 Drawing phase diagrams

We introduce the normalized magnetization and staggered magnetization as $m=M/L^{2}$ and $m_{s}=M_{s}/L^{2}$ respectively, which are the order parameters of our system. After obtaining $g(E,M,M_{s})$ at $H=A=0$ through the simulations, we can use Eq. (10) to get $g(E,M,M_{s})$ at arbitrary points in the phase diagram and calculate different quantities as follows.

We can define the constrained partition function of any macrostate $(m,m_{s})$ as

[TABLE]

The overall partition function of the system is then

[TABLE]

The joint probability of finding the system in a macrostate $(m,m_{s})$ is

[TABLE]

where $\Delta m$, $\Delta m_{s}$ are the step sizes, both chosen to be the same value, $M_{\rm{const}}/L^{2}$. The free energy of macrostate $(m,m_{s})$ is

[TABLE]

As $(m_{A},m_{B})$ has a one-to-one relation with $(m,m_{s})$, we may express these quantities in terms of $(m_{A},m_{B})$, as well. The inset in Fig. 6 shows a free-energy contour diagram $F(m_{A},m_{B})$ close to the critical temperature for $H=A=0$.

We can sum over the contributions of the joint probability (Eq. 17) in one direction, obtaining the marginal probability density as

[TABLE]

With these distributions, we can calculate the expectation values of the order parameters and other quantities. In a complicated phase diagram that involves metastable phase regions, the stable phase will be the phase that has the larger total area in the marginal probability density, rather than the phase that shows the higher peak.

To locate and classify the critical points or lines between ordered and disordered phases in phase diagrams, the fourth-order Binder cumulant is often used Landau and Binder (2015); Binder (1981); Binder and Landau (1984); Challa et al. (1986); Nakada et al. (2011); Chan and Rikvold (2015),

[TABLE]

where $p$ is the order parameter of the system. As an illustration, we consider $A=0$ in the Ising-ASFL model, which is just a purely antiferromagnetic Ising model. It is commonly accepted that this critical line is in the Ising universality class, which (assuming isotropy and periodic boundary conditions) has a cumulant value near $0.61$ Kamieniarz and Blöte (1993). Therefore, using Eqs. (20) and (21) with $p=m_{s}$, we locate the critical line by finding the phase point with cumulant closest to $0.61$. The results are shown in Fig. 6, using $L=32$. The critical line obtained with our method is smooth, and the excellent agreement with the Wu and Wu analytic approximation Wu and Wu (1990) and the simulation results of Lourenço and Dickman Lourenço and Dickman (2016) indicate that our procedure of restarting ‘stuck’ simulation runs (see Secs. IV.7 and VI.2) does not lead to significant numerical inaccuracies.

Another common boundary line that separates different phases in the phase diagram is a first-order phase transition (coexistence) line. We locate it by looking at the order-parameter variance, which is proportional to the susceptibility times the temperature,

[TABLE]

This quantity has a local maximum value when evaluated at a point on the coexistence line, which serves as an accurate tool to determine this line. The line that becomes straight vertical for low $T$ in Fig. 7 is a coexistence line obtained by using Eqs. (19) and (22), with $p=m$ in the Ising-ASFL model with $A=1$.

The metastable phase regions in this phase diagram are bounded by the spinodal lines. We can express the free energy in terms of one order parameter as

[TABLE]

and the spinodal lines can be located by finding the points where the free energy changes from having two local minima to one local minimum. Figure 7 illustrates two spinodal lines, and Fig. 8(a) shows the shape of the free energy at constant $T$ when the system is at and near a spinodal line.

It may happen that a coexistence line or a critical line between two different phases lies in the metastable regions of other phases. In those cases, we have to remove contributions from these non-relevant phases. Figure 8(b) illustrates one example.

V.3 Probability densities and free energies at selected phase points

As we can shift $g(E,M,M_{s})$ to obtain different quantities for any point in the phase diagram, we can analyze the free energy and probability density at any selected phase point in detail.

Figure 9 shows what happens when moving along the temperature axis at $H=0$ and crossing the critical line in the antiferromagnetic Ising model, i.e., $A=0$. At low $T$, the joint probability density $P(m_{A},m_{B})$ has peaks only in the two antiferromagnetic (AFM) phases in the two opposite corners of the $(m_{A},m_{B})$ plane, $m_{s}=\pm 1$. Then, the two peaks connect weakly at the critical temperature. Above the critical temperature, the two peaks join, corresponding to a disordered phase. The free-energy contour diagram for the critical point at $H=0$ is shown as an inset in Fig. 6.

Similarly, Fig. 10 demonstrates what happens when $H$ is increased at constant, low $T$ to cross the critical line. The joint probability density $P(m_{A},m_{B})$ changes from two AFM phase peaks when the system is inside the critical line, to two peaks weakly connected at the critical line, and to a single ferromagnetic peak when the system is outside the critical line. The ferromagnetic peak moves toward the $(m,m_{s})=(1,0)$ corner and becomes more symmetric as $H$ further increases.

VI Generalizing the scheme to bigger systems

VI.1 Scaling effect in choosing $M_{\rm const}$

If the system size is large, in regions where phase transitions occur, the probability densities usually only show peaks in a small region of the order-parameter space. Therefore, the resolution, $M_{\rm const}$, must be small enough to observe these peaks. From finite-size scaling theory, plotting $L^{-\beta/\nu}P_{L}(m_{s})$ vs $L^{\beta/\nu}m_{s}$ for different system sizes $L$ will give curves that coincide (Fig. 11). That means that if we we can get a perfect result with a small system size $L_{s}$, when we consider a bigger system $L_{b}$, the increment has to be chosen such that $M_{{\rm const},b}\leq M_{{\rm const},s}\times(L_{b}/L_{s})^{2-\beta/\nu}$, where $\beta$ and $\nu$ are critical exponents. For critical points that belong to the Ising universality class, $\beta/\nu=1/8$, whereas for the mean-field class, $\beta/\nu=1/2$.

As the number of $(M,M_{s})$ pairs that must be sampled is $L^{4}/M_{\rm const}^{2}$, the minimum number of $(M,M_{s})$ pairs required for a big system is proportional to $(L_{b}/L_{s})^{2\beta/\nu}$. In this sense, as the Ising-ASFL model contains mean-field critical points for large values of $A$ Rikvold et al. (2016), when the system size is doubled, the number of $(M,M_{s})$ pairs must also be doubled. If one is working with a purely short-range ferromagnetic or antiferromagnetic Ising model, all critical points will be in the Ising class, so that the required number of $(M,M_{s})$ pairs one has to consider will be nearly unchanged ($\times 1.189$) if $L$ is doubled.

VI.2 Other problems for large systems, ‘flatness’ criteria, and possible solutions

The occurrence of ‘stuck’ simulation runs has also been observed in other WL based algorithms Yin and Landau (2012); Koh et al. (2013, 2015), particularly in phase regions of extreme energy and/or very low DOS. It is thus not a specific consequence of the macroscopic constraints in the present method.

In the work presented here, simulation runs for larger systems sometimes get ‘stuck,’ particularly at phase points with $M_{s}\approx 0$. These points indeed correspond to extreme energy configurations arranged as perfect strips or droplets (Fig. 3). These states, or states that have similar energies, are also very rare and thus have a very low $g(E|M,M_{s})$. Many of them are found in the initialization process, but they are very hard to reach during the subsequent simulation run. For $L=32$, this problem can still be solved by rejecting and restarting the run without causing significant sampling error. (See the numerical results in Secs. V and VI.1.) However, further increase in system size to $L=64$ causes more simulations to get ‘stuck.’

Using a ‘strict’ flatness criterion as in Wang and Landau (2001a, b), i.e., requiring the histogram of every energy level to not deviate too much from the average histogram, simulations may have great difficulty finishing. Although the presence of these extreme states is known through the initialization process, they are very hard to reach in the random walk process, so that the ‘flat’ histogram may not be attainable. Significant improvement is obtained by using the relaxed flatness criterion Brown (2014) as in Eq. (14), by which the histogram is accepted as ‘flat,’ even though it deviates significantly from the average for these few energy levels. Even without partitioning the energy spectrum into windows, runs that cannot finish in several weeks with the strict flatness criterion can finish in a few days with the relaxed criterion.

However, this method also brings another problem. If a few energy levels are very difficult to reach and are not reached within the first few times a ‘flat’ histogram is accepted, then, if one of these states is reached later, it will have $\ln g(E|M,M_{s})=0$ while other states that previously were visited can have $\ln g(E|M,M_{s})$ of the order of $10^{7}$. The large difference between the two neighboring states causes the walker to stay in this state to increase its $\ln g(E)$ until it reaches the order of $10^{7}$. However, every time a ‘flat’ histogram is obtained, the increment, $\ln f$ is reduced by a factor of two, so that in this run, the increment step size may have already dropped to a very small value, like $10^{-6}$. Then it will take a very long time to raise $\ln g(E|M,M_{s})$ to the order of $10^{7}$, and the program will get stuck. As the difference between two neighboring states (or slope of $g(E|M,M_{s})$ vs $E$) increases with the system size, a simulation under these conditions may easily get ‘stuck’ for large $L$.

To generalize the method to larger systems, one promising solution is dividing the edge windows into smaller windows and applying REWL Vogel et al. (2013); Li et al. (2014); Vogel et al. (2014) in these windows, so that each walker is confined to a smaller energy range. This forces it to sample all these energy levels with low density of states, and at the same time the replica can sample through all possible states in the edge windows. Yet another way that may also alleviate the problem is to let the walker occasionally jump to a previous state stored in a configuration database Vogel and Perez (2015).

VII Conclusion

A macroscopically constrained WL method is proposed, which may be useful in finding DOS with more than one variable, and in obtaining complex phase diagrams. The method converts a multi-dimensional random-walk process into many one-dimensional random walks, with each walker constrained to fixed values of certain macroscopic order parameters. The method is demonstrated and validated on a two-dimensional antiferromagnetic Ising model with ferromagnetic long-range interaction. We obtained the joint DOS $g(E,M,M_{s})$, through simulations at $H=A=0$. The DOS for arbitrary values of $(H,A)$ then follow by a simple transformation of the total system energy, and all the thermodynamic quantities for any point in the phase diagram can then be found. We demonstrate how to use the DOS obtained to efficiently draw phase diagrams and free-energy landscapes for the Ising-ASFL model.

The detailed physics of the Ising-ASFL (spin-crossover material) model, including complex phase diagrams for several values of $A$, will be described in forthcoming papers Chan et al. (a, a).

ACKNOWLEDGMENTS

Chor-Hoi Chan thanks Alexandra Valentim and Ying-Wai Li for helpful discussions of the Wang-Landau and replica-exchange Wang-Landau methods. The simulations were performed at the Florida State University High Performance Computing Center This work was supported in part by NSF Grant No. DMR-1104829.

Appendix A

The exact combinatorial calculation of $g(M,M_{s})$ for the Ising-AFSL model is described below.

Let $N_{A}$ be the number of sites on sublattice $A$ that have spin up, and $N_{B}$ be the number of sites on sublattice $B$ that have spin up. Then, the magnetization ($M$) and staggered magnetization ($M_{s}$) in Eqs. (8) and (9) can be rewritten as

[TABLE]

As the joint DOS $g(M,M_{s})$ is defined as the total number of spin configurations (microstates) that have certain $(M,M_{s})$ values, $g(M,M_{s})$ can be visualized as the total number of ways to allocate $N_{A}$ upspins on sublattice $A$ and $N_{B}$ upspins on sublattice $B$. This can be expressed as the product of two binomial factors,

[TABLE]

where $N/2=L\times L/2$ is the total number of sites on each sublattice. The general binomial recursive formula,

[TABLE]

with boundary values

[TABLE]

is used to speed up the calculation. Overflow problems will be present if the value of $n!$ is too large, so Stirling’s approximation is employed when $\ln g(M,M_{s})$ is greater than $700$:

[TABLE]

Appendix B

The symmetry considerations used to map the DOS from region 0 to regions 1-7 in Fig. 2 are described below.

Region 1: $M_{A}\rightarrow-M_{A}$ and $E\rightarrow-E$.

This is reflection about the $M_{B}$ axis, which means that, if we change $M_{A}$ to $-M_{A}$ by flipping all the spins on sublattice A, the energy changes from $E$ to $-E$, and the new $M$ and $M_{s}$ are related to the original $M_{A}$ and $M_{B}$ through

[TABLE]

Thus,

[TABLE]

Region 2: $M_{B}\rightarrow-M_{B}$ and $E\rightarrow-E$.

This is reflection about the $M_{A}$ axis, which means flipping all the spins on sublattice B. Thus,

[TABLE]

Region 3: $M_{A}\rightarrow-M_{A}$ and $M_{B}\rightarrow-M_{B}$ and $E\rightarrow E$.

This is a combination of reflection about the $M_{A}$ axis and $M_{B}$ axis, which means flipping all the spins on both sublattices. This preserves the energy, and we have

[TABLE]

Region 4: $M_{A}\leftrightarrow M_{B}$ and $E\rightarrow E$.

This is reflection about the $M$ axis, which means exchanging all the spins between the two sublattices. This preserves the energy, and we have

[TABLE]

Region 5: $M_{A}\rightarrow M_{B}$ and $M_{B}\rightarrow-M_{A}$ and $E\rightarrow-E$.

This is a combination of reflection about both the $M_{s}$ and $M_{A}$ axes, which means replacing spins on sublattice A by spins on sublattice B, and spins on sublattice B by the flipped spins on sublattice A. This reverses the energy, and we have

[TABLE]

Region 6: $M_{A}\rightarrow-M_{B}$ and $M_{B}\rightarrow M_{A}$ and $E\rightarrow-E$.

This is a combination of reflection about the $M_{s}$ and $M_{B}$ axes, which means replacing spins on sublattice A by flipped spins on sublattice B, and spins on sublattice B by the spins on sublattice A. This reverses the energy, and we have

[TABLE]

Region 7: $M_{A}\rightarrow-M_{B}$ and $M_{B}\rightarrow-M_{A}$ and $E\rightarrow E$.

This is reflection about the $M_{s}$ axis, which means replacing spins on sublattice A by flipped spins on sublattice B, and spins on sublattice B by the flipped spins on sublattice A. This preserves the energy, and we have

[TABLE]

Appendix C

Details of the implementation scheme for not too large systems are given below.

Step 1: Determine the desired combination of $M$ and $M_{s}$. While $L=6$ and $12$ were simulated with all possible $(M,M_{s})$ pairs (i.e., step size $M_{\rm const}=2$ ), $L=32$ was simulated with $M_{\rm const}=32$, which gives good results. The symmetries described in Sec. IV.4 and Appendix B ensure that we only have to choose data points within one octant of the $(M,M_{s})$ space.

Step 2: Each chosen pair of $(M,M_{s})$ is submitted as a parameter to identical, independent WL programs running on separate processing cores (or sequentially on one core). For each WL process, the values of $M$ and $M_{s}$ are conserved in every time step. Note that fixing an ($M$,$M_{s}$) pair is the same as fixing an ($M_{A}$,$M_{B}$) pair (refer to Eqs. (8) and (9) and Fig. 2).

Step 3: Find the maximum and minimum energies of the system for each chosen $(M,M_{s})$ pair. For each $(M,M_{s})$ pair, divide the energy spectrum into windows if the number of energy levels is large. We keep around $200$ energy levels in each window (around $125$ energy levels in the edge windows), using $50\%$ overlap between neighboring windows.

For a fixed ($M$, $M_{s}$) pair, macrostates that have energies close to the extreme energies are in general either arranged in a configuration close to a strip or a nearly square droplet Leung and Zia (1990) (Fig. 3). We first prepare all the spins pointing up (purely ferromagnetic state), and then flip the spins separately on the two sublattices sequentially one by one until we get a configuration (microstate) that satisfies the chosen ($M_{A}$,$M_{B}$) pair. This locates the highest energy microstate that is arranged close to a strip shape. To locate the highest energy states that are close to a droplet form, we again prepare all the spins pointing up (purely ferromagnetic state), and calculate the approximate size of the droplet if a group of down-spins are arranged in a nearly droplet shape. Then we again flip the spins on each sublattice separately until $M_{A}$ and $M_{B}$ are satisfied. To locate the lowest energy states, we first prepare all spins arranged alternately up and down (purely antiferromagnetic), and then we perform similar flips to get states close to antiferromagnetic strip and droplet shapes.

Step 4: List as many energy levels as possible that satisfy this ($M$,$M_{s}$) pair and initialize $\ln g(E|M,M_{s})=0$ for each energy level found. Using the microstates found in Step 3, which have near-extreme energies and already satisfy the ($M$,$M_{s}$) constraint, we randomly perform spin exchange separately on the two sublattices many times to arrive at most of the energy levels that satisfy the same $M_{A}$ and $M_{B}$. For each energy window, we store one spin configuration as the initial configuration of that energy window. If some energy levels are still missing in this step, they will be visited later during the random walk process. Each newly found state is immediately initialized to $\ln g(E|M,M_{s})=0$.

Step 5: Randomly choose the energy window to start with, and load the spin configuration stored during the initialization stage to use as the starting configuration. Set the histogram $H(E|M,M_{s})=0$ for all the energies $E$ in that energy window.

Step 6: Propose a move in microstates by doing spin-exchange (Kawasaki dynamics) between two sites with different spins on the same sublattice, and decide whether the move is accepted according to

[TABLE]

Step 7: Update the DOS $g(E|M,M_{s})$ for each ($M$,$M_{s}$) pair as $\ln g(E|M,M_{s})\rightarrow\ln g(E|M,M_{s})+\ln f$ whenever the energy level $E$ is visited. At the same time, its histogram is updated as $H(E|M,M_{s})\rightarrow H(E|M,M_{s})+1$

Step 8: Using the root-mean-square ‘flatness’ criterion in Eq. (14), check whether a ‘flat’ histogram $H(E|M,M_{s})$ is obtained for each ($M$,$M_{s}$) pair after every $100\times L^{2}$ time steps, and reduce the modification factor $f$ to its square root value whenever a ‘flat’ histogram is reached.

When the code is found to be stuck and cannot finish after a sufficiently long time for a particular $(M,M_{s})$ pair, we reject that run and re-start the simulation for that $(M,M_{s})$ pair. This reject and rerun process may have to be performed a few times on one $(M,M_{s})$ pair in order to obtain a final converged result. Further discussion of the convergence is given in Secs. IV.7 and VI.2.

$\ln g(E|M,M_{s})$ is shifted to a small value every time we get a ‘flat’ histogram, so that we can avoid $\ln g(E|M,M_{s})$ accumulating to very large values that cause overflow problems. Restart the run at one of the configurations that has nearly minimum histogram in the last run. Repeat the random walk process until the modification factor is less than $10^{-8}$.

Step 9: Join the $\ln g(E|M,M_{s})$ curves obtained in different windows for each ($M$,$M_{s}$) pair at the energy where $\ln g(E|M,M_{s})$ in the two adjacent windows have the closest slopes with respect to $E$. Use formula (27) to obtain the exact number of microstates for each pair of ($M$,$M_{s}$). The overall DOS, $g(E,M,M_{s})$ or $\ln g(E,M,M_{s})$, can then be obtained for the whole system through Eq. (7).

Step 10: The $g(E,M,M_{s})$ obtained still has numerical error. To average out the error, 10 different $g(E,M,M_{s})$ were obtained, and their ensemble average was used as $g(E,M,M_{s})$.