The Wang and Landau algorithm proposed by Fugao Wang and David P. Landau^[1] is an extension of Metropolis Monte Carlo sampling. It 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. 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 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.

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^[2] and the folding of proteins.^[3][4]

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.^[5] 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 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.

Hence, we can use a simple harmonic oscillator potential to test the accuracy of Wang-Landau algorithm because we know already the analytic form of the density of states.

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.

References