A portable and flexible implementation of the Wang–Landau algorithm in order to determine the Density of States

The density of states (DOS), denoted by $\Omega(E)$ and defined as $\Omega(E)=\int d\mathbf{r}d\mathbf{p}\,\delta(E-H(\mathbf{r},\mathbf{p}))$, where $H(\mathbf{r},\mathbf{p})$ is the Hamiltonian of the system, is a fundamental function of a system such that $\Omega(E)dE$ is the number of states that have an energy between $E$ and $E+dE$. This function can be used to compute all of the thermodynamic properties of a system. Examples of such properties are the caloric curve and the heat capacity, both being highly relevant to the study of a wide variety of systems. Therefore, computing the DOS is of fundamental importance, and consequently several methods for the calculation of this function for a given system have been developed [1, 2, 3], some of them based on generalized ensembles and reweighting techniques.

Among the generalized-ensemble methods, Wang-Landau (WL) is a well known Monte Carlo technique for computing the DOS of a system [4]. This method was originally developed for discrete systems and having single processors in mind. However, since then the method has been expanded to off-lattice systems and parallel processing[5]. This has been done in order to take full advantage of the powerful multicore processing units that are the standard today. Although the WL algorithm is widely used and has a large number of variants and improvements [6, 7, 8], there are still practical difficulties in its implementation, ranging from the fine-tuning needed for an optimal operation to the lack of a universal implementation able to be “plugged into” traditional simulation codes.

In this work, we present a state-of-the-art implementation of WL using the flexibility and portability of the Python programming language [9], coupled with the efficiency of the C++ language. This implementation takes advantage of two well-known libraries: The powerful LAMMPS [10] library, which will be used to compute the energy of a complex configuration of atoms, and the easy-to-use Python multiprocessing library for parallel computing [11]. Furthermore, this implementation supports both discrete and continuous systems, while either supporting or being ready to support established techniques to improve the algorithm such as the optimization of flatness criteria [12], parallelization scheme by replica-exchange [5], and more efficient Monte Carlo trial moves [13] among others. Some improvements to the algorithm can be relatively easy to incorporate in the code, such as the WL-$1/t$ variant of Belardinelli et al. [12], which is already included in our code, and or the Stochastic Approximation in Monte Carlo (SAMC) of Liang et al. [14]. Here we are presenting a well-modularized system that allows the scientific community to use and incorporate different aspects of the WL algorithm efficiently.

The applicability of the WL algorithm is quite broad, covering different research areas such as multidimensional integration [15], protein folding [16], thermodynamics of polymers [17], free energy profiles [18], and the study of spin-glass models [19, 20, 21], among others. Here we focus the technique mainly on the calculation of the caloric-curve for different cases, to expand its use towards the condensed matter community. We tested this implementation using a variety of well-established benchmark systems such as the Ising model, the Potts model, classical and quantum harmonic oscillators, the Lennard-Jones model and the embedded atom model (EAM). The goal of these tests was to prove the soundness of our implementation through the comparison with previously known results.

This paper is organized as follows. Following this introduction, in Section 2 we describe the inner workings of the Wang–Landau algorithm. In Section 3 we offer a brief description of this implementation. In Section 4 we provide details of the benchmark simulations performed. A discussion of our findings is given in Section 5.

The WL algorithm is used to find the DOS of a many-particle system with $N$ degrees of freedom. This system can be described by a configurational state $X=(X_{1},X_{2},\ldots,X_{N})$ and a function of $X$ that describes the system energy $E=\mathcal{H}(X)$. This function is called the Hamiltonian.

In order to produce the DOS of a system, the Wang–Landau algorithm takes advantage of the following property. Suppose we sample a large number of possible states according to a distribution that depends only on the Hamiltonian of each state $P(X)=\rho(\mathcal{H}(X))$. If this distribution is chosen to be

then the energies will be distributed according to

Consequently, we would build up a flat histogram. Thereby, we can use this property to test how good a candidate DOS is, and even better, we can determine how to make a small change to it so that it is closer to the true DOS. We keep on making those small changes to the candidate DOS until the histogram is flat enough. The family of algorithms that takes advantage of this property are known as flat histogram methods [22].

The specific steps to reproduce the above procedure for discrete systems is detailed now. At the start of our simulation we define a candidate DOS of the system, $\Omega(E)=1$ for all $E$. We define a walker as an entity that wanders over the configuration space standing on an unique configuration at every simulation step. The energy of a walker is the energy of the configuration where this walker is currently standing on. We start out with a walker on a random state $X$, and then we produce a trial move from $X$ to $X^{\prime}$. The transition acceptance that is in concordance with Eq. 1 and the detailed balance condition is

If the trial move is accepted, the walker moves to the new state $X^{\prime}$; otherwise, it remains on its previous state $X$. Whether the change was accepted or not, we update the value of the DOS at the current energy $E$ according to $\Omega(E)\leftarrow\Omega(E)\cdot f$, where $f>1$. Based on the numeric nature of the algorithm, the values of these constants and their multiplication $\Omega(E)\leftarrow\Omega(E)\cdot f$ must be expressed in terms of their logarithms. Additionally, we update the histogram that counts the number of times certain energy has been visited. This process continues until the histogram is sufficiently flat. The flatness criterion is that the minimum value of the histogram must be at least 80% of its mean value. Once this criterion has been met, the process is interrupted, the value of $f$ is decreased according to $f\leftarrow\sqrt{f}$, the histogram is reset, and the same process is started over again. Finally, the algorithm is stopped when $f$ reaches a value close enough to unity, such as $f\approx 1+10^{-8}$, and the final $\Omega(E)$ is considered to be the true DOS.

For continuous systems, the energies are discretized into small bins of width $\Delta E$. Each bin has its own $\Omega$ and histogram, and they are updated when the energy of the walker lies between $E-\Delta E/2$ and $E+\Delta E/2$. The transition probability is computed after previously performing an interpolation on the discrete $\Omega$, thus obtaining $\Omega(E)$ from the interpolated function. This last step provides far better results than merely choosing the discrete $\Omega$ of the bin that contains $E$. This is because the walker will prefer a transition to a less probable energy, even if this new energy and the current one belong to the same bin. Furthermore, instead of the discrete case where the way in which we produce a change is obvious, we must select an appropriate step size, or alternatively, we can use a variable step size through the definition of a target acceptance. A target acceptance of 30% has been used in the continuous systems.

The algorithm is parallelized in the following way. Firstly, the entire energy range is split up into overlapping windows of energy (Fig. 1). An overlap of 75% has been reported to provide good results [5]. Secondly, every window is given its own walker and then a so-called “sweep” is executed—the performing of a certain number of simulation steps inside of every window. These steps run in parallel and they are completely independent. Thirdly, a replica exchange procedure [5] is tried. This procedure consists in the exchange of walkers between two windows. For this exchange to occur, the windows must be adjacent and their walkers must have an energy that lies over the overlapping interval of their respective windows. If the previous conditions are met for windows $i$ and $j$, the exchange is performed according to the probability

where $w$, $E$ and $\Omega$ are the walker, the walker’s current energy and the current DOS of their respective windows. This procedure allow walkers to travel through the entire energy landscape thus avoiding the presence of energy traps. Fourthly, we repeat the last two steps until every window has ended its own simulation. The windows that have already ended their simulations are required to continue running them without updating their DOS. This is required because we need that those windows keep on performing replica exchanges with the active ones. Finally, we reconstruct the full DOS from the pieces of DOS obtained at every window, as follows.

In order to reconstruct $\Omega$, we must consider that we want to produce the smoothest curve possible, in particular, we want the microcanonical caloric curve $\beta(E)=\partial(\log\Omega(E))/\partial E$ to be smooth. Having this in mind, we select the pieces of $\Omega$ that belong to the first ($\Omega_{1}$) and second ($\Omega_{2}$) windows, and then we join them through the following procedure. Firstly, we compute $\beta_{1}(E)=\partial(\log\Omega_{1}(E))/\partial E$ and $\beta_{2}(E)=\partial(\log\Omega_{2}(E))/\partial E$ and we select the energy $E^{*}$ where the best match between these two functions is. In our second step, we define $\bar{\Omega}_{2}=[\Omega_{1}(E^{*})/\Omega_{2}(E^{*})]\Omega_{2}$ so that $\Omega_{1}$ and $\bar{\Omega}_{2}$ match their values at $E^{*}$. Thirdly, we join $\Omega_{1}$ and $\bar{\Omega}_{2}$ at $E^{*}$, cutting out the leftover “tails” of each function. Finally, the above steps are repeated between the previous pasted $\Omega$ and the next piece of DOS until there are no more pieces left. By this previous procedure, we can now join the DOS from the energy windows into a DOS that spans the whole energy range.

Before discussing the implementation details, we need to address some early issues. Firstly, we need to generalize the definition of the DOS to include discrete systems. For discrete systems $\Omega(E)$ will be simply the number of states with exact energy $E$. Secondly, to avoid treating with very large numbers, the logarithm of the DOS ($\log\Omega(E)$) will be computed instead of the plain DOS [23, 24]. All equations will be transformed accordingly. Thirdly, after each sweep $\log\Omega(E)$ will be normalized such as that $\log\Omega(E_{\text{min}})=0$ or, equivalently, $\Omega(E)$ will be computed in units of $\Omega(E_{\text{min}})$.

This implementation is based on several modules. Firstly, we have a module that loads the simulation data such as the energy range, step size, target acceptance and so on from an input file. Secondly, there is the module that runs the actual simulations. This module spawns the subprocesses parallelization that will run the WL simulation on every window, performs the replica exchanges, merges the DOS pieces and writes the result into a file (Fig. 2). Finally, there are the walker modules, one for each system. These walkers keep track of the state of the system, perform changes, and compute the energy. Furthermore, those walkers are totally independent of the algorithm itself, so they can be used for any other algorithm and they can be coded in another language and wrapped into Python code. For example, a LAMMPS walker has been developed which contains a wrapper for the C++ LAMMPS library. Any other library can be used through a wrapper class. Because of this, this implementation is not affected by the lesser computational efficiency of an interpreted language such as Python, since the time consuming functions can be implemented in a faster language like C++ or Julia.

It is crucial to notice that in this current version, the parallelization is not based on the MPI library but in the Python multiprocessing module, which is designed for multi-core systems. The user has to look out for possible conflicts between this module and the parallelization of the walker in case it is an external process.

A set of different simulations were performed in this work in order to determine the accuracy and practical usefulness of our implementation. For every system, instead of computing the full DOS, we only computed the configurational density of states (CDOS). This means that instead of using the full Hamiltonian $\mathcal{H}=K+\phi$, we only used the configurational part $\phi$. This is done because the contribution of the kinetic part $K$ is always the same and it is easy to compute analytically.

The caloric curve for both microcanonical and canonical ensembles has been computed using the CDOS ($\Omega(\phi)$) obtained in the simulation. If $d$ is the number of degrees of freedom in the system, the microcanonical curve has been computed according to

where

is a function that is proportional to the full DOS [25], and the canonical curve has been computed according to the well known equation

where

is the configurational part of the partition function, and $d/2\beta$ is the kinetic energy term.

The results are presented below. The discrete models and the classical harmonic oscillator are presented as a proof of the soundness of our implementation. The Lennard–Jones and EAM solids are the most important results because they show how our code works seamlessly with the LAMMPS library which enables us to study a great quantity of systems with high precision.

The number of steps per second reported for each system corresponds to standard simulations performed using 88 walkers running on an Intel(R) Xeon(R) CPU E5-2699 v4 @ 2.20GHz with 88 virtual cores on 44 physical cores.

Our work presents a modular implementation of the original Wang–Landau algorithm enhanced with state-of-the-art developments such as off-lattice systems and parallelization. It also integrates seamlessly with the LAMMPS library for fast computing of the energy of atomic configurations, which means that we can, in principle, study any system that can be described by a potential supported by LAMMPS, or, in principle, any other atomistic simulation package for which a wrapper can be written. The results obtained in our simulation of well known systems matched what we expected, including the specially important cases of the Lennard-Jones system and the EAM system because these are the models that run with the LAMMPS library. Consequently, our implementation of the WL algorithm can be used from now on to efficiently study the properties of more realistic models of matter.

SD acknowledges partial support by ANID FONDECYT 1171127 and ANID PIA ACT172101 grants. FM thanks to Universidad Andrés Bello for providing the doctorate scholarship and acknowledges funding from ANID BECAS/DOCTORADO NACIONAL 21212267. Computational work was supported by the supercomputing infrastructures of the NLHPC (ECM-02), and FENIX (UNAB).