# Flat histogram Monte Carlo method

## 1 Introduction

The basic problem in equilibrium statistical mechanics is to compute the canonical average

$$A_T=\frac{\underset{\{\sigma \}}{}A(\sigma )\mathrm{exp}\left(H(\sigma )/kT\right)}{_{\{\sigma \}}\mathrm{exp}\left(H(\sigma )/kT\right)}.$$
(1)   
In addition, free energy

$$F=kT\mathrm{ln}\underset{\{\sigma \}}{}\mathrm{exp}\left(H(\sigma )/kT\right)$$
(2)   
and related quantity entropy are also very important. Standard Monte Carlo method, e.g., Metropolis importance sampling algorithm, is simple and general. However, computation of free energy  is difficult with such method.

Over the last decade, there have been a number of methods addressing this problem . A common theme in these approaches is to evaluate the density of states $`n(E)`$ directly. If this can be done with sufficient accuracy, then the summation over all configuration states can be rewritten as sum over energy only, e.g.,

$$F=kT\mathrm{ln}\underset{E}{}n(E)\mathrm{exp}(E/kT).$$
(3)   
Can we evaluate $`n(E)`$ with a uniform relative accuracy for all $`E`$? Our experience with flat histogram method suggests a positive yes.

## 2 Broad histogram equation

Oliveira et al  showed the validity of following equation relating density of states with the microcanonical average number of moves in a single spin flip dynamics:

$$n(E)N(\sigma ,\mathrm{\Delta }E)_E=n(E+\mathrm{\Delta }E)N(\sigma ^{},\mathrm{\Delta }E)_{E+\mathrm{\Delta }E}.$$
(4)   
This equation is equivalent to a detail balance condition in an infinite temperature transition matrix Monte Carlo simulation . The quantity $`N(\sigma ,\mathrm{\Delta }E)`$ is the number of ways that the system goes to a state with energy $`E+\mathrm{\Delta }E`$ by a single spin flip, given that the current state is $`\sigma `$ with energy $`E`$. The average $`\mathrm{}_E`$ is performed over all the states with a fixed initial energy $`E`$ (i.e., a microcanonical average).

## 3 A flat histogram dynamics

Consider the follow Monte Carlo dynamics 

1. Pick a site at random.
2. Flip the spin with probability $`r(E^{}|E)`$.
3. Sample $`N(\sigma ,\mathrm{\Delta }E)`$, i.e., accumulate the statistics for $`N(\sigma ,\mathrm{\Delta }E)`$.
4. Go to 1.

The flip probability $`r`$ is given by

$$r(E^{}|E)=\mathrm{min}(1,\frac{N(\sigma ^{},\mathrm{\Delta }E)_E^{}}{N(\sigma ,\mathrm{\Delta }E)_E}),$$
(5)   
where the current state $`\sigma `$ has energy $`E`$, and new state $`\sigma ^{}`$ with one spin flipped has energy $`E^{}=E+\mathrm{\Delta }E`$.

With the above choice of flip rate, we can show that detail balance is satisfied

$$r(E^{}|E)P(\sigma )=r(E|E^{})P(\sigma ^{})$$
(6)   
if $`P(\sigma )=\mathrm{const}/n(E(\sigma ))`$, since this equation is equivalent to the broad histogram equation, Eq. (4). The histogram in energy is then $`H(E)P(\sigma )=\mathrm{constant}`$. It turns out that the choice of flip rate is not unique; many other formulas are possible .

Due to Eq. (4), the flip rate is also equal to $`\mathrm{min}[1,n(E)/n(E^{})]`$. This is exactly Lee’s method  of entropy sampling (which is equivalent to multicanonical method ). However, since neither $`n(E)`$ nor $`N(\sigma ,\mathrm{\Delta }E)_E`$ is known before simulation, the way by which simulation gets boot-trapped is quite different. Our method is very efficient in this respect. Another important difference is that we take $`N(\sigma ,\mathrm{\Delta }E)_E`$ as our primary statistics, from which we derive the density of states $`n(E)`$. Apart from the number of iterations needed, the transition matrix results are in general more accurate than that obtained using the energy histogram .

## 4 Simulation procedures

The first approximate scheme we use is to replace the true microcanonical average by a cumulative sample average

$$N(\sigma ,\mathrm{\Delta }E)_E\frac{1}{M}\underset{i=1}{\overset{M}{}}N(\sigma ^i,\mathrm{\Delta }E),$$
(7)   
where the samples $`\sigma ^i`$ are configurations generated during simulation with energy $`E`$. Each state in the simulation contributes to $`N(\mathrm{})`$ for some $`E`$. For those $`E^{}`$ where data are not available, we set $`r=1`$. This biases towards unvisited energy states. This dynamics does not satisfy detail balance exactly since the transition rate is fluctuating. However, test on small systems shows that the sample averages do converge to the exact value with errors proportional to the inverse of square root of Monte Carlo steps.

A two-stage simulation will guarantee detail balance. Stage one is the same as described above. In stage two, we adjust the approximate transition matrix obtained in stage one such that detail balance is satisfied exactly. In this stage, the flip rate is fixed and will not fluctuate with the simulation. The second stage simulation is dynamically equivalent to Berg’s multicanonical or Lee’s entropy sampling dynamics. The stage two can be iterated so that the simulated ensemble approaches the ideal multicanonical ensemble, but we found two-stage or even a single stage simulation already gives excellent results.

The simulation can also be combined with the N-fold way  with little overhead in computer time, since the quantity $`N(\sigma ,\mathrm{\Delta }E)`$ needed in the N-fold way is already computed. In addition, not only we can have equal histogram (multicanonical), we can also generate “equal-hit” ensemble, where each energy of distinct states is visited equally likely .

## 5 Density of states from transition matrix

The density of states is related to the transition matrix $`T_{E,\mathrm{\Delta }E}=N(\sigma ,\mathrm{\Delta }E)_E/N`$ by Eq. (4), where $`N`$ is the total number of possible moves. Since there are more equations than unknown $`n(E)`$, we use a least-square method to obtain “optimal” solution. Let $`S(E)=\mathrm{ln}n(E)`$, we consider

$$\mathrm{minimize}\underset{E,E^{}}{}\frac{1}{\sigma _{E,E^{}}^2}\left(S(E^{})S(E)\mathrm{ln}\frac{T_{E,\mathrm{\Delta }E}}{T_{E^{},\mathrm{\Delta }E}}\right)^2$$
(8)   
subject to all the conditions known. For example, for the Ising model, we have $`n(E_{min})=n(E_{max})=2`$, and $`_En(E)=2^N`$. The variance $`\sigma ^2`$ is the variance of the quantity $`\mathrm{ln}T_{E,\mathrm{\Delta }E}/T_{E^{},\mathrm{\Delta }E}`$ obtained from sets of Monte Carlo data. It is also possible to work with the matrix $`T`$ directly with the conditions that $`T`$ must satisfy.

## 6 Conclusion

We proposed an algorithm which samples energy $`E`$ uniformly. Comparing to multicanonical simulation, the method offers a very easy way of starting the simulation. The dynamic characteristics are similar to well-converged multicanonical Monte Carlo dynamics. For example, the tunneling time for 10-state Potts model in two dimensions is about $`\tau L^{2.6}`$ for $`L\times L`$ system. It is very easy to combine statistics from several simulations, including parallel simulations. It is an efficient method for computing density of states and all the related thermodynamic quantities.

## Acknowledgements

The work presented here is in collaborations with R.H. Swendsen, T.K. Tay, L.W. Lee, and Z.F. Zhen.