On generalized cluster algorithms for frustrated spin models

P.D. Coddington† and L. Han

Physics Department, Syracuse University, Syracuse, NY 13244, USA

May 21, 1993

Standard cluster update algorithms have proven to be very effective for many different spin models, however they fail for frustrated spin systems. Recently a generalized cluster algorithm was introduced which works extremely well for the fully frustrated Ising model on a square lattice. However we have applied this method to other frustrated spin models, including the fully frustrated Ising model on triangular and cubic lattices, and found that it works poorly, producing clusters which are too large or too small.

† Current address: Northeast Parallel Architectures Center, 111 College Place, Syracuse University, Syracuse, NY 13244, USA.

I. INTRODUCTION

The cluster update algorithm of Swendsen and Wang (SW) [1] provides a great improvement in the efficiency of Monte Carlo simulations over the Metropolis algorithm [2, 3] and other local update schemes by dramatically reducing critical slowing down [3, 4] in certain ferromagnetic and anti-ferromagnetic spin models. However the SW algorithm fails for frustrated spin models, which are the systems for which critical slowing down can be the most extreme. The reason for this failure is that at the critical point the algorithm produces a cluster which encompasses a large fraction (if not all) of the sites in the lattice, so that flipping the spins in this cluster is virtually the same as a trivial global spin flip.

Recently Kandel, Ben-Av and Domany (KBD) introduced a generalized cluster algorithm [5, 6] which works very well for the two dimensional fully frustrated Ising model [7]. The only other cluster algorithm which has been shown to work effectively for a frustrated spin model is the Replica Monte Carlo algorithm [8], however this has only been successfully applied to the 2-d Ising spin glass [9].

In section II we introduce the KBD generalized cluster algorithm, and outline the method we used for implementing the algorithm for a given spin model. We then apply this methodology to the fully frustrated Ising model (FFIM) on a square lattice in section ??, the triangular lattice FFIM in (section IV), some other 2-d Ising models in section V, and finally to the 3-d FFIM in section ??.

II. FORMULATING A GENERALIZED CLUSTER ALGORITHM

First we will briefly describe the SW cluster algorithm for the Ising model [3], which has an interaction energy given by

E = − ∑_{<i,j>} J_{ij}σ_{i}σ_{j}, (1)

where the spins σ_{i} can take the values +1 or −1. We will consider the case where the interaction strength J_{ij} takes on values +J or −J, for some constant parameter J. Let us define a link to be the connection between two neighboring sites on the lattice. A link is said to be satisfied if the two sites are in a state of minimum energy, which means that the spins are the same for J_{ij} > 0 or opposite for J_{ij} < 0. Otherwise it is unsatisfied. In the SW algorithm, we introduce bonds between spins on neighboring sites with probability 1 − e^{−2K} if the link is satisfied, and 0 if it is unsatisfied. Here K = J/kT is a dimensionless coupling constant, where T is the temperature and k is Boltzmann’s constant.

This procedure creates clusters of bonded sites. The update part of the algorithm consists of flipping all the spins in a cluster with probability \frac{1}{2}. Since the clusters can be quite large, this clearly produces large, non-local changes in the spin configuration, and thus decorrelates the configurations much faster than local update algorithms.

This algorithm works very well for the ferromagnetic (J_{ij} = +J) and antiferromagnetic (J_{ij} = −J) Ising model, but fails for frustrated systems, which have a mixture of positive and negative couplings in such a way that it is not possible for all links to be satisfied in the ground state. The simple reason for this is that the critical point for such systems usually occurs at low temperature (i.e. large K), where the probability of putting down a bond is large. For example, the two dimensional fully frustrated Ising model has a critical point at zero temperature, where there is one unsatisfied link per plaquette. Hence the SW algorithm bonds each site to 3 of its neighbors, which is more than enough bonds to connect all the sites in the lattice into a single cluster. A more fundamental reason is that the SW algorithm only uses information from links, from which it is not possible to see the frustration in the system. To do that one needs information from the spins over at least a plaquette, which is the basic idea behind the KBD generalized cluster algorithm.

The generalized cluster algorithm of KBD begins by expressing the energy of the system as a sum E = ∑_{l} V_{l}. We will consider l to label subregions of the lattice, such as links, plaquettes, or possibly some other group of sites. For each l we stochastically assign one of n possible operations with a probability P_{i}^{l}(u), which will depend on the spin configuration u and the operation i. These operations will involve placing bonds between certain sites in the region l, which will create clusters of connected sites. The probabilities must of course satisfy

∑_{i} P_{i}^{l}(u) = 1,  0 ≤ P_{i}^{l}(u) ≤ 1,  (2)

for all l and u.

In order to satisfy detailed balance, KBD show that the probabilities must also satisfy

E_{i}(u) − \frac{1}{β} log P_{i}(u) + C_{i} = 0, (3)

E_{i}(u) is the energy for the configuration u and C_{i} is a constant (independent of the configuration), and we have supressed the subscript l labeling the specific collection of sites. The above equation can also be written as

P_{i}(u) = e^{βE_{i}} or 0. (4)

Imposing the normalization condition of equation 2, we have

????? (5)

which can be written as a matrix equation

B\underline{c} = 1, (6)

where the vector c is made up of the C_{i}, and the elements B_{ij} of the matrix B are either e^{βE_{i}} or 0, where j labels ??? (the u’s??).

EXPLAIN HOW TO USE IT, BY SETTING UP A MATRIX EQUATION. OUT-LINE THE CONSTRAINTS ON THIS MATRIX TO GET IT TO WORK, E.G.

MUST BE (OR JUST BETTER TO BE?) UPPER TRIANGULAR, MUST HAVE MORE PARAMETERS THAN CONSTRAINTS, ETC.

III. THE SQUARE LATTICE FULLY FRUSTRATED ISING MODEL

The fully frustrated Ising model (FFIM) was introduced by Villain [7] as a simple regular frustrated system which lacks the extra complication of disorder that is present in other systems such as spin glasses. The critical point of this model is at zero temperature, where the configurations have one unsatisfied link per plaquette. The SW algorithm for this model puts a bond on all satisfied links, so that sites are connected to at least 2, and on average to 3, of their 4 neighbors. This produces a single cluster (consisting of all the sites in the lattice) at the critical point, so the SW algorithm does not work for this model. However the generalized cluster algorithm of KBD, which deletes bonds using information on the state of plaquettes rather than links, works extremely well in this case [5, 10].

A. The KBD Algorithm

In the KBD algorithm there are seven basic freeze/delete operations, which are shown in Fig. 1, where the bold lines indicate frozen bonds. The probabilities for performing these operations for each of the two energy states (+J and −J) are given in Table I. (SHOW HOW WE GET THIS USING THE MATRIX FORMALISM, AS AN EXAMPLE??)

Note that different orientations of the same basic tolopogy of links, such as operations 2 and 3, or operations 4-7, in Fig. 1, count as different operations. In this case, this means we have some extra free parameters, in that there is some flexibility in choosing the probabilities of the different orientations. However to make the matrix algebra easier, we could combine these different orientations into a single operation, and then choose each orientation with equal probability. Both these approaches are correct (i.e. satisfy detailed balance): one is simpler, the other has extra freedom which may allow some improvement in the performance of the algorithm, as we shall see later.

FIX THIS UP!!

Updating these clusters does not produce the large changes in magnetization which occur for the square lattice model, where the average maximum cluster size is 0.6432(5) for this lattice size (64 × 64).

The main reason for this is that the square lattice has two energy conserving operations (numbers 2 and 3 in Fig. 1), which means that there is complete freedom in choosing

and KBD were able to prove that for a particular choice of this ratio the lattice was split into at least two clusters.

B. Dynamic Critical Exponents

KBD measured a dynamic critical exponent of z ≈ 0.55 for their algorithm [10], compared to z ≈ 2 for the Metropolis algorithm (see Refs. [3, 4] for discussion and definitions of autocorrelation times and dynamic critical exponents). However they had very low statistics (100,000 sweeps) and fairly small lattice sizes (up to 128^{2}), and it has been seen for the Ising ferromagnet that it is very difficult to get an accurate determination of the z from such results [11]. We have therefore obtained much better data for the autocorrelations for this algorithm, in order to get a better determination of the dynamic critical exponent.

We have measured the autocorrelation function

ρ_{A}(t) = \frac{< A(0)A(t) > − < A(0) >^{2}}{< A(0)A(0) > − < A(0) >^{2}} (7)

for the magnetization (A = M) as well as the square and the absolute value of the magnetization (A = M^{2} and A = |M|), and used this to extract exponential autocorrelation times τ_{exp}, via

ρ(t) ∼ e^{−t/τ_{exp}}, (8)

as well as integrated autocorrelation times

τ_{int} = \frac{1}{2} + ∑_{t=1}^{∞} ρ_{A}(t), (9)

which are relevant to the increase in the statistical error due to correlated configurations. Details of the methods used for measuring the autocorrelations, doing the fits, and estimating the errors are given in Ref. [11].

In Table IV we give the results for the autocorrelation times of M^{2} for different lattice sizes L up to L = 256. All the results are from at least 5 × 10^{6} sweeps. Fig. 2 shows a log-log plot of these results. Straight line χ^{2} fits give z_{int} = 0.26(3) and z_{exp} = 0.6(1). The integrated autocorrelations are also well fit by a logarithm (the dashed line in Fig. 2), so that z_{int} could also be zero. The same situation also occurs for the autocorrelations of cluster algorithms for the ferromagnetic Ising model [11], since it is very difficult to differentiate between a logarithm and a small power.

Let us also define τ_{init} to be the autocorrelation time obtained from an exponential fit to ρ(t) at only two points, t = 0 and t = 1. This is plotted in Fig. 3 on a log-linear plot. We can see that τ_{init} grows extremely slowly with L (the dashed line shows the function log(log L)). Thus is in marked contrast to τ_{exp}, which grows as a substantial power of L. Unfortunately this means that it is very difficult to get a good asymptotic fit to ρ(t) at large t, since it has already fallen off so much at t = 1 that the signal is very small and noisy in the region where we need to fit to obtain τ_{exp}. The fact that τ_{init} increases so slowly means that the τ_{int} grows much slower than τ_{exp}, and consequently z_{int} appears to be smaller than z_{exp}.

We also measured the autocorrelation times for the magnetization and the absolute value of the magnetization. For the former the results were zero for all lattice sizes, while the latter gave results almost identical to those given above for the square of the magnetization.

In almost every iteration of the standard KBD algorithm, the lattice is split into 2 clusters, and only rarely into 3 or more. The update consists of flipping the spins in a cluster with probability 1/2. If there are two clusters, this means that in half of the iterations either no spins are flipped, or all the spins are flipped, which means there is no change in the configuration apart from a trivial global spin flip. Only half of the iterations change one of the clusters and not the other, and therefore make a non-trivial update of the configuration. If there are 2 clusters, a more effective update scheme is to pick one of them at random and flip its spins. This guarantees that there are no “wasted” iterations.

One would expect that this new update scheme would decrease the error in measurable quantities by a factor of \sqrt{2}, and indeed we found this to be the case. One would also expect this to be due to a decrease in τ_{int} by a factor of 2, however this is not the case. In fact the autocorrelation function for this algorithm is oscillatory, as shown in Fig. 5. Thus using the definition for τ_{int} (equation 9) gives a value close to zero, and much less than half of the value for the standard KBD algorithm.

WHY DOES THIS GIVE NEGATIVE AUTOCORRELATIONS? (LOOKS A BIT LIKE WHAT HAPPENS FOR THE SW REPLICA METHOD). DO THE FOR-MULAS BREAK DOWN IF THE AUTOCORRELATION FUNCTION IS NEGA-TIVE? CHECK SOKAL’S PAPERS.

C. A Variant of the KBD Algorithm

ANGLED OPERATIONS.

We measured the autocorrelations for this alternate version of the KBD algorithm, using only the angle operations (4-7) rather than the operations (2,3) used by KBD. We used the subgroup (4,5) with probability p_{sub} and the subgroup (6,7) with probability 1 − p_{sub}. We tried three different values of p_{sub}: 1.0, 0.8, and 0.5. Integrated autocorrelation times for these three variations of the KBD algorithm are shown in Fig. 4, along with results for the Metropolis algorithm, and the KBD results from Fig 2.

The results are quite surprising. As expected, p_{sub} = 0.5 gives results which are very similar to the Metropolis algorithm (z ≈ 2, since the biggest cluster size is very large and the cluster update does very little. For larger values of p_{sub} we have seen that the biggest cluster size is reduced very slightly, however this seems to be enough to dramatically change the autocorrelations — p_{sub} = 0.8 has z ≈ 0.??, while p_{sub} = 1.0 has z ≈ 0.??, which is very close to that of the standard KBD algorithm! It is amazing that in going from p_{sub} = 0.5 to 1.0 we only reduce the biggest cluster size by about ??%, yet we reduce z from nearly 2 to nearly 0!

This result is very promising, since it implies that critical slowing down can be reduced dramatically even when the biggest cluster size is quite large (??% of the lattice in this case), so for other frustrated spin models we may not have to unfreeze the lattice very much to get good results.

The KBD algorithm works extremely well for the square lattice FFI, almost completely eliminating critical slowing down. However it is not clear that it can be successfully applied to other frustrated models. It is only by a very fortuitous geometrical happenstance that the standard KBD algorithm is able to split the lattice into two large clusters at the critical point, and the algorithm we have investigated which uses only the angled operations only works well because it has a tunable parameter p_{sub} which allows the maximum cluster size to be reduced just enough to give a greatly reduced z. We were therefore interested to see whether the generalized cluster algorithms will work for other frustrated spin models.

IV. THE TRIANGULAR LATTICE ISING ANTIFERROMAGNET

We first attempted to apply the KBD cluster algorithm to an even simpler fully frustrated model, the anti-ferromagnetic Ising model on a triangular lattice. This model is very similar to the square lattice FFIM, being in the same universality class and also having a critical point at zero temperature [12, 13].

The ground state of this model has one unsatisfied link per triangular plaquette. The SW algorithm puts a bond on all satisfied links, so that sites are connected on average to 4 of their 6 neighbors, and it is easy to see that each site must have connections to at least 3 of its neighbors. So as for the square lattice, the SW algorithm freezes the lattice into a single cluster at zero temperature.

We can get some idea of the percentage of bonds which need to be deleted to break up this single cluster by noting that the bond percolation threshold for a triangular lattice is 0.35792 [14]. For the square lattice the bond percolation threshold is 0.5 [14], which is also the ratio of frozen bonds to links in the KBD algorithm for the square lattice FFIM at the zero temperature critical point. This is of course only a rough pointer to what is required, since for bond percolation the bonds are placed on the lattice at random, whereas for the KBD algorithm on the square lattice we must have two bonds for each plaquette at zero temperature. Even with the bond/link ratio being the same as the percolation threshold, the largest cluster for the KBD algorithm percolates at a temperature well above the critical point.

The application of the KBD algorithm to the triangular lattice FFIM is very simple. As with the square lattice, we choose the basic element to be a plaquette (a triangle in this case), which has two possible energy states, +J and −J. For the triangular lattice there are only four possible freeze/delete operations, which are shown in Fig. 6. We can either delete all the bonds, or freeze one of the two satisfied bonds.

The freeze/delete operations 2-4 conserve the energy of the ground state. Solving the constraints as done by KBD [6], or by solving the matrix equation constructed using the method of section II, (DO THIS EXPLICITLY AS AN EXAMPLE??) gives the probabilities given in Table V.

At the zero temperature critical point, we only perform the energy conserving freeze/delete operations. This gives a bond/link ratio of \frac{1}{3}, which is below the percolation threshold, although in this case every plaquette has one bond, so the bonds are not placed randomly. However we expect to at least be able to create multiple clusters, as with the KBD algorithm with angled operations in the square lattice case.

We tested this algorithm at zero temperature, following every cluster update by a Metropolis sweep to ensure ergodicity. We found that the algorithm produces one very large cluster, and a number of very small clusters. For a 64^{2} lattice there are an average of 205.1(1) clusters, with the average size of the largest cluster being 0.9233(1) of the lattice volume.

This is very similar to the results for the square lattice case using the angled operations with equal probability. In that case there is a free parameter (p_{sub}) which allows us to bias the choice of the operations so as to reduce the largest cluster size and greatly improve the performance of the algorithm. However in this case there is no such freedom, and in order to satisfy detailed balance the two possible energy conserving operations for each plaquette must be chosen with equal probability. We checked that biasing the choice of operations does indeed reduce the largest cluster size, but of course it also gives incorrect results due to the violation of detailed balance.

NEED TO MEASURE AUTOCORRS FOR METROP AND CLUSTER TO SHOW THAT CLUSTER IS NO BETTER. TO DO THIS WE NEED TO LOOK AT A CLEVER DEFINITION OF THE MAGNETIZATION, WHICH TAKES THE ABSOLUTE VALUE OF SUMS OVER AXES, RATHER THAN SUMS OVER THE WHOLE LATTICE.

One idea we investigated to try to improve the algorithm is to make the basic object something larger than a plaquette. We tested the algorithm using the double plaquette shown in Fig. 7a. This offers the possibility of reducing the bond/link ratio by sometimes freezing only the link common to both the triangular plaquettes. The problem with this method is that the double plaquette does not cover all the links in the lattice. As shown in Fig. 8, there is an extra link which must be dealt with. If we treat this separately, by for instance using the SW algorithm and freezing it only if it is satisfied, we end up not changing the bond/link ratio, and the cluster sizes are virtually unchanged (IS THIS TRUE??).

Another alternative is to incorporate the extra link into the basic object, giving the shape shown in Fig. 7b. This will now cover all the links when replicated. However both (?? OR JUST THIS ONE??) these algorithms have another problem.

Fundamental to the viability of the generalized cluster algorithm is that there are freeze/delete operations which conserve the energy of the ground state configuration. However because of the shape of the object in Fig. 7b, which is not a closed curve, a ground state configuration can give different energies, and different numbers of frustrated and satisfied links, for this object, depending on which links it covers. In this particular case it is not possible to find a freeze/delete operation which will conserve the energy for all possible ground state configurations. (BUT IT’S NOT JUST THIS - IT’S THAT IT SCREWS UP THE MATRIX, SINCE WE HAVE THE LOWEST ENERGY OPERATION ALSO OCCURING FOR THE SECOND LOWEST ENERGY. RIGHT??)

Thus although it is possible to construct a simple generalized cluster algorithm for the triangular lattice FFIM which does much better than the standard Swendsen-Wang cluster algorithm, it still fares no better than the Metropolis algorithm at the critical point, since the biggest cluster size is too large, and there are no tunable parameters in the algorithm which might enable us to reduce it. Kandel et al. have also tried to apply their method to this model, with the same negative result [15].

V. OTHER 2-D FRUSTRATED MODELS

We also attempted to generalize the KBD algorithm to some other two dimensional frustrated models, the Ising spin glass and the FFIM on a honeycomb lattice.

Let us consider the ±J spin glass, where the couplings J_{ij} in equation 1 are chosen to be ferromagnetic (+J, J > 0) or antiferromagnetic (−J) at random [16]. For the spin glass we can use a very simple generalization of the KBD algorithm. For plaquettes with an odd number of ferromagnetic bonds, we use the KBD algorithm, while for plaquettes with an even number of such bonds, we use the SW algorithm on each of the bonds (i.e. just freeze the satisfied bonds at zero temperature). However in this case there is not the regular geometry which allows the KBD algorithms to split the lattice into two clusters, and at zero temperature we again have the problem of the lattice freezing into a single cluster.

The problems in applying the KBD algorithm to the triangle lattice FFIM were two-fold. Firstly, the simpler model meant that there were less possible freeze/delete operations, and consequently no free parameters to tune; and secondly, the coordination number (and thus the connectivity) of the triangular lattice is very high (6 as opposed to 4 for the square lattice), which means it is easier to produce a very large cluster.

We therefore might expect to do better on a honeycomb lattice, where the coordination number is only 3, and where the fundamental plaquette has 6 bonds instead of 4, giving a much larger number of possible freeze/delete operations. However there is a major problem in constructing a generalized cluster algorithm on a honeycomb lattice – it is not possible to do a black/white partitioning of the lattice so that each bond belongs to only one plaquette, as it is for the square and triangular lattices. It is therefore not obvious how to apply the generalized cluster algorithm formalism in this case.

VI. THE 3-D FULLY FRUSTRATED ISING MODEL

We next attempted to apply the generalized cluster Monte Carlo method to the three dimensional FFIM [17–19]. An efficient cluster algorithm for this model would be very useful, since unlike the 2-d FFIM this model is not exactly solvable, and standard Monte standard Monte Carlo methods are very inefficient. In fact the local Metropolis algorithm does not work at all at zero temperature, since flipping any single spin of a ground state configuration will increase the energy. The model is believed to have a phase transition at low temperature, however due to critical slowing down, simulations are restricted to small lattice sizes at the transition, and even the order of this transition is in doubt [18, 19]. This model also exhibits interesting low temperature behavior, which is even more difficult to investigate using the standard Monte Carlo algorithms.

We first checked the performance of the standard SW cluster algorithm for this model. The bond percolation threshold for a cubic lattice is 0.2492 [14]. We measured the bond/link ratio of the SW algorithm for the 3-d Ising ferromagnet at its critical point to be 0.24078(4), with the maximum cluster containing on average 0.2599(1) of the total number of sites of a 16^{3} lattice (in three dimensions, the SW bonds do not percolate at the critical point, as they do in two dimensions). However at the critical point (K_{c} = 0.738 [19]) of the 3-d FFIM we found that the bond/link ratio was 0.5325(3), with the average maximum cluster size being 0.9942(1) on a 16^{3} lattice. So while the SW algorithm produces good sized clusters and works very well for the ferromagnet, it works extremely poorly for the fully frustrated model – much worse than the Metropolis algorithm, in fact. The problem again is that we need to delete more bonds to break up the largest cluster.

We have generalized the KBD algorithm to three dimensions, with the unit cube replacing the plaquette as the basic element. This means that each freeze/delete operation operates on 12 links, rather than the 4 links of a square plaquette. The number of possible operations is therefore very large, and we have a lot of freedom in choosing them. For the plaquette there are only 2 possible energy states, whereas there are 7 for the cube (E = −6J, −4J, −2J, 0, 2J, 4J, 6J). A generalized cluster algorithm for the 3-d FFIM Ising model is therefore very much more complicated than for the 2-d case.

As in two dimensions, the most important operations are those which conserve the energy. In Fig. 9 we show all possible freeze/delete operations for which flipping the bonded spins will conserve the energy of a ground state configuration (E = −6J). Here the unsatisfied links are dashed lines, and the frozen bonds are bold lines. Note that we only allow satisfied links to be frozen, (?? WHY ??) and it follows that it is not possible to apply these operations to any higher energy configurations. Note that here we show only one of the possible orientations of each operation in the unit cube. In fact each different orientation counts as a different operation, as with operations ?? and ?? in the 2-d algorithm (see Fig. 1. Initially we will clump all the orientations together into a single operation, with each orientation chosen with equal probability. Later we will investigate choosing them with different probabilities, as with the angled operations for the 2-d model in section ??.

We have a large amount of freedom in choosing the other (non-energy conserving) freeze/delete operations. There are 6 of them, since there are 6 energy states above the ground state. If we construct the energy/operation probability matrix with the 7 energy states, and the 6 operations of Fig. 10 plus one more operation representing any of the energy conserving operations of Fig. 9 (or an arbitrary linear combination of these), we will have a square (7×7) matrix which we can solve for the probabilities. In fact the operations have been carefully chosen such that if we tabulate the number of possible ways that the geometry allows us to perform operation i on a configuration of energy j, then we get the upper triangular form shown in Table II, which we have seen in section II is necessary for the probabilities to be positive at all values of the coupling.

We have tried to look at all possible freeze/delete operations which are likely to have a substantial contribution at low temperatures, specifically those which have zero probability except at energies E = −4J and −2J, i.e. the operations in Fig. 10(a) and (b). At the critical point these operations, along with the energy conserving ones, are used in about ??% of the cubes, and are the ones which are play the biggest part in determining the size of the largest cluster.

For the 4 other freeze/delete operations (Fig. 10(c-f)) we chose operations which are simple, and which reduce to the basic plaquette operations of Fig. 1 on the faces of the cube (i.e. plaquettes). This is by no means an exhaustive list of such operations, but we believe that these are among the best possible choices, and that choosing different operations will not change our results in any significant way. Later we will show that this is certainly the case for the lower energy operations of Fig. 9 and Fig. 10, in that varying the choice of these operations makes very little difference to the size of the biggest cluster at the critical point.

We used Mathematica [20] to solve the matrix equation 6, giving the probabilities in Table VII, all of which can be shown to be positive for K ≥ 0. We therefore have a valid algorithm which satisfies detailed balance, and is guaranteed to be ergodic (for T > 0) if we add a Metropolis sweep between each cluster update.

A problem with generalized cluster algorithms of this type is that there is no way of knowing a priori what kind of cluster size distribution they will give. The algorithm has to be implemented and run to obtain these results. We can however get an indication of the possible performance of the algorithm by looking at the zero temperature behavior, since the critical point for this model is at low temperature. At zero temperature we use only the energy conserving operations of Fig. 9, so the bond/link ratio is 0.5, which is unfortunately well above the percolation threshold. We are free to choose any linear combination of these 3 operations, and we have experimented with using the operations in the ratios r_{1}, r_{2}, and r_{3}. For all the choices of these parameters which we tried, the lattice froze into a single cluster at zero temperature.

However the bond/link ratio, and consequently the size of the largest cluster, will decrease as the temperature is increased. It is clearly unlikely that the bond/link ratio at K_{c} will decrease by the required factor of two to get near the percolation threshold, but since the two dimensional algorithm works very well with strongly percolating clusters, this may not be necessary for the algorithm to be effective.

Further runs were done at the infinite volume critical point of the model, again using the measured value K_{c} = 0.738(1) of Ref. [19]. Our results for the average maximum cluster size and average number of clusters at K_{c} for a 16^{3} lattice for various values of the ratios r_{i} are shown in Table VI. We also measured the bond/link ratio at K_{c} to be 0.4009(5) for this lattice size, which is still well above the percolation threshold, and as can be seen, the maximum cluster sizes are depressingly large for all of the choices of the ratios.

We also tried using some different freeze/delete operations which do not conserve energy. For operation 6, which is used on configurations of energy −6J and −4J, there is one other possibility, which is also shown in Fig. 10(a). We tried using this operation instead of the one or the other of these two possible operations, as well as a linear combination of the two, however the results obtained were virtually the same. Changing the other freeze/delete operations should have very little effect, since all of them combined are used on only about ??% of the elementary cubes at the critical point.

In all cases the average maximum cluster size is much too large to make effective changes in the configuration at each iteration. This can also be seen if we look at the autocorrelation times for the magnetization using the Metropolis and generalized cluster algorithms, which are shown in Fig. 7. One must also take into account the fact that one iteration of the cluster algorithm takes approximately ?? times longer to perform than a Metropolis iteration.

ADD THE LATEST WORK WE DID. THERE ARE EXTRA FREE PARAM-ETERS (ESPECIALLY FOR THE LOW ENERGY OPERATIONS) BECAUSE OF ORIENTATION. IN WHAT I SAID ABOVE WE CHOOSE THE ORIENTATION AT RANDOM IF WE HAVE MORE THAN ONE CHOICE. CAN WE BIAS THIS AT ALL?

VII. CONCLUSIONS

One of the main ideas of KBD was to extend the basic element of a cluster algorithm from a link connecting only two sites to a larger geometrical object such as a plaquette. This is clearly necessary for frustrated systems, for which knowledge of the frustration of the system cannot be obtained from looking at the interaction of two neighboring spins. However there are practical problems to this approach.

Finding large, regular collections of sites which will cover the lattice in such a way that each site is in one and only one such group is sometimes not an easy task. For example the honeycomb lattice does not allow a simple covering of this form, and we have seen that trying a larger covering for the triangular lattice than a simple triangular plaquette leads to other problems.

Even when a simple tiling can be found, as for the triangular and cubic lattices, there is no guarantee that a valid algorithm with positive probabilities for each possible freeze/delete operation can be constructed. If such an algorithm does exist, the resulting clusters may be too big, as in both the cases which we have investigated. Unfortunately it appears that although the KBD generalized cluster algorithm seems very promising as an approach to frustrated systems, and works very well for the square lattice fully frustrated Ising model, some new ideas are necessary before it can be successfully applied to general frustrated systems.

ACKNOWLEDGMENTS

We are grateful to John Apostolakis, Ofer Biham, Danny Kandel, Enzo Marinari and Alan Sokal for valuable discussions. This work was sponsored in part by Department of Energy grants DE-FG03-85ER25009 and DE-AC03-81ER40050, and by a grant from the IBM corporation.

REFERENCES

[1] R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 58, 86 (1987); J.-S. Wang and R. H. Swendsen, Physica A 167, 565 (1990).

[2] N. Metropolis et al., J. Chem. Phys. 21, 1087 (1953).

[3] K. Binder ed., Monte Carlo Methods in Statistical Physics, (Springer-Verlag, Berlin, 1986); K. Binder and D.W. Heermann, Monte Carlo Simulation in Statistical Physics, (Springer-Verlag, Berlin, 1988); H. Gould and J. Tobochnik, An Introduction to Computer Simulation Methods, Vol. 2, (Addison-Wesley, Reading, Mass., 1988).

[4] A. D. Sokal, in Computer Simulation Studies in Condensed Matter Physics: Recent Developments, eds. D. P. Landau et al. (Springer-Verlag, Berlin-Heidelberg, 1988); A. D. Sokal, in Proc. of the International Conference on Lattice Field Theory, Tallahassee, October 1990, Nucl. Phys. B (Proc. Suppl.) 20, 55 (1991).

[5] D. Kandel, R. Ben-Av and E. Domany, Phys. Rev. Lett. 65, 941 (1990).

[6] D. Kandel and E. Domany, Phys. Rev. A 43, 8539 (1991).

[7] J. Villain, J. Phys. C 10, 1717 (1977).

[8] R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 57, 2607 (1986).

[9] R. H. Swendsen and J.-S. Wang, Phys. Rev. B 38, 4840 (1988).

[10] D. Kandel, R. Ben-Av and E. Domany, Phys Rev B 45, 4700 (1992).

[11] C. F. Baillie and P. D. Coddington, Phys. Rev. B 43, 10617 (1991).

[12] G. H. Wannier, Phys. Rev. 79, 357 (1950).

[13] J. Stephenson, J. Math. Phys. 11, 420 (1970).

[14] D. Stauffer, Introduction to Percolation Theory, Taylor and Francis, London, 1985).

[15] D. Kandel, private communication.

[16] K. Binder and A. P. Young, Rev. Mod. Phys. 58, 801 (1986); D. Chowdhury, Spin Glasses and Other Frustrated Systems, (Princeton University Press, Princeton, NJ, 1986); M. Mezard and G. Parisi, Spin Glass Theory and Beyond, (World Scientific, Singapore, 1987).

[17] D. Blankschtein, M. Ma and A. N. Berker, Phys. Rev. B 30, 1362 (1984).

[18] G. S. Grest, J. Phys. C 18, 6239 (1985).

[19] H. T. Diep, P. Lallemand and O. Nagai, J. Phys. C 18, 1067 (1985).

[20] S. Wolfram, Mathematica: a system for doing mathematics by computer, (Addison-Wesley, Redwood City, Calif., 1990).

FIGURES

FIG. 1. Freeze/delete operations for square lattice FFIM.

FIG. 2. Log-log plot of autocorrelation times of M^{2} for the standard KBD algorithm for the square lattice FFIM at zero temperature.

FIG. 3. Log-linear plot of τ_{init} for the standard KBD algorithm for the square lattice FFIM at zero temperature. The dashed line shows the function log(log L).

FIG. 4. Log-log plot of autocorrelation times of M^{2} for the Metropolis algorithm, the standard KBD algorithm, and the KBD operation using just angled operations with three different values of the free parameter p_{sub}, for the square lattice FFIM at zero temperature.

FIG. 5. Autocorrelation function for the improved KBD algorithm for the square lattice FFIM at zero temperature.

FIG. 6. Freeze/delete operations for triangular lattice FFIM.

FIG. 7. Possible extensions of basic triangular plaquette: (a) double plaquette, (b) double plaquette plus extra link.

FIG. 8. The triangular lattice cannot be covered using only a double plaquette.

FIG. 9. All possible energy conserving freeze/delete operations for the cubic lattice FFIM.

FIG. 10. Some possible energy non-conserving operations for the cubic lattice FFIM. These correspond to operations which have zero probability for cubes with energies greater than −4J (a), −2J (b), 0 (c), +2J (d), and +4J (e), as well as the operation which deletes all bonds (f). Only (a) and (b) show all possible such operations.

TABLES

TABLE I. Probabilities for the freeze/delete operations of Fig. 1 for the square lattice FFIM.

TABLE II. Number of possible occurrences of the freeze/delete operations of Figs. 9 and 10 for the cubic lattice FFIM.

TABLE III. The average maximum cluster size and the average number of clusters for different values of the free parameter p_{sub} for the KBD algorithm using the angled operations (4-7) for the square lattice FFIM at zero temperature on a 64^{2} lattice. Also shown are the value for the standard KBD algorithm using only operations (1-3).

TABLE IV. Autocorrelation times of M^{2} for the standard KBD algorithm applied to the square lattice FFIM at zero temperature.

TABLE V. Probabilities for the freeze/delete operations of Fig. 6 for the triangular lattice FFIM.

TABLE VI. The average maximum cluster size for different ratios r_{1}, r_{2}, r_{3} of the 3 energy conserving operations at the critical point of the 3-d FFIM, for a 16^{3} lattice.

TABLE VII. Probabilities for the freeze/delete operations of Figs. 9 and 10 for the cubic lattice FFIM at the critical point (β = 0.738). These particular probabilities are obtained using only the first of the possible operations shown in the figures.