Broad histogram: tests for a simple and efficient microcanonical simulator

Abstract

The Broad Histogram Method (BHM) allows one to determine the energy degeneracy g(E), i.e. the energy spectrum of a given system, from the knowledge of the microcanonical averages < Nup (E) > and < Ndn (E) > of two macroscopic quantities Nup and Ndn defined within the method. The fundamental BHM equation relating g(E) to the quoted averages is exact and completely general for any conceivable system. Thus, the only possible source of numerical inaccuracies resides on the measurement of the averages themselves. In this text, we introduce a Monte Carlo recipe to measure microcanonical averages. In order to test its performance, we applied it to the Ising ferromagnet on a 32 × 32 square lattice. The exact values of g(E) are known up to this lattice size, thus it is a good standard to compare our numerical results with. Measuring the deviations relative to the exactly known values, we verified a decay proportional to 1/<img src="http:/img/fbpe/bjp/v30n4/24eq01.gif" border="0">, by increasing the counter (counts) of averaged samples over at least 6 decades. That is why we believe this microcanonical simulator presents no bias besides the normal statistical fluctuations. For counts ~ 10(10), we measured relative deviations near 10(5) for both g(E) and the specific heat peak, obtained through BHM relation.

Broad histogram: tests for a simple and efficient microcanonical simulator

Paulo Murilo Castro de Oliveira

Instituto de Física, Universidade Federal Fluminense

av. Litorânea s/n, Boa Viagem, Niterói RJ, Brazil 24210-340

e-mail [e-mail]

Received on 5 August, 2000

The Broad Histogram Method (BHM) allows one to determine the energy degeneracy g(E), i.e. the energy spectrum of a given system, from the knowledge of the microcanonical averages < N^up (E) > and < N^dn (E) > of two macroscopic quantities N^up and N^dn defined within the method. The fundamental BHM equation relating g(E) to the quoted averages is exact and completely general for any conceivable system. Thus, the only possible source of numerical inaccuracies resides on the measurement of the averages themselves. In this text, we introduce a Monte Carlo recipe to measure microcanonical averages. In order to test its performance, we applied it to the Ising ferromagnet on a 32 × 32 square lattice. The exact values of g(E) are known up to this lattice size, thus it is a good standard to compare our numerical results with. Measuring the deviations relative to the exactly known values, we verified a decay proportional to 1/, by increasing the counter (counts) of averaged samples over at least 6 decades. That is why we believe this microcanonical simulator presents no bias besides the normal statistical fluctuations. For counts » 10¹⁰, we measured relative deviations near 10⁵ for both g(E) and the specific heat peak, obtained through BHM relation.

Monte Carlo methods are applied to statistical physics in order to measure the thermal average

of some macroscopic quantity Q (magnetisation, density, etc). The temperature T is fixed, and the Boltzmann constant is set to unity. Both sums run over all possible microstates S available for the system. The energy (quantity Q) corresponding to S is denoted by E_S (Q_S). Within traditional computer simulations, instead of taking all possible states, one takes only a finite set of them, i.e. a Markovian chain of states randomly tossed according to probabilities dictated by the Boltzmann exponential factors exp(-E_S/T), the so-called importance sampling. Thus, one needs to fix a particular value for the temperature T, before running the computer job. In order to determine the full dependence of < Q > _T upon T, one needs to run the job again and again, for different values of T.

An alternative is to re-write the same average as

where the degeneracy g(E) counts the number of states with energy E, and both sums run over all possible energies. The microcanonical average

of the same quantity Q corresponds to a fixed energy, i.e. the sum in (3) runs only over the states S[E] belonging to energy level E. This microcanonical average (3) is simpler than the canonical counterpart (1) or (2), prescribing exactly the same weight to all averaging states, i.e. it is a uniform averaging process within each energy level separately.

Only the Boltzmann factor exp(-E/T) appearing in equation (2) depends on the temperatute, carrying all thermodynamic information about the environment which continuously exchanges energy with the system under study. Contrary to this, both g(E) and < Q(E) > are independent of the particular way this energy exchange occurs, independent of the environment: they are more fundamental properties of the system alone, defined only by its energy spectrum. They are not thermodynamic quantities, and do not depend upon thermodynamic concerns like temperature, equilibrium, etc. In practical terms, equation (2) allows one to determine < Q > _T for any value of T, from the knowledge of the energy functions g(E) and < Q(E) >.

The Broad Histogram Method (BHM) [1] relates g(E) with the microcanonical averages of two macroscopic quantities N^up and N^dn, measured at the current state S. First, one needs to adopt some protocol of allowed movements which could be performed on S, leading to another possible state S¢. One can adopt any such a protocol, the only restriction being its reversibility (if S ® S¢ is allowed, so is S¢® S). Considering the Ising model, for instance, the protocol could be chosen to be the whole set of single-spin flips (among many other alternative choices). Given such a protocol, N^up(S) counts the number of allowed movements which could be performed on S, leading to a fixed energy increment DE (which must be chosen a priori, although the method is also independent of this choice). Analogously, N^dn(S) counts the number of allowed movements decreasing the energy of S by the same fixed amount DE. The fundamental BHM relation [1] is

where the microcanonical averages of N^up and N^dn are defined by equation (3). By knowing these energy functions, equation (4) allows one to determine g(E) along all the energy axis, in steps of DE - a constant, unimportant pre-factor is cancelled out by performing the average in equation (2). Relation (4) is shown to be valid for any energy spectrum, under completely general grounds [2]. Also, the same reasoning could be applied for another basic quantity q, instead of the energy, by considering the degeneracies g(q) instead of g(E). For simplicity, we will restrict ourselves to the case of the energy.

Thus, the Broad Histogram Method consists in: i) to choose some protocol of allowed movements, as well as an energy jump DE; ii) to measure, by any means, the microcanonical average < Q(E) > as a function of the energy, as well as < N^up(E) > and < N^dn(E) > which determine g(E) through equation (4); iii) to obtain the desired thermal average through equation (2). There is no approximation at all, and the final numerical accuracy depends exclusively upon the microcanonical measuring strategy adopted in step ii. The method is reviewed in [3]. Some references where it is used are [4-16].

Let's briefly analyse hereafter some possible computer strategies one can adopt in order to measure the microcanonical averages < Q(E) > , < N^up(E) > and < N^dn(E) > , as functions of the energy. A Markovian chain of states is obtained by performing random movements transforming the current state into another. These movements are tossed among some previously defined set of possibilities, another protocol which has nothing to do with the BHM protocol of virtual movements. Both protocols could even be chosen to be the same, but not necessarily.

The first direct strategy is to keep always the same fixed value E: by starting from some state corresponding to the desired energy, one simply rejects any tossed movement which changes the current energy. After one has already a large enough number of visited states inside this particular energy level, the same process is repeated for other levels. Depending on the adopted protocol of (real) movements, this strategy could lead to ergodicity problems: due to its high rejection rate, one risks to sample only a biased sub-set of states belonging to energy level E, violating the required uniform visitation. Moreover, it is also an inefficient strategy, in what concerns the computer time, again due to its high rejection rate.

The opposed alternative is to accept also movements leading to energy jumps, storing separated averages for each energy level, in parallel. Obviously, this option is much more efficient in what concerns the computer time. However, one cannot simply accept any tossed movement: energy increments would occur more often than decrements, because g(E) is normally a fast increasing function of the energy. As a result, at the end, only states corresponding to the region near the maximum of g(E) would be sampled. Thus, some movement-rejection prescription must be adopted, and we get back to the uniformity violation problem. One can adopt some already well established movement-rejection prescription based on detailed balance arguments. For instance, canonical, fixed-temperature dynamics could be adopted [4, 6], sampling states inside the narrow energy window corresponding to the fixed value of T. In order to get results on a broader energy window, one can simply superimpose the histograms obtained for different computer runs corresponding to different values of T.

Another possibility is to adopt one of the many multicanonical dynamics [7, 9, 11, 13, 14]. Within this approach, one tunes the E-dependent rejection rate during the computer run, in order to obtain the same visit probability for all energy levels, i.e. a flat distribution along the energy axis, at the end. In this case, the visitation probability to each particular state would be proportional to 1/g(E). The so-called multicanonical methods [17-19] are based just on this feature: by recording the acceptance probability for energy-increasing movements E ® E¢, accumulated during the run and which must be equal to g(E)/g(E¢) at the end, one gets the function g(E) except for an unimportant global factor which cancels in equation (2). BHM is completely distinct from multicanonical approaches in many features. In particular, the infomations extracted from each E-state S belonging to the averaging Markovian chain are the macroscopic values of N^up(S) and N^dn(S), not the mere one-more-visit upgrade V(E) ®V(E) + 1. That is why BHM gives more accurate results than multicanonical approaches, even taking into account the same Markovian chain of averaging states, as shown in [15] where the multicanonical dynamics [18] is adopted in order to measure < N^up(E) > and < N^dn(E) > . At the end, g(E) is determined twice, by following the multicanonical traditional way or, alternatively, by the BHM equation (4), both using data taken from the same set of averaging states. A clear accuracy advantage for BHM is reported [15]. Moreover, due to the macroscopic character of the BHM quantities N^up and N^dn this advantage still increases for larger and larger systems.

Another crucial difference relative to multicanonical approaches is that BHM requires only the uniformity of visits among the states inside each energy level, separately, in order to get the correct microcanonical averages. BHM does not need any detailed balance between visits to different levels E and E¢, nor the multicanonical flat distribution along the energy axis. Any dynamic strategy which is good for multicanonical methods will be also good for BHM (besides the accuracy advantage quoted in the last paragraph), but the reverse is not true. Thus, within BHM, other not-so-restricted dynamic strategies could be used.

Profiting from this feature, we decided to test a very simple dynamic strategy inspired by reference [20] (although within a different method, the dynamic rule introduced in [20] is essentially the same as presented hereafter). The idea is to avoid movement rejections at all, within the energy level currently being sampled for averaging purposes. Rejections will be restricted to other energies, whose states are never included into the averaging statistics. Let's consider that the maximum energy jump allowed by the adopted protocol of (real) movements corresponds to n levels above or below the current energy E. Then, let's take an energy window of 2n + 1 adjacent levels, starting from some state inside it. Any movement which keeps the system still inside this window will be accepted. This is the dynamic strategy we propose here. Averages are taken only for the central level E inside the chosen energy window: the system is allowed to visit the other n levels above it, as well as the n levels below it, nevertheless without measuring anything during these side visits. Note that no tossed movement will be rejected, if the current energy is just E. Thus, the averaging process is completely rejection-free, avoiding any systematic bias due to artificial rejections rules. The same strategy can be easily applied also for continuous energy spectra, by taking averages only within a narrow, rejection-free energy window centered inside another broader, free-visit window: any movement to outside this latter would be rejected.

In what follows, we consider a L × L square lattice Ising ferromagnet, with L = 32 for which the exact function g(E) is known [21]. The energy of the current state S is counted as the total number of its unsatisfied bonds, i.e. the total number of neighbouring pairs of spins pointing one up and the other down. The energy spectrum corresponds to all even numbers between 0 and 2,048, i.e. E = 0, 2, 4, 6, 8 ... 2,048, which can also be represented by energy densities e = E/(2L²). The degeneracies are g(E) = 2, 0, 2,048, 4,096, 1,057,792 ... 2, respectively. This spectrum is symmetric in relation to its center at E_max = 1,024 (or e = 0.5), where g(E_max) » 6.3 × 10³⁰⁶. We need only its first half, corresponding to positive temperatures. The critical energy (at the thermodynamic limit) corresponds to e_c» 0.147 (E_c» 300, for L = 32).

We will adopt single-spin flips as the protocol of movements (both for the real movements tossed during the computer run, as well as for the virtual ones we consider in order to count N^up and N^dn). Starting from the current state S with energy E, 1,024 movements would be available, by tossing one random spin to flip. They can be classified according to the possible energy jumps, E ® E ± DE, where DE could be 0, 2 or 4. Thus, our energy window will have 5 adjacent levels, the central one, E, where the averages will be measured, which is rejection-free, plus two neighbouring levels above it, and two others below it. During the random walk performed inside this window, we measure the values of N^up and N^dn for the current state, every time its energy is E, and accumulate them into two E-histograms H^up(E) and H^dn(E). Also a visit counter V(E) is upgraded to V(E)+1. If the energy of the current state is not E, nothing is measured or accumulated. At the end, we take the averages < N^up(E) > = H^up(E)/V(E) and < N^dn(E) > = H^dn(E)/V(E). Note that this is the only role played by the final V(E): no comparison with the neighbouring values V(E ± DE) is needed, no further information is extracted from these values. They must only be large enough in order to provide a good statistics. For each new E, a new 5-levels energy window centered on it is sampled, and the whole process is repeated.

By following this dynamic rule and by using the BHM relation (4), we measured the quantity ln[g(E+DE)/g(E)] for 15 adjacent levels E = 300, 302 ... 328, at the critical region. The deviations from the exact values were averaged (root mean square) over these 15 levels, and are shown as a function of the number of averaging states sampled inside each energy level (counts), in Fig. 1. The squares corresponds to DE = 4, while the diamonds represent DE = 2. The dashed straight line (1/) indicates that no systematic errors besides the normal statistical fluctuations occur, giving credit to our simple dynamic rule. According to these results, to improve the numerical accuracy is a simple matter of taking more and more averaging states inside each energy level, up to the computer time available.

Figure 1. Deviations between measured g(E) and the exact values, as a function of the number (counts) of sampled averaging states. The dashed line corresponds to 1/. Data for 32 × 32 Ising ferromagnet.

In order to perform thermal averages, one does not need the same accuracy along the whole energy axis. The function [g(E)exp(-E/Tc)]² is displayed by the dotted line in Fig. 2, where T_c = 2.293930 is the exact location of the specific heat peak. This curve displays the (squared) relative contribution of each energy to the partition function. As the number of sampled averaging states inside each energy level is proportional to the squared numerical accuracy, the dotted line in figure 2 shows the ideal profile of visits one needs in order to have equally accurate contributions from each energy level. It is a sharp peaked curve, according to which the computational effort can be concentrated only inside a narrow energy window. Profiting from this feature, we have shaped the profile of visits displayed by the solid line, in figure 2. The possibility of designing this profile of visits, sampling different numbers of averaging states for different energies, according to the relative contribution of each energy region, is a further big advantage of BHM over multicanonical methods. Almost all the computational effort is concentrated near the peak.

Fig. 3 shows a detailed comparison of our simulational results with the exact specific heat curve, near its peak. Being a derivative, which corresponds to a mathematical procedure in which numerical accuracy is strongly compromised, this quantity is a good standard for worst-case comparisons, moreover near its peak. Nevertheless, the relative deviations we obtained are compatible with the number (counts = 8 × 10⁹) of averaged states per energy level, at the peak. Moreover, better yet accuracies were obtained along all the temperature range from 0 to ¥, always by using the same simulational data, i.e. the same averaged values for the BHM quantities N^up and N^dn measured during a single computer run.

In short, we have tested a simple dynamics which is very efficient in measuring microcanonical averages. It is essentially the same dynamics as introduced in [20], for other purposes. Here, the aim is to measure the microcanonical averages of some particular macroscopic quantities defined within the broad histogram method [1-3]. Once one knows these averaged quantities as functions of the energy, the method provides the energy degeneracy function g(E) through an exact relation. During the same computer run, the microcanonical average < Q(E) > of the quantity Q of interest is also measured. Then, once one knows g(E) and < Q(E) > , the canonical thermal average < Q > _T can be determined by equation (2), for any, continuously varying temperature T, without resorting again to further computer simulations. According to our tests, the current dynamic rule does not introduce any systematic averaging bias, besides the normal statistical fluctuations which decay proportionally to 1/. Thus, by applying this dynamic rule to BHM, to improve more and more the numerical accuracy is a simple matter of increasing the computer time. Among the further advantages of the microcanonical dynamics tested here, we can quote: i) its implementation simplicity, without detailed balance and other complications; ii) no movement-rejections at all, within the averaging energy level; iii) the possibility to shape the profile of visits along the energy axis, according to the desired accuracy; iv) no randomness at all is used in order to decide to perform or not the currently tossed movement [22]; v) short and non-periodic waiting time between consecutive averaging states [24].

I am indebted to Kim Doochul who warned me about reference [20], reading which I have had the idea to apply its dynamic rules to the broad histogram method. D. Stauffer was kind enough to perform a critical reading of the manuscript. This work is partially supported by Brazilian agencies CAPES, CNPq and FAPERJ.

References

[1] P.M.C. de Oliveira, T.J.P. Penna and H.J. Herrmann, Braz. J. Phys. 26, 677 (1996) (also in Cond-Mat 9610041).

[2] P.M.C. de Oliveira, Eur. Phys. J. B6, 111 (1998) (also in Cond-Mat 9807354).

[3] P.M.C. de Oliveira, Braz. J. Phys. 30, 195 (2000) (also in Cond-Mat 0003300).

[4] P.M.C. de Oliveira, T.J.P. Penna and H.J. Herrmann, Eur. Phys. J. B1, 205 (1998); P.M.C. de Oliveira, in Computer Simulation Studies in Condensed Matter Physics XI, 169, eds. D.P. Landau and H.-B. Schüttler, Springer, Heidelberg/Berlin (1998).

[5] P.M.C. de Oliveira, Int. J. Mod. Phys. C9, 497 (1998).

[6] J.-S. Wang, T.K. Tay and R.H. Swendsen, Phys. Rev. Lett. 82, 476 (1999).

[7] J.-S. Wang, Eur. Phys. J. B8, 287 (1999) (also in Cond-Mat 9810017).

[8] J.D. Muñoz and H.J. Herrmann, Int. J. Mod. Phys. C10, 95 (1999); Comput. Phys. Comm. 121-122, 13 (1999).

[9] R.H. Swendsen, B. Diggs, J.-S. Wang, S.-T. Li, C. Genovese and J.B. Kadane, Int. J. Mod. Phys. C10, 1563 (1999) (also in Cond-Mat 9908461).

[10] P.M.C. de Oliveira, Comput. Phys. Comm. 121-122, 16 (1999), presented at the APS/EPS Conference on Computational Physics, Granada, Spain (1998).

[11] J.-S. Wang, Comput. Phys. Comm. 121-122, 22 (1999).

[12] A.R. de Lima, P.M.C. de Oliveira and T.J.P. Penna, Solid State Comm. 114, 447 (2000) (also in Cond-Mat 9912152).

[13] J.-S. Wang and L.W. Lee, Comput. Phys. Comm. 127, 131 (2000) (also in Cond-Mat 9903224).

[14] J.-S. Wang, Prog. Theor. Phys. Supp. 138, 454 (2000) (also in Cond-Mat 9909177).

[15] A.R. de Lima, P.M.C. de Oliveira and T.J.P. Penna, J. Stat. Phys. 99, 691 (2000) (also in Cond-Mat 0002176).

[16] M. Kastner, J.D. Muñoz and M. Promberger, Phys. Rev. E, to appear (also in Cond-Mat 9906097).

[17] B.A. Berg and T. Neuhaus, Phys. Lett. B267, 249 (1991); A.P. Lyubartsev, A.A. Martsinovski, S.V. Shevkunov and P.N. Vorontsov-Velyaminov, J. Chem. Phys. 96, 1776 (1992); E. Marinari and G. Parisi, Europhys. Lett. 19, 451 (1992); B.A. Berg, Int. J. Mod. Phys. C4, 249 (1993).

[18] J. Lee, Phys. Rev. Lett. 71, 211 (1993).

[19] B. Hesselbo and R.B. Stinchcombe, Phys. Rev. Lett. 74, 2151 (1995).

[20] K.-C. Lee, J. Phys. A23, 2087 (1990).

[21] P.D. Beale, Phys. Rev. Lett. 76, 78 (1996).

[22] Within single-spin flips, for instance, random numbers are used only to toss which spin would be flipped, a coarse grained decision concerning only integer numbers. One does not need to compare them with precisely defined acceptance ratios, a delicate comparison between real numbers which could introduce undesired biases as in

[23] A.M. Ferrenberg, D.P. Landau and Y.J. Wong, Phys. Rev. Lett. 69, 3382 (1992).

[24] In the present case of a 32 × 32 lattice, a new state was sampled after every ~ 40 spin flips, in average, instead of L² = 1,024.

Publication Dates

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

History

• Received
  05 Aug 2000