A new approach to Monte Carlo simulations in statistical physics

A new approach to Monte Carlo simulations in statistical physics

D. P. Landau^I; F. Wang^II

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

^IIParacer Inc., 1371 Shorebird Way Mountain View, CA 94043, USA

ABSTRACT

We describe a new algorithm that approaches Monte Carlo simulation in statistical physics in a different way. Instead of sampling the probability distribution at a fixed temperature, a random walk is performed in energy space to directly extract an estimate for the density of states. The canonical probability can then be found at any temperature by weighting by the appropriate Boltzmann factor, and thermodynamic properties can be determined from suitable derivatives of the partition function.

1 Introduction

During the last half of the 20th century great strides were made in the theoretical and experimental investigation of condensed matter systems, and consequently our understanding of phase transitions was dramatically enhanced. Certain problems remained, however,in part because imperfections in physical systems limited the resolution of experiment and in part because analytical solutions of all but the simplest theoretical models proved elusive. Computer simulation now plays a significant role in statistical physics [1], particularly for the study of phase transitions and critical phenomena. The ''standard'' Monte Carlo method, which is flexible and easy to implement, has been the Metropolis importance sampling algorithm [2], but near phase transitions this method becomes inefficient because of the long time scales that develop. A number of more efficient algorithms have now been introduced to circumvent this problem, but these are somewhat limited in applicability. Beginning with the seminal work of Swendsen and Wang [3] and extended by Wolff [4], cluster flipping algorithms have been used to reduce critical slowing down near 2nd order transitions. An approach which is quite different in spirit, the multicanonical ensemble method [5-8], was introduced to overcome the tunneling barrier between coexisting phases at 1st order transitions and has general utility for systems with a rough energy landscape [7-17]. In all cases, histogram reweighting techniques [18] can be applied in the analysis to increase the amount of information that can be gleaned from simulational data, but the applicability of reweighting is severely limited in large systems by the statistical quality of the ''wings'' of the histogram. This latter effect is quite limiting in systems with frustration that produces a very complicated energy landscape and limits the efficiency of ''standard'' methods.

The partition function of any model can be expressed in terms of a density of states g(E), i.e. the number of all possible states(or configurations) for an energy level E of the system, but direct estimation of this quantity has not been the goal of simulations. Instead, most conventional Monte Carlo algorithms [1], such as Metropolis importance sampling, Swendsen-Wang cluster flipping, etc. generate a canonical distribution at a given temperature. Such distributions are so narrow that multiple runs are usually needed to describe thermodynamic quantities over a significant range of temperature. In the canonical distribution, g(E) does not depend on the temperature and knowledge of g(E) permits construction of canonical distributions at any temperature. (Once g(E) is known the partition function Z can be calculated:

where b = 1/(k_BT), and the model is essentially ''solved" since most thermodynamic quantities can be extracted from it.)

Even though Monte Carlo methods are already very powerful [1], there was no efficient algorithm to calculate g(E) very accurately for large systems. Even for exactly solvable models such as the 2-dim Ising model, g(E) cannot be calculated exactly for a large system [19]. All methods based on accumulation of histogram entries, such as the histogram method of Ferrenberg and Swendsen [18] , Lee's version of the multicanonical method(entropic sampling) [20], the broad histogram method [21-24] and flat histogram method [25] have the problem of scalability for large systems. Thus, an algorithm is still needed to calculate the density of states for large systems.

We now describe a simple new, general, and efficient Monte Carlo algorithm (generally known as ''Wang-Landau sampling") that offers substantial advantages over existing approaches [29-31]. Unlike conventional Monte Carlo methods that generate a canonical distribution at a given temperature , this approach estimates g(E) directly and accurately via a random walk that produces a ''flat'' histogram in energy space. The estimate for g(E) is improved at each step of the random walk, using a carefully controlled modification factor, to produce a result that converges to the correct value very quickly. The thermodynamic quantities can then be extracted by applying canonical average formulas in statistical physics.

2 The ''Wang-Landau" algorithm

''Wang-Landau sampling" is based on the observation that if we perform a random walk in energy space by flipping spins randomly for a spin system, and the probability to visit a given energy level E is proportional to the reciprocal of the density of states , then a flat histogram is generated for the energy distribution. This is accomplished by modifying the estimated density of states in a systematic way to produce a ''flat'' histogram over the allowed range of energy and simultaneously making g(E) converge to the true value. g(E) is modified constantly during each step of the random walk and use the updated density of states to perform a further random walk in energy space. The modification factor f is controlled carefully, and at the end of simulation it should be very close to 1 which is the ideal case of the random walk with the true density of states.

Initially, g(E) is a priori unknown, so the simplest approach is to set all entries to g(E) = 1 for all possible energies E. Then a random walk in energy space is begun by flipping spins randomly with the probability at a given energy level proportional to . If E₁ and E₂ are energies before and after a spin is flipped, the transition probability from energy level E₁ to E₂ is:

Each time an energy level E is visited, the existing value is multiplied by the modification factor f > 1, i.e. g(E) ® g(E) * f. (It is preferable to work with the logarithm i.e.ln(g(E)) ® ln(g(E)) + ln(f) so that all possible g(E) will fit into double precision numbers.) If the random walk rejects a trial move, g(E) is also modified with the same modification factor. A reasonable, although not necessarily optimal, choice of the initial modification factor is f = f₀= e¹@ 2.71828... which allows g(E) to develop rapidly. If f₀ is too small, the random walk will take too long to reach all possible energies; however, large statistical errors result if f₀ is too large. During the random walk, the histogram H(E) (the number of visits at each energy level E) is accumulated; and when it is approximately flat, the density of states will have converged to the true value with an accuracy proportional to that modification factor ln(f). The modification factor is then reduced e.g. f₁ = , the histogram is re-zeroed, and a new random walk is begun. This iterative procedure continues until the modification factor is smaller than a predefined value (e.g. f_final = exp(10^–8) @ 1.00000001). The modification factor acts as a control parameter for the accuracy of g(E) during the simulation and also determines how many MC sweeps are necessary for the entire simulation. It is impossible to obtain a perfectly flat histogram and the phrase ''flat histogram'' that all histogram entries are not less than x% of the average histogram áH(E)ñ, where x% is chosen according to the size and complexity of system and the desired accuracy of g(E).

Clearly, one essential constraint is that g(E) should converge to the true value. The accuracy of the estimate for g(E) is proportional to ln(f) at that iteration; however, ln(f_final) can not be chosen arbitrarily small or the modified ln(g(E)) will not differ from the unmodified one to within the number of digits in the double precision numbers used in the calculation. If this happens, the algorithm no longer converges to the true value, and the program may run forever. Even f_final is within range, but too small, the calculation might take excessively long.

The method can be further enhanced by performing multiple random walks, each for a different range of energy, either serially or in parallel fashion, with the random walk restricted to remain in the desired range. The resultant pieces of g(E) can then be joined together and used to produce canonical averages for the calculation of thermodynamic quantities at any temperature.

Note that the algorithm does not satisfy the detailed balance condition exactly (especially for early stages of iteration) since g(E) is modified constantly during the random walk. However, after many iterations it quickly converges to the true value as f ® 1. If p(E₁® E₂) is the transition probability from energy level E₁ to level E₂, from Eqn. (2), the ratio of the transition probabilities from E₁ to E₂ and from E₂ to E₁ is:

where g(E) is the density of states. In another words, the random walk algorithm satisfies detailed balance:

where is the probability at the energy level E₁ and p(E₁® E₂) is the transition probability from E₁ to E₂ for the random walk. We conclude that the detailed balance condition is satisfied with accuracy proportional to the modification factor ln(f).

Almost all recursive methods update g(E) by using the histogram data directly only after enough histogram entries are accumulated [5, 8, 11, 13, 14, 15, 16, 27, 32, 33, 34]. Because of the exponential growth of g(E) in energy space, this process is inefficient because the histogram is accumulated linearly. Instead, we modify g(E) at each step of the random walk, and this allows us the approach to the true distribution to be much faster than conventional methods, especially for large systems. (Although histogram entries are accumulated during the random walk, they are only used to check if the histogram is flat enough to go to the next level iteration.)

Although the total number of configurations increases exponentially with the size of the system, the total number of possible energy levels increases linearly with the size of system so it is easy to calculate g(E) with a random walk in energy space for a large system. For example, for a q-state Potts model on a L × L lattice with nearest-neighbor interactions [35] the number of possible energy levels for q > 3 is about 2N, where N = L² is the total number of lattice sites. In contrast, the average number of possible states for each energy level is as large as , where q^Nis the total number of possible configurations of the system. Clearly, current computers are unable to realize all possible states to calculate any thermodynamic quantities even though most models in statistical physics are well defined. This is also why efficient and fast simulational algorithms are required in the numerical investigations.

With the density of states, we not only can calculate most of thermodynamic quantities in all temperature region without multiple simulations but can also access some quantities, such as the free energy and entropy, which are not directly available from conventional Monte Carlo simulations. The free energy is calculated using

All other thermodynamic properties can also be calculated from the density of states, e.g. the internal energy is given by:

and the specific heat can be estimated from the fluctuations in the internal energy:

Thus, a single simulation should be sufficient to provide all interesting information about a system.

3 How well does the method work?

3.1 Application to the Ising square lattice

As a test of the accuracy and convergence of the method, we apply it first to L × L Ising square lattices with nearest neighbor coupling. Not only is the Ising model the ''fruit fly'' of statistical mechanics but so many of its properties can be determined exactly so a truly quantitative check of the method is possible (see e.g. Ref. [37]). The exact solution for the partition function for finite-size systems [38] provides information about thermodynamic properties. Using Mathematica the exact density of states can be obtained for small systems [19], and we extended the calculation up to L = 50 with the Mathematica program provided by Beale. With the algorithm described here, g(E) was estimated for the Ising model for L = 50, with the final modification factor being 1.000000001, and these results are compared with exact values in Fig.3.1. To within the resolution of the figure, simulational data and the exact solution overlap perfectly, so in the inset we show the relative error:

Over most of the region e(log(g)) is smaller than 0.01%. Such errors for low energy levels are directly related to the errors for the thermodynamic quantities calculated from the density of states. The average relative error is 0.019% for L = 50. In comparison, a recent broad histogram study of a 2D Ising model on a 32×32 lattice [24] yielded an average deviation of the microcanonical entropy that was substantially larger than that from this random walk algorithm.

The normalized density of states in Fig. 1 is obtained by normalizing to the condition that the number of ground states is 2 (i.e. all up or all down).

For very large values of L we can divide the desired energy region [–2, 0.2] into multiple energy segments and estimate the density of states on each segment with an independent random walk. The resultant density of states can be joined from adjacent energy segments. A comparison of thermodynamic properties determined in overlapping ranges of energy showed that data are actually only needed over restricted ranges, e.g. results for three independent random walks performed in the ranges E/N = [–1.7, –1.2], E/N = [–1.8, –1.1] and E/N = [–1.9,–1.0] showed no discernable difference to the temperature dependence near the phase transition. Boundary effects at the edges of the energy ranges can be eliminated by simply treating spin flips that are rejected because the new energy would be out of range in the same manner as any rejected spin-flip.

Since the exact g(E) is only available for small systems(L < 50), it is also important to test the accuracy of estimates for thermodynamic quantities of large systems. Using canonical average formulas to calculate internal energy, specific heat, Gibbs free energy and entropy. Ferdinand and Fisher [38] the exact solutions of above quantities for 2D Ising model on finite-size lattices can be compared with simulational results.

In Fig. 2(a), we show the entropy as a function of temperature as determined from the density of states. The exact and simulational data overlap almost perfectly over a wide temperature region from T = 0 to T = 8. Since no difference is visible with the scale used in Fig. 2(a), a more stringent test of the accuracy is provided by the inset which shows the relative errors e(S) For L = 256 the relative error is smaller than 0.1% over almost the entire temperature range.

g(E) can also be used to calculate the specific heat as a function of temperature. Both simulational data and exact results near T_c are shown in Fig. 2(b) for L = 64, 128 and 256 Ising model. To within the error bars there is no difference between the simulational data and exact solutions. In the inset of Fig. 2(b), we show relative errors as a function of temperature for L = 256. The errors for specific heat for a 256 × 256 lattice are smaller than 4.5 % for all temperatures. Recently, Wang et al. [36] estimated the specific heat of the same model on a 64 × 64 lattice by the transition matrix Monte Carlo re-weighting method [39], and for a simulation with 2.5×10⁷ MC sweeps, the maximum error in temperature region T Î [0,8] was about 1%. When we apply this algorithm to the same model on the 64 × 64 lattice, with a final modification factor of 1.000000001 and a total of 2×10⁷ MC sweeps on single processor, the errors in the specific heat are reduced below 0.7% for all temperature. The relatively large errors at low T reflect the small values for the specific heat.

The behavior of the free energy near T_c cannot be extracted from ''standard'' Monte Carlo simulations but can be straightforwardly found from g(E) using Eqn.(5). In Fig. 3(a), the simulational data and exact results are presented in the same figure. To within the resolution of the figure, no difference is visible, and the relative errors are smaller than 0.001% for all temperatures.

3.2 Application to the q = 10 Potts model: A first order phase transition

In this section, we show how the algorithm is well suited to the study of a model with a strong first-order phase transition [40,41]. In such cases, the internal energy, entropy, etc. have discontinuities at the transition and both ordered and disordered states coexist at the transition. We choose the two-dimensional, ten state Potts model since it has a strong temperature-driven first-order phase transitions and a great deal is know about its behavior. We consider the Potts model on L × L square lattices with nearest-neighbor interactions and periodic boundary conditions. The Hamiltonian for this model can be written as:

and q = 1,2,...q. During the simulation, we select lattice sites randomly and choose integers between [1 : q] randomly for new Potts spin values. The modification factor f_i changes from f₀ = e¹ = 2.71828 at the beginning to f_final = exp(10^–8) @ 1.00000001 by the end of the random walk. At the end of the simulations, g(E) is normalized using the fact that the total number of possible states is q^N or that the number of ground states is q, where N = L² is the total number of lattice sites. (Actually one of these two conditions can be used to get the absolute density of states and the other condition to check the accuracy of the result.) To guarantee the accuracy of thermodynamic quantities at low temperatures, we use the condition that the number of the ground states is q to normalize the density of states.

Conventional Monte Carlo simulation (e.g. Metropolis sampling [1, 2]) determines the canonical distribution P(E, T) by generating a random walk Markov chain at a given temperature:

From the simulational result for g(E), we can calculate the canonical distribution by the above formula at essentially any temperature without performing multiple simulations. In Fig. 4 we show the resultant double peaked canonical distribution [41], at the transition temperature T_c for the first-order transition of the q = 10 Potts model. The ''transition temperatures" T_c(L) are determined by the temperatures where the double peaks are of the same height. (Note that the peaks of the distributions are normalized to 1 in this figure.) The valley between two peaks is quite deep. e.g. is 10^–9 for L = 100. The latent heat for this temperature driven first-order phase transition can be estimated from the energy difference between the double peaks. Our results for the locations of the peaks are consistent with the results obtained by multicanonical method [6] and multibondic cluster algorithm [7] for those lattice sizes for which these other methods are able to generate estimates. Wang-Landau sampling also produces results for substantially larger systems than have been studied by these other approaches.

Because of the double peak structure at a first-order phase transition, conventional Monte Carlo simulations are not efficient since an extremely long time is required for the system to travel from one peak to the other in energy space. With ''Wang-Landau sampling" all possible energy levels are visited with equal probability, so it overcomes the tunneling barrier between the coexisting phases in the conventional Monte Carlo simulations.

If we are only interested in a specific temperature range, e.g. near T_c, we could first perform a low precision unrestricted random walk, i.e. over all energies, to estimate the required energy range, and then carry out a very accurate random walk for the corresponding energy region. The inset of Fig. 4 only shows the histograms for the extensive random walks in the energy range between E/N = –1.90 and –0.6. If we need to know the density of states more accurately for some energies, we also can perform separate simulations, one for low energy levels, one for high energy levels, the other for middle energy which includes double peaks of the canonical distribution at T_c. This scheme not only speeds up the simulation, but also increases the probability of accessing the energy levels for which both maximum and minimum values of the distributions occur by performing the random walk in a relatively small energy range. If we perform single random walk over all possible energies, it will take a long time to generate rare spin configurations. Such rare energy levels include the ground energy level or low energy levels with only few spins with different values and high energy levels where all, or most, adjacent Potts spins have different values.

If the system is not larger than 100 × 100, the random walk in important energy regions (i.e. that which includes the two peaks of the canonical distribution at T_c) can be carried out efficiently with a single processor. However, for a larger system it is advisable to accelerate the calculation by performing random walks in different energy regions, each using a different processor. The densities of states for 150 × 150 and 200 × 200 , shown in Fig. 4, were obtained by joining together the estimates obtained from 21 independent random walks, each constrained within a different regions of energy. The histograms from the individual random walks are shown in the inset of Fig. 4 both for 150 × 150 and 200 × 200 lattices. In this case, we only require that the histogram of the random walk in the corresponding energy segment is sufficiently flat without regard to the relative flatness over the entire energy range. The probability distributions for large lattices show pronounced double peaks for the canonical distributions at temperatures T_c(L) = 0.70127 for L = 150 and T_c(L) = 0.701243 for L = 200. The exact result is T_c = 0.701232.... for the infinite system. Considering the minimum for L = 200 is as deep as 10^–9, we can understand why is impossible for conventional Monte Carlo algorithms to overcome the tunneling barrier with available computational resources.

To check the accuracy of the results we can estimate the transition temperature of the q = 10 Potts model for L = ¥ since the exact solution is known. According to finite-size scaling theory, the ''effective" transition temperature for finite systems behaves as:

where T_c(L) and T_c(¥) are the transition temperatures for finite- and infinite-size systems, respectively, L is the linear size of the system and d is dimension of the lattice.

The transition temperature varies as L^–d and the transition temperature extrapolated from our simulational data is consistent with the exact solution for the infinite system. We also compared our results with the existing numerical data such as estimates of transition temperatures and double peak locations obtained with the multicanonical simulational method by Berg and Neuhaus [5] and the Multibondic cluster algorithm by Janke and Kappler [7]. With the random walk simulational algorithm, we can calculate the density of states up to 200 × 200 within 10⁷ visits per energy level to obtain a good estimate of the transition temperature and locations of the double peaks. Using the multicanonical method and a finite scaling guess for the density of states, Berg et. al. only obtained results for lattices as large as 100 × 100 [5], and multibondic cluster algorithm data [7] were not given for systems larger than 50 × 50.

When we calculate the internal energy near the transition, we find that a very sharp ''jump" appears in the internal energy at T_c. The magnitude of this jump equals the latent heat for the (first-order) phase transition and is related to the double peak distribution of the first-order phase transition. When T is slightly different than T_c, one of the double peaks increases dramatically in magnitude and the other decreases. Since we only perform simulations on finite lattices, and use a continuous function to calculate thermodynamic quantities, all our quantities for finite-size systems will appear to be continuous if we use a very small scale. The same density of states can be used to calculate the internal energy for temperatures very close to T_c. On this scale the ''discontinuity'' at the first-order phase transition disappears and a smooth curve can be seen instead of a sharp ''jump''. When the system size goes to infinity, the discontinuity will be real.

The specific heat in the vicinity of the transition temperature T_c shows pronounced finite size behavior. We find it has a finite maximum value for a given lattice size L that, according to finite-size scaling theory for first-order transitions, should vary as:

where c(L,T) = C(L,T)/N is the specific heat per lattice site, L is the linear lattice size, d = 2 is the dimension of the lattice. T(L = ¥) = 0.70123.... is the exact solution for the Q = 10 Potts model [35]. Indeed, the simulational data for systems with L = 60, 100 and 200 can be well fitted by a single scaling function, moreover this function is completely consistent with the one obtained from lattice sizes from L = 18 to L = 50 by standard Monte Carlo [41].

The results for the free energy per lattice site are shown in Fig. 3(b) as a function of temperature. Since the transition is first-order the free energy appears to have a ''discontinuity'' in the first derivative at T_c. This is typical behavior for a first-order phase transition, and even with the fine scale used in the inset of Fig. 3(b), this property is still apparent even though the system is finite. The transition temperature T_c is determined by the point where the first derivative appears to be discontinuous. With a coarse temperature scale we can not distinguish the finite-size behavior of our model; however, we can see a very clear size dependence when we view the free energy on a very fine scale as in the inset.

The entropy is another very important thermodynamic quantity that cannot be calculated directly in conventional Monte Carlo simulations. It can be estimated by integrating over other thermodynamic quantities, such as specific heat, but the result is not always reliable since the specific heat itself is not easy to determine accurately, particularly considering the ''divergence" at the first-order transition. With an accurate density of states estimated by our method, we already know the Gibbs free energy and internal energy for the system, so the entropy can be calculated easily:

It is very clear that the entropy is very small at low temperature and at T = 0 is given by the density of states for the ground states. The entropy has a very sharp ''jump'' at T_c, just as does the internal energy, and this is easily visible in a plot of the entropy near T_c. The change of the entropy at T_c can be obtained by the latent heat divided by the transition temperature, and the latent heat can be obtained by the jump in internal energy at T_c.

With the histogram reweighting method [18], it is possible to use simulational data at specific temperatures to obtain complete thermodynamic information near, or between, those temperatures. Unfortunately it is usually quite hard to get accurate information in the region far away from the simulated temperature due to difficulties in obtaining good statistics, especially for large systems where the canonical distributions are very narrow. With ''Wang-Landau" sampling, the histogram is ''flat" for the random walk and we have essentially the same statistics for all energy levels. Since the output of the simulation is the density of states, which does not depend on the temperature at all, we can then calculate most thermodynamic quantities at any temperature without repeating the simulation. We also believe the algorithm is especially useful for obtaining thermodynamic information at low temperature or at the transition temperature for the systems where the conventional Monte Carlo algorithm is not so efficient.

4 Systems with a Complex Energy Landscape: The 3-dim EA model

It is straightforward to apply Wang-Landau sampling to systems with a rough energy landscape, i.e. with many minima and maxima in the probability distribution. Examples are spin glasses [42] in which the interactions between the magnetic moments produce frustration which in turn produces a complex energy landscape. One of the simplest spin-glass models is the Edwards-Anderson model [43] (EA model); however, because of the rough energy landscape of such disordered systems, the relaxation times of the conventional Monte Carlo simulations are very long [44, 45, 46]. Consequently, simulations can only be performed on rather small systems, and many properties of spin glasses are still unclarified [47-54]. Here, too, Wang-Landau sampling is efficient in estimating the density of states and also very aggressive in finding the ground states. From a random walk in energy space, we can estimate the ground state energy and the density of states very easily. The absolute density of states by the condition that total number of states is 2^N. The entropy at zero temperature can be calculated from either S₀ = ln(g(E₀)) or lim_{T ® 0}, where E₀ is the energy in a ground state. Estimates for s₀ = S₀/N and e₀ = E₀/N per lattice site agree with the corresponding estimates made with the multicanonical method.

If we are only interested in the quantities directly related to the energy, such as free energy, entropy, internal energy and specific heat, one dimensional random walks in energy space will suffice. However for spin glass systems, one of the most important quantities is the order parameter [43]

Here, N = L³ is the total number of the spins in the system, L is the linear size of the system, q(T,t) is the auto-correlation function, which depends on the temperature T and the evolution time t, and q(T,0) = 1. When t ® ¥, q(T,t) becomes the order parameter of the spin glass and

The value at T = 0 can differ from 1 because of the degenerate ground state. There is no temperature introduced during the random walk and it is more efficient to perform a random walk in a single system than two replicas. So the order-parameter can be defined

where {} is one of spin configurations at ground states and {S_i} is any configuration during the random walk. q as defined above is similar to the order-parameter defined by the Edwards and Anderson [43] and was used in the early numerical simulations by Binder et al.[55,56].

First a bond configuration is generated and a one-dimensional random walk in energy space is used to find a spin configuration {} for the ground states. Since the order-parameter is not directly related to the energy, a two-dimensional random walk is needed to get a good estimate of this quantity, i.e. to obtain the density of states g(E,q) with a flat histogram in E-q space. This also allows us to overcome the barriers in parameter space (or configuration space) for such a complex system. The rule for the 2D random walk is the same as for the 1D random walk in energy space.

The canonical distribution as a function of the order-parameter is given by:

In Fig. 5, we show a plot of the canonical distribution at different temperatures for one bond configuration of the L = 8 EA model. At low temperatures, there are multiple peaks, and the depth of the valleys between peaks depends upon temperature. When the temperature is high a single central peak is all that can be seen. At low temperature the relative depth of the minima is as great as 10^–30, and there are several local minima even at higher temperatures. With conventional Monte Carlo simulations, it is impossible to overcome the barriers at the lower temperatures, so the simulation will get trapped in one of the local minima. Therefore, even two decades after the model was proposed, there are disagreements about the existence of a finite phase transition between a glass phase and a disordered phase. For example, using Monte Carlo simulations on systems as large as (64² × 128), Marinari et al. [57] expressed doubt about the existence of the ''well-established" finite-temperature phase transition of the 3D Ising spin glass [44, 47]. Their simulational data could be described either by a finite-temperature transition or by an unusual T = 0 singularity. Kawashima and Young could also not rule out the possibility of T_g = 0 [48]. Thus even the existence of the finite-temperature phase transition is still controversial, and thus the nature of the spin glass state is uncertain. Although the best available simulational results [13, 52, 58] have been interpreted as a mean-field like behavior with replica-symmetry breaking(RSB) [59], Moore et al. found evidence for the droplet picture [60] of spin glasses within the Migdal-Kadanoff approximation. They argued that the failure to see droplet model behavior was because all existing Monte Carlo simulations were done at temperatures so close to the transition that system sizes larger than the correlation length were not used. It is possible to heat the system to increase the possibility of escape from local minima by simulated annealing and the more recent simulated tempering method [61] and parallel tempering method [62,63], but it is still very difficult to reach equilibrium at low temperatures. Recently Hatano and Gubernatis proposed a bivariate multicanonical Monte Carlo method for the 3D ±J spin glass model, and their result also favors the droplet picture [16, 64]. Marinari et al. argued, however, that the data were not thermalized [58]. The nature of spin glasses thus remains controversial [51].

5 Summary and Outlook

We have described an efficient algorithm that allows the calculation of the density of states directly for large systems. By modifying the estimate at each step of the random walk in energy space and carefully controlling the modification factor, we can determine g(E) very accurately and then calculate thermodynamic quantities at essentially any temperature . It is also noteworthy that Wang-Landau sampling allows the direct calculation of the Gibbs free energy and entropy, quantities that can only be obtained with thermodynamic integration using conventional Monte Carlo simulations. The method is applicable to a wide range of systems and is easy to implement. Although we have described its implementation in terms of single spin-flip sampling, it is straightforward to use other types of sampling for those cases in which it will further accelerate the sampling.

The application to the Ising model and the q = 10 Potts model on square lattices shows that the method is effective for systems that show first order or second-order phase transitions. Although the presentation concentrated on the random walk in energy space the idea is very general and can be applied to any parameters [11]. The energy levels of the models treated here are discrete and the total number of possible energy levels is known before simulation, but in a general model such information is not available. Since the histogram of the random walk with our algorithm tends to be flat, it is very easy to probe all possible energies and monitor the histogram entry at each energy level. For some models where all possible energy levels can not be fitted in the computer memory or the energy is continuous, e.g. the Heisenberg model, a discretization of the energy can be implemented.

In this paper, we only applied the Wang-Landau algorithm to simple models; but since the algorithm is very efficient even for large systems it should be very useful in the studies of general, complex systems with rough landscapes. There have already been a significant number of other applications of Wang-Landau sampling including efforts to improve sampling [65], to simulate proteins [66] and polymer films [67], studies of fluids [68], to study antiferromagnetic Potts models [69], reaction coordinates [70], for quantum Monte Carlo [71], for the Kondo problem [72], for the study of biological circuits [73], and for combinatorial number theory [74]. Some understanding of the convergence and performance limitations of the method have already been provided [75]; however, more investigation is needed to better determine under which circumstances the method offers substantial advantage over other approaches and how the method can be further improved.

Acknowledgments

This research has been supported by the National Science Foundation under Grant No. DMR-0094422.

[74] V. Mustonen and R. Rajesh (preprint)

Received on 5 September, 2003

• [1] D. P. Landau and K. Binder, A Guide to Monte Carlo Methods in Statistical Physics, (Cambridge U. Press, Cambridge, 2000).
• [2] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. M. Teller and E. Teller, J. Chem. Phys. 21, 1087, (1953).
• [3] R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 58, 86 (1987).
• [4] U. Wolff, Phys. Rev. Lett. 62, 361 (1989).
• [5] B. A. Berg and T. Neuhaus, Phys. Rev. Lett. 68, 9 (1992).
• [6] B. A. Berg and T. Neuhaus, Phys. Lett. B 267, 249 (1991).
• [7] W. Janke and S. Kappler, Phys. Rev. Lett. 74, 212 (1995).
• [8] B. A. Berg and T. Celik, Phys. Rev. Lett. 69, 2292 (1992).
• [9] N. S. Alves and U. Hansmann, Phys. Rev. Lett. 84, 1836 (2000).
• [10] W. Janke, Int. J. Mod. Phys. C 3, 375 (1992).
• [11] B. A. Berg and U. Hansmann and T. Neuhaus, Phys. Rev. B 47, 497 (1993).
• [12] W. Janke, Physica A 254, 164 (1998).
• [13] B. A. Berg and W. Janke, Phys. Rev. Lett. 80, 4771 (1998).
• [14] B. A. Berg, Nucl. Phys. B, 63, 982 (1998).
• [15] B. A. Berg T. Celik and U. Hansmann, Eur. Phys. Lett. 22, 63 (1993).
• [16] N. Hatano and J. E. Gubernatis in Computer Simulation Studies in Condensed Matter Physics XII, D. P. Landau, S. P. Lewis and H.-B. Schuttler (eds) (Springer, Berlin, Heidelberg, 2000).
• [17] B. A. Berg J. Stat. Phys. 82, 323 (1996).
• [18] A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 61, 2635 (1988), 63, 1195 (1989).
• [19] P. D. Beale, Phys. Rev. Lett. 76, 78 (1996).
• [20] J. Lee, Phys. Rev. Lett. 71, 211 (1993).
• [21] P. M. C. de Oliveira, T. J. P. Penna and H. J. Herrmann, Braz. J. Phys. 26, 677 (1996).
• [22] P. M. C. de Oliveira, T. J. P. Penna and H. J. Herrmann, Eur. Phys. J. B. 1, 205 (1998).
• [23] P. M. C. de Oliveira, Eur. Phys. J. B. 6, 111 (1998).
• [24] A. R. Lima, P. M. C. de Oliveira and T. J. P. Penna, J. Stat. Phys. 99, 691, (2000).
• [25] Jian-Sheng Wang and Lik Wee Lee, Comput. Phys. Commun. 127, 131, (2000).
• [26] Jian-Sheng Wang, Comput. Phys. Commun. 121, 22 (1999).
• [27] B. A. Berg and U. Hansmann, Eur. Phys. J. B 6, 395 (1998).
• [28] Jian-Sheng Wang, Eur. Phys. J. B. 8, 287 (1998).
• [29] F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001).
• [30] F. Wang and D. P. Landau, Phys. Rev. E 64, 056101 (2001).
• [31] D. P. Landau and F. Wang, Comput. Phys. Commun. 147, 674 (2002).
• [32] U. Hansmann, Phys. Rev. B 56, 6200 (1997).
• [33] U. Hansmann and Y. Okamoto, Phys. Rev. E 54, 5863 (1996).
• [34] W. Janke, B. A. Berg and A. Billoire, Comp. Phys. Commun. 121-122, 176 (1999).
• [35] F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982).
• [36] Jian-Sheng Wang, T. K. Tay and R. H. Swendsen, Phys. Rev. Lett. 82, 476 (1999).
• [37] D. P. Landau, Phys. Rev. B13, 2997 (1976).
• [38] A. E. Ferdinand and M. E. Fisher, Phys. Rev. 185, 832 (1969).
• [39] R. H. Swendsen, J. S. Wang, S. T Li, B. Diggs, C. Genovese and J. B. Kadane, Int. J. Mod. Phys. C, 10, 1563 (1999).
• [40] K. Binder, K. Vollmayr, H. P. Deutsch, J. D. Reger M. Scheucher and D. P. Landau, Int. J. of Mod. Phys. C5, 1025, (1992).
• [41] M. S. S. Challa, D. P. Landau and K. Binder Phys. Rev. B 34, 1841, (1986).
• [42] K. Binder and A.P. Young, Rev. Mod. Phys. 58, 801 (1986).
• [43] S.F. Edwards and P.W. Anderson, J Phys. F. Metal Phys. 5, 965 (1975).
• [44] A.T. Ogielski and I. Morgenstern, Phys. Rev. Lett. 54, 928 (1985).
• [45] F. Wang, N. Kawashima and M. Suzuki, Europhys. Lett. 33, 165 (1996).
• [46] F. Wang, N. Kawashima and M. Suzuki, Int. J. Mod. Phys. C. 7, 573 (1996).
• [47] R.N. Bhatt and A.P. Young, Phys. Rev. Lett. 54, 924 (1985).
• [48] N. Kawashima and A.P. Young, Phys. Rev. B53, R484 (1996).
• [49] M. Palassini and A. P. Young, Phys. Rev. Lett. 82, 5128 (1999).
• [50] M. Palassini and A. P. Young, Phys. Rev. Lett. 83, 5126 (1999).
• [51] M. Palassini and A. P. Young, Phys. Rev. Lett. 85, 3017 (2000).
• [52] E. Marinari, Phys. Rev. Lett. 82, 434 (1999).
• [53] M. A. Moore, H. Bokil and B. Drossel, Phys. Rev. Lett. 81, 4252 (1999).
• [54] J. Houdayer and O. C. Martin, Phys. Rev. Lett. 82, 4934 (1999).
• [55] I. Morgenstern and K. Binder, Phys. Rev. B22, 288 (1980).
• [56] K. Binder in Fundamental problems in statistical mechanics V, edited by E. G. D. Cohen, (North-Holland, 1980).
• [57] E. Marinari, G. Parisi and F. Ritort, J. Phys. A27, 2687 (1994).
• [58] E. Marinari, G. Parisi, F. Ricci-Tersenghi and F. Zuliani, J. Phys. A 34, 383 (2001).
• [59] G. Parisi, Phys. Rev. Lett. 43, 1754 (1979), 50, 1946 (1983).
• [60] D . S. Fisher and D. A. Huse, Phys. Rev. B38, 386 (1988).
• [61] E. Marinari and G. Parisi, Europhys. Lett. 19, 451 (1992).
• [62] K. Hukushima and K. Nemoto, J. Phys. Soc. Japan 65, 1604 (1996).
• [63] K. Hukushima, in Computer Simulation Studies in Condensed Matter Physics XIII, D. P. Landau, S. P. Lewis and H.-B. Schuttler (eds) (Springer, Blin, Heidelberg, 2000).
• [64] N. Hatano and J. E. Gubernatis, Phys. Rev. B 65, 140513 (2002).
• [65] B. J. Schulz, K. Binder, and M. Mueller, Int. J. Mod. Phys. C 13, 477 (2002); C. Yamaguchi and N. Kawashima, Phys. Rev. 65, 0556710 (2002); B. J. Schulz, K. Binder, M. Mueller, and D. P. Landau, Phys. Rev. E 67, 067102 (2003).
• [66] 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).
• [67] T. S. Jain and J. J. de Pablo, J. Chem. Phys. 116, 7238 (2002).
• [68] Q. Yan, R. Faller, and J. J. de Pablo, J. Chem Phys. 116, 8745 (2002);
• [69] C. Yamaguchi and Y. Okabe, J. Phys. A 34, 8781 (2001).
• [70] F. Calvo, Mol. Phys. 100, 3421 (2002).
• [71] M. Troyer, S. Wessel, and F. Alet, Phys. Rev. Lett. 90, 120201 (2003)
• [72] W. Koller, A. Prull, and H. G. Evertz, Phys. Rev. B 67, 104432 (2003)
• [73] H.-B. Schttler, et al. (private communication)
• [75] P. Dayal, S. Trebst, S. Wessel, D. Würtz, M. Troyer, S. Sabhapandit, and S. N. Coppersmith, cond-mat/0306108; C. Zhou and R. H. Bhatt, cond-matt/0306711

Publication Dates

• Publication in this collection
  01 Sept 2004
• Date of issue
  June 2004

History

• Received
  05 July 2003
• Accepted
  05 July 2003