Energy degeneracies from microcanonical averages

Computer Physics Communications 121–122 (1999) 16–21 www.elsevier.nl/locate/cpc

Energy degeneracies from microcanonical averages Paulo Murilo Castro de Oliveira 1,2 Laboratoire de Physique et Mécanique des Milieux Hétérogènes, École Supérieure de Physique et de Chimie Industrielles, 10, rue Vauquelin, 75231 Paris Cedex 05, France

Abstract A method for the direct calculation of the energy degeneracy g(E) has been recently introduced. One considers all possible potential modifications one can perform, starting from the current state of the system (for instance all possible spin flips), and counts the number Nup (Ndn ) of such potential movements increasing (decreasing) the energy by a fixed value 1E. The degeneracy g(E) can thus be determined from the microcanonical averages hNup (E)i and hNdn (E)i of these macroscopic quantities, measured as functions of the energy E. Any microcanonical simulator can be used in order to get these averages. Instead of simply measuring the histogram of the number of visits where each new state increases by unity the counter in the corresponding energy channel, here the counters are incremented by macroscopic quantities for each new visited state. Thus, we reach very good accuracy with little computational effort.  1999 Elsevier Science B.V. All rights reserved.

The canonical equilibrium average of some thermodynamic quantity Q reads as P S QS exp(−ES /T ) , (1) hQiT = P S exp(−ES /T ) where the physical system is supposed to be in contact with a heat bath at fixed temperature T , and the Boltzmann constant is set to unity. Both sums run over all possible microstates available for the system, and ES (QS ) is the value of its energy (quantity Q) at the particular microstate S. Most computer simulational methods used in statistical physics deal with the measurement of such an average, where some particular value of T is fixed at the beginning. Instead of taking into account all possible states S, only a set of randomly chosen ones are sampled according to 1 E-mail: [email protected]. 2 Permanent address: Instituto de Física, Universidade Federal Fluminense, av. Litorânea s/n, Boa Viagem, Niterói RJ 24210-340 Brazil.

the energy-dependent probabilities proportional to the Boltzmann exponential factors in (1). After running such a sampling routine on a computer, one has an approximation for hQiT [1]. In order to control the quality of this approximation, many precautions are needed (see [2], or many other books edited by Kurt Binder). In order to obtain hQiT as a function of T , one needs to run the computer program again and again, for different temperatures. In most applications, one needs a continuous curve, and this goal can be achieved by reweighting methods [3] which allow one to extent to an interval T0 − 1T /2 < T < T0 + 1T /2 the measurements actually performed at the fixed value T0 . However, the larger the system, the smaller is the width 1T compatible with good accuracy, and anyway one must perform different runs for different fixed values of T0 . By classifying the summation terms in Eq. (1) according to their energies, one can rewrite it as P hQ(E)ig(E) exp(−E/T ) , (2) hQiT = EP E g(E) exp(−E/T )

0010-4655/99/$ – see front matter  1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 1 0 - 4 6 5 5 ( 9 9 ) 0 0 2 6 9 - 6

P.M.C. de Oliveira / Computer Physics Communications 121–122 (1999) 16–21

where now the sums run over all allowed energy levels E, each one corresponding to a total of g(E) degenerate microstates. The microcanonical average P

 S(E) QS (3) Q(E) = g(E) is now restricted only to the microstates belonging to the energy level E, uniformly weighted. Although both determine the same equilibrium canonical average hQiT , Eqs. (1) and (2) have entirely distinct conceptual interpretations. The degeneracy g(E) has nothing to do with thermodynamics: it is a more fundamental property of the physical system itself, independent of its possible interactions with the environment, independent whether the system is in equilibrium or not, etc. It depends only on the distribution of energies among the system’s possible microstates. It is the signature of the physical system, defining its behaviour under any circumstance, whatever is its current environment. Indeed, hQ(E)i is also a property of the system itself as well. Thus, in Eq. (2), all the interactions of the system with its environment (in this particular case, a fixed-temperature heat bath) are included only into the Boltzmann exponential factors. According to this thermodynamics-free interpretation, a simulational method able to measure directly g(E) and hQ(E)i would be more efficient than to repeat the canonical simulations again and again for different temperatures. Once one has obtained g(E) and hQ(E)i, the effects of the temperature T can be included a posteriori, the computer simulation already being over. Such a general simulational method was introduced two years ago [4]. It is based on two quantities Nup and Ndn measuring the number of potential modifications one can perform into the current system’s microstate, respectively increasing or decreasing its energy by a fixed value 1E. The degeneracy function g(E) can thus be obtained from the microcanonical averages hNup (E)i and hNdn (E)i of these two quantities, through 



(4) g(E) Nup (E) = g(E + 1E) Ndn (E + 1E) , which is the fundamental equation introduced in [4]. In order to apply the method to some particular system, one needs first to adopt a protocol of allowed movements, i.e. allowed modifications which can be performed into the current microstate. Considering an Ising system, for instance, one can adopt singlespin flips. However, the method does not depend on

17

which is the particular adopted protocol. Once such a protocol is chosen, one can measure the quantities Nup and Ndn at the current state. Note that none of these potential movements will be actually performed. What one needs is to measure the microcanonical averages hNup (E)i and hNdn (E)i, no matter the approach adopted in order to perform this task. In Monte Carlo sampling, one constructs a Markovian chain of states, going from the current state to the next by performing movements prescribed by some protocol which, however, does not need to be the same adopted to count Nup and Ndn . One can adopt another completely different protocol, any other dynamics, provided the microcanonical averages hNup (E)i and hNdn (E)i are correctly measured: any good microcanonical simulator can be used. In order to have a fast computer run, for instance, one can follow some cluster updating dynamics [5], while counting Nup and Ndn within the single-spin flip protocol. The broad histogram method [4] was tested for different simulational dynamics [6–8], for the Ising ferromagnet in two and three dimensions, as well as for the Ising spin glass in three dimensions [6]. The fundamental equation (4) was proven to be exact for any system, under completely general conditions [8]. It can also be applied to continuous-energy spectra as the XY model [9]. Combining the idea of measuring the thermodynamics-independent quantities hNup (E)i and hNdn (E)i with the matrix approach of the transition probability method [10], an alternative, elegant formulation is also available [11] allowing one to study dynamical features. Here, we return to the Ising ferromagnet in two and three dimensions, in order to test the numerical performance. All the results are obtained by adopting the single-spin flip protocol, and the following particular dynamics. Let’s assume the two-dimensional case, namely a L × L square lattice, with even L and periodic boundary conditions. The energy E = 0, 4, 6, 8, . . ., 2L2 can be counted as the number of unsatisfied bonds (E = 0 corresponding to complete spin alignment, whereas E = 2L2 is the Néel state energy). The density e = E/2L2 varies from 0 to unity, but the physical (positive-temperature) region is the interval 0 6 e 6 0.5. The critical value being ec ≈ 0.15, we restrict our simulations to the interval 0.05 6 e 6 0.3, dividing it into windows containing 4 energy levels each. We start our simulation inside the first, lowest energy window: a random new state

18

P.M.C. de Oliveira / Computer Physics Communications 121–122 (1999) 16–21

is sampled by flipping one random spin of the current state, provided the new energy does not surpasses the window’s limits. A new averaging state is taken after just L2 spin flips were effectively performed, and so on, until a previously defined number of such whole lattice sweeps were completed. Then, the energy window is shifted one level up, losing its lowest energy level, keeping the other 3, and including now a new adjacent, not-yet-sampled higher-energy level. The simulation process continues inside this new window, and so on. However, the spin flips are not always accepted (as in [7]). Each spin flip can lead to five distinct possible values of the energy jump 1E = −4, −2, 0, +2 or +4. At the current state, we count separately Ndn (−4), Ndn (−2), Nup (+2) and Nup (+4), i.e. the numbers of potential movements corresponding to each possible (non zero) energy jump. Then, we randomly choose some spin to flip: if it fits into one of the 3 first possibilities, i.e. 1E = −4, −2 or 0, the movement is accepted and performed, provided the resulting energy is still inside the current window. The 2 other possibilities are accepted according to probabilities p(+2) = Ndn (−2)/Nup(+2) and p(+4) = Ndn (−4)/Nup(+4), respectively. Accepting a movement independently of the value of its energy jump, depending only upon its sign, corresponds to the energy random walk dynamics [12] designed to find optimal solutions (energy minima) in “zero temperature” complex systems. However, if one wants to measure canonical statistical averages, different energy (jumps) must be treated within different probabilities, namely p(+4) = [p(+2)]2 for the present case (in average), and not taking the simplification p(+4) = p(+2) as in [13], which violates the Boltzmann distribution. Both approaches lead to non-biased random walks along the energy axis, but only the energy-discriminating one [4] obeys the Boltzmann prescriptions not mixing different energies within the same transition probability. Nevertheless, this technical detail has nothing to do with the method [4] designed to measure directly g(E) from Eq. (4): any other (correct) microcanonical simulator could be used. Actually, there is a trick [4] allowing one to include all positive energy jumps into only one histogram Nup , still obeying Boltzmann weights. It consists in first computing the numbers Nup (1E) for all possible positive values of 1E (namely, Nup (+2) and Nup (+4)

Fig. 1. Profile of visits along the energy axis (crosses), where the upper, central plateau corresponds to 32000 averaging states per energy level. Also shown are the magnetization (diamonds) and the dimensionless function τ [14] (squares) for the 32 × 32 square lattice Ising ferromagnet.

for the square lattice Ising model) at the current state, and then taking X  1/1E . (5) Nup (1E) Nup = 1E>0

Analogously, one takes X  1/1E . Ndn (1E) Ndn =

(6)

1E<0

The flipping probability for a positive energy jump 1E can also be replaced by p(1E) = (Ndn /Nup )1E measured at the current state, or alternatively from the values already accumulated at the histograms up to the current stage of the simulation. The present results were obtained within this approach. Fig. 1 shows the energy histograms obtained for 32 × 32 lattices. The profile of visits from e = 0.05 up to e = 0.3 was chosen with the trapezoidal format, in order to get better statistics in the critical region. The largest number of visits per energy level (at the central

P.M.C. de Oliveira / Computer Physics Communications 121–122 (1999) 16–21

plateau) is 100, for each of the 320 distinct lattices (distinct starting microstate) we have sampled. The other two curves correspond to the magnetization and the function τ defined in [14] (the steepest curve). This function has properties similar to the Binder cumulant U [15,2]: (1) it scales as L0 at the critical point; (2) it is a step function in the thermodynamic limit (τ = 1 for T < Tc and τ = 0 for T > Tc ); (3) the crossings of τL (T ) and τL0 (T ) for different lattice sizes L and L0 approach the critical temperature as L, L0 → ∞; (4) the plot of −dτ/dT as a function of T scales as L−1/ν along the horizontal axis, and as L1/ν along the vertical one. Its advantage over the Binder cumulant is that it also obeys all the correct scaling (not critical) behaviours at both T = 0 and T = ∞ (the other two scalinginvariant points, besides Tc ). Maybe because of this property, τ seems to reach the asymptotic region (L → ∞) earlier than U , i.e. the same accuracy should be achieved with smaller lattices. For the particular size 32 × 32, we also completed Fig. 1 by adding two further plateaus of averaging states to the left and to the right. These further plateaus are 100 times lower than the highest one at the central region, and complete the whole energy range 0 < e < 0.5. Fig. 2 shows the averaged energy and specific heat, for this case, compared with the exact curves (circles) available up to this lattice size [16]. Note the broad, continuous range of temperatures obtained, from a single simulation. The inset shows a blow up of the specific heat peak located at Tc (32) = 2.29363(56) with height cmax (32) = 1.9066(27), with the error bars inside parentheses corresponding to the last two digits: the exact values [16] are 2.29393 and 1.9045, respectively. If one wants to sample only the critical region, as shown in the inset, a much narrower energy range would be enough. Fig. 3 shows the data collapsing plots obtained for τ and −dτ/dT , from square lattices with 6 different sizes. The same quality is also obtained in three dimensions. Table 1 shows the critical temperatures obtained by finite size extrapolations of the peak positions of: specific heat, (log of) magnetic susceptibility, −dτ/dT and −dU/dT , for both square and cubic lattices. The simulated sizes were L = 22, 32, 44, 64, 90 and 128 in

19

Fig. 2. Averaged energy and specific heat for the 32 × 32 square lattice Ising ferromagnet. We choose to show data for this particular size, because the corresponding curves (circles) are exactly known [16]. Note the broad temperature range covered by the simulation, on one hand, and also the numerical accuracy obtained at the critical region, on the other hand (inset for the specific heat peak, from the same data). Table 1 Extrapolated values of the infinite lattice critical temperature, using data from simulated finite systems. The exact values correspond to the Onsager one for the square lattice, or to accurate, large lattice simulational results [17] for the cubic one Source

cmax log χmax −dτ/dTmax −dU/dTmax Exact

square lattice 2.2707 2.2693

2.2711

2.2672

2.2692

cubic lattice 4.5101 4.5096

4.5120

4.5062

4.5115

two, and L = 8, 10, 12, 16, 20 and 26 in three dimensions. The total computer time was around 300 CPU hours on DEC Alpha 400 workstations. For the cubic lattice, we grouped 6 energy levels per window, and the energy interval adopted was 0.1 < e < 0.5 for the smaller lattices, and 0.27 < e < 0.37 for L = 20 and 26. An important feature of this method, distinguishing it from generalized ensembles [18], is that it is not based on the profile V (E) of visits along the energy axis, during the sampling dynamics. In Fig. 1, for instance, we choose the trapezoidal profile only

20

P.M.C. de Oliveira / Computer Physics Communications 121–122 (1999) 16–21

Fig. 3. Data collapsing plots of the dimensionless quantity τ [14] and its temperature derivative, for L × L square lattices. Here, Tc = 2.269185 is the Onsager critical temperature for the infinite lattice.

because it is convenient in order to have better statistics near the critical point. However, our measures are not related to this particular distribution of visits, but to the microcanonical averages of Nup and Ndn allowing the direct determination of g(E). The macroscopic character of these two averaged quantities allows one to obtain better accuracies: each new averaging visited state contributes with a macroscopic value to the averaging histograms, instead of merely incrementing V (E) by unit. Although the method itself does not depend on it, the particular dynamics adopted in this work also depends on the macroscopic character of these quantities. Very near the ground state, for instance, the values of Ndn may be not macroscopic, and the currently used dynamics is no longer accurate: one can hardly see some deviations of the specific heat occurring for very low temperatures, in Fig. 2. Generalized ensemble methods [18] introduce controlled biases into the otherwise narrow canonical distribution Vcan (E), obtaining thus a broader distribution Vbias(E). After the simulation is over, Vcan (E) is restored by subtracting the bias from the actually measured distribution Vbias (E). During the simulation, the bias is introduced step by step, following an ad hoc

iterative procedure: from the current sampling distribution, one guesses a further bias and introduce it in order to explore new, not-yet visited regions. At each step, the guess must be recorded, in order to subtract the whole bias at the end. A better way to control the guesses within this iterative process is the transition probability method [10], where the dynamic evolution starts from a low probability microstate, and the inferences about the probability concerning each energy level are made from the record of transitions rather than the record of visited states. A probability transition matrix between distinct energy levels is constructed by counting the movements effectively performed. Better yet is to obtain the same matrix by taking into account all the potential, not implemented movements one could perform at the current microstate. This is done by measuring the microcanonical averages hNup (E)i and hNdn (E)i introduced in [4], a procedure called the transitional Monte Carlo [11]. The broad histogram method [4], however, is based only on the exact Eq. (4), independent of the particular dynamic process followed in order to measure the averages hNup (E)i and hNdn (E)i. Concluding, we have shown numerical data concerning the Ising ferromagnet as a testing model for the broad histogram method [4]. The energy spectrum degeneracy g(E) is directly obtained free from thermodynamic concerns. Thus, by using always the same data measured from a single computer run, one can include the effects of a heat bath after the computer simulation is already over, for any continuously varying value of the temperature. As the method is based on an exact relation between g(E) and the microcanonical averages of two particular quantities (see text), any numerical accuracy can be achieved depending exclusively on the simulational method one chooses to measure the quoted microcanonical averages.

Acknowledgements I am indebted to my collaborators Thadeu and Hans. I thank also Suzana Moss and Dietrich Stauffer for helpful discussions, and Graham R. Smith for drawing my attention to references [10] and for a critical reading of the manuscript. This work is supported by Brazilian agencies CAPES, CNPq, FAPERJ and FINEP.

P.M.C. de Oliveira / Computer Physics Communications 121–122 (1999) 16–21

References [1] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, J. Chem. Phys. 21 (1953) 1087. [2] K. Binder, Monte Carlo Methods in Statistical Physics, Topics Current Phys., Vol. 7 (Springer, Berlin, 1986). [3] Z.W. Salzburg, J.D. Jacobson, W. Fickett, W.W. Wood, J. Chem. Phys. 30 (1959) 65; A.M. Ferrenberg, R.W. Swendsen, Phys. Rev. Lett. 61 (1988) 2635. [4] P.M.C. de Oliveira, T.J.P. Penna, H.J. Herrmann, Braz. J. Physics 26 (1996) 677 (also in Cond-Mat 9610041). [5] R.H. Swendsen, J.-S. Wang, Phys. Rev. Lett. 58 (1987) 86; U. Wolff, Phys. Rev. Lett. 62 (1989) 361. [6] P.M.C. de Oliveira, T.J.P. Penna, H.J. Herrmann, Eur. Phys. J. B1 (1998) 205. [7] P.M.C. de Oliveira, Int. J. Mod. Phys. C 9 (1998) 497. [8] P.M.C. de Oliveira, Eur. Phys. J. B 6 (1998) 111. [9] J.D. Muñoz, H.J. Herrmann, Int. J. Mod. Phys. C 10 (1999) 95; also these Proceedings, Comput. Phys. Commun. 121–122 (1999) 13. [10] G.R. Smith, A.D. Bruce, J. Phys. A 28 (1995) 6623; Phys. Rev. E 53 (1996) 6530; Europhys. Lett. 39 (1996) 91.

21

[11] J.-S. Wang, R.H. Swendsen, Phys. Rev. Lett. 82 (1999) 476; also these Proceedings, Comput. Phys. Commun. 121–122 (1999) 22. [12] B. Berg, Nature 361 (1993) 708. [13] B. Berg, U.H.E. Hansmann, Eur. Phys. J. B 6 (1998) 395. [14] P.M.C. de Oliveira, Europhys. Lett. 20 (1992) 621; J. Monteiro de F. Neto, S. Moss de Oliveira, P.M.C. de Oliveira, Physica A 206 (1994) 463. [15] K. Binder, Z. Phys. B 43 (1981) 119; Phys. Rev. Lett. 47 (1981) 693. [16] P.D. Beale, Phys. Rev. Lett. 76 (1996) 78. [17] H.W.J. Blöte, E. Luijten, J.R. Heringa, J. Phys. A 28 (1995) 6289; A.L. Tapalov, H.W.J. Blöte, J. Phys. A 29 (1996) 5727. [18] G.M. Torrie, J.P. Valleau, Chem. Phys. Lett. 28 (1974) 578; B.A. Berg, T. Neuhaus, Phys. Lett. B 267 (1991) 249; A.P. Lyubartsev, A.A. Martsinovski, S.V. Shevkunov, P.N. Vorontsov-Velyaminov, J. Chem. Phys. 96 (1992) 1776; E. Marinari, G. Parisi, Europhys. Lett. 19 (1992) 451; B.A. Berg, Int. J. Mod. Phys. C 4 (1993) 249; J. Lee, Phys. Rev. Lett. 71 (1993) 211; B. Hesselbo, R.B. Stinchcombe, Phys. Rev. Lett. 74 (1995) 2151.