Wang and Landau algorithm

The Wang and Landau algorithm is used to obtain an estimate for 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.[1] (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.[6]

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

${\displaystyle 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, ${\displaystyle S(E_{i})=0\ \ i=1,2,...,N}$
• initializing ${\displaystyle f=1}$ and
• initializing the system randomly, by putting in a random configuration ${\displaystyle {\boldsymbol {r}}\in \Omega }$.

The algorithm then performs a multicanonical ensemble simulation:[1] a Metropolis–Hastings random walk in the phase space of the system with a probability distribution given by ${\displaystyle P({\boldsymbol {r}})=1/{\hat {\rho }}(E({\boldsymbol {r}}))=\exp(-S(E({\boldsymbol {r}})))}$ and a probability of proposing a new state given by a probability distribution ${\displaystyle g({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')}$. A histogram ${\displaystyle 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 ${\displaystyle {\boldsymbol {r}}'\in \Omega }$ according to the arbitrary proposal distribution ${\displaystyle g({\boldsymbol {r}}\rightarrow {\boldsymbol {r}}')}$
2. accept/refuse the proposed state according to

${\displaystyle 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 ${\displaystyle S=S(E({\boldsymbol {r}}))}$ and ${\displaystyle S'=S(E({\boldsymbol {r}}'))}$.

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

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

This is the crucial step of the algorithm, and it is what makes the 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 ${\displaystyle E_{i}}$, that proposal is now more likely refused; in this sense, the algorithm forces the system to visit all of the spectrum equally.[1] The consequence is that the histogram ${\displaystyle 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.[7] To better and better approximate the exact entropy (and thus histogram's flatness), f is decreased after M proposal-acceptance steps:

${\displaystyle f\leftarrow f/2}$.

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

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

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

The analytical DOS is given by,

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

by performing the last integral we obtain

${\displaystyle \rho (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 estimated density of states ${\displaystyle {\hat {\rho }}(E)}$ obtained by the Wang–Landau algorithm with ${\displaystyle \rho (E)}$.

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

${\displaystyle 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

Molecular dynamics (MD) is usually preferable to Monte Carlo (MC), so it is desirable to have a MD algorithm incorporating the basic WL idea for flat energy sampling. That algorithm is Statistical Temperature Molecular Dynamics (STMD), developed [8] by Jaegil Kim et al at Boston University.

An essential first step was made with the Statistical Temperature Monte Carlo (STMC) algorithm. WLMC requires an extensive increase in the number of energy bins with system size, caused by working directly with the density of states. STMC is centered on an intensive quantity, the statistical temperature, ${\displaystyle T(E)=1/(dS(E)/dE)}$, where E is the potential energy. When combined with the relation, ${\displaystyle \Omega (E)=e^{S(E)}}$, where we set ${\displaystyle k_{B}=1}$, the WL rule for updating the density of states gives the rule for updating the discretized statistical temperature,

${\displaystyle {\tilde {T}}_{j\pm 1}^{\prime }=\alpha _{j\pm 1}{\tilde {T}}_{j\pm 1},}$

where where ${\displaystyle \alpha _{j\pm 1}=1/(1\mp \delta f{\tilde {T}}_{j\pm 1}),\delta f=(lnf/2\Delta E),\Delta E}$ is the energy bin size, and ${\displaystyle {\tilde {T}}}$ denotes the running estimate. We define f as in,[1] a factor >1 that multiplies the estimate of the DOS for the i'th energy bin when the system visits an energy in that bin.

The details are given in Ref.[8] With an initial guess for ${\displaystyle T(E)}$ and the range restricted to lie between ${\displaystyle T_{L}}$ and ${\displaystyle T_{U}}$, the simulation proceeds as in WLMC, with significant numerical differences. An interpolation of ${\displaystyle {\tilde {T}}(E)}$ gives a continuum expression of the estimated ${\displaystyle S(E)}$ upon integration of its inverse, allowing the use of larger energy bins than in WL. Different values of ${\displaystyle S(E)}$ are available within the same energy bin when evaluating the acceptance probability. When histogram fluctuations are less than 20% of the mean, ${\displaystyle f}$ is reduced according to ${\displaystyle f\rightarrow {\sqrt {f}}}$.

STMC was compared with WL for the Ising model and the Lennard-Jones liquid. Upon increasing energy bin size, STMC gets the same results over a considerable range, while the performance of WL deteriorates rapidly. STMD can use smaller initial values of ${\displaystyle f_{d}=f-1}$ for more rapid convergence. In sum, STMC needs fewer steps to obtain the same quality of results.

Now consider the main result, STMD. It is based on the observation that in a standard MD simulation at temperature ${\displaystyle T_{0}}$ with forces derived from the potential energy ${\displaystyle E([x])}$, where ${\displaystyle [x]}$ denotes all the positions, the sampling weight for a configuration is ${\displaystyle e^{-E([x])/T_{0}}}$. Furthermore, if the forces are derived from a function ${\displaystyle W(E)}$, the sampling weight is ${\displaystyle e^{-W(E([x]))/T_{0}}}$.

For flat energy sampling, let the effective potential be ${\displaystyle T_{0}S(E)}$ - entropic molecular dynamics. Then the weight is ${\displaystyle e^{-S(E)}}$. Since the density of states is ${\displaystyle e^{+S(E)}}$, their product gives flat energy sampling.

The forces are calculated as

${\displaystyle F=(-d/dx)T_{0}S(E)=T_{0}(dS/dE)(-d/dx)E([x])=(T_{0}/T(E))F^{0}}$

where ${\displaystyle F^{0}}$ denotes the usual force derived from the potential energy. Scaling the usual forces by the factor ${\displaystyle (T_{0}/T(E))}$ produces flat energy sampling.

STMD starts with an ordinary MD algorithm at constant ${\displaystyle T_{0}}$ and V. The forces are scaled as indicated, and the statistical temperature is updated every time step, using the same procedure as in STMC. As the simulation converges to flat energy sampling, the running estimate ${\displaystyle {\tilde {T}}(E)}$ converges to the true ${\displaystyle T(E)}$. Technical details including steps to speed convergence are described in [8] and.[9]

In STMD ${\displaystyle T_{0}}$ is called the kinetic temperature as it controls the velocities as usual, but does not enter the configurational sampling, which is unusual. Thus STMD can probe low energies with fast particles. Any canonical average can be calculated with reweighting, but the statistical temperature, ${\displaystyle T(E)}$, is immediately available with no additional analysis. It is extremely valuable for studying phase transitions. In finite nanosystems ${\displaystyle T(E)}$ has a feature corresponding to every “subphase transition”. For a sufficiently strong transition, an equal-area construction on an S-loop in ${\displaystyle 1/T(E)}$ gives the transition temperature.

STMD has been refined by the BU group,[9] and applied to several systems by them and others. It was recognized by D. Stelter that despite our emphasis on working with intensive quantities, ${\displaystyle ln(f)}$ is extensive. However ${\displaystyle \delta f=(ln(f)/2\Delta E)}$ is intensive, and the procedure ${\displaystyle f\rightarrow {\sqrt {f}}}$ based on histogram flatness is replaced by cutting ${\displaystyle \delta f}$ in half every fixed number of time steps. This simple change makes STMD entirely intensive and substantially improves performance for large systems.[9] Furthermore, the final value of the intensive ${\displaystyle \delta f}$ is a constant that determines the magnitude of error in the converged ${\displaystyle T(E)}$, and is independent of system size. STMD is implemented in LAMMPS as fix stmd.

STMD is particularly useful for phase transitions. Equilibrium information is impossible to obtain with a canonical simulation, as supercooling or superheating is necessary to cause the transition. However an STMD run obtains flat energy sampling with a natural progression of heating and cooling, without getting trapped in the low energy or high energy state. Most recently it has been applied to the fluid/gel transition [9] in lipid-wrapped nanoparticles.

Replica exchange STMD [10] has also been presented by the BU group.