Wang and Landau algorithm

The Wang and Landau algorithm, proposed by Fugao Wang and David P. Landau, is a Monte Carlo method designed to calculate the density of states of a system. The method performs a non-markovian random walk to build the density of states by quickly visiting all the available energy spectrum. The Wang and Landau algorithm is an important method to obtain the density of states required to perform a multicanonical simulation.

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. The Wang-Landau sampling is related to the metadynamics algorithm.

Overview

The Wang and Landau algorithm is used to obtain the density of states of a system characterized by a cost function. It uses a non-markovian stochastic process which asymptotically converges to a multicanonical ensemble. (I.e. to a Metropolis-Hastings algorithm with sampling distribution inverse to the density of states.) The major consequence is that this sampling distribution leads to a simulation where the energy barriers are invisible. This means that the algorithm visits all the accessible states (favorable and less favorable) much faster than a metropolis algorithm.

Algorithm

Consider a system defined on a phase space $\Omega$, and a cost function, E, (e.g. the energy), bounded on a spectrum $E\in \Gamma =[E_{\min },E_{\max }]$, which has an associated density of states $\rho (E)\equiv \exp(S(E))$, which is to be computed. Because Wang and Landau algorithm works in discrete spectra, the spectrum $\Gamma$ is divided in N discrete values with a difference between them of $\Delta$, such that

$N={\frac {E_{\max }-E_{\min }}{\Delta }},$.

Given this discrete spectrum, the algorithm is initialized by:

• setting all entries of the microcanonical entropy to zero, $S(E_{i})=0\ \ i=1,2,...,N$
• initializing $f=1$ and
• initializing the system randomly, by putting in a random configuration ${\boldsymbol {r}}\in \Omega$.

The algorithm then performs a multicanonical ensemble simulation: a Metropolis-Hastings random walk in the phase space of the system with a probability distribution given by $P({\boldsymbol {r}})=1/\rho (E({\boldsymbol {r}}))=\exp(-S(E({\boldsymbol {r}})))$ and a probability of proposing a new state given by a probability distribution $g({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')$. A histogram $H(E)$ of visited energies is stored. Like in the Metropolis-Hastings algorithm, a proposal-acceptance step is performed, and consists in (see Metropolis–Hastings algorithm overview):

1. proposing a state ${\boldsymbol {r}}'\in \Omega$ according to the arbitrary proposal distribution $g({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')$
2. accept/refuse the proposed state according to

$A({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')=\min \left(1,e^{S-S'}{\frac {g({\boldsymbol {r}}'\rightarrow {\boldsymbol {r}})}{g({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')}}\right)$

where $S=S(E({\boldsymbol {r}}))$ and $S'=S(E({\boldsymbol {r}}'))$.

After each proposal-acceptance step, the system transits to some value $E_{i}$, $H(E_{i})$ is incremented by one and the following update is performed:

$S(E_{i})\leftarrow S(E_{i})+f$.

This is the crucial step of the algorithm, and it is what make this Wang and Landau algorithm non-markovian: the stochastic process now depends on the history of the process. Hence the next time there is a proposal to a state with that particular energy $E_{i}$, that proposal is now more likely refused; in this sense, the algorithm forces the system to visit all the spectrum equally. The consequence is that the histogram $H(E)$ is more and more flat. However, this flatness depends on how well approximated the calculated entropy is to the exact entropy, which naturally depends on the value of f. To better and better approximate the exact entropy (and thus histogram's flatness), f is decreased after M proposal-acceptance steps:

$f\leftarrow f/2$.

It was later shown that updating the f by constantly dividing by two can lead to saturation errors. A small modification to the Wang and Landau method to avoid this problem is to use the f factor proportional to $1/t$, where $t$ is proportional to the number of steps of the simulation.

Test system

We want to obtain the DOS for the harmonic oscillator potential.

$E(x)=x^{2},\,$

The analytical DOS is given by,

$g(E)=\int \delta (E(x)-E)\,dx=\int \delta (x^{2}-E)\,dx,$

by performing the last integral we obtain

$g(E)\propto {\begin{cases}E^{-1/2},{\text{for }}x\in \mathbb {R} ^{1}\\{\text{const}},{\text{for }}x\in \mathbb {R} ^{2}\\E^{1/2},{\text{for }}x\in \mathbb {R} ^{3}\\\end{cases}}$

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. Therefore we compare the density of states $\rho (E)$ obtained by the Wang–Landau algorithm with $g(E)$.

Sample code

The following is a sample code of the Wang–Landau algorithm in Python, where we assume that a symmetric proposal distribution g is used:

$g({\boldsymbol {x}}'\rightarrow {\boldsymbol {x}})=g({\boldsymbol {x}}\rightarrow {\boldsymbol {x}}')$

The code considers a "system" which is the underlying system being studied.

currentEnergy = system.randomConfiguration() # a random initial configuration

while (f > epsilon):
    system.proposeConfiguration() # a proposed configuration is proposed
    proposedEnergy = system.proposedEnergy() # the energy of the proposed configuration computed

    if (random() < exp(entropy[currentEnergy]-entropy[proposedEnergy])):
        # if accepted, update the energy and the system:
        currentEnergy = proposedEnergy
        system.acceptProposedConfiguration()
    else:
        # if rejected
        system.rejectProposedConfiguration()
    
    H[currentEnergy] += 1
    entropy[currentEnergy] += f
    
    if (isFlat(H)): # isFlat tests whether the histogram is flat (e.g. 95% flatness)
        H[:] = 0
        f *= 0.5 # refine the f parameter

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:

$T'(E)\rightarrow {\dfrac {\partial S(E)}{\partial E}}\cdot T(E),$

where $S(E)$ is the entropy of the system, $T(E)$ the micro-canonical temperature and $T'(E)$ is the "scaled" temperature used in the simulation.