Systematic enumeration of configuration classes for entropic sampling of Ising models
Departamento de Física, Instituto de Ciências Exatas, and National Institute of Science and Technology for Complex Systems,
Universidade Federal de Minas Gerais
C.P. 702, 30123-970, Belo Horizonte, MG, Brazil
We describe a systematic method for complete enumeration of configuration classes (CCs) of the spin-1/2 Ising model in the energy-magnetization plane. This technique is applied to the antiferromagnetic Ising model in an external magnetic field on the square lattice, which is simulated using the tomographic entropic sampling algorithm. We estimate the number of configurations, , and related microcanonical averages, for all allowed energies and magnetizations for to 30, with . With prior knowledge of the CCs, we can be sure that all allowed classes are sampled during the simulation. Complete enumeration of CCs also enables us to use the final estimate of to obtain good initial estimates, , for successive system sizes () through a two-dimensional interpolation. Using these results we calculate canonical averages of the thermodynamic quantities of interest as continuous functions of temperature and external field . In addition, we determine the critical line in the - plane using finite-size scaling analysis, and compare these results with several approximate theoretical expressions.
Keywords: Ising model; antiferromagnet; configuration classes; Monte Carlo simulation; phase diagram.
The most important task of entropic sampling algorithms  is to visit the full configuration space (CS) to obtain good estimates of the number of configurations, , as functions of the energy E and other quantities of interest. In studies of the Ising model in an external field, for example, we require with the magnetization; each allowed pair defines a class of configurations (CC). Only if we know beforehand the possible values of for a given system size, can we be sure that all CCs are sampled during the simulation.
Figure 1 shows the CCs in the plane for the spin-1/2 Ising model with nearest-neighbor (NN) interactions, on a square lattice of sites with periodic boundaries. [We use to denote the number of NN pairs of spins with the same orientation; the interaction energy of the antiferromagnetic (AF) Ising model is .] Although the CCs tend to fill in a triangular region, some gaps are evident near the lower apex and along the upper edge. Knowing just which values are allowed for a given lattice size is important if we are to implement entropic sampling with confidence. In this paper we present a method for systematically enumerating all CCs of Ising models in the plane.
We study the spin-1/2 antiferromagnetic Ising model in an external magnetic field, whose energy is given by
where , indicates a sum over NN pairs of spins, is the external field, and is the number of spins; the model is defined on a square lattice of sites, with periodic boundary conditions. Unlike the ferromagnetic Ising model (), which exhibits a unique critical point in the plane and has an exact solution , the AF model () possesses a critical line, which is not completely understood.
Various approximate methods have been applied to determine the critical line of the AF Ising model on the square lattice ; these results, however, do not agree altogether. Binder and Landau  estimated the critical line via Monte Carlo simulation, obtaining very good agreement with the approximate closed-form expression of Müller-Hartmann and Zittartz , raising the possibility that the latter expression was in fact exact. (It was later shown that this is not so .) Hwang et al.  studied the AF Ising model on the square and triangular lattices using a microcanonical transfer matrix method and the Wang-Landau algorithm . They performed an exact enumeration of the number of configurations, , and found CSs in the plane similar to that shown in Fig. 1.
Using the tomographic entropic sampling (TES) algorithm  we estimate , and associated microcanonical averages, for lattice sizes to 30, with an increment . We then calculate the canonical averages of the thermodynamic quantities of interest. Using these data we map out the critical line in the plane, and compare our results with several theoretical expressions.
Prior determination of the set of allowed CCs is an important tool to verify the quality of the sampling: we want to be sure that all CCs are visited during the simulation. Since this algorithm uses an initial guess, , to begin the study, it is convenient to use the final estimate (after the -th iteration) to obtain the initial guess of the next system size to be studied . As we will show, good initial estimates, , can be obtained using a two-dimensional interpolation because we know a priori the set of allowed CCs.
This paper is organized as follows. In Sec. 2 we define the basic CCs and the respective allowed values of for the spin-1/2 Ising model on the square lattice; the main goal of this section is to find all gaps in the plane. In Sec. 3, this information is used in simulations of the AF Ising model via the TES algorithm. There we describe the method used to determine via two-dimensional interpolation of the final estimate, , of the previous system size studied. Simulation results are reported in Sec. 4, for the order parameter and the staggered susceptibility as functions of and . Points along the critical line in the plane are obtained using finite size scaling analysis, and the results compared with several theoretical expressions. We summarize our findings in Sec. 5.
2 Configuration Classes
As pointed out above, one of the main problems in entropic sampling methods is the prior determination of the complete set of configuration classes for a given system size. Let us denote by and the number of up and down spins, respectively, on a square lattice with spins. The number of pairs of opposite spins, , and are related to and , respectively, via
Thus, once the possible values of are determined, so are those of . Note that , and can only take even values. The gaps in CS fall near the maximum and minimum values of ( and , respectively) for a given . Therefore, we will identify the possible values of near its maximum, , and minimum, , for a given ; note that the number of configurations is symmetric under interchange of and .
2.1 Determining and nearby classes
2.1.1 Compact configurations
Compact configurations consist of a square or rectangular cluster with up spins. To begin, consider the case of a square cluster of size with , as is illustrated in Fig. 2. It is evident that this configuration corresponds to the minimum value of , . A configuration with the same number of up spins, and is obtained by transferring an up spin from one of the corners of the square to an edge, as shown in Fig. 3. From this configuration, further rearrangements leading to , etc., are possible. When is not a square number, the most compact configuration (i.e., with the smallest perimeter) is a rectangle, or a square or rectangle with an incomplete layer of sites along one edge. For we have , while for , . In all cases, moving a corner site to an edge, one constructs a configuration with , and further arrangements yield additional increases in .
2.1.2 Extended configurations
Thus there are no gaps in the large- region due to compact configurations. Such compact configurations, however, are not necessarily the ones that minimize for a given value of . Consider for example the case with an integer. The cluster of up spins can be arranged to wrap around the lattice, yielding , independent of . We call such a configuration extended. For greater than a certain value, on the order of , the configurations that minimize are extended rather than compact. Suppose now that for a given , is obtained with an extended configuration, and that all compact configurations have strictly greater than . If lies between and , then the configuration that minimizes has at least two corners, and by the previous construction, configurations with exist. But if , the minimizing configuration has no corners and any modification yields a configuration with . This is how the gaps near the maximum values of arise.
Summarizing, if and are such that is obtained with a compact configuration, then there are configurations with , , …, etc., and no gap exists. If, on the other hand, is realized only for an extended configuration, and is an integer multiple of , then there are no configurations with .
2.2 Determining and nearby classes
The largest possible value of , , occurs in a configuration such that with spins arranged in a chess board (CB) configuration, such that all up spins have down spins as NNs and vice versa. One readily verifies that for , the maximum number of unlike NN pairs is . Starting from the CB configuration, we can reduce by exchanging an up and a down spin. If the exchanged spins are NNs, the resulting configuration has ; otherwise one has . Thus, for , there are no configurations with or . One readily verifies that configurations with , , etc., can be obtained via further exchanges of spins.
For , we have . Such configurations can be constructed by flipping one up spin in the CB configuration. Starting from this configuration, exchanging a pair of spins, one can reduce by 4, 6 or 8, but there is no rearrangement which reduces by just two. Thus for , there is no configuration having ; configurations with , , etc. do exist.
Finally, we note that for , there are no gaps in the neighborhood of . To verify this, consider a configuration obtained by flipping up spins in the CB configuration, so that , where . By hypothesis there are now at least four more sites with down spins than with up spins. Starting from a configuration with , in which each up spin is surrounded by down spins, we can create a single NN pair of up spins, with all six neighbors down, and with all remaining up spins completely surrounded by down spins. In this manner, is reduced by two. Configurations with , , etc. can be obtained by further exchanges of up and down spins.
Using the facts summarized above, it is straightforward to construct an algorithm that determines which values of are possible, for a given and , and thereby which values are accessible for a given system size. In the simulations reported below, we have verified that our entropic sampling scheme converges to visit all allowed classes.
Using tomographic sampling, we study the antiferromagnetic Ising model in an external field on the square lattice; we consider periodic boundary conditions and NN interactions. The CCs of the systems are defined in the energy-magnetization space .
The TES method is applied in order to generate estimates of . For the smallest system size we begin with a guess, , obtained using a mean-field approximation. For subsequent system sizes, however, we use a two-dimensional interpolation of (the final estimate of for the smaller system size, ) to construct , where . For most studies we use five iterations, each one with initial configurations, which are simulated for lattice updates or Monte Carlo steps.
Let us denote by and the CCs that contain configurations and , respectively. The simulation uses a single spin-flip dynamics, so that the possible variations of and are and . At iteration , the acceptance probability for the transition is
These probabilities are stored in a table. For each configuration generated, be it a new one (if it is accepted) or the same (if it is rejected), we update the sums used to calculate the microcanonical and canonical averages of , , , and , where
is the order parameter (staggered magnetization); are the magnetizations of the two sublattices. At the end of each iteration the estimate of is refined according to
where is the histogram containing the number of times the CC is visited during the sampling, and is the average of over all accessible CCs; the acceptance probability [Eq. (4)] is updated using the new estimate of and the histogram, , is set to zero.
A single iteration of our method consists of ten independent simulations, each involving lattice updates, and each beginning from a different initial configuration. (By a lattice update we mean one attempted flip per spin; the initial configurations include, high, low, and intermediate interaction energies and both signs of the magnetization.)
3.1 Determining : mean field approach
As noted above, for the smallest system size we use an estimate of obtained via a mean-field approximation, specifically
where , and
To obtain this expression, we first note that , and use Stirling’s approximation. The dependence on is then obtained by estimating and for a given and , using a random-mixing approximation, and supposing that follows a Gaussian distribution.
In Fig. 4 we plot the estimate of given by Eq. (7). We present in Fig. 5 the final estimated value of after the 5 iteration of the simulation; this result is quite similar to that presented by Hwang et al.  for the square lattice, which was obtained using Wang-Landau sampling . We note that the differences between the initial estimate obtained via mean field approximation and the final simulation result of are more evident near the edges of the CS, specially near to the maximum values of . To have a better idea of how close this initial estimate is to the final simulation result, we plot in Fig. 6 the difference between that estimate and the final result of the simulation for . Analyzing this figure, it is again clear that the differences are larger near the edges of the CS.
3.2 Determining : interpolation
We expect that the closer the estimate is to the (unknown) exact value, , the faster the simulation will converge, and the more accurate our final estimate will be. Since the mean-field estimate worsens as the system size grows, we only use it for the smallest system size studied; initial estimates of subsequent system sizes are obtained by interpolating the final estimate, , of the previous system size studied. Our procedure is based on the existence of the limiting microcanonical entropy density as a function of the intensive parameters and (the bond and magnetization densities, respectively),
where and . (Since and are restricted to even integers we have and similarly for .) The idea is then to write
where , , and the r.h.s. is evaluated by extending to noninteger and via extrapolation and interpolation. Using this approach, we obtain better estimates, as is shown in Fig. 7. (Note that the largest differences continue to fall along the edges of the CS.)
3.3 Extrapolation and interpolation
Suppose we wish to construct the initial estimate on the basis of the simulation results for a smaller system, . We could do this via interpolation if every CC of the larger system were surrounded by four points of the smaller one in the - plane. Fig. 8 shows, however, that along the edges of the CS, the points corresponding to classes of the larger system are not surrounded by four points of the smaller one. For those points, one could in principle use extrapolation rather than interpolation. We found, however, that direct extrapolation yields poor estimates for .
We obtain better estimates by first extrapolating the points along the edges of the CS for the smaller system, such that every point of the larger system is surrounded by four points of the smaller. Figure 9 shows the CS with accessible and extrapolated CCs for . Following this extrapolation we perform a linear two-dimensional interpolation as per Eq. (10).
In this section we present results of the AF Ising model on the square lattice in an external field; periodic boundary conditions are employed. We use TES to simulate systems of sizes to 30, with . To calculate the uncertainties we perform five independent studies for each system size. We plot in Fig. 10 the staggered magnetization (order parameter) per site, , as a function of at . We can see that decreases considerably between and ; this behavior suggests a critical point, , marking a phase transition from the AF to the paramagnetic state. Figure 11 shows the staggered susceptibility per site,
as a function of at . As grows the peaks tend to the critical point, . The specific heat per site, , (not shown) has a similar behavior.
4.1 Phase diagram
Using finite size scaling analysis  we estimate the critical line, , or, equivalently, , via the relations:
where is the field at which (the specific heat or the staggered susceptibility) takes its maximum for a given temperature and system size; is defined in an analogous manner. The estimates obtained using the maximum of and are averaged to yield and ; Figure 12 illustrates the procedure for estimating , for . Our estimated points for the phase boundary are shown in Fig. 13 along with the approximate expression derived by Müller-Hartmann and Zittartz :
Our simulation results are in good agreement with their expression. For the critical field, the greatest relative difference between theory and simulation is about 0.9%, which occurs at .
In Table 1 we compare our simulation estimates for the critical magnetic field with some theoretical approximations. For temperature we find good agreement with the estimates of Monroe  (whose analysis involves a free parameter ), Wu and Wu , and Blöte and Wu , whereas there are significant discrepancies in relation to the other approximations. At higher temperatures, differences appear between simulation and the predictions of Monroe, Wu and Wu, and Blöte and Wu. These may reflect a small systematic error or an underestimate of simulation uncertainties. We intend to examine this issue in greater detail in future work.
|TES: Tomographic entropic sampling |
|MHZ: Müller-Hartmann and Zittartz |
|WW: Wu and Wu |
|BW: Blöte and Wu |
|WK: Wang and Kim |
The complete enumeration of CCs is of fundamental importance to study the antiferromagnetic Ising model using tomographic entropic sampling. The determination of CCs for entropic sampling of Ising models also enables us to obtain good initial estimates for the configuration numbers , using two-dimensional linear interpolation; the initial estimate is reasonably close to the final estimate . Despite the relatively small system sizes used in this study, we obtain a good estimate for the critical line in the temperature-magnetic field plane. Further details on critical behavior will be published elsewhere .
We are grateful to CNPq and CAPES, Brazil, for financial support.
-  B. A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 (1992).
-  J. Lee, Phys. Rev. Lett. 71, 211 (1993).
-  F. Wang and D.P. Landau, Phys. Rev. Lett., 84, 10 (2001).
-  F. Wang and D.P. Landau, Phys. Rev. E, 63, 056101 (2001).
-  M.E.J Newman and G.T. Barkema. Monte Carlo Methods in Statistical Physics (Oxford University Press, New York, 2001), p. 161 - 169.
-  R. Dickman and A.G. Cunha-Netto, Phis. Rev. E, 84, 026701 (2011).
-  L. Onsager, Phys. Rev., 65, 117 (1944).
-  E. Müller-Hartmann and J. Zittartz, Z. Phys., 73, 261 (1977).
-  L. Sneddon, J. Phys. C, 12, 3051 (1979).
-  R.R. Santos, J. Phys. C, 18, L1067 (1985).
-  M. Kaufman, Phys. Rev. B, 36, 3697 (1987).
-  X.N. Wu and F.Y. Wu, Phys. Lett. A, 144, 123 (1990).
-  H.W.J. Blöte and X-N. Wu, J. Phys. A, 23, L627 (1990).
-  M. Badehdah, A. Benyoussef and M. Touzani, J. Magn. Magn. Mater., 172, 254 (1997).
-  X-Z. Wang and J.S. Kim, Phys. Rev. Lett., 78, 413 (1997).
-  X-Z. Wang and J.S. Kim, Phys. Rev. E, 56, 2793 (1997).
-  J.L. Monroe, Phys. Rev. E, 64, 016126 (2003).
-  S.J. Penney, V.K. Cumyn and D.D. Betts, Physica A, 330, 507 (2003).
-  K. Binder and D.P. Landau, Phys. Rev. B, 21, 5 (1980).
-  C-O. Hwang, S-Y. Kim, D. Kang and J.M Kim, J. Stat. Mech: Theory Exp., 2007, L05001 (2007).
-  V. Privman. Finite Size Scaling Analysis and Numerical Simulations of Statistical Systems (World Scientific, London, 1990).
-  B.J. Lourenço and R. Dickman, in preparation.