# Cluster-Flip Algorithm

The central idea comes from the following observation. In the vicinity of the transition, the correlation length is very long (see previous sections). Large clusters of parallel spins are formed just above Tc (we are in the case of a ferromagnetic crystal). The

Metropolis algorithm which updates spin after spin will take a long time to consider unnecessarily all spins in a cluster because, except for a small fraction at spins at the boundary, the spins in a cluster do not need to change orientation (they are ready for the ferromagnetic state at Tc). Therefore, to perform single-spin flips for all spins in a cluster near Tc is a loss of time. In addition, due to the existence of large clusters near Tc, the relaxation time is very long [see Eq. (6.40)]. This phenomenon is called "critical slowing-down.” Wolff [364] and Swendsen and Wang [337] have proposed to update simultaneously all spins of a cluster by flipping the entire cluster (we have in mind the case of the Ising model). The method is simple to implement. To simplify the description, we consider a system of Ising spins with a ferromagnetic interaction J between nearest neighbors. The principle of the algorithm is the following:

• (1) For a given spin configuration, we consider a spin 5, and we "construct" a "cluster" around S, as follows. We examine the neighboring spins: If a neighbor in one direction is parallel to Sj, then it belongs to the cluster with a probability p = 1 - exp(-2fij) where ft = (kBT)_1. We consider the spin next to that spin in that direction, and we continue the cluster construction. The limit of the cluster in the considered direction is where the cluster encounters an antiparallel spin or if a random number taken between 0 and 1 is larger than p. We have to go to all directions to determine the boundary of the cluster. Note that there is a very efficient algorithm for the cluster construction proposed by Hoshen and Kopelman [156].
• (2) We flip the cluster and we calculate Д£ the difference in energy with the previous state: AE = 2/[C(+-F) - C(—1-)] where C(++) the number of broken parallel links along the boundary of the cluster, and C(—1-) the number of antiparallel links along the boundary. An example is shown in Fig. 6.6. The new orientation of the cluster is accepted or not, following the Metropolis criterion as for a single spin.
• (3) We take another spin outside the above cluster and we begin again another cluster construction, etc.

The difference between the method by Wolff and that by Swendsen-Wang is the following. In the former, we flip only large

Figure 6.6 Example of a cluster, limited by the discontinued contour, constructed by cluster-construction algorithm. Black circles = t spins, white circles = | spins. The number of broken parallel links C(+-b) is 1, that of antiparallel spin links C(—h) is 11.

clusters, while in the latter one we flip all clusters. It is true that near Tc, the two methods are equivalent because most of clusters are large. However, a little bit further from Tc where there are many small clusters, we spend much time to flip small clusters in the Swendsen-Wang method, this does not significantly improve the result.

In practice, we can combine the Metropolis algorithm and the cluster-flip algorithm in one simulation: We use the latter from time to time because the cluster construction takes time. We need it only very near Tc to get rid of the critical slowing-down. The result is striking: The dynamic exponent z ~ 2 as obtained by the Metropolis algorithm alone becomes z ~ 0.5 using the cluster-flip method near Tc.

# Histogram Method

In standard Metropolis Monte Carlo simulations, we calculate average physical quantities at discrete temperatures. We extrapolate results between discrete temperatures. However, near the transition temperature, extrapolation is impossible because physical quantities diverge at a precise single temperature value. Practically, we cannot find the exact location of the transition temperature using discrete temperatures. If the chosen temperature is not exactly Tc, we will not have the precision on the critical exponents calculated with the heights of Cv and у which are not at Tc.

To avoid this difficulty, Ferrenberg and Swendsen [110] have proposed the histogram method which consists of making a simulation at a temperature T0 as close as possible to Tc and recording the instantaneous system energy E during the simulation as long as possible to establish a histogram H[E). Using this histogram we can calculate the canonical probabilities P[T, E) at neighboring temperatures T around T0 at as many points as we wish. Using these probabilities P[T, E), we can calculate average values of physical quantities as continuous functions of T around T0, by the formulas of the canonical description (see Appendix A), not by simulations. It suffices to choose T0 close to Tc (not necessarily at Tc since Tc is not known before the simulation), we can find the exact location of the maximum of Cv and y, for example. How to choose T0? The answer is we have to make a preliminary run with many discrete temperatures using the Metropolis algorithm. From these preliminary results, we take T0 at the maximum of Cv or у for the histogram run.

The method is described in the following. We consider the partition function at T0 (see Appendix A):

where the sum is taken over all energies of microscopic states, W[E) denotes the degeneracy of the energy level E independent of T0. The probability of the level E at T0 is then

Now, by simulation we obtain the energy histogram H[E) at T0. The probability P(T0, E) is nothing but

where NMc is the number of Monte Carlo steps used to establish H (E). By comparison of (6.50) to (6.49) we have

We consider now a temperature T near T0. The probability P[T, £) is written by

where we have used (6.51) to replace W[E) and the following relation to replace Z(T):

The histogram H[E) established for T0 is thus used to calculate the probability at another temperature T by (6.52). Using this probability we can calculate without simulation the average value of a quantity A by the formula

Remarks:

• (1) To get a precise H(E), we have to use a veiy large number of Monte Carlo steps N^c in order to include as many as possible microscopic states in the sum.
• (2) Since H(E) is obtained with an importance sampling at T0, if T is rather far from T0, the probability P[T, E) calculated by (6.52) using H{E) is not precise. It is therefore very important to verify the form of each P[T, E) by looking at its plot before using it. In general, we have a Gaussian form as shown in

Figure 6.7 Probability P(E) obtained by simulation in the case of a face-centered cubic antiferromagnet with interaction J between nearest neighbors. Left: Energy histogram taken at temperature T0 = 1.76] (ks = 1), just above the transition temperature. Right: Histogram taken at the transition temperature To = Tc = 1.755/ E is the total energy of the system.

Fig. 6.7 (left). If T is far from T0, P(T, E) calculated has an irregular, asymmetric form. We should not use it [110].

• (3) To determine with precision the transition temperature TC(L], we have to choose T0 in the critical region, as close as possible to the presumed TC(L). To have a good choice, we have to make several simulations with Metropolis algorithm to detect a good value of To.
• (4) When the transition is of first order, Р(Г0, E) presents a double peak if T0 coincides or very close to TC(L). This is shown in Fig. 6.7 (right) where the energies at the peaks E and £2 correspond to energies of the ordered and disordered phases, which coexist at T0. The histogram between two peaks is almost zero, indicating an energy discontinuity, namely the latent heat of the system.

The determination of the critical exponents with the histogram method is very precise (see, for example, the original papers of Ferrenberg and Swendsen [110]).

# Multiple-Histogram Technique

The multiple-histogram technique is known to reproduce with very high accuracy the critical exponents of second-order phase transitions [50, 111], It is more complicated to be implemented because we have to realize many histograms in independent simulations, but it gives much better results in difficult cases.

The principle consists of the following steps [111]:

• (1) First, we realize independent simulations at n temperatures Г,
• (/ = 1, • • • , n). For each temperature Tit the number of Monte Carlo steps is Л/,. The histogram taken during the simulation at that temperature is H[E, Г,): One has Tt) = N,.
• (2) Second, we calculate the density of states p(E) by

where the partition function Z(Г,) is

We see here that p(E) and Z(T,) should be calculated self-consistently. The choice of neighboring temperatures T, T2, ■ ■ ■ , T„ should be guided as the choice of T0 discussed in the single histogram technique shown above.

(3) Once p(E) and Z(T) are obtained, we can calculate the thermal average of any physical quantity A at T by

Thermal averages of physical quantities are thus calculated as continuous functions of T. The results are valid over a much wider range of temperature than for any single histogram. The calculation of the critical exponents is more precise than with a single histogram technique.

# Wang–Landau Flat-Histogram Method

Wang and Landau [351] have proposed a Monte Carlo algorithm for classical statistical models which allowed us to study systems with difficultly accessed microscopic states. In particular, it permits to detect with efficiency weak first-order transitions [250, 251]. The algorithm uses a random walk in the energy space in order to obtain an accurate estimate for the density of states g[E) which is defined as the number of spin configurations for any given E. This method is based on the fact that a flat energy histogram H[E) is produced if the probability for the transition to a state of energy E is proportional to giEy1.

We summarize how this algorithm is implemented here. At the beginning of the simulation, the density of states is set equal to one for all possible energies, g[E) = 1. We begin a random walk in energy space (E) by choosing a site randomly and flipping its spin with a probability proportional to the inverse of the temporary density of states. In general, if E and E' are the energies before and after a spin is flipped, the transition probability from E to Ef is

Each time an energy level E is visited, the density of states is modified by a modification factor / > 0 whether the spin is flipped or not, i.e., g(E) -> g(E) f. At the beginning of the random walk, the modification factor / can be as large as e1 ~ 2.7182818. A histogram H(E) records the number of times a state of energy E is visited. Each time the energy histogram satisfies a certain "flatness” criterion, / is reduced according to f -*■ J~J and H[E) is reset to zero for all energies. The reduction process of the modification factor / is repeated several times until a final value /fjnai which is close enough to one. The histogram is considered as flat if

for all energies, where x% is chosen between 70% and 95% and (tf (£)) is the average histogram.

Thermodynamic quantities [49, 351] can be evaluated using g[E). For example,

where Z is the partition function defined by

The canonical distribution at any temperature can be calculated simply by

In practice, we have to choose an energy range of interest [223, 311] [Emin, Emax). We divide this energy range into R subintervals, the minimum energy of each subinterval is E,'nin for / = 1, 2, • • ■ , R, and the maximum of the subinterval; is E'max = E'^ + 2 Д E, where ДЕ can be chosen large enough for a smooth boundaiy between two subintervals. The Wang-Landau algorithm is used to calculate the relative density of states of each subinterval (E,'nin, E'max) with the modification factor /f,nal = exp(10-9) and flatness criterion x% = 95%. We reject the suggested spin flip and do not update g[E) and the energy histogram Я(Е) of the current energy level E if the spin-flip trial would result in an energy outside the energy segment. The density of states of the whole range is obtained by joining the density of states of each subinterval [E,'„in + ДЕ, E'max — AE). Numerous examples of applications of this method can be found in the literature and in part II of this book.