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

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 multivariable density of states is obtained for one set of values of fieldlike 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 multivariable 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 fieldlike parameters in this model are an effective magnetic field and the strength of the long-range interaction.

Figure 1

The two-dimensional joint DOS, $g(M,M_s)$, for $L=12$, obtained by (a) exact combinatorial calculation as described in Appendix pp1, and by the two-stage method, using different values of $H$ and $A$, (b) $(H,A)=(0,0)$, (c) $(H,A)=(3,0)$, and (d) $(H,A)=(0,7)$. The $g(M,M_s)$ should be independent of $H$ and $A$. The deviations from (a) of the results shown in (b)–(d) imply that the two-stage method does not give correct results for the Ising-ASFL model, even for this small system size.

Figure 2

Uniform sampling of ($M,M_s$) pairs is performed in region 0. Some data points are illustrated as black crosses in the figure. After the sampling in region 0 is finished, symmetries are used to obtain data for all the other seven regions.

Figure 3

Energy states that are close to the extreme energy states. The spins are aligned close to either a strip or a droplet shape [63].

Figure 4

$lng(E,M,M_s)$ vs $E/4, M/2$, and $M_s/2$ for (a) $(H,A)=(0,0)$, and (b) $(H,A)=(3,7)$, both using $L=6$ for improved visibility. The results for larger systems are similar. $g(E,M,M_s)$ at $(H,A)=(3,7)$ is obtained by shifting the $g(E,M,M_s)$ at $H=A=0$ through Eq. (10). The color of the data points shows the relative magnitude of the natural logarithm of the DOS, ranging from red (smallest) to magenta (largest). Only results for $M_s\le 0$ are shown as there is reflection symmetry about the $M_s=0$ plane.

Figure 5

$lng(E,M)$ vs $E/4$ and $M/2$ for (a) $(H,A)=(0,0)$, and (c) $(H,A)=(3,7)$. $lng(E,M_s)$ vs $E/4$ and $M_s/2$ for (b) $(H,A)=(0,0)$, and (d) $(H,A)=(3,7)$. All using $L=12$. All the data are obtained by summing over the contribution of different directions in $g(E,M,M_s)$. $g(E,M,M_s)$ at $(H,A)=(3,7)$ is obtained by shifting the $g(E,M,M_s)$ at $H=A=0$ through Eq. (10), i.e., all data are obtained from $g(E,M,M_s)$ at $H=A=0$.

Figure 6

Critical line for an antiferromagnetic Ising system. The blue critical line is obtained by increasing $H$ from 0 in steps of 0.01 or 0.02, then performing a temperature ($T$) scan, choosing $\mathrm{\Delta }T$ to be 0.001 to 0.005, and locating the critical line by choosing the point that gives the cumulant value [Eq. (21) with $p=m_s]$ closest to 0.61 [68]. Data points on the negative $H$ side are obtained by reflection. The analytically approximated critical line (red), obtained by the method of Ref. [69], is also plotted. The two results coincide at this resolution. The inset shows the free-energy contour diagram at the critical temperature for $H=0$ [refer to Eq. (18)].

Figure 7

The low-temperature portion of the phase diagram for $A=1$, near where the critical line ends at a tricritical point, using $L=32$. Metastable phases appear, leading to the appearance of a coexistence line flanked by spinodal lines. The critical line is located by the cumulant method [Eq. (21)] in the same way as in Fig. 6 with $p=m_s$. The coexistence line is located by finding the maximum susceptibility [Eq. (22)] with $p=m$ when $H$ is changed at constant $T$. Spinodal lines are located by finding the value of $H$ where the free energy $F\left(m\right)$ [Eq. (23)] changes from having two local minima to only one local minimum. The positions of the lines at $T=0$ agree with exact ground-state calculations [24].

Figure 8

(a) Free energy $F\left(m\right)$ [Eq. (23)] for $A=1$ and $T=0.11$, when the system is in the metastable region with two local minima at $H=3.689$ (red dashed), at the spinodal point at $H=3.789$ (black solid), and outside the metastable region at $H=3.889$ (blue dotted). (b) Three different phases are shown in this free-energy diagram at $(H,T,A)=(2.267,0.033,8)$. The system is lying on the coexistence line between the phases in the middle and on the right-hand side. When locating the coexistence line through the susceptibility [Eq. (22)], one has to remove all the contributions from the metastable phase on the left-hand side of the black point.

Figure 9

(a)–(d) Joint probability density $P(m_A,m_B)$ [Eq. (17) with a change of variables]; (e)–(h) corresponding free energy $F\left(m_s\right)$ [Eq. (24)] when the system ($L=32$) is moved from a low temperature to a high temperature, crossing the critical line at $H=0$. (a),(e) $T=2<T_c$, where system is in one of the AFM phases with equal probability. (b),(f) $T=2.275=T_c$ for this system size, where the two peaks are connected by a “bridge” and the cumulant is approximately equal to 0.61. (c),(g) $T=2.4>T_c$, where system is in the disordered phase with large AFM fluctuations. (d),(h) $T=4\gg T_c$, where the AFM fluctuations are much less pronounced. Note the different scales in the free-energy plots (e), (f), (g), (h).

Figure 10

Joint probability density $P(m_A,m_B)$ [Eq. (17) with a change of variables] when the system ($L=32$) is moved parallel to the $H$ axis at a low temperature, $T=0.06$, and crosses the critical line. (a) $H=3.95<H_c$, where system is close to the critical line but still in the AFM phases. (b) $H=3.96=H_c$ for this system size and temperature, where the two peaks are connected by a bridge and the cumulant is approximately equal to 0.61. (c) $H=3.97>H_c$, where system is in the $\mathrm{FM}+$ phase. (d) $H=4>H_c,$ where the $\mathrm{FM}+$ phase peak is more symmetric and closer to the corner at $m=1$.

Figure 11

(a) Marginal probability density [Eq. (20)] $P_L\left(m_s\right)$ vs $m_s$ for $L=6$ (green), 12 (red), and 32 (blue) at $H=0$, and at their corresponding critical temperatures as determined from the cumulant value. (b) Shows the scaled plot, i.e., $L^{-\beta /\nu }{P}_L\left(m_s\right)$ vs $L^{\beta /\nu }{m}_s$. The excellent data collapse indicates that even these small systems are fully in the asymptotic scaling regime, at least in this region of the phase diagram. It also indicates that any errors that might be caused by our procedure of restarting stuck simulation runs (see Secs. 4g and 6b) are insignificant.