Density of states for systems with multiple order parameters: a constrained Wang-Landau method
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] D ec Density of states for systems with multiple order parameters: aconstrained Wang-Landau method
C. H. Chan, ∗ G. Brown, † and P. A. Rikvold ‡ Department of Physics, Florida State University,Tallahassee, Florida 32306-4350, USA
Abstract
A macroscopically constrained Wang-Landau Monte Carlo method was recently proposed tocalculate the joint density of states (DOS) for systems with multiple order parameters. Here wedemonstrate results for a nearest-neighbor Ising antiferromagnet with ferromagnetic long-rangeinteractions (a model spin-crossover material). Its two relevant order parameters are the magneti-zation M and the staggered magnetization M s . The joint DOS, g ( E, M, M s ) where E is the totalsystem energy, is calculated for zero external field and long-range interaction strength, and thenobtained for arbitrary values of these two field-like model parameters by a simple transformationof E . Illustrations are shown for several parameter sets. ∗ [email protected] † [email protected] ‡ [email protected] . MOTIVATION, METHOD AND MODEL To obtain phase diagrams for many-particle systems with multiple order parameters andfield-like system parameters by standard importance sampling Monte Carlo (MC) simula-tions, one typically needs to perform simulations for a range of temperatures and field-likeparameters. This is a quite tedious and computationally intensive task. Much effort canoften be saved by calculating the density of states (DOS) in energy space, g ( E ), using theWang Landau (WL) method [1, 2] or one of its many later refinements. However, if a joint DOS involving the energy and several order parameters, g ( E, O , O , ... ) is needed,the random walks that lie at the heart of the WL method become multi-dimensional andconsequently quite slow [3].Motivated by a computational problem of this nature (see specifics below), we recentlydeveloped an algorithm based on the WL method, in which the multidimensional randomwalk in energy and order parameters is replaced by many independent, one-dimensionalwalks in energy space, each constrained to fixed values of the order parameters [3]. Forthis reason, we call it the macroscopically constrained Wang-Landau (MCWL) method. (Asimilar, but less general, method was recently proposed by Ren et al. [4].) Our implemen-tation of the algorithm is designed to execute in trivially parallel fashion on a large set ofindependent processors. Through further, symmetry based simplifications [3], the methodprovides an accurate estimate of a joint DOS with three variables in a relatively short time.All the simulations necessary to obtain the main numerical results presented here and inRefs. [3] and [5] were finished within one week, running independently on about two hundredseparate computing cores. Here, we use the method to obtain the densities of states of amodel spin-crossover material with two macroscopic order parameters, needed to producephase diagrams under arbitrary external conditions. However, it has the potential to beapplied to other complicated systems and phase spaces of higher dimension. Details of thealgorithm and its implementation are given in Ref. [3].The problem that inspired our development of the MCWL method is a model spin-crossover material. These materials are molecular crystals composed of magnetic ions sur-rounded by organic ligands. The molecules can be in two different magnetic states, a highlydegenerate, high-spin (HS) and a less degenerate, low-spin (LS) state. The HS state hashigher energy and larger molecular volume than the LS state, and the volume difference leads2o elastic distortions of crystals with a mixture of HS and LS molecules. A simple model ofphase transitions in some spin-crossover materials can be constructed as an antiferromag-netic Ising-like model with elastic distortions [6]. The elastic interactions can be furtherapproximated as a ferromagnetic mean-field (equivalent neighbor or Husimi-Temperley) in-teraction, leading to the following Hamiltonian on an L × L square lattice with periodicboundary conditions [3, 7, 8], H = J X h i,j i s i s j − HM − A L M , (1)with J >
0. Each spin, s i is either up (+1, HS) or down ( −
1, LS). J P h i,j i s i s j representsthe short-range interaction between neighboring spins, which favors adjacent spins aligningalternately up and down in an antiferromagnetic (AFM) configuration. M = P i s i is themagnetization and H is the strength of the applied field (which is actually an effective fieldin the spin-crossover material [8, 9]). The term − HM favors all spins aligned with H inthe same direction. L is the total number of sites, and A ≥ M , suchthat − A L M has its most negative value if all the spins are aligned in the same direction.Therefore, this amounts to a ferromagnetic (FM) long-range interaction. We will refer tothis model as the two-dimensional Ising Antiferromagnetic Short-range and FerromagneticLong-range (2D Ising-ASFL) model. Its two field-like system parameters are H and A .Throughout the paper, the total energy ( E ), magnetic field ( H ), and long-range inter-action ( A ), will be expressed in units of J . The lattice is divided into two interpenetratingsublattices, A and B in a checkerboard fashion, with sublattice magnetizations M A and M B ,respectively. The two order parameters of the model are the magnetization ( M = M A + M B )and the staggered magnetization ( M s = M A − M B ).The conditional DOS, g ( E | M, M s ), is determined through independent, one-dimensionalWL simulations for fixed M and M s . The order parameters are kept constant by performingmicroscopic spin-exchange (Kawaski) updates separately on each sublattice. The joint DOSis then obtained by an analogue of Bayes’ theorem as g ( E, M, M s ) = g ( E | M, M s ) P E g ( E | M, M s ) g ( M, M s ) , (2)where g ( M, M s ) is determined by exact combinatorial calculation. It is only necessary tocalculate g ( E | M, M s ) for one single set of H and A . Here we use H = A = 0, which3orresponds to a simple square-lattice Ising antiferromagnet [10]. The joint DOS for anyarbitrary value of ( H, A ) can be obtained by the simple transformation of E , g ( E, M, M s ) → g ( E − HM − AM L , M, M s ) . (3)This result is based on the fact that all microstates are equally shifted in energy when afield-like parameter, such as H or A , couples to a function of a global property, such as M (or M s ), as shown in Eq. (1). II. RESULTS
Figure 1 shows the joint DOS, g ( E, M, M s ), for different magnetic fields H and long-rangeinteraction strengths A . Simulations were only performed for the case H = 0 and A = 0,shown in (a). Only 1 / M, M s )plane as described in [3]. Data points in (b) to (f) are all obtained through transformation ofthe energy axis in (a), using Eq. (3). When the field H is applied to the system as in (d) and(e), the term − HM in the Hamiltonian pulls down the states with positive magnetization( M >
0) linearly, while it pulls up the states with negative magnetization (
M < A ) present, due to the term − AM / L , stateswith non-zero magnetization ( M ) are pulled down in a quadratic manner. If both H and A are applied, some states are shifted up, and some are shifted down as shown in (f). Figure1 shows results for the largest system size we studied, L = 32. While Fig. 1 shows the DOSversus three variables ( E, M, M s ), Figs. 2 and 3 show the DOS after summing the resultsalong one axis, i.e., when the DOS is plotted against two variables only, E and M , and E and M s , respectively.Figure 2 shows g ( E, M ). A similar effect as discussed for Fig. 1 is seen. I.e., when H isapplied, the energies for states with positive M are pulled down, and the energies for stateswith negative M are pulled up. When A is applied, the energies for all states with non-zero M are pulled down quadratically.Figure 3 shows g ( E, M s ). Parts (a), (e), and (i) show the effect of increasing H with A = 0. The energy range of states with M s near zero increases as states with negative M are pulled up in energy, and states with positive M are pulled down. The latter effect4 IG. 1. Natural logarithm of the joint DOS, ln g ( E, M, M s ), of the 2D Ising-ASFL model forsystem size L = 32 for different field strengths ( H ) and long-range interaction strengths ( A ). Thethree axes are M/ M s /
64, and E/
4. The colors of the data points show the relative magnitudeof the natural logarithm of the DOS, ranging from red (lowest) to magenta (highest). Only resultsfor M s ≤ M s = 0 plane. M and M s are bothsampled with an increment of 64. All possible energy levels are included, giving an increment of 4for E at fixed M and M s . also causes the energy corresponding to the highest DOS to shift to a lower value. Theother graphs in Fig. 3 show that increasing the long-range interaction strength A lowers theminimum energy of states with M s near zero, as the energies of the corresponding stateswith large positive and negative M are pulled down quadratically.5 IG. 2. Plots of the logarithmic DOS vs energy ( E ) and magnetization ( M ), ln g ( E, M ), fordifferent fields H and long-range interactions A . The smaller system size, L = 12, is used forbetter graphical resolution. All possible values of E and M for this system size are considered. III. SUMMARY
As an example of the application of the recently proposed macroscopically constrainedWang-Landau Monte Carlo method for systems with multiple order parameters, we havepresented the joint DOS of a square-lattice Ising model with antiferromagnetic short-rangeinteractions and ferromagnetic long-range interactions (a model spin-crossover material [8]).Based on data from simulation for one single parameter set, results are shown for various val-ues of the model’s field-like parameters, applied field H and long-range interaction strength, A . Further results on the phase diagrams of this model spin-crossover material are given inRef. [5]. 6 IG. 3. Plots of the logarithmic DOS vs energy ( E ) and staggered magnetization ( M s ), ln g ( E, M s ),for different fields H and long-range interactions A , again using L = 12. All possible values of E and M s for this system size are considered. ACKNOWLEDGMENTS
Chan thanks Alexandra Valentim and Ying Wai Li for helpful discussions of the WLmethod. The Ising-ASFL model was first proposed by Seiji Miyashita, and we thank himfor useful discussions. The simulations were performed at the Florida State University HighPerformance Computing Center. This work was supported in part by NSF grant No. DMR-1104829. 7
EFERENCES [1] Wang F and Landau DP 2001
Phys. Rev. Lett . Phys. Rev. B Phys. Rev. E Phys. Rev. E Phys. Rev. B Phys. Rev. B Phys. Proc. Phys. Rev. B J. Phys. (Paris), Colloq. C1-91[10] Louren¸co BJ and Dickman R 2016
J. Stat. Mech.033107