Uncovering the secrets of unusual phase diagrams: applications of two-dimensional Wang-Landau sampling

Uncovering the secrets of unusual phase diagrams: applications of two-dimensional Wang-Landau sampling

Shan-Ho Tsai^{I, II}; Fugao Wang^I,*; D. P. Landau^I

^ICenter for Simulational Physics, University of Georgia, Athens, GA 30602, USA

^IIEnterprise Information Technology Services, University of Georgia, Athens, GA 30602, USA

ABSTRACT

We use a two-dimensional Wang-Landau sampling algorithm to calculate the density of states for two discrete spin models and then extract their phase diagrams. The first system is an asymmetric Ising model on a triangular lattice with two- and three-body interactions in an external field. An accurate density of states allows us to locate the critical endpoint accurately in a two-dimensional parameter space. We observe a divergence of the spectator phase boundary and of the temperature derivative of the magnetization coexistence diameter at the critical endpoint in quantitative agreement with theoretical predictions. The second model is a Q-state Potts model in an external field H. We map the phase diagram of this model for Q > 8 and observe a first-order phase transition line that starts at the H = 0 phase transition point and ends at a critical point (T_c,H_c), which must be located in a two-dimensional parameter space. The critical field H_c(Q) is positive and increases with Q, in qualitative agreement with previous theoretical predictions.

Keywords: Ising and Potts models; Critical endpoints; Two-dimensional Wang-Landau

I. INTRODUCTION

Conventional Monte Carlo methods[1] generate a canonical distribution at a fixed temperature T and with other relevant parameters, such as the external field H, fixed as well. Through histogram reweighting[2] this canonical distribution can then be extended to a region of T and H near the respective values used in the initial simulation. For large system sizes, however, the region in T and H for which the results can be reweighted is very narrow. Therefore, numerous simulations have to be carried out for different values of T (and H) in order to locate phase transitions in an arbitrary case. In contrast, the Wang-Landau sampling method estimates the density of states directly and accurately[3–5] and allows determination of most thermodynamical quantities for all values of T and H with a single simulation. Wang-Landau sampling has proven very useful and efficient in many different applications, including studies of complex systems with rough energy landscapes.

In this paper we show how the Wang-Landau sampling method can be used to study discrete spin models in two-dimensional parameter spaces. The first model we consider is an asymmetric Ising model on a triangular lattice with two- and three-body interactions in a uniform external field. For a range of coupling parameters[6] this model has a critical endpoint (CE), i.e. a point in the phase diagram where a critical line meets and is truncated by a first-order transition line. The bulk critical exponents are believed[7] to be unchanged at the CE, but this has not been checked beyond phenomenological theory and renormalization calculations[8]. In 1990 Fisher and Upton argued that a new universal singularity in the first-order phase transition line should arise at a CE[7, 9]. This prediction was confirmed by Fisher and Barbosa's phenomenological studies for an exactly solvable spherical model[9]. Later, Wilding[10] studied critical endpoint behavior in a symmetrical binary fluid mixture using a large scale Monte Carlo simulation with multicanonical[11] and histogram reweighting[2] methods. Wilding predicted and observed a singularity in the derivative of the diameter of the liquid-gas coexistence curve at the CE. He also provided the first numerical evidence of the divergence in the curvature of the spectator phase boundary[10], in accordance with previous theoretical prediction [7, 9]; however, a quantitative analysis of singularities at the CE was not carried out because much larger system sizes would have been required. Our current work provides a quantitative test of the theoretical predictions[7, 9, 10] for these two divergences at the CE in an asymmetric model.

The second model that we study with the Wang-Landau sampling method is a Q-state Potts model in the presence of an external field H. An accurate estimate of the density of states allows us to map out the phase diagram in the (T,H) plane and test the theoretical prediction by Goldschmidt[12] that this model exhibits a first-order phase transition line starting at the H = 0 phase transition point and ending at a critical point (T_c,H_c) with H_c(Q) > 0. The value of H_c(Q) is predicted to increase with Q.

This paper is organized as follows: Sec.II describes the Wang-Landau sampling method, Sec.III defines the asymmetric Ising model and presents its results, and Sec.IV is dedicated to the Potts model in an external field. Some concluding remarks are given in Sec.V.

II. WANG-LANDAU SAMPLING

The Hamiltonian of the two models studied in this paper has the form

where E is an interaction energy term and M is the magnetization of the system. The density of states g(E, M) is the number of spin configurations with energy E and magnetization M. To overcome the barriers in both E and M spaces, we have to perform a two-dimensional (2D) random walk as in the spin glass problem [4].

We provide a succinct description of the sampling method here, and more details can be found elsewhere[3–5]. At the start of the simulation, the density of states is unknown, so we simply set g(E, M) = 1 for all possible (E, M). Then we begin a random walk in (E,M) space by choosing a site randomly and flipping its spin with a probability proportional to 1/g(E, M). If we denote A¢(E¢, M¢) and A¢¢(E¢¢, M¢¢) as the points before and after a spin is flipped, respectively, the transition probability from A¢ to A¢¢ is:

If the point A¢¢(E¢¢, M¢¢) is accepted we modify the existing density of states by a modification factor f > 1, that is, g(E¢¢, M¢¢) ® g(E¢¢, M¢¢) * f, and we update the histogram h of visited states, that is h(E¢¢,M¢¢) ® h(E¢¢,M¢¢)+1. If A¢¢(E¢¢, M¢¢) is not accepted, we update g(E¢,M¢) ® g(E¢,M¢) * f and h(E¢,M¢) ® h(E¢,M¢)+1.

We continue performing the random walk until the histogram h(E,M) is "flat" in (E,M) space. The modification factor f introduces a systematic error whose magnitude for ln[g(E,M)] is proportional to ln(f). To reduce this source of error, we systematically reduce the modification factor using a function like f_i+1 = f_i^1/n (n > 1) where i here denotes the level of iteration in the algorithm (in each iteration a separate 2D random walk is performed in (E,M) space). To be explicit, when iteration i generates a "flat" histogram, we reset the histogram to zero, and begin the next iteration with a new factor f_i+1 and the last g(E,M) as the starting value. We end the simulation when the modification factor is smaller than a predefined value (such as ln f_final = 10^-8 used in Sec.III and 10^-6 used in Sec.IV). To speed up the convergence of the density of states to the true value, the initial modification factor was as large as f = f₀ = e ~ 2.71828..., and n = 4 in our simulations.

We should point out that it is impossible to obtain a perfectly flat histogram. In this work the histogram is considered "flat" when the number of entries larger than or equal to 2000 remains unchanged for L²×10⁶ spin flip trials. For small lattices (up to L = 24) our tests showed that the results so obtained agree well with those found when using a more stringent flatness criterion, namely the histogram is considered "flat" when h(E,M) for all possible (E,M) is not less than 80% of the average histogram áh(E,M)ñ. For large lattices we may need to divide the two-dimensional parameter space into regions and perform separate random walks[3, 4, 13], because with our current computational resources it is very difficult to obtain a flat histogram for a single random walk over the entire (E,M) parameter space. However, problems arise from different choices of the parameter space decomposition, as we will discuss in Sec.IV.

With an accurate estimate of the density of states g(E,M), we can calculate thermodynamic quantities at any temperature T, and external magnetic field H for the system with Hamiltonian given in Eq.(1).

A mathematical proof of the convergence of the Wang-Landau method was given by Zhou and Bhatt [14] and applications of this method already include simulations of fluids [15] (continuum systems), protein folding [16], polymer films [17], polymer collapse [18–20], binary Lennard-Jones glass [21], liquid crystals [22], random spin systems [23, 24], atomic clusters [25], optimization problems [26], combinatorial number theory [27], Blume-Capel model [28], and 3D Potts model [29, 30]. This method has been improved by using N-fold way [31] and multibondic sampling [32, 33], and it also has been generalized to perform quantum Monte Carlo simulations [34, 35] and sampling along the reaction coordinates for a molecular system [36]. By suitable adaptation of the problem, multidimensional numerical integration has also been successfully treated using Wang-landau sampling [37]. A rigorous derivation for off-lattice implementation of this algorithm was given in Ref. 38. Further generalizations and studies of this sampling technique have also been carried out [13, 39–44].

III. ASYMMETRIC ISING MODEL

The model considered here is described by the Hamiltonian

where S_i = ±1 is an Ising spin on site i of a triangular lattice with linear size L, áijñ denotes nearest-neighbor spin-pairs, and áijkñ represents the three spins on the elementary triangles. The parameters J and J₃ are two- and three-body nearest-neighbor spin couplings, respectively, and H is an external field. Previous Monte Carlo simulations by Chin and Landau [6] showed that this model exhibits critical endpoints for a range of values of J and J₃. In particular, for J = 1 and J₃ = 2 (parameters used here) there is a CE that is well-separated from a critical point at the end of the first-order phase transition line.

We use Wang-Landau sampling[3–5] described in Sec.II with E = å_áijñS_iS_j + 2å_áijkñS_iS_jS_k and M = å_iS_i to determine the density of states g(E,M), and from it we obtain thermodynamical quantities such as the magnetization, specific heat, susceptibility, order parameter[45, 46], etc. Multiple, independent random walks were performed to look for unexpected behavior resulting from different random number streams as well as to allow the determination of error bars.

Figure 1 shows the phase diagram in the (T,H) plane extrapolated to L = ¥ using finite size scaling relations. For finite lattices, the phase boundaries were determined by the locations of the maximum of the susceptibility [46]. There are three distinct phases with different values of magnetization at T = 0: phase A(+ + +) with M(0,H) = 1 (all spins are up), phase B(- - -) with M(0,H) = -1 (all spins are down) and a ferrimagnetic phase C(- - +) with M(0,H) = -1/3 (the spins on two sublattices are down and on the other are up).

The critical endpoint for L = ¥ is estimated as (T,H)_CE = (2.443(10),-2.934(10)) and the Ising-like critical point at the end of the first-order line is at (T,H)_c = (6.125(3), -1.879(2)). Our tests[45, 46] showed that the critical line T_c(H) for this asymmetric Ising model is in the same universality class as the two-dimensional Q = 3 Potts model[47], as expected due to symmetry considerations. We recall that the conjectured values[47] for the critical exponents of the 2D Q = 3 Potts model are a = 1/3, b = 1/9, g = 13/9, and n = 5/6.

In Fig.2(left), the curvatures d²H_s(T)/dT² of the spectator phase boundary show a very clear singularity at the CE. According to a general finite-size scaling argument [7, 10], this curvature should diverge at the CE with a specific heat-like form

where a and n are critical exponents defined on the critical line T_c(H) and a₀ is a background. Using the conjectured exponents for the 2D Q = 3 Potts model[47], the predicted scaling exponent (see Eq. (4)) is a/n = 2/5.

Figure 2(right) shows a plot of the maximum of d²H_s/dT², denoted as f(L), as a function of L^a/n = L^0.4. The scaling relation is obtained by a linear fit to the data for L = 18 to 42 to Eq. (4), with a₀ and a₁ as fitting parameters and with the exponent fixed at a/n = 2/5 = 0.4. Because we obtained quite a good fit to the scaling function f(L) = a₀+a₁ L^0.4 [solid line in Fig. 2(right)] we can conclude that our data are in good agreement with the predicted scaling exponent (see Eq. (4)), but the background term is not negligible. Our data also indicate that there are small correction terms to the finite-size scaling; however, the resolution of our data and the lattice sizes used here are not adequate to estimate these correction terms.

The magnetization coexistence diameter M_d is defined as the average magnetization along the first-order transition line H_s(T). Using generalized scaling arguments, Wilding [10] predicted that the coexistence diameter derivative diverges as

where a, b, and n are critical exponents defined on the critical line T_c(H). We observe a clear divergence of the derivative of the magnetization coexistence diameter (see Fig. 3(left)) near the critical endpoint. Figure 3(right) shows the scaling of the maximum of -dM_d/dT with lattice size L. The solid line in Fig. 3(right) is a least squares fitting of the maximum of -dM_d/dT using Eq. (5) with L = 15 to 42 data. For the symmetric binary fluid studied by Wilding [10], the coefficient c₂ in Eq. (5) is predicted to vanish; therefore dM_d/dT should diverge with L^a/n at the CE. However, for the asymmetric model studied here, we obtained a better fit when all terms in the right-hand side of Eq. (5) were included.

IV. Q-STATE POTTS MODEL IN AN EXTERNAL FIELD

The second model considered here is described by the Hamiltonian

where s_i = 1,2,¼Q on site i of a two-dimensional L × L square lattice, áijñ denotes pairs of nearest-neighbor spins, J > 0 is a ferromagnetic coupling, and H is an external field. We use the Wang-Landau sampling[3–5] described in Sec.II with E = -å_áijñ,s_j and M = å_i,₁ to determine the density of states g(E,M) and from it we obtain the energy and magnetization probability distributions at different values of T and H.

Figure 4 shows the phase diagram for the Q = 10 Potts model, for lattices with L = 24,32,48 and the extrapolation to L = ¥ using finite size scaling. For finite lattices, three "1st-order" phase boundaries meet at a triple point which occurs for non-zero field; however, this point extrapolates to H = 0 in the thermodynamical limit. This means that the "ordinary", 1st order Potts model transition is actually a triple point. The phase transition lines are determined from equal height peaks in the magnetization probability distribution P(M,T,H). Near the triple point, P(M,T,H) has three peaks, corresponding to (i) a phase rich in state 1, (ii) a phase poor in state 1, and (iii) a disordered state, as illustrated in Fig.5. For each set of parameters (T,H) the right-most peak in Fig.5 corresponds to phase (i), the left-most peak to phase (ii), and the middle peak to (iii). For L = ¥ we observe a first-order phase transition line starting at the H = 0 phase transition point and ending at a critical point (T_c,H_c), shown as a solid circle in Fig.4. Note that for finite lattices P(M,T,H) continues to have two equal height peaks for H > H_c; however, the logarithm of the value of the peak height over the value of the minimum of the distribution between the two peaks decreases[48] with L for H above H_c. For each L the phase transition line separating phases (i) and (ii) goes to (T,H) = (0,0). The H = 0 phase transition temperature obtained in our simulation is T_pt = 0.701228(10), in good agreement with the exact result = 0.701231... .

For Q = 8, 9 we have also located a first-order phase transition line starting at the respective H = 0 phase transition points and ending at a critical point. The critical field H_c(Q) increases with Q, as shown in Fig.6, in qualitative agreement with previous theoretical predictions[12]. For smaller Q the critical field becomes tiny, and determining its location as Q approaches 4 will be challenging.

The phase diagrams for Q = 8,9,10 were obtained using L < 48, and a single random walk in the entire (E,M) space. For smaller Q values (Q = 5,6,7), the finite size effects are stronger and larger lattices are needed in order to locate the critical point accurately. However, with our current computational resources it is very difficult to carry out a single 2D random walk over the entire (E,M) parameter space for L much larger than 48. Parallel, one-dimensional Wang-Landau sampling, in which independent random walks are carried out in restricted phase space regions, has been implemented with success in previous work[3, 4, 13]. Therefore, we have also tried to parallelize the 2D algorithm by performing independent random walks on overlapping (E,M) regions. We divided the parameter space into strips with: (i) restricted values of E and unrestricted M, (ii) restricted M and unrestricted E, and (iii) both E and M restricted. The independent random walks were carried out in parallel, and each generated a density of states for the respective random walk regions. The densities of states from two adjacent random walk regions were then normalized with one overlapping (E,M) point. However, we found that the resulting density of states over the entire (E,M) range had some small discontinuities, which carried over to the probability distributions of energy and magnetization, as illustrated in Fig.7. For the division in case (i), the distribution P(M,T,H) agreed well with the result from a single random walk over the entire (E,M) space; however, the energy distribution P(E,T,H) had a number of discontinuities. For case (ii), P(E,T,H) was correct, but P(M,T,H) had a number of small jumps. The parallelization in (iii) also generated jumps in P(M,T,H). The general source of these difficulties seems to be due to the difficulty in matching surfaces at the boundaries rather than curves as in one-dimensional random walks. Thus, the simple parallelization described above cannot be used to obtain accurate phase diagrams. In the parallel Wang-Landau sampling of g(E) performed in Ref.4 the random walk within each energy region was restarted periodically from independent spin configurations, to ensure that all spin configurations with energies in the restricted range can be equally accessed. This approach might be helpful too for the parallel Wang-Landau sampling of g(E,M). Another viable simulational approach to study the Q = 5,6,7 models is to perform 2D Wang-Landau sampling on systems with L < 48 and use these results to guide more standard hybrid Monte Carlo algorithms for larger lattices at a few values of (T,H).

V. CONCLUSION

We used Wang-Landau sampling with a two-dimensional random walk to determine the density of states g(E,M) for an asymmetric Ising model with two- and three-body interactions on a triangular lattice, in the presence of an external field. An accurate g(E,M) allowed us to map out the phase diagram accurately and to observe a clear divergence of the spectator phase boundary and of the magnetization coexistence diameter derivative at the critical endpoint. The exponents for both divergences agree well with the predicted values[7, 10].

We also applied the 2D Wang-Landau sampling method to obtain g(E,M) for the Q-state ferromagnetic Potts model in the presence of an external field. We mapped out the phase diagram for Q = 8,9,10 and observed a first-order phase transition line starting at the H = 0 phase transition temperature and ending at a critical point (T_c,H_c) with H_c > 0. The critical field increases with Q, in qualitative agreement with previous large-Q series expansion results[12]. Our attempts to parallelize the algorithm demonstrated that further development is needed to avoid subtle errors that arise from the process.

Acknowledgments

We thank K. Binder, H. Meyer, and J. A. Plascak for helpful discussions, and P. Brunk for technical assistance with computer resources. Simulations were performed on an Itanium cluster at SDSC and on computer resources at the Research Computing Center at UGA. This work is partially supported by the NSF under Grants DMR-0341874 and DMR-0307082.

[1] D. P. Landau and K. Binder, A Guide to Monte Carlo Methods in Statistical Physics, second edition (Cambridge U. Press, Cambridge, 2005).

[2] A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 61, 2635 (1988); 63, 1195 (1989).

[3] F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001).

[4] F. Wang and D. P. Landau, Phys. Rev. E. 64, 056101 (2001).

[5] D. P. Landau and F. Wang, Comput. Phys. Commun. 147, 674 (2002); D. P. Landau, S.-H. Tsai, and M. Exler, Amer. J. Phys. 72, 1294 (2004); D. P. Landau and F. Wang, Braz. J. Phys. 34, 354 (2004). Y Y

[6] K. K. Chin and D. P. Landau, Phys. Rev. B 36, 275 (1987).

[7] M. E. Fisher and P. J. Upton, Phys. Rev. Lett. 65, 2402 (1990);65, 3405 (1990); M. E. Fisher, Physica A 172, 77 (1991). Y

[8] T. A. L. Ziman, D. J. Amit, G. Grinstein and C. Jayaprakash, Phys. Rev. B 25, 319 (1982).

[9] M. E. Fisher and M. C. Barbosa, Phys. Rev. B 43, 11177 (1991); M. C. Barbosa and M. E. Fisher, Phys. Rev. B 43, 10635 (1991); M. C. Barbosa, Phys. Rev. B 45, 5199 (1992). Y Y

[10] N. B. Wilding, Phys. Rev. Lett. 78, 1488 (1997); Phys. Rev. E 55, 6624 (1997). Y

[11] B. A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 (1992).

[12] Y. Y. Goldschmidt, Phys. Rev. B 24, 1374 (1981).

[13] B.J. Schulz, K. Binder, M. Müller, and D. P. Landau, Phys. Rev. E 67, 067102 (2003).

[14] C. Zhou and R. N. Bhatt, Phys. Rev. E 72, 025701(R) (2005).

[15] Q. Yan, R Faller, and J. J. de Pablo, J. Chem. Phys. 116, 8745 (2002); T. S. Jain and J. J. de Pablo, J. Chem. Phys. 118, 4226 (2003). Y

[16] N. Rathore and J. J. de Pablo, J. Chem. Phys. 116, 7225 (2002); N. Rathore, T. A. Knotts, and J. J. de Pablo, J. Chem. Phys. 118, 4285 (2003). Y

[17] T. S. Jain and J. J. de Pablo, J. Chem. Phys. 116, 7238 (2002).

[18] F. Rampf, K. Binder, and W. Paul, J. Polym. Sci. Part B: Polym. Phys. 44, 2542 (2006).

[19] D. F. Parsons and D. R. M. Williams, Phys. Rev. E 74, 041804 (2006); J. Chem. Phys. 124, 221103 (2006). Y

[20] D. T. Seaton, S. J. Mitchell, and D. P. Landau, Braz. J. Phys. 36, 623 (2006).

[21] R. Faller and J. J. de Pablo, J. Chem. Phys. 119, 4405 (2003).

[22] E. B. Kim, R. Faller, Q. Yan, N. L. Abbott, and J. J. de Pablo, J. Chem. Phys. 117, 7781 (2002).

[23] Y. Okabe, Y. Tomita, and C. Yamaguchi, Comput. Phys. Commun. 146, 63 (2002).

[24] Y. Wu and J. Machta, Phys. Rev. B 74, 064418 (2006).

[25] F. Calvo and P. Parneix, J. Chem. Phys. 119, 256 (2003).

[26] M. A. de Menezes and A. R. Lima, Phys. A 323, 428 (2003).

[27] V. Mustonen and R. Rajesh, J. Phys. A 36, 6651 (2003).

[28] C. J. Silva, A. A. Caparica, and J. A. Plascak, Phys. Rev. E 73, 036702 (2006).

[29] C. Yamaguchi and Y. Okabe, J. Phys. A. 34, 8781 (2001).

[30] Y. Okabe and H. Otsuka, J. Phys. A 39, 9093 (2006).

[31] B. J. Schulz, K. Binder, and M. Müller, Int. J. Mod. Phys. C 13, 477 (2002).

[32] C. Yamaguchi and N. Kawashima, Phys. Rev E 65, 056710 (2002).

[33] B. A. Berg and W. Janke, Phys. Rev. Lett. 98, 040602 (2007).

[34] M. Troyer, S. Wessel, and F. Alet, Phys. Rev. Lett. 90, 120201 (2003).

[35] P. N. Vorontsov-Velyaminov and A. P. Lyubartsev, J. Phys. A 36, 685 (2003).

[36] F. Calvo, Mol. Phys. 100, 3421 (2002).

[37] Y. W. Li, T. Wüst, D. P. Landau, and H. Q. Lin, Comput. Phys. Commun. 177, 524 (2007).

[38] M. S. Shell, P. G. Debenedetti, and A. Z. Panagiotopoulos, Phys. Rev. E 66, 056703 (2002).

[39] Q. Yan and J. J. de Pablo, Phys. Rev. Lett. 90, 035701 (2003).

[40] N. Rathore, Q. Yan, and J. J. de Pablo, J. Chem. Phys. 120, 5781 (2004).

[41] E. A. Mastny and J. J. de Pablo, J. Chem. Phys. 122, 124109 (2005).

[42] A. Malakis, A. Peratzakis, and N. G. Fytas, Phys. Rev. E 70, 066128 (2004); A. Malakis, S. S. Martinos, I. A. Hadjiagapiou, N. G. Fytas, and P. Kalozoumis, Phys. Rev. E 72, 066120 (2005); A. Malakis and N. G. Fytas, Phys. Rev. E 73, 056114 (2006); 73, 016109 (2006). Y Y

[43] C. Zhou, T. C. Schulthess, S. Torbrügge, and D. P. Landau, Phys. Rev. Lett. 96, 120201 (2006).

[44] H. K. Lee, Y. Okabe, and D. P. Landau, Comput. Phys. Commun. 175, 36 (2006).

[45] S.-H. Tsai, F. Wang, and D. P. Landau, Braz. J. Phys. 36, 635 (2006).

[46] S.-H. Tsai, F. Wang, and D. P. Landau, Phys. Rev. E 75, 061108 (2007).

[47] F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982).

[48] J. Lee and J. M. Kosterlitz, Phys. Rev. Lett. 65, 137 (1990).

Received on 28 September, 2007

*

Current address: Intel Corporation, PC1-105, 44235 Nobel Drive, Fremont, CA 94538

Publication Dates

• Publication in this collection
  27 Mar 2008
• Date of issue
  Mar 2008

History

• Received
  28 Sept 2007