A new approach to Monte Carlo simulations in statistical physics: Wang-

Landau sampling
D. P. Landau, Shan-Ho Tsai, and M. Exler

Citation: Am. J. Phys. 72, 1294 (2004); doi: 10.1119/1.1707017

View online: http://dx.doi.org/10.1119/1.1707017
View Table of Contents: http://ajp.aapt.org/resource/1/AJPIAS/v72/i10
Published by the American Association of Physics Teachers

Additional information on Am. J. Phys.

Journal Homepage: http://ajp.aapt.org/
Journal Information: http://ajp.aapt.org/about/about_the_journal
Top downloads: http://ajp.aapt.org/most_downloaded
Information for Authors: http://ajp.dickinson.edu/Contributors/contGenInfo.html

Downloaded 26 Sep 2012 to 136.159.235.223. Redistribution subject to AAPT license or copyright; see http://ajp.aapt.org/authors/copyright_permission
A new approach to Monte Carlo simulations in statistical physics:
Wang-Landau sampling
D. P. Landau, Shan-Ho Tsai,a) and M. Exlerb)
Center for Simulational Physics, The University of Georgia, Athens, Georgia 30602
共Received 15 December 2003; accepted 20 February 2004兲
We describe a Monte Carlo algorithm for doing simulations in classical 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 extract an estimate for the density of states. The probability can be
computed at any temperature by weighting the density of states by the appropriate Boltzmann factor.
Thermodynamic properties can be determined from suitable derivatives of the partition function
and, unlike ‘‘standard’’ methods, the free energy and entropy can also be computed directly. To
demonstrate the simplicity and power of the algorithm, we apply it to models exhibiting first-order
or second-order phase transitions. © 2004 American Association of Physics Teachers.
关DOI: 10.1119/1.1707017兴

I. INTRODUCTION such as Metropolis importance sampling and Swendsen-

Wang cluster flipping, generate an unnormalized canonical
Computer simulation now plays a major role in statistical distribution
physics,1 particularly for the study of phase transitions and P 共 E,T 兲 ⫽g 共 E 兲 e ⫺E/k B T 共1兲
critical phenomena, and is an important tool for teaching and
understanding thermodynamics and statistical mechanics.2 at a given temperature T (k B denotes the Boltzmann con-
The reason for its importance is that all but the simplest stant兲. Such distributions are so narrow that multiple runs are
models are theoretically intractable, and only approximate usually needed to describe thermodynamic quantities over a
methods can generally be used. In particular, stochastic tech- significant range of temperatures. Because g(E) does not
niques known as Monte Carlo 共MC兲 simulations have proven depend on the temperature, we can construct canonical dis-
to be very powerful. The standard MC method, developed a tributions at any temperature if we can estimate g(E) with
half-century ago, is the Metropolis importance sampling high accuracy for all energies. Once g(E) is known, we can
algorithm,3 but more recently new, more efficient algorithms calculate the partition function as
have begun to play a role in allowing simulations to achieve
the resolution that is needed to accurately locate and charac- Z⫽ 兺
兵 configurations其
e ⫺E/k B T ⫽ 兺E g 共 E 兲 e ⫺E/k T ,
B 共2兲
terize phase transitions.1
The motivation for these new developments is that tradi- and the model is essentially ‘‘solved,’’ because most thermo-
tional methods exhibit long time scales, thus requiring long dynamic quantities can be calculated from Z.
simulations. At first-order phase transitions metastable states Although MC methods are already very powerful,1 there
appear, and critical slowing down becomes a problem at con- has been no efficient algorithm to calculate g(E) very accu-
tinuous transitions. For spin systems, beginning with the rately for large systems. Even for exactly solvable models,
seminal work of Swendsen and Wang,4 and extended by such as the two-dimensional 共2D兲 Ising model, g(E) cannot
Wolff,5 cluster algorithms have been used to reduce critical be calculated exactly for large systems.14 All methods based
slowing down. The multicanonical ensemble method6 –10 was on accumulation of histogram entries13,15–18 have the prob-
introduced to overcome the tunneling barrier between coex- lem of scalability for large systems.
isting phases at first-order transitions, and has general utility In this paper we describe a new, general, and efficient MC
for systems with a rough energy landscape.7,11,12 In both situ- algorithm 共generally known as the ‘‘Wang-Landau algo-
ations, histogram reweighting techniques13 can be applied in rithm’’兲 that offers substantial advantages over existing
the analysis to increase the amount of information that can approaches.19 We will explain the algorithm in detail and
describe its application for first- and second-order phase tran-
be gleaned from simulational data, but the applicability of
sitions. Unlike conventional MC methods that directly gen-
reweighting is severely limited in large systems by the sta-
erate a canonical distribution g(E)e ⫺E/k B T at a given tem-
tistical quality of the wings of the histogram. This limitation
is important in systems with competing interactions for perature T, this approach estimates g(E) directly and
which short-range order effects might occur over very broad accurately via a random walk that produces a flat histogram
temperature ranges or give rise to frustration and produce a in energy space. The estimate for g(E) is improved at each
very complicated energy landscape, thus reducing the effi- step of the random walk, using a carefully controlled modi-
fication factor, to produce a result that converges to the real
ciency of standard methods.
value quickly.
The partition function can be expressed in terms of a den-
Wang-Landau sampling19 has proven to be very useful and
sity of states g(E), the number of all possible states 共or efficient in many different applications, including studies of
configurations兲 for an energy level E of the system, but di- complex systems with rough energy landscapes. For ex-
rect estimation of g(E) has not usually been the goal of ample, the method has been used in studies of a Potts
simulations. Instead, most conventional MC algorithms,1 antiferromagnet,20 random spin systems,21 quantum

1294 Am. J. Phys. 72 共10兲, October 2004 http://aapt.org/ajp © 2004 American Association of Physics Teachers 1294

Downloaded 26 Sep 2012 to 136.159.235.223. Redistribution subject to AAPT license or copyright; see http://ajp.aapt.org/authors/copyright_permission
systems,22–24 fluids,25,26 binary Lennard-Jones glass,27 liquid 1 is smaller than or equal to the ratio g(E 1 )/g(E 2 )]. If the
crystals,28 polymers,25,29 proteins,30,31 other molecular trial state with energy E 2 is accepted, we multiply the exist-
systems,32,33 atomic clusters,34 optimization problems,35 and ing value of g(E 2 ) by a modification factor f ⬎1, that is,
combinatorial number theory.36 Generalizations and further g(E 2 )→ f ⫻g(E 2 ), and we update the existing entry for
studies of this sampling technique have been carried out by
H(E 2 ) in the energy histogram, that is, H(E 2 )→H(E 2 )
several authors.37– 41
Although the Wang-Landau method can be applied to ⫹1. If the random walk rejects the trial move and remains at
many different types of systems, we will describe it here the same energy level E 1 , we modify the existing density of
only in the context of classical spin systems with discrete states g(E 1 ) by the same modification factor; that is, g(E 1 )
energy values. Therefore, when we refer to the density of → f ⫻g(E 1 ), and we update the existing entry for H(E 1 );
states g(E), we do not mean an actual density, but the num- that is, H(E 1 )→H(E 1 )⫹1. Because g(E) becomes very
ber of states for a given energy E. The two simple models of large, in practice it is preferable to work with the logarithm
interest are the Ising model,42 which has a second-order of the density of states, so that all possible ln关g(E)兴 will fit
phase transition, and the Q-state Potts model43 with Q⫽8, into double precision numbers. Therefore, each update of the
which undergoes a first-order phase transition. density of states is implemented as ln关g(E)兴→ln关g(E)兴
⫹ln(f ), and the ratio of density of states in Eq. 共3兲 is com-
puted as exp兵ln关g(E1)兴⫺ln关g(E2)兴其.
II. THE WANG-LANDAU ALGORITHM
A reasonable, although not necessarily optimal choice of
If we perform an unbiased random walk in energy space the initial modification factor is f ⫽ f 0 ⫽e 1 ⯝2.71828, which
by changing the states of the spins at random and accepting allows us to reach all possible energy levels quickly even for
all energy values thus obtained, the histogram of the energy a large system. If f 0 is too small, the random walk will spend
distribution should converge to the density of states g(E) in a very long time to reach all possible energies; however, too
the limit of a very long random walk that visits all possible large a choice of f 0 will lead to large statistical errors. We
spin configurations of the system. In practice it is forbid- proceed with the random walk in energy space until we ob-
dingly difficult to realize such a long random walk with our tain a ‘‘flat’’ histogram H(E). We typically check whether
current computer resources, given the extremely large num- the histogram is flat after every 10 000 MC sweeps, where
ber of spin configurations. For example, the Ising model on a one MC sweep corresponds to randomly picking N spins and
10⫻10 square lattice already has 2 100⬇1.3⫻1030 spin con- thus generating N trial states (N denotes the total number of
figurations! spins on the lattice兲. When the histogram is flat, all the pos-
The Wang-Landau sampling method performs random sible energies have been roughly visited an equal number of
walks in energy space by changing the states of spins ran- times, and the density of states converges to the true value
domly, but the energy E associated with each spin configu- with an accuracy proportional to the modification factor
ration is only accepted with a probability that is proportional ln(f ). We then reduce the modification factor by using a
to the reciprocal of the density of states. During the random function such as f 1 ⫽ 冑 f 0 , reset the histogram to H(E)⫽0
walk, we also accumulate the histogram H(E) in energy for all values of E, and begin the next level random walk
space, a quantity that keeps track of the number of visits at during which we modify the density of states with the
each energy level E 关each time an energy E is visited, the
smaller modification factor f 1 for each step. Each level ran-
corresponding entry in H(E) is incremented by 1兴. The al- dom walk is referred to as one iteration in the algorithm.
gorithm modifies the estimate of the density of states by a Note that the spin configuration and the density of states are
multiplicative factor f , and uses the updated density of states never reset during the simulation. We continue performing
to perform a further random walk in energy space. With this the random walk until the histogram H(E) is flat again, and
choice of acceptance probability, each random walk gener-
then we reduce the modification factor f i⫹1 ⫽ 冑 f i , reset the
ates a flat histogram for the energy distribution. The modifi-
cation factor f is carefully controlled, and at the end of the histogram to H(E)⫽0 for all values of E, and restart the
simulation, it should be very close to 1, which is the ideal random walk. We stop the simulation when the modification
case of the random walk with the true density of states. factor is smaller than a predefined value 共such as f final
At the beginning of the simulation, g(E) is unknown, and ⫽exp(10⫺8 )⯝1.000 000 01). The modification factor acts as
we make an initial guess for it. The simplest approach is to a control parameter for the accuracy of the density of states
set g(E)⫽1 for all possible energies E. The initial spin con- during the simulation and also determines how many MC
figuration for the entire lattice can be chosen arbitrarily. sweeps are necessary for the whole simulation.
Then, a random walk in energy space is begun by forming It is impossible to obtain a perfectly flat histogram and the
trial states, each of which is produced by randomly picking a phrase ‘‘flat histogram’’ in this paper means that the histo-
spin and randomly changing its state. In general, if E 1 and gram H(E) for all possible E is not less than x% of the
E 2 are energies before and after a spin value is changed, the average histogram 具 H(E) 典 , where x% is chosen according
transition probability from energy E 1 to E 2 is to the size and complexity of the system and the desired

冉 冊
accuracy of the density of states. For the 2D Ising model
g共 E1兲 with only nearest-neighbor couplings on small lattices, this
p 共 E 1 →E 2 兲 ⫽min ,1 . 共3兲 percentage can be chosen as high as 95%, but for large sys-
g共 E2兲
tems the criterion for ‘‘flatness’’ may never be satisfied if we
Equation 共3兲 implies that if g(E 2 )⭐g(E 1 ), the state with choose too high a percentage, and the program might run
energy E 2 is accepted; otherwise it is accepted with a prob- forever.
ability g(E 1 )/g(E 2 ) 关that is, the state with energy E 2 is ac- Clearly, one essential constraint is that g(E) should con-
cepted if a random number picked uniformly between 0 and verge to the true value. The accuracy of the estimate for

1295 Am. J. Phys., Vol. 72, No. 10, October 2004 D. P. Landau, Shan-Ho Tsai, and M. Exler 1295

Downloaded 26 Sep 2012 to 136.159.235.223. Redistribution subject to AAPT license or copyright; see http://ajp.aapt.org/authors/copyright_permission
g(E) is proportional to ln(f ) at that iteration. However, exponential growth of the density of states in energy space,
ln(f final) cannot be chosen arbitrarily small or the modified this process is inefficient because the histogram is accumu-
ln关g(E)兴 will not differ from the unmodified one to within the lated linearly. Instead, in Wang-Landau sampling we modify
number of digits in the double precision numbers used in the g(E) at each step of the random walk, and this modification
simulation. If this happens, the algorithm no longer con- allows us to approach the true distribution much faster than
verges to the true value, and the program may run forever. If conventional methods, especially for large systems. 共We also
f final is within the double precision range but is too small, the accumulate histogram entries during the random walk, but
we only use them to check whether the histogram is flat
calculation might take excessively long to finish.
enough to go to the next level random walk.兲
A simple recipe for reducing the modification factor is to
Although the total number of configurations increases ex-
take a square-root function, and f approaches 1 as the num-
ponentially with the size of the system, the total number of
ber of iterations approaches infinity. 共There is no reason why possible energies increases linearly with the size of system,
any function cannot be used as long as it decreases f mono- so it is easy to calculate g(E) with a random walk in energy
tonically to 1. A simple and efficient formula is f i⫹1 ⫽ f 1/n
i , space for a large system. Consider, for example, a Q-state
where n⬎1. The value of n can be chosen according to the Potts model on a L⫻L lattice with nearest-neighbor
available CPU time and the expected accuracy of the simu- interactions.43 For Q⭓3, the number of possible energies is
lation. For the systems that have been studied, the choice of
about 2N, where N⫽L 2 is the total number of the lattice
n⫽2 yields good accuracy in a relatively short time, even for
sites. However, the average number of possible states for
large systems.兲
each energy level is as large as Q N /2N, where Q N is the total
For the initial modification factor ln(f 0)⫽1 and the final
number of possible configurations of the system. This large
factor ln(f final)⫽10⫺8 , the total number of iterations is 27. number is why we cannot simply use a computer to realize
We do not set a predetermined number of MC sweeps for all possible states and why efficient and fast algorithms are
each iteration, but rather let the program check periodically required.
whether the established criterion for a flat histogram is sat- At the end of the simulation, the Wang-Landau algorithm
isfied. Generally, the number of MC sweeps needed to satisfy provides only a relative density of states for different ener-
the criterion increases as we reduce the modification factor, gies. To extract the correct density of states g n (E) for the
but we cannot predict the exact number of MC sweeps
Q-state Potts model, we can either use the fact that the total
needed for each iteration before the simulation. It is prefer-
able to allow the program to decide how much simulational number of possible states is 兺 E g n (E)⫽Q N , or that the num-
effort is needed for a given modification factor f i . Nonethe- ber of ground states 共where E⫽⫺2N) is Q. By using the
less, we need to perform some test runs to make sure that the former rescaling condition, the correct normalized density of
program will finish within a given time. states g n (E) can be obtained from the simulation data g(E),
The simulation method can be further enhanced by per- by the relation ln关gn(E)兴⫽ln关g(E)兴⫺ln关 兺Eg(E)兴⫹N ln(Q),
forming multiple random walks, each for a different range of whereas the latter condition leads us to use ln关gn(E)兴
energy, either serially or in parallel. We can restrict the ran- ⫽ln关g(E)兴⫺ln关g(E⫽⫺2N)兴⫹ln(Q). For simplicity, we will
dom walk to remain in the range by rejecting any move out denote the normalized density of states simply as g(E) in the
of that range.19,41 The resultant parts of the density of states following. The latter normalization guarantees the accuracy
can then be joined together. of the density of states at low energy levels, which is impor-
During the random walk 共especially in the early itera- tant in the calculation of thermodynamic quantities at low
tions兲, the algorithm does not satisfy the detailed balance temperature. With this normalization, when T⫽0, we can
condition exactly, because g(E) is modified constantly dur- obtain exact solutions for the internal energy, entropy, and
ing the random walk. After many iterations, however, g(E) free energy when we calculate these quantities from the den-
converges to the true value as the modification factor ap- sity of states. If we apply the normalization that the total
proaches 1. If p(E 1 →E 2 ) is the transition probability from number of states is Q N , we cannot guarantee the accuracy of
energy E 1 to energy E 2 , the ratio of the transition probabili- g(E) for energies at or near the ground state, because the
ties from E 1 to E 2 and from E 2 to E 1 can be calculated very rescaling factor is dominated by the maximum density of
easily as states. We can use one of these two normalizations to obtain
the absolute density of states, and use the other normaliza-
p 共 E 1 →E 2 兲 g 共 E 1 兲 tion to check the accuracy of the result.
⫽ , 共4兲
p 共 E 2 →E 1 兲 g 共 E 2 兲 One of the advantages of the Wang-Landau method is that
the density of states does not depend on the temperature. For
where we have used Eq. 共3兲. In other words, the random example, the internal energy U(T) can be calculated by
walk algorithm satisfies the detailed balance:
兺 E Eg 共 E 兲 e ⫺E/k B T
1 1 U共 T 兲⫽ ⬅具E典, 共6兲
p 共 E 1 →E 2 兲 ⫽ p 共 E 2 →E 1 兲 , 共5兲 兺 E g 共 E 兲 e ⫺E/k B T
g共 E1兲 g共 E2兲
and the specific heat C(T) can be determined from the fluc-
where 1/g(E 1 ) is the probability at the energy E 1 and tuations in the internal energy
p(E 1 →E 2 ) is the transition probability from E 1 to E 2 . We
conclude that the detailed balance condition is satisfied with ⳵ U 共 T 兲 具 E 2典 ⫺ 具 E 典 2
accuracy proportional to the modification factor ln(f ). C共 T 兲⫽ ⫽ . 共7兲
⳵T k BT 2
Almost all recursive methods update the density of states
by using the histogram data directly, and only after enough We can also access some quantities, such as the Helmholtz
histogram entries are accumulated.6,9,11,44 –51 Because of the free energy and entropy, that are not directly available from

1296 Am. J. Phys., Vol. 72, No. 10, October 2004 D. P. Landau, Shan-Ho Tsai, and M. Exler 1296

Downloaded 26 Sep 2012 to 136.159.235.223. Redistribution subject to AAPT license or copyright; see http://ajp.aapt.org/authors/copyright_permission
conventional MC simulations. For example, by using con-
ventional MC methods the entropy can be estimated by inte-
grating over other thermodynamic quantities, such as the
specific heat, but the result is not always reliable because the
specific heat itself is not easy to accurately determine, par-
ticularly considering its divergence at a phase transition.
However, the free energy F(T) can be calculated directly
from the partition function Z using

F 共 T 兲 ⫽⫺k B T ln共 Z 兲 ⫽⫺k B T ln 冉兺 E

冊
g 共 E 兲 e ⫺E/k B T , 共8兲

and the entropy can then be easily computed by

U 共 T 兲 ⫺F 共 T 兲
S共 T 兲⫽ . 共9兲
T
We point out that even for relatively small lattices, the
Fig. 1. Logarithm of the density of states, ln关g(E)兴, of the 2D Ising model
partition function may be too large to fit into a double pre- for L⫽16. The relative errors of the simulational densities of states are
cision number, in which case it cannot be easily computed in shown in the inset.
practice. Nevertheless, the thermodynamic quantities can
still be readily computed if we note that

兺E X 共 E 兲 g 共 E 兲 e ⫺E/k T ⫽e ␭ 兺E X 共 E 兲 e ln[g(E)]⫺E/k T⫺␭ ,

B B netic Ising model42 with nearest-neighbor interactions on a
L⫻L square lattice with periodic boundary conditions. Each
共10兲 of the N⫽L 2 lattice sites i has a spin denoted as ␴ i , which
where X(E) is a general function of E and ␭ is the largest can assume the values ␴ i ⫽⫹1 for spin up and ␴ i ⫽⫺1 for
exponent, ln关g(E)兴⫺E/kBT. The summation on the right-hand spin down. The interaction Hamiltonian is given by
side of Eq. 共10兲 can be computed and the factor e ␭ does not
have to be evaluated. Because this factor appears in the nu-
merator and the denominator of Eq. 共6兲, it cancels. This can-
H⫽⫺ 兺
具 i, j 典
␴ i␴ j , 共11兲

cellation also occurs in the evaluation of the specific heat. where 具 i, j 典 denotes distinct pairs of nearest-neighbor sites;
The free energy, which is proportional to the logarithm of the
the number of energies for this system is N⫺1 for even L.
partition function, can also be computed without evaluating
This model provides an ideal benchmark for new
e ␭ explicitly.
algorithms,13,52 and is also an ideal laboratory for testing
Statistical errors in the thermodynamic quantities can be
theory,14,53 because this model can be solved exactly.
estimated by repeating the simulation several times using
With the exact solution for the partition function on finite-
different random number sequences, and then computing the
size systems,54 and the expansion of the expression by Math-
averages and fluctuations in these quantities.
ematica, the density of states for the Ising model on a square
With the histogram reweighting method,13 it is possible to
lattice can be obtained exactly.14 Beale14 obtained the exact
use simulational data at specific temperatures to obtain com-
density of states up to L⫽32, and using Beale’s program,
plete thermodynamic information near, or between, those
temperatures. Unfortunately, it is usually quite difficult to Wang and Landau19 were able to compute g(E) for L⫽50.
obtain accurate information in the region far away from the In this paper we show results for L⫽16, but Wang-Landau
simulated temperature due to difficulties in obtaining good sampling has been used to determine g(E) for lattices up to
statistics, especially for large systems where the canonical L⫽256 for which there is currently no exact solution.19
distributions are very narrow. With Wang-Landau sampling, The estimate of the density of states for L⫽16 using
the histogram is ‘‘flat,’’ and we have essentially the same Wang-Landau sampling is shown in Fig. 1, along with the
statistics for all energy levels. Because the output of the exact results by Beale.14 The initial and final modification
simulation is the density of states, which does not depend on factors for the random walks were ln(f 0)⫽1 and ln(f final)
the temperature, we can then calculate most thermodynamic ⫽10⫺8 . The histogram H(E) was considered flat when all
quantities at any temperature without repeating the simula- entries were not less than 80% of the average 具 H(E) 典 . The
tion. The algorithm is especially useful for obtaining thermo- absolute density of states in Fig. 1 is obtained by the condi-
dynamic information at low temperatures, or at the transition tion that the number of ground states is 2 for the 2D Ising
temperature where the conventional MC algorithm is not so model. With the logarithmic scale used in Fig. 1, the simu-
efficient. lational data and exact solution overlap perfectly with each
other. In the inset of Fig. 1, we show the relative error ␧,
III. APPLICATION TO A SECOND-ORDER PHASE which is defined by the ratio between the error of the simu-
TRANSITION lational data and exact values for any quantity X as
Wang-Landau sampling is very efficient for the study of 兩 X sim⫺X exact兩
second-order phase transitions, because it sidesteps critical ␧共 X 兲⬅ . 共12兲
X exact
slowing down at the critical temperature T c and the slow
dynamics at low temperature. To check the accuracy and We see that ␧ 关 ln(g)兴 is smaller than 0.2% for most of the
convergence of the method, we apply it to the 2D ferromag- region.

1297 Am. J. Phys., Vol. 72, No. 10, October 2004 D. P. Landau, Shan-Ho Tsai, and M. Exler 1297

Downloaded 26 Sep 2012 to 136.159.235.223. Redistribution subject to AAPT license or copyright; see http://ajp.aapt.org/authors/copyright_permission
 the insets, which show the relative errors for the respective
thermodynamic quantities. The relative errors are quite small
for the entire temperature region from k B T⫽0 – 8.
Note that because the system has a second-order phase
transition, the first derivative of the free energy is a continu-
ous function of temperature. There are no jumps in either the
internal energy or the entropy even in the limit as the system
size goes to infinity.
The random number generator used in our simulation was
a shift-register algorithm denoted as R1279.1 The average
number of visits to each energy for the entire duration of the
simulation 共adding the average number of visits for all itera-
tions兲 was roughly 106 . The CPU time of the simulation to
obtain the density of states shown in Fig. 1 was less than 3
min using a GNU compiler on a Pentium 4 共1.3 GHz兲 pro-
Fig. 2. The canonical distribution at the transition temperature P(E,T c ) cessor.
⫽g(E)e ⫺E/k B T c for the L⫽16 Ising model. The inset shows the canonical
distribution at a temperature slightly above and below T c for the same
system.
IV. APPLICATION TO A FIRST-ORDER PHASE
TRANSITION
In this section, we apply the algorithm to a model with a
We can calculate the canonical distribution using Eq. 共1兲 at
first-order phase transition.55,56 In such cases, the internal
essentially any temperature without performing multiple
energy and the entropy have discontinuities at the transition,
simulations. In Fig. 2 we show the resultant canonical distri-
at which both ordered and disordered states coexist. We con-
bution at the critical temperature T c , which exhibits a single
sider the 2D Q⫽8 Potts model43 on L⫻L square lattices
peak. The distributions at temperatures above and below T c
with nearest-neighbor interactions and periodic boundary
are also single peaked, as illustrated in the inset of Fig. 2.
conditions. The total number of spins is N⫽L 2 , and the
It is also important to study the influence of the errors in
Hamiltonian can be written as
the density of states on the calculated thermodynamic quan-
tities. In Fig. 3 we show the internal energy, the specific heat,
the Helmholtz free energy, and the entropy as a function of H⫽⫺ 兺
具 i, j 典
␦ 共 q i ,q j 兲 , 共13兲
temperature for L⫽16. Both the simulational results com-
puted with Eqs. 共6兲–共9兲, and the exact solutions are plotted where q i ⫽1,2,...Q denotes the Potts spin at site i and
and overlap almost perfectly over a wide temperature region ␦ (q i ,q j ) is a Kronecker delta. During the simulation, we
from k B T⫽0 – 8. Because no difference is visible in these select lattice sites randomly and choose integers between
figures, more stringent tests of the accuracy are provided by 关 1,Q 兴 randomly for new Potts spin values. The modification

Fig. 3. Thermodynamic quantities for the L⫽16 2D Ising model calculated from the density of states. The relative errors with respect to the exact solutions
by Ferdinand and Fisher54 are shown in the insets: 共a兲 internal energy, 共b兲 specific heat, 共c兲 Helmholtz free energy, and 共d兲 entropy.

1298 Am. J. Phys., Vol. 72, No. 10, October 2004 D. P. Landau, Shan-Ho Tsai, and M. Exler 1298

Downloaded 26 Sep 2012 to 136.159.235.223. Redistribution subject to AAPT license or copyright; see http://ajp.aapt.org/authors/copyright_permission
Fig. 4. Logarithm of the density of states g(E) for the 2D Q⫽8 Potts model Fig. 6. Final histogram of energy for the last iteration of random walks to
as a function of energy per lattice site, E/N, for L⫽8, 12, and 16. With the estimate the density of states of a Q⫽8 Potts model for L⫽16.
scale in the figure, the errors of the simulational data are within the width of
the lines.
ference between the double peaks. When T is slightly away
from T c , one of the double peaks increases dramatically in
factor ln(f i) changes from ln(f 0)⫽1 at the beginning, to magnitude and the other decreases as shown in the inset of
ln(f final)⫽10⫺8 by the end of the random walks. The histo- Fig. 5.
Because of the double-peaked structure at a first-order
gram of energy H(E) is considered flat when all entries are
phase transition, conventional MC simulations are not effi-
not less than 80% of the average 具 H(E) 典 . To guarantee the cient because an extremely long time is required for the sys-
accuracy of thermodynamic quantities at low temperatures, tem to travel from one peak to the other in energy space.
we use the condition that the number of the ground states is With the Wang-Landau algorithm, all possible energy levels
Q⫽8 to normalize the density of states. The densities of are visited with equal probability, so it overcomes the barrier
states for L⫽8, 12, and 16 lattices are shown in Fig. 4. We between the coexisting phases in the conventional MC simu-
see that the maximum density of states from our data for L lations. The final flat histogram for L⫽16 is shown in Fig. 6,
⫽16 is very close to e 530, which is about 1.5⫻10230. and it describes the total number of visits to each energy
In Fig. 5 we show the double-peaked canonical probability level for the random walk of the last iteration.
distribution56 at the transition temperature T c for the first- Figure 7 illustrates some thermodynamic quantities calcu-
order transition, computed from the simulational data using lated from the density of states using Eqs. 共6兲–共9兲. Near the
Eq. 共1兲. The ‘‘transition temperature’’ k B T c (L) is approxi- transition temperature T c , the internal energy, shown in Fig.
mately 0.7519 for L⫽16 and is the temperature where the 7共a兲, has a steplike change that becomes sharper as the lattice
double peaks are of the same height. The transition tempera- size increases and transforms into a discontinuous jump
ture for the infinite lattice is known exactly to be k B T c when the system size goes to infinity. The magnitude of this
⫽1/ln(1⫹冑Q)⬇0.7449. 43 The valley between the two peaks jump 关shown in Fig. 7共a兲 for an infinite lattice兴 equals the
latent heat for the phase transition.
is approximately 0.37 for L⫽16, and becomes deeper as L
The specific heat, shown in Fig. 7共b兲, has a peak in the
increases. The latent heat for this temperature-driven first-
order phase transition can be estimated from the energy dif- vicinity of T c , and both the maximum value and the position
of the peak depend on the finite size of the lattice. As L
increases, the peak in the specific heat becomes narrower and
goes to a delta function in the thermodynamic limit.
Our results for the Helmholtz free energy per lattice site
are shown in Fig. 7共c兲 as a function of temperature. Because
the transition is of first-order, the first derivative of the free
energy has a discontinuity at T c . 共The location of this dis-
continuity can be used as an estimate of T c .)
Like the internal energy, the entropy shown in Fig. 7共d兲
has a steplike change near T c . This change becomes sharper
as L increases and becomes a discontinuous jump when L
→⬁. Because the jump in the internal energy equals the
latent heat of the phase transition, and the free energy is
continuous at T c , the magnitude of the jump in the entropy
equals the latent heat divided by T c 关see Eq. 共9兲兴.

V. RANDOM WALK IN ENERGY AND ORDER

Fig. 5. The canonical distribution at the transition temperature P(E,T c ) PARAMETER SPACE
⫽g(E)e ⫺E/k B T c for the Q⫽8 Potts model for L⫽16. The inset shows the
canonical distribution at a temperature slightly above and below T c for the To study the effect of an applied magnetic field on the
same system. Ising and Potts models, we have to perform a random walk

1299 Am. J. Phys., Vol. 72, No. 10, October 2004 D. P. Landau, Shan-Ho Tsai, and M. Exler 1299

Downloaded 26 Sep 2012 to 136.159.235.223. Redistribution subject to AAPT license or copyright; see http://ajp.aapt.org/authors/copyright_permission
Fig. 7. Thermodynamic quantities calculated from the density of states for the Q⫽8 Potts model for L⫽4, 8, and 16: 共a兲 internal energy, 共b兲 specific heat, 共c兲
Helmholtz free energy, and 共d兲 entropy.

in both the energy and order parameter space. 2D random near first-order phase transitions. This metastability hinders
walks may also be required to study systems with more com- studies of first-order phase transitions using standard MC
plex orders, such as a three-dimensional spin glass model, methods.
even in the absence of external fields.
To illustrate a simple case where a 2D random walk is VI. DISCUSSION AND CONCLUSION
required, we consider the 2D Ising model in the presence of
an external magnetic field h. The Hamiltonian is given by We have described an efficient algorithm to calculate the
N
density of states directly for large systems. By modifying the
estimate at each step of the random walk in energy space and
H⫽⫺ 兺
具 i, j 典
␴ i ␴ j ⫺h 兺 ␴i .
i⫽1
共14兲 carefully controlling the modification factor, we can deter-
mine the density of states very accurately. Using the density
The order parameter is the magnetization, defined as M ⬘ of states, we can then calculate thermodynamic quantities at
⫽ 兺 i⫽1
N
␴ i , and we denote the exchange energy as E ⬘ essentially any temperature. An important advantage of this
⫽⫺ 兺 具 i, j 典 ␴ i ␴ j . The algorithm works as before, except that approach is that we can also calculate the Helmholtz free
the random walk is now performed in both the energy E ⬘ and energy and entropy, quantities that are not directly available
the order parameter M ⬘ , and a 2D histogram H(E ⬘ ,M ⬘ ) is
accumulated.
With the estimate of the density of states g(E ⬘ ,M ⬘ ), the
partition function can be computed as

Z 共 T,h 兲 ⫽ 兺 g 共 E ⬘ ,M ⬘ 兲 e ⫺(E ⬘ ⫺hM ⬘ )/k B T . 共15兲

E ,M⬘ ⬘
From the partition function we can obtain thermodynamic
quantities for all values of the temperature and magnetic
field. For example, the mean magnetization of the system can
be computed as
兺 E ⬘ ,M ⬘ M ⬘ g 共 E ⬘ ,M ⬘ 兲 e ⫺(E ⬘ ⫺hM ⬘ )/k B T
M 共 T,h 兲 ⫽ . 共16兲
兺 E ⬘ ,M ⬘ g 共 E ⬘ ,M ⬘ 兲 e ⫺(E ⬘ ⫺hM ⬘ )/k B T
M (T,h) is shown in Fig. 8 as a function of the external
magnetic field h, for different values of the temperature T.
Note that for fixed T⬍T c the Ising model has a first-order Fig. 8. Magnetization per site for the L⫽16 Ising model as a function of
applied field for T⫽1.0/k B ⬍T c 共circle兲, T⫽2.3/k B ⬇T c 共upright triangle兲,
phase transition at h⫽0. Standard MC methods generate T⫽3.0/k B ⬎T c 共inverted triangle兲 obtained with Wang-Landau sampling,
hysteresis in the magnetization curve at low T 共shown as the and the hysteresis curve obtained with the Metropolis algorithm at T
solid line in Fig. 8兲, because of metastable states that appear ⫽1.0/k B 共solid line兲. The dotted line is a guide to the eye.

1300 Am. J. Phys., Vol. 72, No. 10, October 2004 D. P. Landau, Shan-Ho Tsai, and M. Exler 1300

Downloaded 26 Sep 2012 to 136.159.235.223. Redistribution subject to AAPT license or copyright; see http://ajp.aapt.org/authors/copyright_permission
 11
from conventional MC simulations. The method is applicable B. A. Berg and T. Celik, ‘‘New approach to spin-glass simulations,’’ Phys.
to a wide range of systems and is easy to implement. Al- Rev. Lett. 69, 2292–2295 共1992兲.
12
N. A. Alves and U. H. E. Hansmann, ‘‘Partition function zeros and finite
though we have described its implementation in terms of
size scaling of helix-coil transitions in a polypeptide,’’ Phys. Rev. Lett. 84,
single spin-flip sampling, it is straightforward to use other 1836 –1839 共2000兲.
types of sampling for cases in which it will further accelerate 13
A. M. Ferrenberg and R. H. Swendsen, ‘‘New Monte Carlo technique for
the simulation. studying phase transitions,’’ Phys. Rev. Lett. 61, 2635–2638 共1988兲; ‘‘Op-
Applications to the 2D Q⫽8 Potts model and to the 2D timized Monte Carlo data analysis,’’ ibid. 63, 1195–1198 共1989兲.
14
Ising model show that the method is effective for systems P. D. Beale, ‘‘Exact distribution of energies in the two-dimensional Ising
that exhibit first-order or second-order phase transitions. Our Model,’’ Phys. Rev. Lett. 76, 78 – 81 共1996兲.
15
J. Lee, ‘‘New Monte Carlo algorithm: Entropic sampling,’’ Phys. Rev.
presentation concentrated on the random walk in energy Lett. 71, 211–214 共1993兲.
space 共and order-parameter space兲; however, the idea is very 16
P. M. C. de Oliveira, T. J. P. Penna, and H. J. Herrmann, ‘‘Broad histogram
general and can be applied to any parameters. The energy method,’’ Braz. J. Phys. 26, 677– 683 共1996兲; ‘‘Broad histogram Monte
levels of the models treated here are discrete, and the total Carlo,’’ Eur. Phys. J. B 1, 205–208 共1998兲.
17
number of possible energies is known before the simulation, J.-S. Wang and L. W. Lee, ‘‘Monte Carlo algorithms based on the number
but in general such information is not available. For models of potential moves,’’ Comput. Phys. Commun. 127, 131–136 共2000兲.
18
where all the possible energy levels cannot be fitted in the A. R. Lima, P. M. C. de Oliveira, and T. J. P. Penna, ‘‘A comparison
between broad histogram and multicanonical methods,’’ J. Stat. Phys. 99,
computer memory or the energy is continuous, for example, 691–705 共2000兲.
the Heisenberg model, we must bin the energy. Statistical 19
F. Wang and D. P. Landau, ‘‘Efficient, multiple-range random walk algo-
and systematic errors in the density of states, and thus of the rithm to calculate the density of states,’’ Phys. Rev. Lett. 86, 2050–2053
thermodynamic quantities, are controlled by the flatness of 共2001兲; ‘‘Determining the density of states for classical statistical models:
the histogram at the end of each iteration and the final modi- A random walk algorithm to produce a flat histogram,’’ Phys. Rev. E 64,
fication factor f final . These errors can be decreased by requir- 056101 共2001兲.
20
C. Yamaguchi and Y. Okabe, ‘‘Three-dimensional antiferromagnetic
ing a more strict condition for a flat histogram and by using q-state Potts models: application of the Wang-Landau algorithm,’’ J. Phys.
a ln(f final) that is closer to zero. A 34, 8781– 8794 共2001兲.
In this paper, we only applied the Wang-Landau algorithm 21
Y. Okabe, Y. Tomita, and C. Yamaguchi, ‘‘Application of new Monte Carlo
to simple models on small lattices, but the method is also algorithms to random spin systems,’’ Comput. Phys. Commun. 146, 63– 68
efficient for large systems and has proven to be useful in the 共2002兲.
22
studies of general, complex systems with rough landscapes M. Troyer, S. Wessel, and F. Alet, ‘‘Flat histogram methods for quantum
systems: Algorithms to overcome tunneling problems and calculate the
共see references given in Sec. I兲. However, more investigation free energy,’’ Phys. Rev. Lett. 90, 120201 共2003兲.
is needed to better determine under which circumstances 23
P. N. Vorontsov-Velyaminov and A. P. Lyubartsev, ‘‘Entropic sampling in
the method offers substantial advantages over other the path integral Monte Carlo method,’’ J. Phys. A 36, 685– 693 共2003兲.
approaches.57 24
W. Koller, A. Prüll, H. G. Evertz, and W. von der Linden, ‘‘Uniform
hopping approach to the ferromagnetic Kondo model at finite tempera-
ture,’’ Phys. Rev. B 67, 104432 共2003兲.
ACKNOWLEDGMENTS 25
T. S. Jain and J. J. de Pablo, ‘‘Calculation of interfacial tension from
density of states,’’ J. Chem. Phys. 118, 4226 – 4229 共2003兲.
We would like to thank Laura Nurminen and Stefan Tor- 26
Q. L. Yan, R. Faller, and J. J. de Pablo, ‘‘Density-of-states Monte Carlo
brügge for their critical reading of the manuscript. The re- method for simulation of fluids,’’ J. Chem. Phys. 116, 8745– 8749 共2002兲.
search project was supported by the National Science Foun- 27
R. Faller and J. J. de Pablo, ‘‘Density of states of a binary Lennard-Jones
dation under Grant No. DMR-0094422. glass,’’ J. Chem. Phys. 119, 4405– 4408 共2003兲.
28
E. B. Kim, R. Faller, Q. Yan, N. L. Abbott, and J. J. de Pablo, ‘‘Potential
a兲
Author to whom correspondence should be addressed; electronic mail: of mean force between a spherical particle suspended in a nematic liquid
tsai@hal.physast.uga.edu crystal and a substrate,’’ J. Chem. Phys. 117, 7781–7787 共2002兲.
29
b兲
Present address: Fachbereich Physik, Universität Osnabrück, T. S. Jain and J. J. de Pablo, ‘‘A biased Monte Carlo technique for calcu-
D-49069 Osnabrück, Germany. lation of the density of states of polymer films,’’ J. Chem. Phys. 116,
1
D. P. Landau and K. Binder, A Guide to Monte Carlo Methods in Statis- 7238 –7243 共2002兲.
30
tical Physics 共Cambridge U. P., Cambridge, 2000兲. N. Rathore and J. J. de Pablo, ‘‘Monte Carlo simulation of proteins
2
D. P. Landau and R. Alben, ‘‘Monte Carlo calculations as an aid in teach- through a random walk in energy space,’’ J. Chem. Phys. 116, 7225–7230
ing statistical mechanics,’’ Am. J. Phys. 41, 394 – 400 共1973兲. 共2002兲.
3 31
N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. N. Rathore, T. A. Knotts, and J. J. de Pablo, ‘‘Density of states simulations
Teller, ‘‘Equation of state calculations by fast computing machines,’’ J. of proteins,’’ J. Chem. Phys. 118, 4285– 4290 共2003兲.
32
Chem. Phys. 21, 1087–1092 共1953兲. T. J. H. Vlugt, ‘‘Measurement of chemical potentials of systems with
4 strong excluded volume interactions by computing the density of states,’’
R. H. Swendsen and J.-S. Wang, ‘‘Nonuniversal critical dynamics in
Monte Carlo simulations,’’ Phys. Rev. Lett. 58, 86 – 88 共1987兲. Mol. Phys. 100, 2763–2771 共2002兲.
5 33
U. Wolff, ‘‘Collective Monte Carlo updating for spin systems,’’ Phys. Rev. F. Calvo, ‘‘Sampling along reaction coordinates with the Wang-Landau
Lett. 62, 361–364 共1989兲. method,’’ Mol. Phys. 100, 3421–3427 共2002兲.
6 34
B. A. Berg and T. Neuhaus, ‘‘Multicanonical ensemble: A new approach to F. Calvo and P. Parneix, ‘‘Statistical evaporation of rotating clusters. I.
simulate first-order phase transitions,’’ Phys. Rev. Lett. 68, 9–12 共1992兲. Kinetic energy released,’’ J. Chem. Phys. 119, 256 –264 共2003兲.
7 35
W. Janke and S. Kappler, ‘‘Multibondic cluster algorithm for Monte Carlo M. A. de Menezes and A. R. Lima, ‘‘Using entropy-based methods to
simulations of first-order phase transitions,’’ Phys. Rev. Lett. 74, 212–215 study general constrained parameter optimization problems,’’ Physica A
共1995兲. 323, 428 – 434 共2003兲.
8 36
W. Janke, ‘‘Multicanonical simulation of the two-dimensional 7-state Potts V. Mustonen and R. Rajesh, ‘‘Numerical estimation of the asymptotic
model,’’ Int. J. Mod. Phys. C 3, 1137–1146 共1992兲. behaviour of solid partitions of an integer,’’ J. Phys. A 36, 6651– 6659
9
B. A. Berg, U. Hansmann, and T. Neuhaus, ‘‘Simulation of an ensemble 共2003兲.
37
with varying magnetic field: A numerical determination of the order-order C. Yamaguchi and N. Kawashima, ‘‘Combination of improved multibon-
interface tension in the D⫽2 Ising model,’’ Phys. Rev. B 47, 497–500 dic method and the Wang-Landau method,’’ Phys. Rev. E 65, 056710
共1993兲. 共2002兲.
10 38
W. Janke, ‘‘Multicanonical Monte Carlo simulations,’’ Physica A 254, B. J. Schulz, K. Binder, and M. Müller, ‘‘Flat histogram method of Wang-
164 –178 共1998兲. Landau and N-fold way,’’ Int. J. Mod. Phys. C 13, 477– 494 共2002兲.

1301 Am. J. Phys., Vol. 72, No. 10, October 2004 D. P. Landau, Shan-Ho Tsai, and M. Exler 1301

Downloaded 26 Sep 2012 to 136.159.235.223. Redistribution subject to AAPT license or copyright; see http://ajp.aapt.org/authors/copyright_permission
39 50
A. Hüller and M. Pleimling, ‘‘Microcanonical determination of the order U. H. E. Hansmann and Y. Okamoto, ‘‘Monte Carlo simulations in gener-
parameter critical exponent,’’ Int. J. Mod. Phys. C 13, 947–956 共2002兲. alized ensemble: Multicanonical algorithm versus simulated tempering,’’
40
M. S. Shell, P. G. Debenedetti, and A. Z. Panagiotopoulos, ‘‘Generaliza- Phys. Rev. E 54, 5863–5865 共1996兲.
51
tion of the Wang-Landau method for off-lattice simulations,’’ Phys. Rev. E W. Janke, B. A. Berg, and A. Billoire, ‘‘Multi-overlap simulations of free-
66, 056703 共2002兲. energy barriers in the 3D Edwards-Anderson Ising spin glass,’’ Comput.
41
B. J. Schulz, K. Binder, M. Müller, and D. P. Landau, ‘‘Avoiding boundary Phys. Commun. 121Õ122, 176 –179 共1999兲.
effects in Wang-Landau sampling,’’ Phys. Rev. E 67, 067102 共2003兲. 52
J.-S. Wang, T. K. Tay, and R. H. Swendsen, ‘‘Transition matrix Monte
42
See, for example, L. Onsager, ‘‘Crystal statistics. I. A two-dimensional Carlo reweighting and dynamics,’’ Phys. Rev. Lett. 82, 476 – 479 共1999兲.
model with an order-disorder transition,’’ Phys. Rev. 65, 117–149 共1944兲. 53
D. P. Landau, ‘‘Finite-size behavior of the Ising square lattice,’’ Phys. Rev.
43
F. Y. Wu, ‘‘The Potts model,’’ Rev. Mod. Phys. 54, 235–268 共1982兲. B 13, 2997–3011 共1976兲.
44 54
B. A. Berg and W. Janke, ‘‘Multioverlap simulations of the 3D Edwards- A. E. Ferdinand and M. E. Fisher, ‘‘Bounded and inhomogeneous Ising
Anderson Ising spin glass,’’ Phys. Rev. Lett. 80, 4771– 4774 共1998兲. models. I. Specific-heat anomaly of a finite lattice,’’ Phys. Rev. 185, 832–
45
B. A. Berg, ‘‘Algorithmic aspects of multicanonical simulations,’’ Nucl. 846 共1969兲.
Phys. B 63, 982–984 共1998兲. 55
K. Binder, K. Vollmayr, H. P. Deutsch, J. D. Reger, M. Scheucher, and D.
46
B. A. Berg, T. Celik, and U. Hansmann, ‘‘Multicanonical study of the 3D P. Landau, ‘‘Monte Carlo methods for first-order phase transitions: some
Ising spin-glass,’’ Europhys. Lett. 22, 63– 68 共1993兲. recent progress,’’ Int. J. Mod. Phys. C 3, 1025–1058 共1992兲.
47 56
N. Hatano and J. E. Gubernatis, ‘‘Bivariate multicanonical Monte Carlo of M. S. S. Challa, D. P. Landau, and K. Binder, ‘‘Finite-size effects at
the 3D ⫾J spin glass,’’ in Computer Simulation Studies in Condensed temperature-driven first-order transitions,’’ Phys. Rev. B 34, 1841–1852
Matter Physics XII, edited by D. P. Landau, S. P. Lewis, and H.-B. Schüt- 共1986兲.
tler 共Springer, Berlin, 2000兲, pp. 149–161. 57
See EPAPS Document No. E-AJPIAS-72-006406 for a sample code of the
48
B. A. Berg and U. H. E. Hansmann, ‘‘Configuration space for random walk Wang-Landau algorithm for the 2D Ising model. This document may also
dynamics,’’ Eur. Phys. J. B 6, 395–398 共1998兲. be retrieved via the EPAPS homepage 共http://www.aip.org/pubservs/
49
U. H. E. Hansmann, ‘‘Effective way for determination of multicanonical epaps.html兲 or from ftp.aip.org in the directory /epaps. See the EPAPS
weights,’’ Phys. Rev. E 56, 6200– 6203 共1997兲. homepage for more information.

1302 Am. J. Phys., Vol. 72, No. 10, October 2004 D. P. Landau, Shan-Ho Tsai, and M. Exler 1302

Downloaded 26 Sep 2012 to 136.159.235.223. Redistribution subject to AAPT license or copyright; see http://ajp.aapt.org/authors/copyright_permission