# Monte Carlo Algorithms based on the Number of Potential Moves

## 1 Introduction

The traditional Monte Carlo method applied to statistical physics  is mostly a sampling method to generate standard statistical ensembles, e.g., the canonical ensemble or microcanonical ensemble. In recent years, other ensembles have been used which do not correspond to thermodynamically meaningful ensembles, but used only as a vehicle for computing quantities of interests by Monte Carlo method. The earliest such method is the umbrella sampling . Other important recent developments are the multi-canonical simulations , $`1/k`$-sampling , and broad histogram methods . According to one definition , the multi-canonical ensemble is an ensemble that the probability $`P(E)`$ of having energy $`E`$ is a constant. This may be realized in a piecewise fashion. In the $`1/k`$-sampling method, the probability for a state having energy $`E`$ is given by $`1/_{E^{}<E}n(E^{})`$, where $`n(E)`$ is density of states. The broad histogram dynamics does not have a well characterized distribution, but $`P(E)`$ is much broader than the canonical ensemble. The canonical distribution is well approximated by a Gaussian function. It was also pointed out that this dynamics is not entirely correct . The ultimate goal of generating these distributions is usually to compute thermodynamic averages, which are mostly canonical averages. Some form of reweighting is then used to obtain the desired distribution.

Recently, we proposed a dynamics  which can generate a flat histogram $`P(E)=const`$. This method is exact when a “self-consistency” is achieved. The meaning of which will be made clear later. Similar to the broad histogram method, the central quantity is $`N(\sigma ,\mathrm{\Delta }E)_E`$ , the microcanonical average of the number of ways to move from one energy level $`E`$, to a nearby energy level, $`E+\mathrm{\Delta }E`$. We can then construct either the density of states or the canonical distribution at any temperature. The canonical distribution is determined from an artificial dynamics which we call it transition matrix Monte Carlo. We discuss these methods and present some preliminary results in the later section.

## 2 Sampling the Inverse Density of States

We illustrate the method using a two-dimensional Ising model on a square lattice as an example. First of all, we choose a type of permissible moves. For purpose of connection with standard single-spin-flip Glauber (or Metropolis) dynamics, we take the set of moves to be all single-spin flips. For a given state $`\sigma `$, we can obtain $`N=L^2`$ new states through flipping each of the spins in an $`L\times L`$ system. If the original state has energy $`E=E(\sigma )`$, the new state may have energy $`E+\mathrm{\Delta }E`$. Since the energy spectrum is discrete, we have only a finite number of possibilities for the new energies; for the two-dimensional Ising model, we have five possible energy changes, $`\mathrm{\Delta }E=0,\pm 4J,\pm 8J`$. Let the counts of moves for each energy changes be $`N(\sigma ,\mathrm{\Delta }E)`$. Hence the total number of moves is $`_{\mathrm{\Delta }E}N(\sigma ,\mathrm{\Delta }E)=N`$.

Following the argument of Oliveira , we consider two energy levels $`E`$ and $`E^{}=E+\mathrm{\Delta }E`$. Each move from the state $`\sigma `$ of energy $`E`$ to the state $`\sigma ^{}`$ of energy $`E^{}`$ is through a single spin flip and the reverse move is also allowed. Thus, the total number of moves from all the states with energy $`E`$ to $`E^{}`$ is the same as from $`E^{}`$ to $`E`$:

$$\underset{E(\sigma )=E}{}N(\sigma ,\mathrm{\Delta }E)=\underset{E(\sigma ^{})=E+\mathrm{\Delta }E}{}N(\sigma ^{},\mathrm{\Delta }E).$$
(1)   
The microcanonical average of a quantity $`A(\sigma )`$ is defined as

$$A_E=\frac{1}{n(E)}\underset{E(\sigma )=E}{}A(\sigma ),$$
(2)   
where the summation is over all the states with a fixed energy $`E`$, and $`n(E)`$ is the number of such states. In terms of microcanonical average, we can re-write Eq. (1) as

$$n(E)N(\sigma ,\mathrm{\Delta }E)_E=n(E+\mathrm{\Delta }E)N(\sigma ^{},\mathrm{\Delta }E)_{E+\mathrm{\Delta }E}.$$
(3)   
This is the basic equation of the broad histogram method  and is also our starting point for a flat histogram sampling algorithm.

Consider the following flip rate for a single-spin-flip move from state $`\sigma `$ to $`\sigma ^{}`$ with energy $`E`$ and $`E^{}=E+\mathrm{\Delta }E`$, respectively:

$$r(E^{}|E)=\mathrm{min}(1,\frac{N(\sigma ^{},\mathrm{\Delta }E)_E^{}}{N(\sigma ,\mathrm{\Delta }E)_E}).$$
(4)   
The site of the spin flip is chosen at random. Then the detailed balance condition for this rate

$$r(E^{}|E)P(\sigma )=r(E|E^{})P(\sigma ^{})$$
(5)   
is satisfied for $`P(\sigma )1/n(E(\sigma ))`$. Thus energy histogram is flat,

$$P(E)=\underset{E(\sigma )=E}{}P(\sigma )n(E)\frac{1}{n(E)}=const.$$
(6)   
Suppose that such samples are generated, then in some sense, it is the optimal ensemble for evaluation of $`N(\sigma ,\mathrm{\Delta }E)_E`$. This is because for different $`E`$, we take samples uniformly in $`E`$, and thus the relative errors in $`N(\sigma ,\mathrm{\Delta }E)_E`$ are about the same for all $`E`$.

Since $`N(\sigma ,\mathrm{\Delta }E)_E`$ is not known in general, we cannot start the simulation unless an approximation scheme is used. We can think of the process as finding fixed point value of the system $`x=f(x)`$, where vector $`x`$ represents the whole set of $`N(\sigma ,\mathrm{\Delta }E)_E`$ values. While the function $`f`$ can be evaluated, its explicit form is not known. Some iterative scheme may be useful to speed-up the convergence. To start the iterative process, we use a cumulative average for the true microcanonical average. For those $`E`$ which we do not have any sample yet, we simply set $`r(E^{}|E)`$ to 1. This simple scheme is very good for small systems even without iteration.

## 3 The transition matrix Monte Carlo dynamics

We can construct a Monte Carlo dynamics, in the space of energy, with the average number of moves, $`N(\sigma ,\mathrm{\Delta }E)_E`$, . Let us look at a single-spin-flip Glauber dynamics. Suppose we do not care about the spin states and only want to know the change of energy. The rate of a spin flip is given by the Glauber rate,

$$w(\mathrm{\Delta }E)=\frac{1}{2}\left[1\mathrm{tanh}\frac{\mathrm{\Delta }E}{2kT}\right].$$
(7)   
Since there are (on average) $`N(\sigma ,\mathrm{\Delta }E)_E`$ different ways of going from $`E`$ to $`E^{}=E+\mathrm{\Delta }E`$, the total probability for transition from $`E`$ to $`E^{}`$ is

$$W(E^{}|E)=w(\mathrm{\Delta }E)N(\sigma ,\mathrm{\Delta }E)_E,EE^{}.$$
(8)   
The diagonal elements are fixed by the requirement that $`W(E|E^{})`$ is a stochastic matrix. This transition matrix satisfies detailed balance with respect to the canonical distribution, $`P_T(E)n(E)\mathrm{exp}(E/kT)`$. Thus the stationary distribution is the canonical distribution. This new dynamics in the space of energy $`E`$ is related to the single-spin-flip dynamics by 

$$W(E^{}|E)=\frac{1}{n(E)}\underset{E(\sigma )=E}{}\underset{E(\sigma ^{})=E^{}}{}\mathrm{\Gamma }(\sigma ^{}|\sigma ),$$
(9)   
where $`\mathrm{\Gamma }(\sigma ^{}|\sigma )`$ is the transition matrix of the single-spin-flip dynamics.

An interesting aspect of this dynamics is that it has a much reduced critical slowing down. In fact, one can show that the relaxation time at the critical point $`T_c`$ is proportional to the specific heat. Thus for the two-dimensional Ising model, the divergence of the relaxation time is only logarithmic. In one dimension, the dynamics has a curious dynamical critical exponent of $`z=1`$ as oppose to 2 for the local dynamics and 0 for the Swendsen-Wang dynamics . Since the dynamics can not be realized without first knowing the values $`N(\sigma ,\mathrm{\Delta }E)_E`$, the real usefulness is in the construction of canonical distribution from the samples obtained by flat histogram or any other algorithms that can compute $`N(\sigma ,\mathrm{\Delta }E)_E`$ accurately.

## 4 Results

In Fig. 1, we show the energy histograms for three different types of dynamics of the $`32\times 32`$ two-dimensional Ising model, (a) the Gaussian-like peak for the standard canonical ensemble at the critical temperature $`T_c`$; (b) the broad histogram dynamics with a sharp peak near $`E=0`$; (c) the flat histogram method with an insert showing the fluctuation on a fine scale.

Given the estimates for $`N(\sigma ,\mathrm{\Delta }E)_E`$, there are a number of ways to determine the canonical distribution, $`P_T(E)n(E)\mathrm{exp}(E/kT)`$. For example, we can use Eq. (3) to determine the density of states. We can also determine $`P_T(E)`$ directly from the detailed balance of the transition matrix Monte Carlo dynamics :

$$w(\mathrm{\Delta }E)N(\sigma ,\mathrm{\Delta }E)_EP_T(E)=w(\mathrm{\Delta }E)N(\sigma ^{},\mathrm{\Delta }E)_{E+\mathrm{\Delta }E}P_T(E+\mathrm{\Delta }E),$$
(10)   
where $`w(\mathrm{\Delta }E)`$ is given by Eq. (7). Since there are more equations than unknowns, it is natural to solve these over-determined equations with least-square method. However, a more direct iterative scheme is also quite accurate and more efficient.

In Fig. 2, we show the specific heat (upper part) and relative errors as compared with exact results . The dash lines are for the broad histogram method and solid lines are from the flat histogram sampling. The broad histogram method shows an anomalous peak around $`T=1.3`$, while the flat histogram result agrees with exact values with errors of $`10^2`$ or less.

Since we can compute the density of states $`n(E)`$ easily, we can also compute free energy and entropy with ease. These quantities are more difficult to compute by the traditional methods. Fig. 3 shows the entropy and errors. The flat histogram is again better than the broad histogram method.

All approaches that use reweighting technique, such as the histogram methods of Ferrenberg and Swendsen , Lee’s version of multicanonical method , or the broad histogram method , have the problem of scalability for large systems. Our flat histogram method also suffers from this. While the simple method without an iterative process and without requirement for self-consistency seems to work well for systems $`L32`$, systematic errors are observed for large systems. Substantial deviations (extra anomalous peaks in the specific heat, for example) are present for the $`L=64`$ systems. Such systematic deviations can be measured quantitatively by what we called detailed balance violation :

$$v(E)=\left|1\frac{g(E,E^{\prime \prime })g(E^{\prime \prime },E^{})g(E^{},E)}{g(E,E^{})g(E^{},E^{\prime \prime })g(E^{\prime \prime },E)}\right|$$
(11)   
where $`g`$ is generally a transition rate—for our problem here, we’ll take $`g(E,E^{})=N(\sigma ,E^{}E)_E`$, with $`E^{}=E+4J`$ and $`E^{\prime \prime }=E+8J`$. The quantity $`v(E)`$ should be zero, up to the usual Monte Carlo statistical errors, if the estimates are not systematically biased. In Fig. 4, we show this quantity as a function of $`E`$ for the $`L=64`$ system. The largest violation occurs at the two ends of the distribution. This systematic trend is also present for small systems.

There are a number of ways to fix this problem. One is to do a number of canonical simulations at lower temperatures where the violation of detailed balance is biggest. This indeed proves to be effective for the Ising model. However, this solution is not very satisfactory, as such simulations may be very difficult, for example, for spin glasses. Thus, a more systematic approach is to use an iterative scheme, which can hopefully converge to the true value without any systematic bias.

## 5 Conclusion

We study a recently proposed Monte Carlo dynamics in which the energy histogram is exactly flat in principle. We demonstrated that such method is capable of giving highly accurate results for the thermodynamic quantities in a single or few simulations for the whole temperature region. While some systematic errors are present in our current simple implementation, there are ways to improve the naive algorithm. We expect that this method will be a useful alternative for thermodynamic calculations, especially for free energy and entropy calculations.