Monte Carlo algorithm for the random cluster model
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b Monte Carlo algorithm for the random clustermodel
Christian R. Scullard Lawrence Livermore National Laboratory, Livermore, California, USA
Abstract.
I present a generalization of the Newman-Ziff algorithm, a standardtechnique in numerical percolation, to the random cluster model. The new methodis straightfoward to implement, works for real q >
0, and is completely free of thecritical slowdown problem. By using the method to sample the critical polynomial,accurate estimates of Potts critical points can be found using relatively small latticesizes. Furthermore, results for all values of q can be found at once within a singlesimulation. I demonstrate its utility here by computing critical points for non-integervalues of q on the square lattice to compare with the exact solution, and on the non-planar square matching lattice. onte Carlo algorithm for the random cluster model q . While methods exist that handle real q [7, 8, 9, 10], thereare often trade-offs involved either in efficiency or ease of use compared to SW. Here,I present a Monte Carlo scheme, based not on any of these previous techniques but onthe Newman-Ziff approach to percolation, that has no critical slowdown and works forreal q >
0. In addition, results for an arbitrary number of different q can be obtainedfrom a single simulation. The price is a comparatively steep memory requirement forlarge lattices. There very likely exist strategies for working around this issue, whichI will discuss later, but if we focus for now on the problem of finding critical pointswe need not worry about it. Because the critical polynomial [11, 12, 13] provides veryaccurate estimates of transition points, even on relatively small graphs, it can easily besampled with the present Monte Carlo scheme. In the standard polynomial method,the polynomial, or more often just its root [14], is computed via an analytic technique[15] and, for two-dimensional planar lattices, provides accuracy far surpassing what iscurrently possible with Monte Carlo methods [16]. However, relying as it does on atransfer matrix calculation it is somewhat challenging to generalize to non-planar, notto mention higher-dimensional, lattices. Additionally, even in the planar case, it israther more complex to implement than a typical Monte Carlo algorithm and one mayprefer to trade some accuracy for simplicity.The partition function for the random cluster model with edge probability p is givenby [17] Z = X { ω } p n (1 − p ) N − n q C ( ω ) (1)where n is the number of edges present in the configuration ω , C ( ω ) is the numberof connected components (including isolated vertices), and the sum is over allconfigurations. It is convenient to write Z in the form Z = N X n =0 L X C =1 p n (1 − p ) N − n q C Γ n ( C ) , (2)for a lattice with L vertices and N edges, and where Γ n ( C ) is the number ofconfigurations consisting of n edges and C connected components. This function isobviously not known analytically in general, but it is a simple matter to estimate itduring the simulation. The probability of some event A can then be written P ( A ) = 1 Z N X n =0 L X C =1 p n (1 − p ) N − n q C Γ n ( C ) P ( A | n, C ) (3)where P ( A | n, C ) is probability of A for fixed n and C . Therefore to compute theprobability of an event A , the quantities P ( A | n, C ) and Γ n ( C ) can be estimated for all n and C and the results combined in Equation (3). All configurations with a given n and onte Carlo algorithm for the random cluster model C have the same probability, and thus they might be said to constitute a microcanonicalensemble. More concretely, we have P ( A | n, C ) = Ω n ( A, C )Γ n ( C ) (4)where Ω n ( A, C ) is the number of configurations consisting of n edges and containing C clusters for which the event A occurs. Random sampling is therefore used onlyto count configurations. The numerical approach encapsulated in Equation (3) isthe generalization to arbitrary q of that of Newman and Ziff [18, 19] (NZ), to whichthis method reduces when q = 1. Their sampling procedure will be used here as itis straightforward to implement and hardly needs any generalization for the presentpurposes. Although I describe it below, I recommend the reader consult their paper formore details. I point out here that, also in the spirit of generalizing NZ, Wang, Kozanand Swendsen [20, 21, 22] developed an algorithm they call “binary tree summation”.However, it is significantly different from the present method, which I believe is a moredirect descendant of NZ.The sampling algorithm for a single run works as follows. First, the lattice is empty,with all edges closed. We then add edges one by one in a completely random order untilthey are all open. At step n , we find the number of clusters, C , which is simple to do;adding an edge can only reduce C by one or leave it unchanged, and when n = 0 we have C = L . Performing multiple runs we then find f n ( A, C ), the fraction of runs for whichthe event A occurs at a given n and C . The estimate of the number of configurations isthen Ω n ( A, C ) = (cid:18) Nn (cid:19) f n ( A, C ) . (5)The only difference between runs is the order in which the edges are added to the system.Permutations of the edge order are handled in exactly the same way as in NZ, and in factthey provide a function written in C called permutation() which will work unchangedhere. Their algorithm is as follows. Start with a list of the edges in the system organizedin an arbitrary order. Set i = 1 and then,(i) Choose a number j uniformly and at random from the range i ≤ j ≤ N .(ii) Exchange the edges that are currently in positions i and j (and do nothing if i = j ).(iii) Set i = i + 1 and return to step 1 until i = N .Adding all N edges while tracking clusters takes a time O( N ) [19]. Although we do nothave independence between step n and step n + 1 within a single run, all that matters isthat there be independence between runs and that all configurations with a fixed n and C should appear with equal probability, which is plainly the case here. We can now seehow the critical slowdown and the restriction to integer q are evaded by this algorithm.The configuration sampling knows nothing about the temperature (the probability p inour case) or q , but rather these appear as parameters in a convolution, Eq. (3), thatoccurs after the sampling has been done. For the same reason, we get results for everyvalue of q in a single simulation. onte Carlo algorithm for the random cluster model ptr[i] , which for vertex i either gives its parent sitein the cluster or, if i is a root, its value is the negative of the size of the cluster. Thus,every site is connected by a path of parents to the root. So if an edge joins the two sites i and j , we take the following steps(i) Find the roots of i and j by traversing their trees.(ii) If i and j have the same root, they are in the same cluster. This situation will beimportant below, but as far as the cluster configuration goes we need do no more.(iii) If i and j are different, determine which cluster is smaller. If, say, j is in the smallercluster, we make the root of j point to the root of i and add the size of the smallercluster to the larger one by adjusting ptr[i] . Obviously, if both clusters are thesame size we can choose arbitrarily which cluster to add to which.The algorithm is most efficient if all sites point to their roots. Thus, after we have foundthe roots with step 1, it is worthwhile to backtrack over the path and set all pointersto point directly to the root we have just found. See NZ for a fuller discussion of thispath compression.To locate the critical point, we estimate the critical polynomial. This is defined (upto a sign) on a doubly periodic finite lattice, or basis, B by [13] P B ( p, q ) ≡ P (2D; p, q ) − qP (0D; p, q ) (6)where P (2D) is the probability that there is a cluster that wraps both dimensions and P (0D) is the probability that no cluster wraps either. A third possibility, which wewill encounter later, is referred to as 1D and is the event that a wrapping spans onlyone dimension. See Figure 1 for examples of these events. That the root in [0 , P B ( p ) provides an estimate for the transition point is well established [25, 12, 26]and is a consequence of universality. Estimates are obtained by computing averages of P (2D | n, C ), P (0D | n, C ) and Γ n ( C ) and combining them with Eq. (3). As a function of p , P B ( p, q ) is monotonic in p , taking the value − q at p = 0 and 1 at p = 1. A simpleroot-finding algorithm suffices to find the estimate for p c , and here I use the bisectionmethod. By making good initial guesses, perhaps derived from smaller-sized lattices,one can minimize the number of evaluations needed by this root-finding step.The next piece of the algorithm is the identification of the different wrapping events.This is done by using a variation of an approach first suggested by Machta et al. [5].For each site i , in addition to the pointers ptr[i] , we also associate x[i] and y[i] , the x and y lattice displacements to the root site of the vertex i . Upon joining two verticesthat are already in the same cluster, if at least one vertex’ displacement vector pointsthrough a boundary then this operation might produce a wrapping cluster. That is, theroot of the cluster can now be connected to itself via a path that winds through one or onte Carlo algorithm for the random cluster model Figure 1.
The three possible wrapping events on a periodic lattice, or basis. both of the periodic directions. The exact manner of this winding can be encoded in aconnection vector, as shown in Figure 2. It is a simple matter to compute this vectorfor each join we perform. If we are joining vertices, say 1 and 2, for which the vectorconnecting 1 to 2 is u , and their root displacement vectors are d and d then theconnection vector is given by c = d − d + u mod L (7)(note that the overall sign is irrelevant). If the connection vector is zero, then we havecreated no new wrapping. Upon forming a 1D wrapping, we must store the connectionvector associated with that cluster ‡ . Call this vector c , and if upon later joining twofurther vertices of this cluster we form a new connection vector, c , such that the crossproduct c × c = (8)then we know we have just formed a 2D wrapping.Note that if we were to restrict ourselves only to planar lattices, the state of thewrapping would be a global property. For example, if we formed a wrapping with c = (1 ,
0) after already having encountered a cluster with c = (0 , B ( N, n, C ; p, q ) ≡ (cid:18) Nn (cid:19) p n (1 − p ) N − n q C . (9)In the summation, the maximum C is L and N is of the order L so naively this sumcontains O( L ) terms. In reality, Γ n ( C ) is non-zero only within a range of C boundedby the limits C min ( n ) and C max ( n ). These are lattice-dependent and unknown exactly ‡ as a practical matter it can be stored in the same array as the displacement vectors, as we know thedisplacement vector for the root is zero and its entry can be used for something else onte Carlo algorithm for the random cluster model Figure 2.
The 1D wrapping event of Figure 1 tiled to demonstrate the meaning ofthe connection vector, indicated by the red arrow. The circle denotes the root of thecluster. The connection vector is (2 , −
1) in this case. but it is a simple matter to keep track of them during the calculation. One then findsthat there are really more like L terms. We can obviously not evaluate (9) directly, asindividual terms may be very large or very small even though their product might beof order unity or a few tens. We follow NZ here by setting the largest term, located at n max , to 1 and calculating all other terms relative to that with recursion. When q = 1,the largest term occurs at n max = N p . However, for general q , we must also account forthe n -dependence in q C max ( n ) . Approximating C max ( n ) as a linear function of n , we findthe estimate for n max , n max = N p (cid:20) − (1 − p ) L N ln q (cid:21) (10)which is not perfect § but performs well enough for the present purposes. We then setthe term with n = n max and C = C max ( n ) equal to 1 and use recursion to fill in theother terms in C down to C min ( n max ), B ( N, n, C ; p, q ) = B ( N, n, C + 1; p, q ) /q (11)for n = n max . For n > n max , we compute B ( N, n, C max ( n ); p, q )= B ( N, n − , C max ( n ); p, q ) N − n + 1 n p − p (12)and then fill in the rest of this n down to C = C min ( n ) with (11), repeating the processuntil we reach n = N . For n < n max , we start the recursion in C at C min ( n ) B ( N, n, C min ( n ); p, q )= B ( N, n + 1 , C min ( n ); p, q ) n + 1 N − n − pp (13)and carry this up to C max ( n ) with B ( N, n, C ; p, q ) = B ( N, n, C − p, q ) q. (14) § Linearity is a poor model for C max ; the estimates from this formula are sometimes quite far from thetrue n max but are close enough to allow the recursion to work. onte Carlo algorithm for the random cluster model Z in Equation (2),Equation (3) will be normalized correctly.In Table 1, I present results on the square lattice for various q , where the criticalpoints are known to be [27]: p c ( q ) = √ q √ q . (15)If the critical point of a lattice is known exactly, the critical polynomial will find theexact solution for any size of basis. To illustrate this, calculations are presented for both L = 3 using 10 runs and L = 16 using 10 . Although both are in agreement with theexact solution it is clear that for L = 16 and q >
1, accuracy suffers somewhat because ofestimates that rely on some portion of ( n, C ) space that is relatively difficult to resolve.We would expect this situation to improve eventually as q becomes large because theprobabilities become more weighted toward the least-occupied configurations.For an unsolved problem, the critical polynomial does not give the exact answer forany size of basis. However, the estimates it provides are generally very accurate evenwhen the basis is of a size normally considered small for Monte Carlo. For example, forthe kagome lattice, a basis of size L = 4 already produces estimates of the percolationthreshold accurate to six digits [12]. By the time L ∼
20 a polynomial prediction for agiven system is generally accurate enough that it cannot be checked with Monte Carlo.Thus, using these relatively small lattices is perfectly acceptable in the present approach.As an example of an unsolved problem, I calculate critical points for the matchinggraph of the square lattice, the L = 3 basis of which is shown in Figure 3. This is anon-planar lattice, with the diagonals in each face crossing but not directly connectedto each other. Its site percolation threshold is 1 − p c (square site) ≈ . q including that for q = 1, the bond percolation threshold, which has previously beenfound numerically by several authors [28, 29] and very recently by Xu et al. [30] usinga Monte Carlo-critical polynomial technique. Their value is 0 .
250 368 40(4). Becausenone of these cases is exact, we expect the estimates to depend on L , and I used afew different sizes around L = 25 just to confirm that the same predictions were made.These results were found with 10 total runs, far short of what one might normally usebut my goal is only to demonstrate that the algorithm works, not to present values ofdefinitive accuracy. Note that one may instead want to do computations on smallerlattices and then extrapolate to infinite L , as this may be the best way to maximizethe accuracy of the estimates. There is, however, not yet a good theory governing thisextrapolation, and scaling exponents, which are known to be lattice-dependent [14, 16],must be calculated empirically, so I do not explore this here.There are two ways in which the efficiency of the original NZ algorithm suffersunder the generalization to arbitrary q . The first is the fact that we must now continueevery run until we have added the final edge. In the percolation algorithm, once the2D wrapping, say, has occured we can stop adding edges because we know that all onte Carlo algorithm for the random cluster model Figure 3.
The L = 3 basis for matching graph of the square lattice. q L = 3 L = 16 exact1.0 0.5000002(7) 0.500000(2) 1/21.5 0.5505103(6) 0.550511(3) 0.5505102572...2.5 0.6125740(5) 0.61252(5) 0.6125741133...3.5 0.6516683(4) 0.6515(3) 0.6516685226...9.5 0.7550342(2) 0.7552(2) 0.7550344704...10.0 0.7597467(2) 0.7599(2) 0.7597469266... Table 1.
Potts critical points, p c , on the square lattice computed with the Monte Carloalgorithm for bases of size 3 and 16. These both agree with the known exact solution, asis standard for critical polynomials. These estimates are found by calculating averagesof P (2D | n, C ), P (0D | n, C ), and Γ n ( C ), combining these with Eq. (3) to get P (2D)and P (0D) and numerically finding the root of equation (6). Results for all q are foundin a single simulation. q p c q p c Table 2.
Critical points, p c , for the random cluster on the square matching lattice. onte Carlo algorithm for the random cluster model P (2D | n ) will be 1. But here we are also resolving C , which means every runmust go to completion. This will result in a longer run time by a factor of roughly 2[19]. The other issue is that we must now converge estimates not just for given n butfor given n and C . The loss here is somewhat more difficult to quantify, but it manifestsas the aforementioned dropoff in accuracy after q = 1, as can be seen in Tables 1 and 2.Other questions we can consider are whether this sampling algorithm might beused to find quantities other than critical points and whether it is feasible to work inhigher dimensions. The key barrier in both cases is the memory limitation; we muststore an array in both n and C for all quantities of interest. As argued earlier, this isa memory requirement that is roughly L k . For this reason, as well as the fact thatit also becomes increasingly difficult to evaluate the convolution (3), we may struggleto compute quantities that require large lattices. However, these problems are likelysurmountable in many cases. For example, if one is interested in some function F ( p )only at a particular value of, or within a small range of, p , one does not need to keepevery n and C as most of the terms in the convolution will prove to be irrelevant.Keeping fewer terms also makes the evaluation of the convolution easier. It is clearthat thorough study will be needed to determine the true limitations and the range ofapplicability of this algorithm.I thank Robert Ziff for a stimulating and helpful discussion of this problem. Thiswork was performed under the auspices of the U.S. Department of Energy at theLawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344and was supported by the LLNL-LDRD Program under Project No. 19-DR-013. References [1] Potts R B 1952
Proc. Camb. Phil. Soc. Rev. Mod. Phys. Phys. Rev. Lett. Phys. Rev. Lett. Phys. Rev. E Phys. Rev. Lett. Phys. Rev. B Phys. Rev. E Physica A
Phys. Rev. E Phys. Rev. E J. Phys. A: Math. Theor. J. Phys. A: Math. Theor. J. Phys. A: Math. Theor. J. Phys. A: Math. Theor. Phys. Rev. Research Physica Phys. Rev. Lett. k It is also not really trivial to realize this L behavior because it requires dynamically managing thearray as more C appear onte Carlo algorithm for the random cluster model [19] Newman M E J and Ziff R M 2001 Phys. Rev. E arXiv:cond-mat/0203264 [21] Wang J S 2003 Physica A [22] Wang J S, Kozan O and Swendsen R H 2003
Computer Simulation Studies in Condensed-MatterPhysics XV: Proceedings of the Fifteenth Workshop Athens, GA, USA, March 11-15, 2002 ( Springer Proceedings in Physics vol 90) (Springer, Berlin, Heidelberg) p 189[23] Galler B A and Fisher M J 1964
Commun. ACM Algorithms
Phys. Rev. Lett.
Phys. Rev. E J. Phys. C: Solid State Phys. L445[28] Feng X, Deng Y and Bl¨ote H W J 2008
Phys. Rev. E Phys. Rev. E062101[30] Xu W, Wang J, Hu H and Deng Y 2020