An Efficient Simulation Protocol for Determining the Density of States: Combination of Replica-Exchange Wang-Landau Method and Multicanonical Replica-Exchange Method

TL;DR

This paper introduces a new simulation protocol combining two advanced algorithms to accurately determine the density of states, demonstrated on the Ising model, outperforming individual methods in precision.

Contribution

The paper presents the REWL-MUCAREM protocol, integrating REWL and MUCAREM for improved density of states estimation in statistical models.

Findings

REWL-MUCAREM yields more accurate density of states than individual methods.

The combined protocol outperforms separate REWL and MUCAREM in simulations.

Validation performed on the square-lattice Ising model.

Abstract

By combining two generalized-ensemble algorithms, Replica-Exchange Wang-Landau (REWL) method and Multicanonical Replica-Exchange Method (MUCAREM), we propose an effective simulation protocol to determine the density of states with high accuracy. The new protocol is referred to as REWL-MUCAREM, and REWL is first performed and then MUCAREM is performed next. In order to verify the effectiveness of our protocol, we performed simulations of a square-lattice Ising model by the three methods, namely, REWL, MUCAREM, and REWL-MUCAREM. The results showed that the density of states obtained by the REWL-MUCAREM is more accurate than that is estimated by the two methods separately.

Tables

Table 1. Table 1: Conditions of the present simulations.

			Frequency (in MC sweeps)
Methods	Number of spins	Number of replicas	of flatness evaluation	Total MC sweeps
	$N$	$M$	in Eq. (5)	per replica
MUCAREM	64	4	NA	200 000
	256	8		200 000
	1024	16		200 000
	4096	32		500 000
	16384	64		3 000 000
REWL	64	4	1000	200 000
	256	8		200 000
	1024	16		200 000
	4096	32		500 000
	16384	64		3 000 000
REWL$-$MUCAREM	64	4	1000$-$NA	100 000$-$100 000
	256	8		100 000$-$100 000
	1024	16		100 000$-$100 000
	4096	32		250 000$-$250 000
	16384	64		1 500 000$-$1 500 000

Table 2. Table A-I: Conditions of MUCAREM simulations.

Methods	Number of spins	MUCAREM	Number of replicas		Number of MC sweeps		Total MC sweeps
	$N$	iterations			per replica
			REM	MUCAREM	REM	MUCAREM
MUCAREM1	64	1	4	4	100 000	100 000	$\mathrm{200\hspace{0.17em}000}\times 4$
	256		8	8	100 000	100 000	$\mathrm{200\hspace{0.17em}000}\times 8$
	1024		16	16	100 000	100 000	$\mathrm{200\hspace{0.17em}000}\times 16$
	4096		32	32	250 000	250 000	$\mathrm{500\hspace{0.17em}000}\times 32$
	16384		64	64	1 500 000	1 500 000	$\mathrm{3\hspace{0.17em}000\hspace{0.17em}000}\times 64$
MUCAREM2	64	1	8	4	10 000	180 000	$\mathrm{200\hspace{0.17em}000}\times 4$
	256		16	8	10 000	180 000	$\mathrm{200\hspace{0.17em}000}\times 8$
	1024		32	16	10 000	180 000	$\mathrm{200\hspace{0.17em}000}\times 16$
	4096		64	32	25 000	450 000	$\mathrm{500\hspace{0.17em}000}\times 32$
	16384		128	64	150 000	2 700 000	$\mathrm{3\hspace{0.17em}000\hspace{0.17em}000}\times 64$
MUCAREM3	64	2	8	4	10 000	90 000	$\mathrm{200\hspace{0.17em}000}\times 4$
	256		16	8	10 000	90 000	$\mathrm{200\hspace{0.17em}000}\times 8$
	1024		32	16	10 000	90 000	$\mathrm{200\hspace{0.17em}000}\times 16$
	4096		64	32	25 000	225 000	$\mathrm{500\hspace{0.17em}000}\times 32$
	16384		128	64	150 000	1 350 000	$\mathrm{3\hspace{0.17em}000\hspace{0.17em}000}\times 64$

Equations

$$P_{\mathrm{MUCA}}(E)\propto g(E)W_{\mathrm{MUCA}}(E)\equiv \mathrm{const},$$

$$W_{\mathrm{MUCA}}(E)$$

$$w\left(E\to E^{\prime }\right)=\min \text{sl}\left[1,\frac{W_{\mathrm{MUCA}}\text{sl}(E^{\prime })}{W_{\mathrm{MUCA}}\text{sl}(E)}\right]=\min \text{sl}\left[1,\frac{g(E)}{g(E^{\prime })}\right].$$

$$g(E)=\frac{H(E)}{W_{\mathrm{MUCA}}\text{sl}(E)},$$

$$\frac{H_{\min }}{H_{\max }}>0.5,$$

$$W_{\mathrm{MUCA}}^{\{m\}}(E)=\{\begin{array}{ll}{e}^{-{\beta }_{\mathrm{L}}^{\{m\}}E},& \mathrm{for}E<E_{\min }^{\{m\}},\\ \frac{1}{g_m(E)},& \mathrm{for}{E}_{\min }^{\{m\}}\le E\le E_{\max }^{\{m\}},\\ e^{-{\beta }_{\mathrm{H}}^{\{m\}}E},& \mathrm{for}E>E_{\max }^{\{m\}},\end{array}$$

$$W_{\mathrm{MUCAREM}}(E)=$$

$$w_{\mathrm{MUCAREM}}$$

$$\displaystyle\left\{\begin{array}[]{l}g(E)=\frac{\displaystyle{\sum_{m=1}^{M}}H_{m}(E)}{\displaystyle{\sum_{m=1}^{M}}n_{m}\exp\left(f_{m}\right)W_{{\rm MUCA}}^{\{m\}}(E)}~{},\\ \\ \exp\left(-f_{m}\right)=\displaystyle{\sum_{E}g(E)W_{{\rm MUCA}}^{\{m\}}(E)}~{}.\end{array}\right.$$

$$E=-J\sum _{⟨i,j⟩}{S}_i{S}_j,$$

$$C(T)=\frac{⟨E^2{⟩}_T-⟨E{⟩}_T^2}{T^2},$$

$$⟨A{⟩}_T=\frac{\sum _{E}A(E)g(E)e^{-\beta E}}{\sum _{E}g(E)e^{-\beta E}},$$

$${\epsilon }_A=\sqrt{\frac{\sum _{i=1}^n{\left(A^{\left\{i\right\}}-\overline{A}\right)}^2}{n\left(n-1\right)}},\overline{A}=\frac{\sum _{i=1}^n{A}^{\left\{i\right\}}}{n}.$$

$$\{\begin{array}{ll}G(E)& =\frac{1}{n}\sum _{i=1}^n{G}^{\left\{i\right\}}(E),\\ G^{\left\{i\right\}}(E)& =\frac{g_{\mathrm{sim}}^{\left\{i\right\}}(E)}{g_{\mathrm{EXACT}}(E)}.\text{par}\end{array}$$

$$F\equiv \frac{G_{\min }}{G_{\max }},$$

$$\{\begin{array}{ll}{F}_{\min }& =\frac{G_{\min }-\frac{1}{2}\epsilon (E_{\mathrm{glmin}})}{G_{\max }+\frac{1}{2}\epsilon (E_{\mathrm{glmax}})},\\ F_{\max }& =\frac{G_{\min }+\frac{1}{2}\epsilon (E_{\mathrm{glmin}})}{G_{\max }-\frac{1}{2}\epsilon (E_{\mathrm{glmax}})}.\end{array}$$

$$S_1\left(\frac{1}{N}\right)=k_{\mathrm{B}}\log W_1\left(\frac{1}{N}\right).$$

$$W_1\left(\frac{1}{N}\right)=W_0+a_1{\left(\frac{1}{N}\right)}^{\theta }$$

$$W_0$$

$$S_0=0.815627\pm 0.000039\left[cal/deg\mathrm{mole}\right].$$

Full text

I Introduction

The statistical mechanical expectation value of a physical quantity can be accurately calculated if the density of states (DOS) is given. However, in many cases, we do not know DOS *a priori *and it is often difficult to obtain it theoretically or experimentally. In recent decades, many methods were developed for the determination of DOS by using Monte Carlo (MC) and/or molecular dynamics (MD) simulations (e.g., see Refs. UMBRELLA –R12c ). One of the earliest such methods may be Umbrella Sampling UMBRELLA . Muticanonical Algorithm R1 –R3 , Simulated Tempering ST1 –ST3 , Replica-Exchage Method R6 –R7 , Wang-Landau method R4a ; R4b , and Metadynamics META –META3 were then developed, and generalizations and extensions of these methods were further proposed R7b –R12c . These methods are closely related. For example, it has been shown that Statistical Temperature Molecular Dyanamics STMD is equivalent to Metadynamics RESTMD_META . We also remark that Metadynamics can be considered to be Wang-Landau method in reaction coordinate space (rather than energy space) RevYO . These methods have been successfully applied to a wide range of problems in condensed matter and statistical physics including spin glasses, liquid crystals, polymers, and proteins. Nevertheless, the problem still remains that the exact calculation of DOS cannot be achieved when the systems become large and complex.

In this article, we propose an efficient simulation protocol to obtain the most precise DOS by combining the Replica-Exchange Wang-Landau method (REWL)R11 ; R12 and the Multicanonical Replica-Exchange Method (MUCAREM)R8 ; R9 ; R10 .

This article is organized as follows. In Sec. II, we explain the methods. In Sec. III, the computational details are given. In Sec. IV, the results and discussion are presented. Sec. V is devoted to conclusions.

II Computational methods

We first introduce three basic generalized-ensemble algorithms, Multicanonical Algorithm, Wang-Landau method, and Replica-Exchange Method. The Multicanonical Algorithm (MUCA)R1 ; R2 is one of the representative methods. A simulation in multicanonical ensemble is based on a non-Boltzmann weight factor, which we refer to as the multicanonical weight factor. This is inversely proportional to DOS of the system, and a free random walk in potential energy space is realized so that a wide configurational space may be sampled. The DOS is often not known a priori. The multicanonical weight factor is usually determined by iterations of short trial simulationsR3 ; R20 . After a production run with the determined MUCA weight factor, the single-histogram reweighting techniquesR14a are employed to obtain an accurate DOS. However, the weight factor determination process can be tedious and difficult. The Wang-Landau method (WL)R4a ; R4b solved this problem drastically. In the WL sampling, the weight factor, which is also inversely proportional to DOS, is updated during the simulation by multiplying a constant to the weight factor. This procedure leads to a uniform histogram in potential energy space, and the modified weight factor converges to the inverse of the DOS. Another powerful algorithm is the Replica-Exchange Method (REM)R6 ; R6a (it is also referred to as parallel temperingMPRL ). Closely related method was independently developed in SW . In this method, several copies (replicas) of the original system at different temperatures are simulated independently and simultaneously by conventional canonical MC or MD. Every few steps, pairs of replicas are exchanged with a specified transition probability. This exchange process realizes a random walk in temperature space, which in turn induces a random walk in potential energy space. After a production simulation, the multiple histogram reweighting techniquesR14b ; R15 (an extension of which is also referred to as the Weighted Histogram Analysis Method (WHAM)R15 ) are used in order to determine the most accurate DOS from all the histograms of sampled potential energy at different temperatures.

These basic simulation methods can be combined for more effective sampling. One method is referred to as the Multicanonical Replica-Exchange Method (MUCAREM)R8 ; R9 ; R10 . In this method, the total energy range where we want to calculate the DOS is divided into smaller regions, each corresponding to a replica, and MUCA simulations are performed independently and simultaneously in each replica. Every few steps, a pair of neighboring replicas are exchanged like REM. The configurations can be sampled more effectively than ordinary MUCA because of replica exchange. The final, most accurate estimation of DOS is obtained by the multiple-histogram reweighting techniques again R8 ; R9 ; R10 . A similar method is the Replica-Exchange Wang-Landau (REWL) methodR11 ; R12 . The idea is almost the same as in MUCAREM except for using WL instead of MUCA for each replica. After a REWL simulation, DOS pieces are obtained for different energy regions. Connecting these pieces at the point where the slope of DOS is coincident, we can obtain the final estimation of DOS over the entire energy range.

We found that the DOS with the highest accuracy can be obtained by combining these two methods. The REWL is employed in the first half of the total number of MC (or MD) steps in order to get a rough estimate of MUCAREM weight factor and the MUCAREM is performed in the second half in order to refine the DOS. We refer to this protocol as REWL-MUCAREM. The DOS thus obtained has higher accuracy than that is estimated by the two methods separately.

A brief explanation of MUCA is now given here. The multicanonical probability distribution of potential energy $P_{{\rm MUCA}}(E)$ is defined by

[TABLE]

where $W_{{\rm MUCA}}(E)$ is the multicanonical weight factor and the function $g(E)$ is the DOS. $E$ is the total potential energy of a system. By omitting a constant factor, we have

[TABLE]

In MUCA MC simulations, the trial moves are accepted with the following Metropolis transition probability $w\left(E\rightarrow E^{\prime}\right)$:

[TABLE]

Here, $E$ is the potential energy of the original configuration and $E^{\prime}$ is that of a proposed one. After a long production run, the best estimate of DOS can be obtained by the single-histogram reweighting techniques:

[TABLE]

where $H(E)$ is the histogram of sampled potential energy. Practically, the $W_{{\rm MUCA}}(E)$ is set $\exp[-\beta E]$ at first and modified by repeating sampling and reweighting. Here, $\beta$ is the inverse of temperature $T$ ($\beta=1/k_{{\rm B}}T$ with $k_{{\rm B}}$ being the Boltzmann constant).

The WL also uses $1/g(E)$ as the weight factor and the Metropolis criterion is the same as in Eq.(3). However, $g(E)$ is updated dynamically as $g(E)\rightarrow f\times g(E)$ during the simulation when the simulation visits a certain energy value $E$. $f$ is a modification factor. We continue the updating until the energy histogram becomes flat. If $H(E)$ is flat enough, a next simulation begins after resetting the histogram to zero and reducing the modification factor (usually, $f\rightarrow\sqrt{f}$). The flatness evaluation can be done in various ways. In this article, we considered that the histogram is sufficiently flat when

[TABLE]

where $H_{\rm min}$ and $H_{\rm max}$ are the least number and the largest number of nonzero entries in the histogram, respectivelyR16 . This process is terminated when the modification factor attains a predetermined value $f_{{\rm final}}$ and $\exp(10^{-8})\simeq 1.000\,000\,01$ is often used as $f_{{\rm final}}$. Hence, the estimated $g(E)$ tends to converge to the true DOS of the system within this much accuracy set by $f_{{\rm final}}$. (Several reports argue that there is a possibility that conventional WL algorithm has a systematic error which does not decrease any more. See, e.g., Ref.TWL .)

In MUCAREM, the entire energy range of interest $\left[E_{{\rm min}},E_{{\rm max}}\right]$ is divided into $M$ sub-regions, $E_{{\rm min}}^{\{m\}}\leq E\leq E_{{\rm max}}^{\{m\}}$ $(m=1,2,\cdots,M)$, where $E_{{\rm min}}^{\{1\}}=E_{\rm min}$ and $E_{{\rm max}}^{\{M\}}=E_{\rm max}$. There should be some overlaps between the adjacent regions. MUCAREM uses $M$ replicas of the original system. The weight factor for sub-region $m$ is defined byR8 ; R9 ; R10 :

[TABLE]

where $g_{m}(E)$ is the DOS for $E^{\{m\}}_{{\rm min}}\leq E\leq E^{\{m\}}_{{\rm max}}$ in sub-region $m$, $\beta_{\rm L}^{\{m\}}=d\log\left[g_{m}(E)\right]/dE~{}(E=E^{\{m\}}_{\rm min})$ and, $\beta_{\rm H}^{\{m\}}=d\log\left[\it g_{m}(E)\right]/dE~{}(E=E^{\{m\}}_{\rm max})$. The MUCAREM weight factor $W_{{\rm MUCAREM}}(E)$ for the entire energy range is expressed by the following formula:

[TABLE]

After a certain number of independent MC steps, replica exchange is proposed between two replicas, $i$ and $j$, in neighboring sub-regions, $m$ and $m+1$, respectively. The transition probability, $w_{{\rm MUCAREM}}$, of this replica exchange is given by

[TABLE]

where $E_{i}$ and $E_{j}$ are the energy of replicas $i$ and $j$ before the replica exchange, respectively. If replica exchange is accepted, the two replicas exchange their weight factors $W^{\{m\}}_{{\rm MUCA}}(E)$ and $W^{\{m+1\}}_{{\rm MUCA}}(E)$ and energy histogram $H_{m}(E)$ and $H_{m+1}(E)$. The final estimation of DOS can be obtained from $H_{m}(E)$ after a simulation by the multiple-histogram reweighting techniques or WHAM. Let $n_{m}$ be the total number of samples for the $m$-th energy sub-region. The final estimation of DOS, $g(E)$, is obtained by solving the following WHAM equations self-consistently by iterationR9 :

[TABLE]

These MUCAREM sampling and WHAM reweighting processes can, in principle, be repeated to obtain more accurate DOSR10 . We remark that REM is often used to obtain the first estimate of DOS in the MUCAREM iterations. We also remark that when REM instead of MUCAREM is performed, the best estimate of DOS can be obtained by solving Eq. (12), where $W^{\{m\}}_{{\rm MUCA}}(E)$ is replaced by $\exp(-\beta_{m}E)$ with temperature $T_{m}~{}(\beta_{m}=1/k_{\rm B}T_{m})$ for ($m=1,2,\cdots M$).

The REWL method is essentially based on the same weight factors as in MUCAREM, while the WL simulations replace the MUCA simulations for each replica. This simulation is terminated when the modification factors on all sub-regions attain a certain minimum value $f_{{\rm final}}$. After a REWL simulation, $M$ pieces of DOS fragments with overlapping energy intervals are obtained. The fragments need to be connected in order to determine the final DOS in the entire energy range $\left[E_{{\rm min}},E_{{\rm max}}\right]$. The joining point for any two overlapping DOS pieces is chosen where the inverse microcanonical temperature $\beta=d\log\left[\it g(E)\right]/d\it E$ coincides bestR11 ; R12 .

III Computational details

In order to compare the effectiveness of the REWL-MUCAREM with other methods, we performed simulations of a 2-dimensional Ising model with periodic boundary conditions.

In a square-lattice Ising model, the total energy $E$ is defined by

[TABLE]

where $i$ and $j$ are labels for lattice points. $J$ is the magnitude of interaction between neighboring spins. In this article, $J$ and $k_{{\rm B}}$ are set to one for simplicity. $\langle i,j\rangle$ represents pairs of nearest-neighbor spins. $S_{i}$ is the state of spin on a lattice point $i$ and takes on values of $\pm 1$. Beale calculated the exact DOS of the model of finite sizesR17 ; R18 .

Table I lists the conditions of our simulations. The total number of spins $N$ is $L^{2}$, where $L$ is the length of a side of the square lattice, The total number of spins considered was $N=64,256,1024,4096$, and $16384$. One MC sweep is defined as an evaluation of Metropolis criteria $N$ times. The cost of computations (for example, the total number of MC sweeps) was set equal. However, we should point out that while the ordinary REWL algorithm is terminated when the recursion factor $f$ converged to $f_{\rm final}$, but our REWL simulations were finished after a certain fixed number of flatness evaluations had been made.

A Marsaglia random number generator was employed and we used the program code on open sourceR19 ; R20 . The number of replicas was set equal to $L/2$. Each replica performed a MUCA simulation in MUCAREM or a WL simulation in REWL within their energy sub-regions, which had an overlap of about $80$ percent between neighboring sub-regions. In the cases of REWL and REWL-MUCAREM simulations, WL flatness criterion was tested every $1000$ MC sweeps. If the histogram of energy distribution is sufficiently flat at this time, the WL recursion factor was reduced. Replica exchange was tried every $100$ MC sweeps. The cost of calculation in our simulations was measured by the total number of MC sweeps because we spend most of computational time to perform MC simulations. With the conditions in Table I, we made $n=25$ independent runs with different initial random number seeds. In this work, we did not iterate the DOS evaluation during the MUCAREM simulations for simplicity. In the present MUCAREM simulation, the first half of the total MC sweeps was run with REM and the remaining of the simulation was MUCAREM with the DOS obtained from the REM simulation. We evaluated the effectiveness of iterations of MUCAREM and WHAM (see Appendix A). In the REM simulation, $M$ temperature values were evenly distributed between $\beta_{1}=1.0$ and $\beta_{M}=0.01$.

IV Results and discussion

The four figures in Figs. 1(A) show the specific heat which was calculated from the estimated DOS by using the following equation:

[TABLE]

where

[TABLE]

and $A(E)$ is any physical quantity that depends on $E$. The errors were obtained by the standard error estimation:

[TABLE]

Here, $A^{\left\{i\right\}}$ is obtained from the $i$-th simulation $(i=1,2,\cdots,n)$. Although the exact values of specific heat in finite sizes were obtained by Ferdinand and FisherR21 , we calculated the exact specific heat in Fig. 1(A)(a) from the exact DOS, $g_{{\rm EXACT}}(E)$, of BealeR17 ; R18 by using Eqs. $(\ref{EQ11})$ and $(\ref{EQ12})$. We used the Mathematica code, which is given in R18 , for the calculations of $g_{{\rm EXACT}}(E)$. All the algorithms could reproduce the exact solutions very well.

The differences between exact values and simulation results are shown in Fig. 1(B). Note that $|\Delta C(T)|$ takes maximum values around the phase transition temperature $T_{c}=2/\log(1+\sqrt{2})\simeq 2.269$ in each method. The results imply that the results of the three methods agree with the exact ones in the order of REWL-MUCAREM, REWL, and MUCAREM. It means that REWL-MUCAREM could get more accurate DOS than the other two methods.

In order to directly compare the accuracy of DOS among the three methods, we show the mean local flatness $G(E)$ in Fig. 2 for the system of $N=32\times 32$, where

[TABLE]

Here, $g_{{\rm sim}}^{\left\{i\right\}}$ is the DOS estimated from the $i$-th simulation $(i=1,2,\cdots,n)$. We matched $\log g_{{\rm EXACT}}(E)$ and $\log g_{{\rm sim}}^{\left\{i\right\}}(E)$ at $E/N=-0.5$. If $g_{{\rm sim}}^{\left\{i\right\}}(E)$ is equal to $g_{{\rm EXACT}}(E)$, $G^{\left\{i\right\}}(E)$ becomes flat ($=1$) ideally in the entire energy range. The red curves and the yellow vertical bars in Fig. 2 are the values of $G(E)$ and the error bars, respectively. The errors were also estimated by the standard error estimation. The difference between red curve and black base line became large at lower energy region. The tendency became stronger in larger systems, where the phase transition became stronger (see Fig 1). Because the error bars are the smallest at lower energy region among the three methods, REWL-MUCAREM could obtain more precise DOS than REWL and MUCAREM.

In order to examine the accuracy of DOS further, we define the degree of global flatness $F$ by the following formula:

[TABLE]

where $G_{\rm min}$ is the minimum value of $G(E)$ over the entire energy range and $G_{\rm max}$ is the maximum one. $F$ takes on values between 0 and 1. The closer the calculated $g_{\rm sim}^{\left\{i\right\}}(E)$ is to $g_{\rm EXACT}(E)$ globally, the closer $F$ is to $1$. Fig. $3$ shows the measured flatness $F$. As a measure of errors, we define the minimum value of $F$ and maximum one by

[TABLE]

Here, $\varepsilon(E)$ is the standard error of $G(E)$, and $E_{{\rm glmin}}$ and $E_{{\rm glmax}}$ are the energy values where $G_{{\rm min}}$ and $G_{{\rm max}}$ are obtained, respectively. It is obvious that the value deteriorates as the size of system gets larger. This means that it was difficult to estimate the DOS of large systems because of the large degrees of freedom. We needed more samples in order to obtain DOS for larger systems. We could not find much differences among the three methods up to $N=16\times 16$. One bad data was founded in MUCAREM for the size of $N=32\times 32$, and it made the error bar larger than the other methods. If we consider the system larger than $N=32\times 32$, REWL-MUCAREM gave the best results among the three methods.

In the present implementation of MUCAREM, we performed two multiple-histogram reweighting (WHAM) operations: one after the REM simulation in the first half of the run and second after the MUCAREM simulation in the second half of the run. Although a second WHAM operation converges quickly, the calculation cost of the first WHAM operation can become non-negligible in large systems. The REWL-MUCAREM uses only the second WHAM and this WHAM converges even more quickly than in MUCAREM, because a good estimate of DOS is already prepared by the preceding REWL simulation. Hence, REWL first and MUCAREM second is the order that we want to adopt in REWL-MUCAREM. Note also that in REWL-MUCAREM, we do not need the piece-connecting process of DOS required in REWL, because WHAM automatically gives DOS in the entire energy range of interest.

V Conclusions

In this article, we investigated an effective simulation protocol to estimate the density of states with highest accuracy. We proposed REWL-MUCAREM that combines the advantages of REWL and MUCAREM, where REWL is performed first and MUCAREM is performed next. This protocol was compared with existing two methods in a square-lattice Ising model, and the results showed that REWL-MUCAREM gave the most accurate density of states.

REWL-MUCAREM is effective with other systems and it can be extended to the MD simulation, because MUCA MD R2b ; R24 and WL MD R25 have already been developed. We have already calculated the residual entropy of Ice Ih R22 and we applied the protocol to helix-coil transitions of homo-polymers R23 . The pre-views of these results are given in Appendix B. REWL-MUCAREM MD simulations of protein folding are now under way. We remark that there is an article which says that improvements can be transfered between MC and MD broad histogram methods RESTMD_META . The protocol of REWL-MUCAREM is easy to implement and it can give more reliable results.

Acknowledgements:

Some of the computations were performed on the supercomputers at the Supercomputer Center, Institute for Solid State Physics, University of Tokyo.

Appendix A: Optimization of Conditions in MUCAREM

We would like to discuss the optimization conditions in MUCAREM in order to obtain more accurate DOS in REWL-MUCAREM. We performed two more MUCAREM simulations with different conditions. Table A-I lists the conditions of the additional simulations (MUCAREM$2$ and MUCAREM$3$) together with the first MUCAREM in Table I (which is now referred to as MUCAREM$1$). The major differences of the additional MUCAREM simulations from the previous MUCAREM simulation lies in the following: number of MC sweeps for REM and MUCAREM, number of replicas used for REM, and number of iterations of MUCAREM. In the additional MUCAREM simulations, the $10~{}\%$ of the total MC sweeps was run with REM and the remaining $90~{}\%$ of the simulation was MUCAREM with the DOS obtained from the preceding REM simulation. In REM simulations in MUCAREM$2$ and MUCAREM$3$, replica exchange was proposed every $10$ MC sweeps. On the other hand, in MUCAREM in MUCAREM$2$ and MUCAREM$3$, replica exchange was proposed every $100$ MC sweeps.

It is often said that the DOS obtained by MUCAREM simulations becomes better by iterating MUCAREM and WHAM reweightingR10 . We iterated MUCAREM simulations once for MUCAREM$3$. It should be mentioned that because we obtained clearly wrong DOS, we performed one extra run for each system $N=32\times 32$ and $N=128\times 128$ in MUCAREM$2$, and simply discarded these apparently bad runs. (We did not find a bad run in REWL and REWL-MUCAREM simulations.) We think the problem came from the difficulty of uniformly sampling over a wide energy range in REM. The rough DOS obtained from the first REM was not good and the inaccuracy had a bad influence to the sampling of the following MUCAREM. Giving apparently wrong results suggests that MUCAREM simulations are unstable comparing to REWL and REWL-MUCAREM simulations, which implies that the total number of MC sweeps and/or the number of runs $n$ should be longer for MUCAREM simulations to conclude with confidence.

Fig. A1 shows the global flatness $F$ in Eq. $(\ref{EQ18})$. Although there are little differences in $F$ among the three MUCAREM simulations up to the system $N=64\times 64$, we found a large difference in the system $N=128\times 128$. MUCAREM$3$ could obtain a good estimate of DOS compared with MUCAREM$1$ and MUCAREM$2$. It implies that the DOS became better by iterating MUCAREM simulations rather than by using more sampling obtained from a single, long run of MUCAREM. The error bar for the system $N=64\times 64$ is large if it was compared to other size of systems in MUCAREM$3$. We found a bad result from one run out of the $25$ runs, which made the error bar larger. These again suggest that MUCAREM simulations are unstable compared to REWL and REWL-MUCAREM simulations. Note that the estimation of DOS under the conditions of MUCAREM$3$ for $N=128\times 128$ is even better than that of REWL-MUCAREM in Fig. $3$, although MUCAREM seems to be unstable. We expect that REWL-MUCAREM could remove the instability of MUCAREM and it would give better DOS if its MUCAREM simulation in REWL-MUCAREM was also iterated.

We can obtain good estimate of DOS under appropriate conditions and the DOS becomes more accurate by iterating the MUCAREM and WHAM reweighting operations. Although the most suitable conditions will depend on systems and methods, the combination of REWL and MUCAREM can give accurate DOS without worrying about bad data.

Appendix B: Applications of REWL-MUCAREM

We applied the REWL-MUCAREM MC protocol to two characteristic systems, Ice Ih and biopolymer. In this Appendix, we show some preliminary results.

The first example is the estimation of the residual entropy of Ice Ih by REWL-MUCAREM R22 . According to Pauling’s theory, the residual entropy is obtained from the degrees of freedom of orientations of water molecules which are observed in the groundstate ICE_PAU . Two simple Potts-like models, which are referred to as the 2-site model and 6-state model, with nearest-neighbor interactions on three-dimensional hexagonal lattice were introduced and the residual entropy was estimated by a MUCA simulation R16 . We applied our protocol to the 2-site model for obtaining the residual entropy with high accuracy. The final estimation of the residual entropy is estimated by extrapolation, taking the thermodynamic limit.

The estimations of the degrees of freedom of orientations of one water molecule $W_{1}(1/N)$ are shown in Fig. B1(a). Here, $N$ stands for the total number of water molecules of the system. The relationship between the degrees of freedom $W_{1}(1/N)$ and the residual entropy $S_{1}(1/N)$ is given by

[TABLE]

In Fig. B1(a), three data points ($N=1600,2880,4704$) are plotted. It was not possible to obtain the value for $N=4704$ by the previous MUCA simulations R16 ; ICE2012 , and we needed the REWL-MUCAREM to obtain this new value. A fit (the green curve in Fig B1(a)) for the data to the form

[TABLE]

is shown. Here, $W_{0}$, $a_{1}$ and $\theta$ are fitting parameters and $W_{0}$ converts to the final estimation of the residual entropy $S_{0}$ in the thermodynamical limit ($N\rightarrow\infty$).

Our final estimation of $W_{0}$ is R22

[TABLE]

and the final residual entropy $S_{0}$ is

[TABLE]

These results agree well with the other results in ICE_J ; ICE_V , which were estimated by other simulation methods.

Another example is a folding simulation of a simple biopolymer R23 . In order to examine the effectiveness of our protocol for protein folding simulation, we studied the helix-coil transition of a deca-alanine (which is a helix former) with AMBER99/GBSA force field. REWL-MUCAREM MC protocol was employed and the dihedral angles between residues were updated by the Metropolis criterion during simulations.

The temperature dependence of the average helix length of deca-alanine is shown in Fig. B1(b). Because the structures of the terminal residues are frayed, the maximum helix length is 8. The deca-alanine is in a coil state above the transition temperature of $T_{c}~{}\approx~{}350$ K and in a helix state below $T_{c}$. Folding simulations of larger and more complex proteins are under way R23 .