The Wang and Landau algorithm proposed by Fugao Wang and David P. Landau

David P. Landau

Not to be confused with Lev Landau.David P. Landau is Distinguished Research Professor of Physics and founding Director of the Center for Simulational Physics at the University of Georgia. He is a Fellow of the American Physical Society. He won the Aneesur Rahman Prize for Computational Physics...

 is an extension of Metropolis Monte Carlo sampling. It is designed to calculate the density of states

Density of states

In solid-state and condensed matter physics, the density of states of a system describes the number of states per interval of energy at each energy level that are available to be occupied by electrons. Unlike isolated systems, like atoms or molecules in gas phase, the density distributions are not...

 of a computer-simulated system, such as an Ising model

Ising model

The Ising model is a mathematical model of ferromagnetism in statistical mechanics. The model consists of discrete variables called spins that can be in one of two states . The spins are arranged in a graph , and each spin interacts with its nearest neighbors...

 of spin glasses, or model atoms in a molecular force field. The idea of the method is to build the density of states using the fact that a system, when being simulated, is generating its density of states. It is the general method to obtain the DOS in a multicanonical simulation.

Initially designed for physical systems, the Metropolis Monte Carlo method was based on Boltzmann's insight that at a given temperature molecule

Molecule

A molecule is an electrically neutral group of at least two atoms held together by covalent chemical bonds. Molecules are distinguished from ions by their electrical charge...

s are distributed between high energy, or unfavorable states, and low energy, or favorable states, with a probability given by both the energy difference and the density of states. Thus if there is an extremely large number of less-favored states, the less-favored states will dominate, if the temperature is sufficiently high to overcome the energy difference. This is seen, for instance, when wax

Wax

thumb|right|[[Cetyl palmitate]], a typical wax ester.Wax refers to a class of chemical compounds that are plastic near ambient temperatures. Characteristically, they melt above 45 °C to give a low viscosity liquid. Waxes are insoluble in water but soluble in organic, nonpolar solvents...

 melts

Melting point

The melting point of a solid is the temperature at which it changes state from solid to liquid. At the melting point the solid and liquid phase exist in equilibrium. The melting point of a substance depends on pressure and is usually specified at standard atmospheric pressure...

. At a low temperature only the favorable, solid states are reachable. As temperature rises, a huge number of less favored states becomes possible, and the solid turns into a liquid.

In principle the Wang-Landau algorithm can be applied to any system which is characterized by a cost (or energy) function. For instance,
it has been applied to the solution of numerical integrals
and the folding of proteins.

Algorithm

Firstly, unless dynamical allocation schemes are employed, the minimum
and maximum possible states of the system need to be calculated. For
instance, for the Ising model in the most favorable state all spins
point in the same direction, and in the least favorable state they
assume a checkerboard pattern. The range is then divided into a given
number of discrete histogram entries. If the maximum and minimum values
of the cost function are and , respectively, and a precision is
required, then the number of bins would be


and an array of entries is needed to store the values of
the density of states (DOS). Care needs to be taken with choice of since
the simulation time increases very fast with .

Initially, the DOS is unknown, so all bins in the
array are set to unity. Since densities typically range over multiple
orders of magnitude, it is common to store the logarithm of the DOS,
i.e., , as discussed below. In addition, a visits
histogram is maintained. Initially, all bins have zero
visits for both and . The bins
are then filled over the course of a Monte-Carlo-like simulation, in
which the probability of acceptance is given by the density of states
instead of the usual Metropolis criterion. Moves are accepted if


where is a uniform random number on [0, 1),
is the energy of the current state, and
the energy of the proposed move. After the move is accepted or rejected, the visits histogram is
incremented by one and the density of states histogram
is multiplied by a constant factor,


where usually the initial choice is
(Euler's number). As soon as the visits histogram
is flat, contains an accurate estimate of the
density of states within the limits of the modification factor. At this
step the visits histogram is set to zero and the modification factor is reduced as


A simple repetition of the above simulation scheme suffers from the
shortcoming that very large entries need to be stored in
. As mentioned before, in order to avoid this problem,
the quantity


is evaluated. Now, the previous choice of turns out
advantageous since . The modification factor is
then updated as


The algorithm works because the move function is, by definition,
sampling from the true density of states. Therefore, if acceptance by
the reciprocal of the density of states produces an even histogram, the
estimate must be accurate. In fact, can never reach
perfect flatness, so some criterion must be used—for instance,
less than 10% difference between the highest and lowest entry.

Updating the modification factor by the root square can lead to
saturation errors. A method to avoid this problem is to use a modification factor
proportional to , where is the simulation
time. In simple models such as the Ising model, this represents a minor
challenge, but it is a serious problem in more complex systems such as
proteins, membranes, and high-dimensional integrals because of the
computing time.

The density of states is independent of temperature. However, it can be
used to determine what the state of the system will be in for any given
temperature.

Test system

We want to obtain the DOS for the harmonic oscillator

Harmonic oscillator

In classical mechanics, a harmonic oscillator is a system that, when displaced from its equilibrium position, experiences a restoring force, F, proportional to the displacement, x: \vec F = -k \vec x \, where k is a positive constant....

 potential.


The analytical DOS is given by,


by performing the last integral we obtain,


in general, the DOS for a multidimensional harmonic oscillator will be given by some
power of E, the exponent will be a function of the dimension of the system.

Sample code

The following is an (unoptimized) sample code of the Wang-Landau algorithm in D programming language:

real[real] log_g; // this associative array stores the log of the density of states
uint[real] H; // this stores the histogram

real E_old = the_system.energy; // this stores the old energy.
real F = 1; // the modification factor.

while(F > epsilon) {
the_system.evolve; // propose a new configuration of the system.
real E_new = the_system.energy; // compute the current system energy.

if (log_g[E_new] < log_g[E_old]) // check if the proposal is accepted.
E_old = E_new;
else if (random < exp(log_g[E_old] - log_g[E_new]))
E_old = E_new;
else
the_system.reject_and_undo;

H[E_old]++; // update the histogram and density of states.
log_g[E_old] += F;

if (is_flat(H)) { // when the histogram is flat,
H[] = 0; // reset it and reduce the modification factor.
F *= 0.5;
}
}

return log_g;

Wang and Landau Molecular Dynamics

It should be noted that the Wang and Landau algorithm can be implemented not only in a Monte Carlo simulation but also in a molecular dynamics simulation. To do this would require an escalation of the temperature of the system as follows:


where is the entropy of the system, the micro-canonical temperature and
is the "scaled" temperature used in the simulation.