Wang-Landau simulations with non-flat distributions
aa r X i v : . [ phy s i c s . c o m p - ph ] F e b Wang-Landau simulations with non-flat distributions
Stefan Schnabel ∗ and Wolfhard Janke † Institut f¨ur Theoretische Physik, Universit¨at Leipzig, IPF 231101, 04081 Leipzig, Germany (Dated: February 2, 2021)We show how the well-known Wang-Landau method can be modified to produce non-flat distri-butions. Through the choice of a suitable profile this can lead to an increase in efficiency for somesystems. Examples for such an enhancement are provided.
INTRODUCTION
Generalized ensemble Monte Carlo methods likereplica exchange [1] or the multicanonical method [2, 3]have been introduced some time ago to improve upon thestandard Metropolis algorithm. Since then it has alsobeen discussed how these generalized ensembles can beespecially designed to increase the performance further.There are for instance strategies to select suitable tem-peratures for the replica exchange method [4, 5]. It is alsounderstood [6, 7] that the multicanonical method can beimproved if one aims at a distribution in energy that isproportional to a non-trivial profile instead of constant.Recently we have shown how simulations of spin glassescan be improved when a specially designed profile is used[8]. However, when it comes to the closely related andwidely applied Wang-Landau method [9] changing to aprofile is not as straightforward as for the other methodsand we are not aware of any application of such a mod-ified Wang-Landau method. In this study we demon-strate how a non-flat Wang-Landau method can be de-vised through a small modification of the original algo-rithm and show with two examples how the performancecan be improved.
WANG-LANDAU ALGORITHM
The Wang-Landau (WL) algorithm employs a function g ( E ) to accept or reject proposed moves from microstate µ to microstate ν with probability P acc ( µ, ν ) = min (cid:18) , g ( E µ ) g ( E ν ) (cid:19) . (1)Here, we assume that the energy E does either assumeonly a finite number of discrete values E ∈ { E , E , . . . } or that a continuous interval E ∈ [ E min , E max ] is dividedinto a finite number of equally wide subintervals (bin-ning) on which g is constant such that there is alwaysa finite number of values for g . Usually g ( E ) is set tounity everywhere in the beginning and after each stepits value for the energy of the currently occupied state µ t is multiplied by a factor f > g ′ ( E µ t ) = f · g ( E µ t ).Depending on the strategy for reducing f the function g ( E ) will converge towards or at least become very simi- lar to the density of states Ω( E ) up to a constant factor.Here we will use a different but equivalent notation. Wewill use a weight function W ( E ) = 1 /g ( E ). If detailedbalance holds – which requires f = 1 among other con-ditions – it is proportional to the probability with whicha microstate with energy E is visited. Since during thecourse of a simulation g ( E ) and equally W ( E ) can ex-tend over many orders of magnitude one often stores anduses its logarithmic values and consequently P acc ( µ, ν ) = min (cid:16) , e ln W ( E ν ) − ln W ( E µ ) (cid:17) (2)and ln W ′ ( E µ t ) = ln W ( E µ t ) − ln f. (3)The concept of the WL method can be expressed asfollows: Reduce the probability of the energy that thewalker is currently at such that over long enough peri-ods of time states that are over-represented in the en-semble are inhibited and a quasi-steady-state is reachedwhere the logarithmic weights ln W ( E ) of all energiesare on average reduced equally. Then the differencesbetween the logarithmic weights of any two energiesln W ( E ) − ln W ( E ) and consequently the ratios of theirweights W ( E ) /W ( E ) remain constant even though theabsolute values change. It is clear that for the standardWL method this steady state is reached if all energies arehit with equal frequency since only then the number ofsubtractions of ln f that each logarithmic weight ln W ( E )experiences is the same for all energies. Therefore, mea-suring this frequency by means of a histogram h ( E ) andtesting for its flatness is a reliable way to evaluate theprogress of the algorithm. NON-FLAT WANG-LANDAU ALGORITHM
The goal is now to alter the method such that thesteady state is reached while the frequency is not con-stant but is proportional to a given profile p ( E ). If wewant to keep the basic procedure of the WL method andto change ln W ( E ) only at the energy of the currentlyoccupied state µ t it is then clear that we have to modifyEq. (3) to ln W ′ ( E µ t ) = ln W ( E µ t ) − ln fp ( E µ t ) . (4)If now each energy is hit with a frequency proportionalto the profile function p ( E ), the accumulated changes ofln W ( E ) will again be independent of E since p ( E ) willcancel out. The rule for accepting updates in Eq. (2) re-mains unaffected. Ignoring any error saturation, the finalresult of this procedure, approached when f is sufficientlyclose to unity, is a weight function W ( E ) = p ( E )Ω( E ) (5)where Ω( E ) is again the density of states. This modifi-cation of the WL method can on the one hand be seenas the introduction of energy-dependent weight modifi-cation factors: ˜ f ( E ) = f /p ( E ) . (6)On the other hand this means that if energy-dependentweight modification factors ˜ f ( E ) are used the non-flatWL algorithm will produce histograms h ( E ) that are pro-portional to 1 / ln( ˜ f ( E )) ∝ p ( E ).Previously we have introduced the profile p ( E ) as afunction proportional to the desired frequency leaving itnot fully determined. However, at this point its mag-nitude becomes relevant since it directly influences themagnitude of the changes to ln W ( E ). A simple choice is p ( E ) ≥
1; it ensures that the weights experience changesequal to or smaller than those that would be imposed bythe original WL algorithm.One important aspect of the WL method is the way themodification factor f is reduced over time. The originalstrategy is to sample with a constant f until a sufficientlyflat histogram h ( E ) has been produced and to reduce f by taking its square root. Since our stated goal is to pro-duce non-flat distributions, flatness of the histogram isnot to be expected. Now, it is the ratio of histogram andprofile h ( E ) /p ( E ) that will approach a constant functionand can be used instead of the histogram alone.Liang et al. [10] as well as Belardinelli and Pereyra [11]have suggested an alternative way of reducing f . Theypropose to use distinct values f t at any time t and oneexample out of a family of possibilities for the history of f is ln f t = t max( t , t ) (7)for some t >
0. With our modification this method canbe used unchanged.To conclude the dicussion of the modification of theWL method, we come back to the original notation withthe direct approximation of the density of states using g ( E ). Using Eq. (5) to replace W ( E ) by p ( E ) /g ( E ) weobtain for the acceptance probability P acc ( µ, ν ) = min (cid:18) , g ( E µ ) p ( E ν ) g ( E ν ) p ( E µ ) (cid:19) (8a)= min (cid:18) , e ln g ( E µ ) − ln g ( E ν ) p ( E ν ) p ( E µ ) (cid:19) (8b) while g ( E ) is modified according toln g ′ ( E µ t ) = ln g ( E µ t ) + ln fp ( E µ t ) . (9)It should be pointed out that a variation of the WLmethod that allows any desired profile in E instead of theconstant distribution that comes with the standard WLalgorithm has already been introduced [10]. However,the proposed procedure is somewhat cumbersome sinceit requires the modification of the weights W ( E ) for allvalues (or intervals) of E at every time step [21]. Themethod we propose here is simpler and also more in tunewith the basic principle of WL sampling. APPLICATIONSProof-of-concept: Ising and Potts Models
In order to demonstrate that the method is working,i.e., that it is able to produce histograms in accordancewith desired profiles we performed as basic test simula-tions of a L = 32 Ising model and a L = 64 , q = 10 Pottsmodel on square lattices with periodic boundary condi-tions. For the reduction of f we modified the originalrecipe from [9] for our algorithm through replacing h ( E )by h ( E ) /p ( E ): The modification factor is reduced ac-cording to f ′ = f / and the histogram reset to h ( E ) := 0if the minimum of h ( E ) /p ( E ) is larger than a certainfraction of its meanmin ( h ( E i ) /p ( E i )) ≥ ρB B X i =1 h ( E i ) /p ( E i ) , (10)where B is the number of values (subintervals) of theenergy and we chose ρ = 0 .
9. The results displayed inFig. 1 show that the method works as expected and thatthe histograms reproduce the desired profiles very well.
Ising Spin Glass
Next, we consider an Edwards-Anderson spin glass [13]on a cubic lattice with the Ising Hamiltonian H = − X h ij i J ij s i s j , s i , J ij ∈ {− , } , (11)where the sum goes over all pairs of adjacent spins andthe bonds are randomly chosen for any new disorder re-alization. Due to their rough energy landscape spin-glasssystems pose an interesting challenge and serve as bench-mark cases for Monte Carlo methods. In a recent study[8] it has been demonstrated that for equilibrium ( f = 1)simulations a power-law profile with a strong emphasis of . . . − − − − h ( E ) / P h ( E ) p ( E ) Ef = ef = exp { − } f = exp { − } f = exp { − } p ( E )00 . . . . − − − −
500 0 500 1000 1500 2000012345678 h ( E ) / P h ( E ) p ( E ) E f = ef = exp { − } f = exp { − } f = exp { − } p ( E ) FIG. 1. Histograms for different values of f from simulationsof the L = 64 , q = 10 Potts model with a sequence of linesegments as profile (top) and the L = 32 Ising model with acurved profile taken from [12] (bottom) on square lattices. t fl a t t prof FIG. 2. Scatter plot of the simulation times required for in-dividual spin-glass samples. low-energy states is superior to a flat histogram and al-lows a more rapid exploration of state space. We expectthat a similar acceleration can be achieved for f close tobut larger than unity and that the performance of the WLalgorithm can thus be improved. Note that the experi-ment described here is kept rather simple. It is intendedto be a proof-of-concept and further improvements to themethod are possible. We generate 100 disorder realiza-tions { J ij } with N = 10 spins and run WL simulationsfor all of them with a constant profile p ( E ) = 1 as wellas with p ( E ) = (cid:18) E (cid:19) − . , (12) which is based on the profile that was used for N ≤ in [8]. In the beginning f = e and the simulation isstopped if f < exp (cid:8) − (cid:9) . After an initial phase of10 N attempted spin flips (steps) every 10 N steps itis tested whether according to Eq. (10) a flat histogramwith ρ = 0 . f is reduced: f ′ = f / . If during the simulationa new lowest energy is found and if f < exp (cid:8) − (cid:9) wereset f = exp (cid:8) − (cid:9) to allow the simulation to adjustthe weight of the states at the new energy and prevent itfrom getting trapped there.The required Monte Carlo time in units of N stepsfor 86 disorder realizations is shown in Fig. 2. For theremaining 14 samples the two methods did not find thesame lowest energy and can, therefore, not be compared.In thirteen cases the WL simulation with the power-lawprofile reached a lower energy while the flat-histogramversion reached a lower energy once. Where comparisonsare possible we see that for the hard samples, i.e., thedisorder realizations that require long simulations, thepower-law profile in Eq. (12) is more than three timesfaster while it leads to less efficient simulations for thevery easy samples. The aggregated running times for all86 samples shown are 1 . × for the power-law profilevs 2 . × for the flat distribution. Lennard-Jones Polymer
As a second example for a useful application we ap-ply the method to a Lennard-Jones polymer. We investi-gated this system some time ago and details of the modeland the results can be found in [14–17]. For our cur-rent purpose it is sufficient to say that it is an off-latticebead-spring polymer model that for the considered sizeof N = 147 beads possesses a very stable state of icosahe-dral geometry (Fig. 3) at low temperatures. At mediumtemperatures one observes an unstructured dense glob-ular droplet and at high temperatures beyond the so-called Θ-point we find extended conformations that re-semble self-avoiding random walks. Therefore, there aretwo transitions [22], one at energies around E ≈ − E ≈ − p ( E ) = (cid:26) Π , if | E + 670 | ≤ , else , (13)in order to enhance the simulation. We perform WLsimulations for different values of Π. We require minimalflatness of the histogram, i.e., we proceed to the nextiteration when all energies have been visited at least once.Then the modification factor is reduced according to f ′ = FIG. 3. Icosahedral low-energy conformation of the polymer.Colors indicate the different layers containing 1,12,42, and 92beads. f / until ln f ≤ − (starting with f = e ). For theMonte Carlo updates we use an elaborate set of movesthat is discussed in detail in [18].At extreme low energies the system undergoes a finalenergy optimization which does not significantly affectthe position of the beads but rearranges the bonds suchthat unfavorable distances are avoided if possible. Thishas for instance the effect that at the global energy mini-mum conformation – the lowest microstate we found with E = − .
161 is shown in Fig. 3 – only exactly one bondconnects the different layers. This process slows downthe simulation in the proximity of the ground state andone might try to counter this by an additional increaseof the profile in this region. However, here we just wantto look on the effect of the ‘solid-liquid’ transition andhence exclude the ground state by restricting the energyrange to − ≤ E ≤ E . Theincrease of p ( E ) enhances this frequency and thus accel-erates the simulation. Although the length of individualWL simulations is to some extent subject to chance andchanges with the seed of the random number generator,the general trend is obvious and Π ∈ [50 ,
75] is clearly asuperior choice to a flat (Π = 1) distribution.This method should work in all cases where such awell localized single bottleneck is hampering the randomwalker and we expect that also studies of other systemswith first-order-like phase transitions can benefit from it. E − − − − − − − − − (ln f ) Π = 1 t − − − − − − − − − (ln f ) Π = 10Π = 20Π = 30Π = 50Π = 75Π = 100Π = 5 FIG. 4. Time series for different levels of enhancement Π ofthe transition region using the profile (13). All time seriescover the same vertical energy range E ∈ [ − , t is measured in units of 10 N = 1 . × updates which is also the interval between individual pointsin the lower plots. For the time series with Π = 1 that intervalis 5 × N . CONCLUSION
We have shown how the Wang-Landau method canwith minimal effort be adapted to produce non-flat his-tograms with a desired profile for any value of the mod-ification factor. As expected the advantages of balancedsimulations with a profile carry over to Wang-Landausampling. Spin glasses can be simulated more efficientlywith the here proposed non-flat Wang-Landau algorithmif the profile is high at low energies and it is likely thatsimilar gains can be achieved for other glassy systems aswell. As shown in the case of a polymer, transitions be-tween different macrostates occur more often if the profileis enhanced in the transition region and the performanceof the simulation method can thus be increased.To introduce the concept of our non-flat Wang-Landaumethod, we have in this Letter focused on an implemen-tation based on Monte Carlo simulations as in the orig-inal publication [9]. It is, however, also easily possibleand straightforward to boost standard flat Wang-Landaumolecular dynamics simulations [19, 20] by employing anon-flat generalization along the same lines as discussedhere.
ACKNOWLEDGEMENT
This project was funded by the Deutsche Forschungs-gemeinschaft (DFG, German Research Foundation) un-der Project No. 189 853 844–SFB/TRR 102 (projectB04). It was further supported by the Deutsch-Franz¨osische Hochschule (DFH-UFA) through the Doc-toral College “L4” under Grant No. CDFA-02-07. ∗ [email protected] † [email protected][1] K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. ,1604 (1996).[2] B. A. Berg and T. Neuhaus, Phys. Lett. B , 249(1991).[3] B. A. Berg and T. Neuhaus, Phys. Rev. Lett. , 9(1992).[4] H. G. Katzgraber, S. Trebst, D. A. Huse, and M. Troyer,J. Stat. Mech., P03018 (2006).[5] E. Bittner, A. Nußbaumer, and W. Janke, Phys. Rev.Lett. , 130603 (2008).[6] S. Trebst, D. A. Huse, and M. Troyer, Phys. Rev. E ,046701 (2004).[7] B. Hesselbo and R. B. Stinchcombe, Phys. Rev. Lett. ,2151 (1995).[8] F. M¨uller, S. Schnabel, and W. Janke, Phys. Rev. E ,053303 (2020). [9] F. Wang and D. P. Landau, Phys. Rev. Lett. , 2050(2001).[10] F. Liang, C. Liu, and R. Carrol, J. Am. Stat. Assoc. ,305 (2007).[11] R. E. Belardinelli and V. D. Pereyra, J. Chem. Phys. , 184105 (2007).[12] A. de Saint-Exup´ery, The Little Prince (Reynal & Hitch-cock, New York, 1943).[13] S. F. Edwards and P. W. Anderson, J. Phys. F , 965(1975).[14] S. Schnabel, T. Vogel, M. Bachmann, and W. Janke,Chem. Phys. Lett. , 201 (2009).[15] S. Schnabel, M. Bachmann, and W. Janke, J. Chem.Phys. , 124904 (2009).[16] S. Schnabel, D. T. Seaton, D. P. Landau, and M. Bach-mann, Phys. Rev. E , 011127 (2011).[17] J. C. S. Rocha, S. Schnabel, D. P. Landau, and M. Bach-mann, Phys. Rev. E , 022601 (2014).[18] S. Schnabel, W. Janke, and M. Bachmann, J. Comput.Phys. , 4454 (2011).[19] T. Nagasima, A. R. Kinjo, T. Mitsui, and K. Nishikawa,Phys. Rev. E , 066706 (2007).[20] C. Junghans, D. Perez, and T. Vogel, J. Chem. TheoryComput.10