# Monte Carlo Algorithm for Free Energy Calculation

## Abstract

We propose a new Monte Carlo algorithm for the free energy calculation based on configuration space sampling. Upward or downward temperature scan can be used to produce . We implement this algorithm for Ising model on square lattice and on triangular lattice. Comparison with the exact free energy shows an excellent agreement. We analyse the properties of this algorithm and compare it with Wang-Landau algorithm which samples in energy space. This method is applicable to general classical statistical models. The possibility of extending it to quantum systems is discussed.

###### pacs:

05.10.Ln, 75.10.Hk, 02.70.Rr## I Introduction

Monte Carlo (MC) simulation is one of the most important numerical methods for solving statistical problems in fields such as chemistry, biology and physics. In condensed matter physics, MC is extensively used to study the properties of many statistical models, phase transitions, and quantum many-body systems [(1)]. Often, besides the expectation values of certain physical quantities, one also needs the free energy of the system in thermal equilibrium. Calculation of free energy is a difficult problem for traditional Metropolis algorithm since partition function plays the role of the normalization constant in the thermal probability density distribution of an ensemble.

In the past three decades, various MC algorithms have been proposed to do the free energy calculation for statistical Hamiltonians. For classical systems, the frequently used methods are the histogram reweighting method [(2)], multiple histogram [(3); (4); (5)], transition matrix MC [(6)], entropic sampling[(7)], flat histogram method [(8)], and the Wang-Landau method [(9)]. Recently Wang-Landau method has been improved in different aspects [(10),(11)]. All these methods have their respective advantages and disadvantages. For examples, the histogram reweighting method produces the energy histogram at a given temperature and employs the reweighting method to recover the distribution at a different temperature. As is shown below, usually the canonical distribution is sharply peaked around which is -dependent. Thus the error of free energy becomes large when is large. For both the entropic sampling and the Wang-Landau method, the ensemble is created in the energy space instead of in the configuration space. This helps to obtain the density of states efficiently, but it is less convenient if one hopes to study other physical quantities in a single simulation, especially when these quantities are not simple functions of energy , such as the correlation function.

For quantum systems, the calculation of free energy is also a pertinent but harder problem. In this regard, the Wang-Landau algorithm has been combined with statistical series expansion method to calculate the free energy of quantum Hamiltonians such as the Heisenberg model [(12)]. The idea of flat histogram has been applied to the diagrammatic MC method to improve the long imaginary-time results [(13)], and to calculate the grand potential of a cluster system with electron bath [(14); (15)]. However, considering that the flat-histogram method or Wang-Landau algorithm have not been implemented under the path integral(PI) quantum Monte Carlo (QMC) methods such as the determinantal QMC [(16); (17)] and the continuous-time PI QMC method [(18)], which are based on the Metropolis sampling in configuration space, it is still desirable to develop a free energy calculation method which can calculate free energy by the configuration space sampling.

It is the purpose of this paper to propose such a new MC algorithm that can calculate the free energy using the configuration-based sampling algorithm. The price that we have to pay is a sequential scan from low to high temperatures, or vice versa. We demonstrate the implementation of this algorithm using the Ising model on square as well as on triangular lattice, for which the exact solutions are known. Comparison with Wang-Landau method shows that both efficiency and accuracy of this method are satisfactory. The additional advantage of this method is that it is based on Metropolis algorithm and hence in principle it can be extended to quantum systems within determinantal or path integral methods.

## Ii Method and Results

In this part, we demonstrate the implementation of our method and analyse its features using the two dimensional Ising model on a square lattice. For comparison purposes, here we use the equivalent Hamiltonian of the two-state Potts model

(1) |

Here, is the spin degrees of freedom on site and it takes integer values from to . if and if . The summation is for pairs of nearest neighbour sites on a square lattice with geometry. This model has been studied extensively as a basic statistical model [(19)]. Its exact critical transition temperature is . In the following we use as the energy unit and set the Boltzmann constant for convenience.

One of the widely used MC algorithms for studying classical statistical models is the Metropolis sampling in configuration space [(20)]. In this algorithm, one begins by choosing a random spin configuration (here we use capital letter to denote the spin configuration of the whole lattice), and update the configuration according to a given proposal probability and an accepting probability . The transition probability must satisfy the detailed balance condition to guarantee that the resulting Markovian chain has the target distribution in the equilibrium limit. For the thermodynamical calculation, we use the Bolzmann distribution as our target distribution,

(2) |

is the partition function . Here is the inverse temperature. After the Markovian chain reaches equilibrium, sampling on this chain can produce the required statistical averages.

However, the free energy cannot be calculated directly, because the partition function is a normalization factor of the probability distribution and hence cannot be treated as a statistical average. To calculate , usually one either employs the concept of energy histogram [(2)] within MC method, the maximum entropy method [(21)], or by numerical integration over the derivative of free energies [(22)]. In the following, we propose a new method which combines the idea of energy histogram and the configuration-space sampling to calculate the free energy over full temperature range.

In the Metropolis algorithm, the energy probability distribution produced by the Markovian chain is

(3) |

where is the degeneracy of energy level of the given Hamiltonian. can be estimated approximately from the energy histogram of the Markovian chain, that is,

(4) |

Here is the number of spin configurations with energy and is the total number of sampled configurations in the Markovian chain. The precision of the above estimation increases as increases. In the limit , we get , and hence

(5) |

In principle, can be calculated from the values of and at any given energy . For a large variety of classical Hamiltonians, the ground state degeneracy is easy to obtain. For the Ising model on square lattice, for an example, . is thus in principle obtainable from the ground state energy histogram .

In practice, however, the above simple scheme does not work at arbitrary because the energy distribution is sharply peaked at an energy (or at many different ’s for some models) which is an increasing function of . In Fig.1, we show for the Ising model on the square lattice with , at different temperatures. Here is the energy per site and . is plotted in logarithmic scale and it decays very fast away from the peak position. At high temperatures, is peaked at high energy and is so small that it is impossible to calculate it accurately from , because at high temperatures the latter is practically zero for finite .

One way to overcome the above difficulty is to use Eq.(3) at a different temperature ,

(6) |

Using the fact that is -independent, the ratio of Eq.(6) and Eq.(3) gives

(7) |

which can be used to evaluate the ratio of partition function at and as

(8) |

This strategy works as long as for some energy , can be measured accurately both at and . It is the basis for the histogram reweighting method [(2)] and the multiple histogram method [(3); (4); (5)] to produce the free energy up to a constant factor. This constant factor can be further fixed by using the infinite temperature partition function [(4)] or the integral of density of states (Ref. (5)), both of which are equal to the total number of configurations.

Here, we will take a different strategy, which is based on Eq.(3) at the same temperature but at a different energy , . From the ratio of it and Eq.(3) one obtains

(9) | |||||

This means that the knowledge of can be transferred from one energy to all , where is the energy window in which is reasonably large. As shown in Fig.1, the center of this window increases with temperature .

For those Hamiltonians where the ground state degeneracy is known, to obtain at higher energies, we need to increase from a small value in such small steps that the energy windows and of adjacent temperatures and have significant overlap. Suppose one knows for a certain energy . one does the MC calculation at to produce the histogram for each . Using the known value, is obtained from Eq.(5) and for each is obtained from Eq.(9). For the next temperature , one calculate the histogram for each . We choose an energy . Using value obtained in the calculation and , is obtained from Eq.(5) and for each is obtained from Eq.(9). This process goes on until the desired high temperature is reached. In this way, for larger and larger energies and at successively higher temperatures can be obtained. We call this scheme the upward temperature scan.

One could also start from the high temperature limit , and employ the direct sampling without Boltzmann factor to produce for energies . By decreasing from step by step, one could reach low and calculate for all temperatures. This scheme is called downward temperature scan.

In the implementation of the above algorithm, an important technical issue is how to select the common energy point . In our calculation, we use the crossing energy , which is determined by . It is the energy where the curves of adjacent temperatures and crosses. The value chosen in this way has the largest value for both and , hence guarantees the optimal precision.

After MC calculation at () is carried out, the rest calculations are done at equal distance temperatures () () for upward (downward) scan. For the upward scan, we use and for downward scan, . We carefully control the interval to reach the optimal precision. A too large will lead to a small and relatively large error in . If is too small, the number of temperatures required in scan from to the final will increase. It leads to longer calculation time and larger accumulated error in the transfer. Therefore, a suitable should be found by testing calculations. This issue is discussed in below (see Fig.6).

### ii.1 upward temperature scan

In our benchmark calculation for the Ising model on a square lattice, we use the cluster update scheme of Wolff [(24)]. It has a relatively high updating speed and weak critical slowing down near the critical temperature. The free energy calculation algorithm described above can also be used with local update algorithm without modification. In Fig.1, we show the normalized energy probability distribution obtained from the Markovian chain, for (in Fig.1(a)) and for relatively high temperatures (in Fig.1(b)). Here is the energy per site. At low , has a peak at and its width broadens with increasing . While for , the peak position begins to increase with temperature and its width gets saturated. The peak position as a function of is shown in the inset of Fig.1(a). It moves from at to in large limit. Note that the sharp increase occurs near , close to the of order-disorder phase transition of this system in the thermodynamic limit [(23)]. Significant overlap in the peak energy windows for adjacent temperatures and is crucial for our algorithm to work. In this work, we use a uniform temperature mesh and choose such that the overlap of the peaks are large enough to guarantee the precision.

In Fig.2(a), the free energy per site obtained from the upward scan for is shown. It is compared with the exact result [(23)]. Except stated otherwise, we use samples to calculate for each temperature . With , producing in the range , requires temperature values . In Fig.2, our results are indistinguishable from the exact curve on the scale of the figure. In the inset, the standard variance is shown as functions of . For each temperature, is measured from independent data of , each of which is produced by sampling. The actual error has same order of magnitude as and it is always less than . Both quantities are small in and increase linearly with temperature in . is close to the order-disorder transition occurring at in the thermodynamic limit. The linear increase of error with in the disordered phase is related to the constant relative standard variance of the partition function, as discussed in Fig.3 below. For the highest temperature that we study , the actual error is smaller than . The relative error first increases with and then saturates to in .

The same comparison is made for larger lattice size in Fig.2(b). The results are similar to and the agreement between MC and the exact result is excellent. The main difference from the case is that the standard variance is larger than the actual error . As temperature increases, an abrupt increase of and occurs at and the linear behavior occurs at . These features are same as case. For the relative error, a similar saturation in is observed in high temperature to about .

### ii.2 downward temperature scan

In Fig.3(a), we show per site as function of obtained from the downward scan method, for the same Ising model on the square lattice. It is noted that the numerical error obtained from downward scan is of the same order as the upward scan in Fig.2. Compared to Fig.2, the and have a sharp peak around , instead of a kink. In order to understand the linear -dependence of the errors in Fig.2 and its difference with Fig.3(a), in Fig.3(b), we show the relative standard variance of the partition function , which is related to the variance of by . It is seen that for the upward scan, surges at by one magnitude, from smaller than below to about above . Further increasing temperature it does not change much. For the downward scan, is small for and surges at to about . It stays at this value to . Clearly, in both cases, the transfer error of from one to another is largest around . As shown in the inset of Fig.1(a), the peak position in curve moves fastest near . For a uniform temperature mesh and constant peak width, the overlap of is therefore smallest for and close to , leading to the largest accumulation of error in . Such behavior of as a function of translates into the the linear behavior of for away from and the kink or the sharp peak for close to , as shown in the insets of Fig.2 and Fig.3(a).

### ii.3 Ising model on triangular lattice

To show that our method works also for other systems, especially the frustrated system, we apply the downward temperature scan to the Ising model on a triangular lattice. This system has been solved exactly [(25); (26)]. The frustration induces huge low energy degeneracy and makes this system particular interesting.

In our calculation, we consider a square lattice with the nearest neighbour coupling as well as the next nearest neighbour coupling in the direction. Both couplings are antiferromagnetic. The Hamiltonian is the same as Eq.(1), but with . Open boundary condition is used. In Fig.4(a), per site is shown as a function of . Compared to the square lattice case, except for the similarity in the shape of curve, there are two important differences. One is that in the low temperature limit, approaches instead of as in the square lattice case. This shows that frustration increase the free energy of the spins. The second difference is that is more smooth and there is no transition point associated with a finite temperature phase transition. This is consistent with the fact that there is no finite temperature transition in the triangular lattice Ising model [(25); (26)]. In the inset of Fig.4(a), the variance of decreases smoothly as is lowered and increases sharply at very low temperatures. It remains to be elucidated whether this behavior is related to the singular ground state correlation of the Ising model on a triangular lattice [(25)].

Due to the difference in Hamiltonian definition, it is difficult to compare our MC free energy with the exact result in Ref. (25). In Fig.4(b), we show the entropy per site as functions of for various lattice size and . It is calculated from the energy and via . In the high temperature limit, approaches (not shown). In the zero temperature limit, our results for system and from exact enumeration shows that , being consistent with the observed two fold degeneracy in the ground state. However, for a low but finite temperature, entropy increases with , as shown in Fig.4(b). For , entropy approaches at , the lowest temperature that we study. This value is already very close to the exact in the thermodynamic limit [(25)]. This supports the reliability of our calculation of free energy for this system. Note that for a given , the error bar increases with lowering temperature decreases mainly due to the sharp increase of error in in downward T-scan. For a fixed temperature, the error bar decreases with increasing because of the factor, showing that our method is robust for large systems.

## Iii Discussions

In this section, we discuss the calculation error of the free energy and compare it with the Wang-Landau’s algorithm. To see the dependence of error on system size , sampling number , and the temperature mesh , we take the data from upward scan for square lattice Ising model with . First, we present the -dependences of and in Fig.5. It is seen that for all the calculated size , and are on the same order, all smaller than . This shows that the efficiency of our algorithm does not deteriorate with increasing , at least for . The relative magnitude of and may vary for different ’s, but the actual error is always within . In the inset of Fig.5, we show the full width at half maximum of as a function of at . It is observed that the peak width scales as . This reminds us that for very large , the curves will be very sharp and we have to use a denser -mesh. At and , however, the actual effective peak width is of the same order as the peak distance even in the worst case of vicinity of . Therefore, we can still reach similar precision as for smaller systems. For even larger size, needs to be scaled as to maintain the balance of peak width and peak distance, and thus maintain the calculation precision.

In Fig.6, we show the standard variance as functions of sampling number for various values. Our results on square lattice shows that with the average , close to the expected from the central limit theorem. For fixed , has a weak dependence on in the range . As shown in the inset of Fig.6, is approximately a parabolic function of , with the minimum reached at a -dependent . For all values we used, we find that the smallest error is reached at values around .

We make quantitative comparison with Wang and Landau’s results [(9)], although this method may not be the most accurate one of free energy calculation [(6)]. We find that when scaled to the same values, our result of has an error about one magnitude larger than that in Ref.(9). Simple optimization of our method can reduced the error significantly, such as using denser T-mesh near and sparser one away , or using upward T-scan for and downward T-scan for .

Unlike the Wang-Landau algorithm which generates random walks in energy space, our algorithm samples directly in the configuration space. As a result, this method could be fitted into certain MC simulation of quantum systems. In the PI QMC[(18)] and the determinantal QMC [(17)] methods, the partition function of a given quantum Hamiltonian is expressed by the summation over configurations of classical auxiliary fields. The proposed algorithm can then be used to calculate , using the same Markovian chain as used for evaluating general expectation values.

One example of the application of appears in the study of the Mott metal-insulator transition [(27)] in the half-filled Hubbard model using the dynamical mean-field theory (DMFT) [(28); (29)], which is exact in infinite spatial dimensions. The transition from the Fermi liquid state in small regime into the Mott insulator in large regime was found to be a special second order phase transition at and a first order one at . To determine the actual transition line at , one needs to compare the free energy of the two coexisting phases within the two spinodal lines. Within DMFT, this task is reduced to the evaluation of free energy of the effective Anderson impurity model [(29)]. This proves to be a difficult problem for QMC methods such as the Hirsch-Fye algorithm [(30); (31)]. Recently, the grand potential of the cluster problem is calculated by Wang-Landau method combined with the continuous time QMC [(14)], and it is used to calculate the grand potential of lattice model within cluster dynamical mean-field theory. It is an interesting topic to apply our algorithm in various QMC methods to handle similar problems. Work in this direction is under progress.

## Iv Acknowledgements

We are grateful to the helpful discussions with Prof. You-Jin Deng. This work is supported by 973 Program of China (2012CB921704), NSFC grant (11374362), Fundamental Research Funds for the Central Universities, and the Research Funds of Renmin University of China.

### References

- D. P. Landau and K. Binder, A Guid to Monte Carlo Simulations in Statistical Physics, 4th ed., Cambridage University Press, 2014.
- A. M. Ferrenberg and R. H. Swendsen, Phys. Rev.Lett. 61, 2635 (1988).
- A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 63, 1195 (1989).
- J. P. Valleau and D. N. Card, J. Chem. Phys. 57, 5457 (1972).
- N. A. Alves, B. A. Berg, and R. Villanova, Phys. Rev. B 41, 383 (1990).
- J. S. Wang and R. H. Swendsen, J. Stat. Phys. 106, 245 (2002).
- J. Lee, Phys. Rev. Lett. 71, 211 (1993).
- B. A. Berg and T. Neuhaus, Phys. Lett. B 267, 249 (1991).
- F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001); Phys. Rev. E 64, 056101 (2001).
- For examples, see R. E. Belardinelli and V. D. Pereyra, Phys. Rev. E 75, 046701 (2007); R. E. Belardinelli, S. Manzi, and V. D. Pereyra, Phys. Rev. E 78, 067701 (2008).
- T. Vogel, Y. W. Li, T. Wüst, and D. P. Landau, Phys. Rev. Lett. 110, 210603 (2013).
- M. Troyer, S. Wessel, and F. Alet, Phys. Rev. Lett. 90, 120201 (2003).
- N. G. Diamantis and E. Manousakis, Phys. Rev. E 88, 043302 (2013).
- G. Li, W. Hanke, A. N. Rubtsov, S. Bse, and M. Potthoff, Phys. Rev. B 80, 195118 (2009).
- E. Gull et al., Rev. Mod. Phys. 83, 349 (2011).
- M. Suzuki, Prog. Theor. Phys. 65, 1454 (1976).
- J. E. Hirsch, Phys. Rev. B 31, 4403 (1985).
- B. B. Beard and U. J. Wiese, Phys. Rev. Lett. 77, 5130 (1996).
- F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982).
- N. Metropolis et al., J. Chem. Phys. 21, 1087 (1953).
- C. Huscroft, R. Gass, and M. Jarrell, Phys. Rev. B 61, 9300 (2000).
- See, for an example, N. H. Tong, S. Q. Shen, and F. C. Pu, Phys. Rev. B 64, 235109 (2001).
- A. E. Ferdinand and M. E. Fisher, Phys. Rev. 185, 832 (1969).
- U. Wolff, Phys. Rev. Lett. 62, 361 (1989).
- G. H. Wannier, Phys. Rev. 79, 357 (1950).
- R. M. F. Houtappel, Physica 16, 425 (1950).
- F. Gebhard, The Mott Metal-Insulator Transition, Springer-Verlag Berlin Heidelberg, 1997.
- W. Metzner and D. Vollhardt, Phys. Rev. Lett 62, 324 (1989).
- A. Georges, G. Kotliar, W. Krauth, M. Rozenberg, Rev. Mod. Phys. 68, 13, (1996).
- J. E. Hirsch and R. M. Fye, Phys. Rev. Lett. 56, 2521 (1986).
- N. Bluemer, Mott-Hubbard Metal-Insulator Transition and Optical Conductivity in High Dimensions, Shaker Verlag, Aachen, 2003.