Wang and Landau algorithm

The Wang and Landau algorithm is an extension of Metropolis Monte Carlo sampling.

Initially designed for physical systems, the Metropolis Monte Carlo method was based on Boltzmann's insight that at a given temperature molecules 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 melts. 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.

The Wang and Landau algorithm is designed to calculate the density of states of a computer-simulated system, such as an Ising model of spin glasses, or model atoms in a molecular force field.

Algorithm

Firstly the minimum and maximum possible states of the system are calculated. For instance for an Ising model the most favorable state is when all cells are spinning in the same direction, and the least favorable state is a checkerboard pattern. The range is then divided into a given number of discrete histogram entries.

Initially the density of states "g"("E") is unknown so all densities are set to unity. In fact densities typically range over such a huge interval that we use log space. A visits histogram "H"("E") is also maintained. Initially all bins have zero visits. The algorithm then runs a Monte Carlo-like simulation, filling the visits histogram. However instead of using the Metropolis criterion for acceptance, the probability of acceptance is given by the density of states. Call the density of states histogram "g". Moves are accepted if

:$$p < minleft{ 1, frac{g(E)}{g(E')} ight},

where "p" is a uniform random number on [0, 1), "E" is the energy of the current state, "E"&prime; the energy of the proposed move.

If the move is rejected the current histogram bin is incremented. The visits histogram is incremented by 1 on each move, the density of states histogram "g"("E") is multiplied by a constant factor, initially usually "e" = 2.72 (Euler's number). When the visits histogram, call it "H", is flat "g"("E") contains an accurate estimate of the density of states, within the limits of the modification factor.

The visits histogram "H"("E") is then reset to zero, and the modification factor reduced, typically to the square root of the previous factor, to produce a finer estimate of "g"("E").

The algorithm works because the move function is, by definition, sampling from the true density of states. So if acceptance by the reciprocal of the density of states produces and even histogram, the estimate must be accurate. In fact "H"("E") can never reach perfect flatness, so some criterion must be used &mdash; typically 10% difference or less between the highest and lowest entry.

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

ample 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 statesuint [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;

References

* F. Wang and D.P. Landau. Efficient Multiple Range Random Walk Algorithm to Calculate Density of States. "Phys Rev Lett" 86, 2050 (2001)