Control of accuracy in the Wang-Landau algorithm

TL;DR

This paper analyzes the Wang-Landau algorithm, revealing how the eigenvalues of its transition matrix can be used to control the accuracy of density of states estimation, supported by analytic and numerical results.

Contribution

It introduces a method to control the accuracy of the Wang-Landau algorithm based on eigenvalue analysis, with analytic expressions and numerical validation.

Findings

Eigenvalue difference indicates estimation accuracy

Analytic expressions derived for 1D Ising model

Numerical validation on Ising and Potts models

Abstract

The Wang-Landau (WL) algorithm has been widely used for simulations in many areas of physics. Our analysis of the WL algorithm explains its properties and shows that the difference of the largest eigenvalue of the transition matrix in the energy space from unity can be used to control the accuracy of estimating the density of states. Analytic expressions for the matrix elements are given in the case of the one-dimensional Ising model. The proposed method is further confirmed by numerical results for the one-dimensional and two-dimensional Ising models and also the two-dimensional Potts model.

Equations

$$g(E_k)P(E_k,E_m)=g(E_m)P(E_m,E_k),$$

$$T(E_k,E_m)=\min \left(1,\frac{g(E_k)}{g(E_m)}\right)P(E_k,E_m),$$

$$Z=\sum _{k=1}^{N_E}g(E_k)e^{-E_k/k_{B}T},$$

$$\tilde{T}(E_k,E_m)=\min \left(1,\frac{\tilde{g}(E_k)}{\tilde{g}(E_m)}\right)P(E_k,E_m),$$

$$T(E_k,E_m)=\min \left(1,\frac{g(E_k)}{g(E_m)}\right)\sum _{i=0}^{2k}\frac{N_i{Q}_i^{E_k\to E_m}}{g(E_k)},$$

$$Z_L=2\sum _{k=0}^{L/2}{C}_L^{2k}{e}^{(L-4k)/(k_{B}T)}.$$

$$T(E_k,E_{k+1})=T(E_{k+1},E_k)=\frac{C_{L-2}^{2k}}{\max \left(C_L^{2k},C_L^{2k+2}\right)}.$$

$$\Delta =\frac{1}{N_E}\sum _E\left|\frac{\tilde{S}(E,t)-S_{\text{exact}}(E)}{S_{\text{exact}}(E)}\right|$$

$$\frac{\tilde{T}(E_k,E_m)}{\tilde{T}(E_m,E_k)}=\frac{\tilde{g}(E_k)}{\tilde{g}(E_m)}\frac{P(E_k,E_m)}{P(E_m,E_k)},$$

$$\frac{\tilde{T}(E_k,E_m)}{\tilde{T}(E_m,E_k)}=\frac{{\eta }_m}{{\eta }_k},$$

$${\eta }_m={\eta }_m\sum _k\tilde{T}(E_m,E_k)=\sum _k\tilde{T}(E_k,E_m){\eta }_k.$$

$$N_i$$

$$N_{2k}$$

$$Q_i^{E_k\to E_{k-1}}$$

$$Q_i^{E_k\to E_{k+1}}$$

Full text

I Introduction

The Wang–Landau (WL) algorithm Wang-Landau ; Wang-Landau-PRE has been shown to be a very powerful tool for directly determining the density of states (DOS) and is also quite widely applicable. It overcomes some difficulties existing in other Monte Carlo algorithms (such as critical slowing down) and allows calculating thermodynamic observables, including free energy, over a wide temperature range in a single simulation.

A number of papers investigated statistical errors of the DOS estimation, and it was found in Yan2003 that errors reach an asymptotic value beyond which additional calculations fail to improve the accuracy of the results. Yet it was established in Zhou2005 ; Lee2006 that the statistical error scales as the square root of the logarithm of the modification factor, if the factor is kept constant.

It follows from the results in Yan2003 that there is a systematic error of DOS estimation by the WL algorithm 111In fact, in the very early presentation of the algorithm in the Rahman Prize Lecture in 2002, David Landau already mentioned the systematic error in the DOS estimation.. It was also confirmed in the case of the two-dimensional Ising model that the deviation of the DOS obtained with the WL algorithm from the exact DOS does not tend to zero 1overt ; 1overt-a . Several improvements of the behavior of the modification factor in the algorithm, which were shown to overcome the problem of systematic error in selected applications, have been suggested 1overt ; 1overt-a ; SAMC ; SAMC2 ; Eisenbach2011 .

There are about fifteen hundred papers that apply the WL algorithm and its improvements to particular problems (e.g., to the statistics of polymers Binder09 ; Ivanov2016 and to the diluted systems Malakis04 ; Fytas2013 , among many others).

In this paper, we address the question of the accuracy of the DOS estimation. We report a method for obtaining information on both the convergence of simulations and the accuracy of the DOS estimation. We numerically apply our algorithm to the one-dimensional and the two-dimensional Ising models, where the exact DOS is known Beale , and to the two-dimensional 8-state Potts model, which undergoes a first-order phase transition. We also present analytic expressions for the transition matrix in the energy spectrum for the one-dimensional Ising model.

Our approach is based on introducing the transition matrix in the energy space (TMES), whose elements show the frequency of transitions between energy levels during the WL random walk in the energy space. Its elements are influenced by both the random process of choosing a new configurational state and the WL probability of accepting the new state.

We consider a chain of random updates (e.g., flips of randomly chosen spins for the Ising model) of a system configuration. Each of the updates is accepted with unitary probability. This random walk in the configurational space is a Markov chain. Its invariant distribution is uniform, i.e., the probabilities of all states of the physical system are equal to each other. For any pair $\Omega_{A}$ and $\Omega_{B}$ of configurations, the probability of an update from $\Omega_{A}$ to $\Omega_{B}$ is equal to the probability of an update from $\Omega_{B}$ to $\Omega_{A}$. Hence, the detailed balance condition is satisfied. Therefore,

[TABLE]

where $g(E)$ is the true DOS and $P(E_{k},E_{m})$ is a probability of one step of the random walk to move from a configuration with the energy $E_{k}$ to any configuration with the energy $E_{m}$. We introduce the notation

[TABLE]

which represents nondiagonal elements of the TMES of the WL random walk on the true DOS. Relation (1) can be rewritten as $T(E_{k},E_{m})=T(E_{m},E_{k})$. Therefore, the TMES of the WL random walk on the true DOS is a symmetric matrix. Because the matrix is both symmetric and right stochastic, it is also left stochastic. This means that the rates of visiting of all energy levels are equal to each other.

In simulations with a reasonable modification of the WL algorithm, the systematic error of determining the DOS can be made arbitrarily small. In this case, we find that the computed TMES approaches a stochastic matrix as the computed DOS approaches the true value. There are several interesting conclusions. First, this explains the criterion of histogram flatness, which is one of the main features of the original WL algorithm Wang-Landau . Because the histogram elements are equal to sums of columns in the TMES, histogram flatness is related to the closeness of the TMES to a stochastic matrix. Second, it gives a criterion for the proximity of the simulated DOS to the true value. We introduce the difference of the largest eigenvalue of the calculated TMES from unity as a parameter. We show that the parameter is closely connected with the deviation of the DOS from the true value. We confirm numerically that the deviation of the DOS from the true value decays in time in the same manner as our parameter decays.

We are not aware of any other method for determining the accuracy of a WL simulation without knowing the exact value of the DOS.

The paper is organized as follows. In Sec. II we describe the variants of the WL algorithm. In Sec. III we introduce the TMES and, in particular, we describe the behavior of the TMES for the one-dimensional Ising model. In Sec. IV we present our main results and discussion, including discussion of properties of the TMES, description of the method and numerical results for the one-dimensional and two-dimensional Ising models and for the two-dimensional Potts model.

II The algorithms

Directly estimating the DOS with the WL algorithm allows calculating the free energy as the logarithm of the partition function

[TABLE]

where $g(E_{k})$ is the number of states (density of states) with the energy $E_{k}$, $N_{E}$ is the number of energy levels, $k_{B}$ is the Boltzmann constant, and $T$ is the temperature.

The main idea of the WL algorithm is to organize a random walk in the energy space. We take a configuration of the system with the energy $E_{k}$, randomly choose an update to a new configuration with the energy $E_{m}$, and accept this configuration with the WL probability $\min\left(1,\tilde{g}(E_{k})/\tilde{g}(E_{m})\right)$, where $\tilde{g}(E)$ is the DOS approximation. The approximation is obtained recursively by multiplying $\tilde{g}(E_{m})$ by a factor $f$ at each step of the random walk in the energy space 222If the new configuration is not accepted, then the configuration is left unchanged, and the step is counted as the move to the energy $E_{k}$, i.e., $\tilde{g}(E_{k})$ is multiplied by the factor $f$.. Each time that the auxiliary histogram $H(E)$ becomes sufficiently flat, the parameter $f$ is modified by taking the square root, $f:=\sqrt{f}$. Each histogram value $H(E_{m})$ contains the number of moves to the energy level $E_{m}$. The histogram is filled with zeros after each modification of the refinement parameter $f$. It is convenient to work with the logarithms of the values $S(E_{k}):=\ln\tilde{g}(E_{k})$ and $F:=\ln f$ (to fit the large numbers into double precision variables) and to replace the multiplication $\tilde{g}(E_{m}):=f\cdot\tilde{g}(E_{m})$ with the addition $S(E_{m}):=S(E_{m})+F$.

At the end of the simulation, the algorithm provides only a relative DOS. Either the total number of states or the number of ground states can be used to determine the normalized DOS.

It is natural to ask the following three questions:

1. Q1

Which condition for the flatness check is optimal? 2. Q2

How does the histogram flatness influence the convergence of the DOS estimation? 3. Q3

Is the choice of the square root rule to modify the parameter $f$ optimal?

A practical answer to question Q1 was given in the original algorithm Wang-Landau : keep the flatness within the accuracy of about 20%. Choosing an accuracy between 1% and 20% is sometimes useful Wust2012 but can result in a substantial increase of the simulation time Wang-Landau-PRE . An answer to question Q3 was obtained in two independent works 1overt and SAMC , which introduced modifications of the WL algorithm, the WL-1/t algorithm and the stochastic approximation Monte Carlo (SAMC) algorithm, respectively.

There are two phases of the WL-1/t algorithm 1overt . The first phase is similar to the WL algorithm except that every test of the histogram flatness is replaced with a simpler check: Is $H(E)\neq 0$ for all $E$? The algorithm enters its second phase if $F\leq N_{E}/t$, where $t$ is the simulation time measured as the number of attempted spin flips. For $t>t_{s}$, the histogram is no longer checked and $F$ is updated as $F=N_{E}/t$ at each step. Here $t_{s}$ is the simulation time when the WL-1/t algorithm enters the second phase.

Both modified WL algorithms exhibit the same long-range behavior of the refinement parameter $F$ proportional to $1/t$ for long simulation times SAMC ; SAMC2 . This is natural due to the following conditions of the convergence: $\sum_{t=1}^{\infty}F(t)=\infty$ and $\sum_{t=1}^{\infty}F(t)^{\zeta}<\infty$ for some $\zeta\in(1,2)$ SAMC ; SAMC2 . The SAMC algorithm has an additional parameter $t_{0}$, which is the simulation time when the algorithm enters its second phase. Obtaining the appropriate value of $t_{0}$ can be quite cumbersome because the rule of thumb for choosing $t_{0}$ given in SAMC is violated even by the $128\times 128$ Ising model Janke2017 . The WL-1/t algorithm and its further improvements Zhou2008 ; Swetnam2011 ; Wust2009 seem to perform more reliably. Here, we use the WL-1/t algorithm, although the main obtained results are qualitatively independent of the modification choice.

III Transition matrix in the energy space

We calculate the TMES for the WL random walk as follows. The elements of the TMES $\tilde{T}(E_{k},E_{m})$ are probabilities for the WL random walk to move from a configuration with the energy $E_{k}$ to a configuration with the energy $E_{m}$. For simplicity, we consider the case of the Ising model with periodic boundary conditions and the energy $E=-\sum_{<i,j>}\sigma_{i}\sigma_{j}$, where the sum ranges pairs of neighboring spins and $\sigma_{i}=\pm 1$. The number of energy levels accessible for the WL random walk is $N_{E}=L/2+1$ for $d=1$ and $N_{E}=L^{2}-1$ for $d=2$, where the even integer $L$ is the linear size of the hypercubic lattice and $d$ is the lattice dimension. A WL random move cannot increase or decrease the energy of the configuration by more than $d$ energy levels, and every column and every row of the TMES therefore contains no more than $1{+}2d$ nonzero elements. The nondiagonal elements of $\tilde{T}(E_{k},E_{m})$ can be represented as

[TABLE]

where $k\neq m$. In general, the structure of the probability $P(E_{k},E_{m})$ depends on both the system dimension and the local lattice properties and is rather complicated.

In the case of the one-dimensional Ising chain of $L$ spins with periodic boundary conditions, the probability to change energy from $E_{k}$ to $E_{m}$ in a WL random move is

[TABLE]

where $k\neq m$. Here $k$ is the number of couples of domains walls in the configuration, which determines the energy level $E_{k}=-\sum_{j=1}^{L}\sigma_{j}\sigma_{j+1}=-L+4k$, $N_{i}(k,L)$ is the number of configurations where $i$ domains consist of only one spin and $2k{-}i$ domains consist of more than one spin, and $Q_{i}^{E_{k}\to E_{m}}(L)$ is the probability that a single spin flip moves the system to the energy $E_{m}$ from such configurations. Occupations of the energy levels of the chain are expressed in terms of binomial coefficients as $g(E_{k})=2C_{L}^{2k}$ because there are exactly $C_{L}^{2k}$ ways to arrange the $2k$ domain walls. Therefore, partition function (3) is

[TABLE]

The detailed analytic expressions for $N_{i}$ and $Q_{i}$ are presented in Appendix B. It follows that

[TABLE]

Equation (7) can be understood as follows. The probability of the system to change energy from $E_{k}$ to $E_{k+1}$ due to a spin flip is equal to the probability that there are no domain walls adjacent to the spin. Therefore, $P(E_{k},E_{k+1})=C^{2k}_{L-2}/C^{2k}_{L}$. Similarly, $P(E_{k+1},E_{k})=C^{2k}_{L-2}/C^{2k+2}_{L}$. We hence obtain (7).

IV Results and discussion

IV.1 TMES and the accuracy of the DOS estimation

The convergence of the WL-1/t algorithm follows from the arguments presented in Zhou2008 . Therefore, there is a final stage of each simulation, where the normalized DOS remains almost the same and is close to the limiting one.

We note that the condition that $F(t)$ is much smaller than one in itself does not guarantee that the algorithm is already in its final stage, because it follows from $\sum_{t=1}^{\infty}F(t)=\infty$ that a substantial cumulative change of the DOS due to a long simulation time is possible. At the same time, a large value of $F(t)$, resulting in a rapid increase of the calculated DOS, does not guarantee a rapid increase of the normalized DOS.

The normalized DOS remains almost the same during a long simulation time of the final stage. Therefore, the rate of increase of the logarithm of the nonnormalized DOS is nearly the same for all energies. The behavior of the algorithm is close to a Markov chain in the final stage, and the TMES remains almost the same. The invariant distribution of the Markov chain has the property that all energy levels are almost equiprobable, while different configurations having the same energy may have different probabilities. Therefore, the TMES is close to a stochastic matrix in the final simulation stage. The following proposition also holds: if the TMES is close to a stochastic matrix, then the obtained normalized DOS is close to the true DOS (see details in Appendix A).

The first phase of the WL-1/t algorithm aims to obtain the first crude approximation for the DOS, while the aim of the second phase (in which the factor $F$ is updated as $F(t)=N_{E}/t$ at each step) is to converge to the true DOS. Both the histogram flatness test in the original WL algorithm and the test whether all energies have been visited in the WL-1/t modification are quickly passed in the final stage of the calculation because all energies are almost equally probable. A much longer simulation time is required to satisfy these tests in the early calculation stage, when the probabilities of energy levels differ substantially.

IV.2 The control parameter

The largest eigenvalue of any stochastic matrix is equal to one, and we therefore propose to use the difference of the largest eigenvalue of the TMES from unity computed during the final stage of the WL simulation as a criterion for the proximity of the DOS to the true value.

We estimate the elements of the TMES in simulations as follows. The auxiliary matrix $U(E_{k},E_{m})$ is initially filled with zeros. The element $U(E_{k},E_{m})$ is increased by unity after every WL move from a configuration with the energy $E_{k}$ to a configuration with the energy $E_{m}$. During the simulations, we compute the normalized matrix $\tilde{T}(E_{k},E_{m})=U(E_{k},E_{m})/\tilde{H}$, where $\tilde{H}=\sum_{k,m}U(E_{k},E_{m})/N_{E}$. The obtained matrix $\tilde{T}$ approaches the stochastic matrix $T$ in the final stage of calculation. The difference of the largest eigenvalue $\lambda_{1}$ of $\tilde{T}$ from unity gives the control parameter $\delta=\left|1-\lambda_{1}\right|$.

There are many algorithms for computing the largest eigenvalue of a matrix, and almost all are suitable for calculating $\delta$. We used the power method, also known as power iteration or Von Mises iteration poweriteration . The algorithm does not compute a matrix decomposition, so it is quite efficient for large sparse matrices. It is terminated when a desired accuracy of the eigenvector approximation is achieved; the eigenvalue estimate is then found by applying the Rayleigh quotient to the resulting eigenvector. The method can be used if $\lambda_{1}$ is the eigenvalue of largest absolute value and $|\lambda_{1}/\lambda_{2}|\neq 1$, where $\lambda_{1},\dots,\lambda_{n}$ is the list of the matrix eigenvalues ordered so that $|\lambda_{1}|\geq|\lambda_{2}|\geq|\lambda_{3}|\geq...\geq|\lambda_{n}|$. The absolute value of any eigenvalue of any stochastic matrix is less than or equal to unity, therefore, the power method is applicable for estimating $\delta$ in the final stage of the WL-1/t algorithm. It is known that $|\lambda^{(k)}-\lambda_{1}|=O(|\lambda_{2}/\lambda_{1}|^{2k})$, where $\lambda^{(k)}$ is the approximation for $\lambda_{1}$ obtained after $k$ iterations Mehl , so the error asymptotically decreases by a factor of $|\lambda_{1}/\lambda_{2}|^{2}$ at each iteration.

The TMES is typically a sparse matrix, and its storage usually requires only $O(N_{E})$ of memory. The matrix-vector multiplications are performed very efficiently if the matrix is sparse, so each iteration of the power method requires only $O(N_{E})$ operations in this case. Software libraries such as ViennaCL ViennaCL contain the implementation of the power method for sparse matrices. The power method may require many iterations if $|\lambda_{1}/\lambda_{2}|\approx 1$. However, we note that the eigenvalue needs to be calculated only occasionally. For example, in our simulations, we calculate $\delta$ only once for each integer $n$, where $n\leq 100\log t<n+1$. Such a simulation applies the power method only several thousands of times during a WL-1/t calculation with $10^{13}$ spin flips, so the computing time used for the eigenvalue calculation is negligible.

IV.3 The histogram flatness

We can calculate the normalized histogram ${\cal H}=H(E_{m})/\sum_{m}H(E_{m})$ as ${\cal H}=\sum_{k}\tilde{T}(E_{k},E_{m})$. Hence, the histogram flatness condition is equivalent to the property that the matrix $\tilde{T}$ is close to stochastic. Thus, the histogram flatness is closely connected at the final simulation stage of the WL-1/t algorithm with the proximity to the true DOS.

For the original WL algorithm, there is no guarantee that the rate of increase of the logarithm of the nonnormalized DOS is the same for all energies in the final stage of the calculation because the parameter modification rule $F:=F/2$ results in a rapid decay of $F$, and the algorithm hence converges because the value of $F$ is negligible. The histogram flatness check is performed with a finite accuracy such as several percent, which results in a finite accuracy of the calculated DOS. The choice of high accuracy in the flatness criterion can result in a slow convergence and a very long simulation time Wang-Landau-PRE .

IV.4 Normalizing the DOS

Normalizing the DOS only at the end of the simulation was suggested in the original papers Wang-Landau ; 1overt ; SAMC . We note that this can limit the accuracy of the estimated DOS. For example, we consider the one-dimensional Ising model with $L=512$, where the transition to the second phase of the WL-1/t algorithm occurs at $t\sim t_{s}=2\cdot 10^{10}$, where $S(E,t_{s})\sim 10^{7}$. After only several hours of the calculation, we have $t=5\cdot 10^{11}$ and $F=N_{E}/t=5\cdot 10^{-10}$. The operation $S(E):=S(E)+F$ is then beyond the capabilities of double-precision floating-point variables because there is already a $17$ orders of magnitude difference between $S(E)$ and $F$. Hence, the operation is in fact not performed and the DOS is not updated after that. Therefore, we recommend normalizing the calculated DOS more frequently during the simulation. For the simulation corresponding to Fig. 1, the calculated DOS is normalized every time the values of $\delta$ and $\Delta$ are calculated.

IV.5 Behavior of the control parameter for the WL-1/t algorithm

The parameter

[TABLE]

estimates the deviation of the computed DOS $\tilde{g}(E_{k})$ from the exact DOS $g(E_{k})$. Figure 1 shows the behavior of $\overline{\Delta}$ and $\overline{\delta}$ as a function of simulation time $t$. The overline means that the data were obtained by averaging over $M$ independent runs of the algorithm to reduce statistical noise, where $M=60$ in Fig. 1.

We note that $\tilde{S}(E,t)$ in Eq. (8) corresponds to the normalized DOS. Here, we use the normalization $\tilde{S}(E,t)=S(E,t)-\Delta S$, where $\Delta S=S(E_{j},t)-S_{\text{exact}}(E_{j})$ and $j$ is chosen such that $S(E_{j})=\max_{j}S(E_{j})$. Both the abovementioned normalization to the total number of states and the normalization to the number of ground states turn out to give values of $\Delta$ close to those presented in Fig. 1. The vertical dashed line marks the average value of $t_{s}$.

Figure 1 demonstrates the monotonic power-law decrease of both the parameters $\delta$ and $\Delta$ during the second phase of the WL-1/t algorithm. We use the logarithmic scale in both axes. A stable power-law decay of the parameter $\delta$ reveals the convergence of $\tilde{T}$ to a stochastic matrix and can be used as a criterion for the convergence of the simulated DOS to the exact DOS.

The fluctuations of the parameters $\delta$ and $\Delta$ are shown in Fig. 2 for the simulations described in Fig. 1. Figure 2 shows $\sigma(\overline{\delta})/\overline{\delta}$ and $\sigma(\overline{\Delta})/\overline{\Delta}$ as functions of $t$. The relative standard deviations were obtained using 60 independent runs of the algorithm. Therefore, the values in Fig. 2 represent the relative magnitudes of the error bars in Fig. 1. It follows from Fig. 2 that $\sigma(\delta)=\sqrt{M}\sigma(\overline{\delta})$ and $\sigma(\Delta)=\sqrt{M}\sigma(\overline{\Delta})$ are of the order of $\overline{\delta}$ and $\overline{\Delta}$, respectively.

The condition $\delta(t_{2})\ll\delta(t_{1})$ observed during the second algorithm phase should result in satisfying the condition $\Delta(t_{2})\ll\Delta(t_{1})$, which allows approximating the value of $\Delta(t_{1})$ as the deviation between the DOS computed at $t=t_{1}$ and $t=t_{2}$. This allows estimating the simulation accuracy in the case where the DOS of the simulated system is not known exactly. In Fig. 3, as an example of such a case, we present the results of simulating the two-dimensional Potts model with $q{=}8$ spin states. The dependence of the parameters $\delta$ and $\tilde{\Delta}$ on $t$ are qualitatively similar to those calculated for the Ising model (Fig. 1). Because we do not have an analytic expression for the DOS in this case, we calculate the deviation of $\tilde{g}(E)$ using the expression $\tilde{\Delta}=1/N_{E}\cdot\sum_{E}\left|(\tilde{S}(E,t)-S_{0}(E))/S_{0}(E)\right|$ and taking $S_{0}(E)=\tilde{S}(E,t_{f})$ for a large value of $t_{f}$ ($t_{f}=2.6\cdot 10^{12}$ in Fig. 3). The control parameter $\delta$ can thus be used to estimate the accuracy of the obtained DOS.

Very similar results to those shown in Fig. 1 were obtained for various values of the lattice size. The calculations were performed with $L$ up to 1024 for the one-dimensional Ising model and up to 64 for the two-dimensional Ising model. Figures 4 and 5 show $\overline{\delta}(t)$ and $\overline{\Delta}(t)$ for several different values of the Ising model lattice size $L$, where $M=40$. Figures 1, 4 and 5 also demonstrate different values of $t_{s}$, which grows with the system size.

IV.6 Behavior of the control parameter for the original WL algorithm

Figure 6 shows $\overline{\delta}(t)$ and $\overline{\Delta}(t)$ for the original WL algorithm described in Wang-Landau . The algorithm was applied to the one-dimensional and two-dimensional Ising models with $L=32$. The data in the left panel were obtained by applying the WL algorithm to the one-dimensional Ising model and averaging over 40 independent runs. The right panel corresponds to a single run of the WL algorithm applied to the two-dimensional Ising model.

Therefore, both $\Delta$ and $\delta$ saturate for the original WL algorithm (see also Sec. IV.3). Using the control parameter $\delta$ thus confirms the systematic error of the original WL algorithm previously reported in [3, 6, 7, 24].

V Conclusion

We have analyzed properties of the algorithms and of the TMES. TMES of the WL random walk on the true DOS is stochastic and symmetric. We present analytic expressions for the TMES in the case of one-dimensional Ising model. We improve the WL algorithm based on the WL-1/t modification of the original algorithm 1overt and propose a method for examining the convergence of simulations to the true DOS and for controlling the accuracy of the DOS calculation. The monotonic power-law decrease of the control parameter $\delta$ during the second phase of the algorithm reveals the convergence of the algorithm, and the values of the control parameter can be used to estimate the accuracy of the DOS calculations.

This approach can be generalized to systems with an intitially unknown discrete spectrum, where the general procedure can be applied for the dynamic change of the TMES. It would be interesting to check its applicability to systems with a continuous energy spectrum.

Acknowledgements.

This work is supported by the grant 14-21-00158 from the Russian Science Foundation.

Appendix A Convergence of the WL-1/t algorithm to the true DOS

We have shown that the TMES $T$ of the WL random walk on the true DOS is stochastic, and also that the TMES $\tilde{T}$ is close to a stochastic matrix in the final stage of the WL-1/t algorithm.

Here we demonstrate that the obtained normalized DOS is close to the true DOS if the TMES $\tilde{T}$ is a stochastic matrix.

It follows from (4) that

[TABLE]

where $\tilde{g}(E)$ is the obtained normalized DOS. Using (1), we hence obtain

[TABLE]

where $\eta_{i}=g(E_{i})/\tilde{g}(E_{i})$ and $g(E)$ is the true DOS. It follows from (10) and the stochasticity of $\tilde{T}$ that

[TABLE]

Because the TMES is a stochastic matrix, the rates of visiting all energy levels are equal to each other. The values of $\tilde{g}(E)$ therefore remain almost the same, and the behavior of the algorithm is close to a Markov chain. Moreover, the invariant distribution of the Markov chain has the property that all energy levels are equiprobable. It follows from (11) that the values $\eta_{i}/\sum_{k}\eta_{k}$ represent the invariant distribution of the Markov chain. Therefore, $\eta_{i}$ is independent of $i$, and the obtained normalized DOS is hence close to the true DOS.

Appendix B Expressions for $N_{i}$ and $Q_{i}$.

We have the relations

[TABLE]

where $\delta_{L,2k}$ is the Kronecker delta.

Expression (12) is derived as follows. We consider the circular chain of $L{-}2k$ spins. We place the first domain wall in front of the first spin. We add another $2k{-}i{-}1$ domain walls in the remaining space between the spins; there are $C_{L-2k-1}^{2k-i-1}$ ways to do this. Therefore, we have $L{-}2k$ spins and $2k{-}i$ domain walls, where the first spin of the first domain is the first spin of the chain.

We then add one more spin in every domain. We also add $i$ domains consisting of only one spin. There are exactly $C_{2k}^{i}$ ways to choose $i$ domains among the $2k$ domains. Each of these choices unambiguously defines how to add $i$ domains, each consisting of only one spin, to the available $2k{-}i$ domains of the chain.

We have thus calculated the number of configurations of the circular chain of $L$ spins containing $2k$ domains such that $i$ domains consist of only one spin, $2k{-}i$ domains consist of more than one spin, and there is a domain wall in front of the first spin. This number is $M_{i}=2C_{2k}^{i}C_{L-2k-1}^{2k-i-1}$.

When $2k$ domain walls are placed among the $L$ spins, the probability that there is a domain wall in front of the first spin is equal to $p=2k/L$. Hence, $N_{i}=M_{i}/p$, i.e., we have obtained Eq. (12).

The justification of Eqs. (13) is as follows. We have $2k$ domains, where $i$ domains consist of only one spin and $2k{-}i$ domains consist of more than one spin. To remove a couple of domains with just a single spin flip, we must choose one of $i$ spins from the domains consisting of only one spin. Therefore, $Q_{i}^{E_{k}\to E_{k-1}}=i/L$.

To add a couple of domains with just a single spin flip, we must choose a spin that is not a boundary spin of a domain. There are $L{-}4k{+}i$ spins satisfying this condition because there are $2k$ spins located to the right of a domain wall, $2k$ spins located to the left of a domain wall, and $i$ spins which are located with a domain wall on both the right and the left. Therefore, $Q_{i}^{E_{k}\to E_{k+1}}=(L-4k+i)/L$. Finally, $Q_{i}^{E_{k}\to E_{k}}=1-Q_{i}^{E_{k}\to E_{k-1}}-Q_{i}^{E_{k}\to E_{k+1}}=(4k-2i)/L$.