Bond counting Monte Carlo for the 2-D Ising model in an external field

Bond counting Monte Carlo for the 2-D Ising model in an external field

H. G. Dias and J. Florencio

Departamento de Física, Universidade Federal de Minas Gerais,

30.161-970 Belo Horizonte, MG, Brazil

Received on 15 August, 2000

We propose an improvement of a Monte Carlo method designed to treat the Ising model in a field [C. Lieu and J. Florencio, J. Low Temp. Phys. 89 , 565 (1992)]. The method involves the counting of bonds linking neighboring like-spins and yields the degeneracy of the system's energy states, hence the partition function. There is no acceptance-rejection procedure and all the randomly generated configurations are kept. The sampling depends on geometry only, so results of a given run can be used for all temperatures and energy parameters. In order to understand the virtues and inadequacies of the method, we obtained exact results for small lattices. We find that a Monte Carlo run must be followed by a Gaussian fit in order to account properly for the rare events not recorded in the sampling. Finally, we also established bounds for the location of the peak for the specifc heat of the Ising model in a magnetic field in two dimensions for several values of the field in the thermodynamic limit.

The thermodynamics of the Ising model in one and two dimensions has been known for quite some time [1, 2]. In 1D, there is no long-range order except at zero temperature. On the other hand, in 2D the model shows a transiton from a paramagnetic to a ferromagnetic state as the temperature is lowered through T_c, the critical temperature [3]. The model can help the understanding of many phenomena observed in some magnetic systems. In addition, it constitutes the starting point for many other systems found in nature.

Monte Carlo computer simulation has been used in the last few decades to study the model in dimensions D > 2 [4]. In spite of the inherent limitations of the lattices sizes used in computer simulations, results can often be extended to the thermodynamic limit. Most Monte Carlo simulations are based on the Metropolis algorithm, which employs a rejection-acceptance procedure to enforce detailed balance [5]. Despite its successes, those simulations face some difficulties. For instance, near the transition temperature it takes a long computer time to equilibrate the system due to critical slowing down. Another point to consider is that distinct runs are needed to obtain the temperature dependence of a given thermodynamic quantity. The determination of the location of narrow peaks of thermodynamic quantities is not very accurate, which is the case in the vicinity of the phase transition. Even though several methods have been proposed to deal with these problems, we present in this work an improvement of the Monte Carlo method of Lieu and Florencio [6] which does not involve the rejection-acceptance feature of the Metropolis method, and that might prove useful in the study of the Ising model.

The Monte Carlo method considered here involves random sampling of spin configurations, with all them taken into account. The method allows for an estimation of the degeneracies of the energy states, leading directly to the partition function and the ensuing thermodynamic functions. Its main advantage is that the sampling depends on geometry only. A single run is valid for all the energy parameters J and B, as well as the temperature T. The method yields analytical expressions, so there is no difficulty in dealing with narrow peaks. The configurations are generated independently, thus critical slowing down is not a factor. We shall present here results for the 2D Ising model in an external field. However the method can be used to deal with higher dimensional lattices.

Consider the Ising model in a magnetic field on a square lattice of N = L × L sites,

with periodic boundary conditions s_i+L,j+L = s_i,j, where s_i,j can take one of the two values, 1 (up-spin) or -1 (down-spin), B is the external magnetic field and J (FM or AFM) is the coupling parameter between neighboring spins.

Consider a state of the system which has p up-spins (and N - p down-spins). Suppose the up-spins form m nearest neighbor pairs (bonds) between themselves. The energy of such state is then E₀ + pD₁ + mD₂, where E₀ = N (B - 2 J) is the energy of the state where all the spins are pointing down, D₁ = 8J - 2B is the energy needed to flip one down-spin originally surrounded by neighboring down-spins, and D₂ = -4J is the correction in case there exists a single up-spin located at a site neighboring the flipped spin. Hence, the partition fuction can be expressed as

where h = e^-bD₁ and x = e^-bD₂. The quantity gives the number of states containing p up-spins forming m bonds. The determination of , the only unknown in the partition function, is a very difficult combinatorics problem. The knowledge of for any lattice size completely solves the thermodynamics of the 2D Ising model in a field.

The core of the method is to determine an average of these coefficients by means of sampling. The a priori probability of picking up any given state with p up-spins is ()^-1, regardless the energy of that state. Therefore the determination of is purely geometric. Since each configuration has the same probability to occur and they are statistically independent, the variance of decreases as (Np^MC)^-1, where Np^MC is the number of elements in the sample of states with p up-spins.

We ask the computer to generate a configuration with p up-spins distributed randomly on the lattice, count the number of bonds m formed in such configuration, and store the result. The process is repeated a number of times and, by counting how many configurations yielded m bonds, we obtain the configurational average . The proper number of states is ensured by the normalization condition å_m = (). Thus, the determination of is simply geometrical in nature, hence the numerical values of the model parameters and of the temperature are not needed in a Monte Carlo run.

By exactly analyzing the m-dependence of for several values of p for small lattice sizes, we find that it can be fitted with a Gaussian. As an illustration, Fig. 1 depicts a Gaussian fit (dashed line) to the the exact values of (circles) on a 7 × 7-lattice in the case where p = 21. We also find that a Gaussian fit becomes better as the lattice size increases. (It is not so good for very small lattice sizes, such as 2 × 2, 3 × 3, etc.) Given that the overwhelming majority of the states have values of m distributed around the Gaussian peak, random sampling will pick up mostly those states. The states in the Gaussian tails will be undercounted, if not neglected entirely. That poses a problem for the method, since the contributions from those states and similar states from Gaussian tails of different p's may go uncounted. In short, simple random sampling favors counting energy states from around the peaks of , while neglecting states from its tails. In order to account for those states, which are unlikely to be chosen by random sampling, we propose to improve the method by using a Gaussian fit to the results obtained from random sampling of the configurations. Another point to consider is that some states located in the tail of the Gaussian fit never occur. For example, in the case of p = 4 and L = 5 or greater, configurations with m = 5, 6, 7 and 8 never occur, because we can have at most four bonds with just four up-spins in a lattice size where L ³ 5. Since the partition function, given by Eq. (2), involves a summation over m up to 2p, we have to exclude such impossible states. We have found a way to track these configurations and it is actually a simple matter to get rid of them. Hence, after the simulation is done, we perform a Gaussian fit and normalize after the exclusion of those impossible states (we set = 0 for such states).

Figure 1.
Exact number of spin configurations on a 7 × 7 lattice versus the number the bonds formed by neighboring up-spins. In the figure shown, the number of up-spins is p = 21.

These improvements of the method turn out to produce a more accurate results. We tested it for the 6 × 6 and the 8 × 8 lattices with J = 1 and B = 0. For the 6 × 6 lattice, we performed 96 Monte Carlo runs, each with 100000 samples per value of p. The exact value for the temperature which maximizes the specific heat is 2.389709. We obtain a value of 2.33 ± 0.06 for the peak temperature. We find equally satisfactory results for the 8 × 8 lattice. Fig. 2 shows our results for the specific heats of those lattice sizes.

Figure 2.
Specific heat using our Monte Carlo method for the 6 × 6 and 8 × 8 lattices, where J = 1 and B = 0.

Finally, we determined bounds on the location of the peaks of the specific heat of the model in the presence of a field, in the thermodynamic limit. It is well known that unlike the case where B = 0, where the specific heat diverges at the critical temperature, there is no divergence of the specific heat when B ¹ 0. Nevertheless, the specific heat has a peak centered at a given temperature, which we denote by T_p, the peak temperature. First, we determine T_p for small lattice sizes, up to the size N = 6 × 6. Fig. 3 shows the peak temperature as a function of 1/N, the reciprocal of the lattice size. The cases (J, B) = (1, 0) and (1, 0.5) are displayed in the figure. By extrapolation, we can infer bounds for the location of the peak temperatures in the thermodynamic limit. The results, for several values of B are given in Table I. We also find that as the field increases the upper and lower bounds approach each other.

Table I.
Bounds for T_p, the location of the specific heat peak, for several values of B.

Figure 3.
Temperature at which the specific heat is peaked vs the reciprocal of the system size, for B = 0 and B = 0.5. The coupling energy J = 1. The dashed lines are just a guide to the eye.

We need, of course, to test the method on larger lattices in order to make more precise predictions on the influence of the field B on the thermal properties of the Ising model. The extension of the method to higher dimensions is straightforward.

Acknowledgments

This work was partially supported by CNPq, FAPEMIG, and MCT (Brazilian agencies).

• [1] E. Ising, Z. Physik 31, 253 (1925).
• [2] L. Onsager, Phys. Rev. 64, 117 (1944).
• [3] B.M. McCoy and T.T. Wu, The two dimensional Ising model (Harvard Univ. Press, Cambridge, Mass., 1973).
• [4] See, e.g. , Computer Simulation in Condensed Matter Physics , Eds. D.P. Landau, S.P. Lewis, and H.-B. Schlütter (Springer-Verlag, New York, 2000).
• [5] N. Metropolis, W.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953).
• [6] C. Lieu and J. Florencio, J. Low Temp. Phys. 89, 565 (1992).

Publication Dates

• Publication in this collection
  09 Jan 2002
• Date of issue
  Dec 2000

History

• Received
  15 Aug 2000