%Filename: paper.tex
\documentstyle[preprint,revtex]{aps}
\begin{document}
\draft
% \preprint{SCCS-???}
%
\begin{title}
On generalized cluster algorithms for frustrated spin models
\end{title}
%
\author{P.D. Coddington$^{\dag}$ and L. Han}
\begin{instit}
Physics Department, Syracuse University, Syracuse, NY 13244, USA
\end{instit}
%
\begin{center}
\bigskip
%
\today
%
\bigskip
\end{center}
\begin{abstract}
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.
\end{abstract}
% \pacs{PACS numbers 05.50.+q, 75.40.Mg}

\vfill
\noindent
{\tenrm 
$^{\dag}$ Current address: Northeast Parallel Architectures Center,
111 College Place, Syracuse University, Syracuse, NY 13244, USA.}

\newpage
\narrowtext

\section{INTRODUCTION}
\label{sec:intro}

The cluster update algorithm of Swendsen and Wang (SW)~\cite{SW} 
provides a great improvement in the efficiency of Monte Carlo simulations 
over the Metropolis algorithm~\cite{metrop,monte} and other local update
schemes by dramatically reducing critical slowing down~\cite{monte,sokalrev} 
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~\cite{kandel90,kandel91a} which works very well for the two 
dimensional fully frustrated Ising model~\cite{villain}. 
The only other cluster algorithm which has been shown to work effectively for 
a frustrated spin model is the Replica Monte Carlo algorithm~\cite{replica},
however this has only been successfully applied to the 2-$d$ Ising spin
glass~\cite{repglass}.
% SOME MORE WAFFLE HERE?

In section~\ref{sec:algorithm} 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~\ref{sec:square}, 
the triangular lattice FFIM in (section~\ref{sec:triangle}), 
some other 2-$d$ Ising models in section~\ref{sec:others},
and finally to the 3-$d$ FFIM in section~\ref{sec:cube}.


\section{FORMULATING A GENERALIZED CLUSTER ALGORITHM}
\label{sec:algorithm}

First we will briefly describe the SW cluster algorithm for the Ising 
model~\cite{monte}, which has an interaction energy given by
%
\begin{equation}
E = - \sum_{<i,j>} J_{ij} \sigma_{i} \sigma_{j},            \label{eq:energy}
\end{equation}
%
where the spins $\sigma_{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 
$1 \over 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 anti-ferromagnetic ($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.

% *************** ?????????????????????????????? ***************************
% EXPLAIN KBD'S GENERALIZED CLUSTER ALGORITHM METHOD.

The generalized cluster algorithm of KBD begins by expressing the energy of the
system as a sum $E = \sum_{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
%
\begin{equation}
\sum_{i} P_{i}^{l}(u) = 1, \qquad 0 \le P_{i}^{l}(u) \le 1,  \label{eq:norm}
\end{equation}
%
for all $l$ and $u$.

In order to satisfy detailed balance, KBD show that the probabilities must
also satisfy
%
\begin{equation}
E_i(u) - {1 \over \beta} \log P_{i}(u) + C_i = 0,          \label{eq:detbal}
\end{equation}
%
$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
%
\begin{equation}
P_{i}(u) = e^{\beta E_i} \quad {\rm or} \quad 0.           \label{eq:probs}
\end{equation}
%
Imposing the normalization condition of equation~\ref{eq:norm}, we have
\begin{equation}
?????						           \label{eq:normprob}
\end{equation}
which can be written as a matrix equation
%
\begin{equation}
{\bf B} \underline{\bf c} = {\bf 1},                       \label{eq:matrix}
\end{equation}
%
where the vector {\bf c} is made up of the $C_i$, and the elements 
$B_{ij}$ of the matrix {\bf B} are either $e^{\beta E_i}$ or 0,
where $j$ labels ??? (the u's??).


% *************** ?????????????????????????????? ***************************
EXPLAIN HOW TO USE IT, BY SETTING UP A MATRIX EQUATION. OUTLINE 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.
% *************** ?????????????????????????????? ***************************



\section{THE SQUARE LATTICE FULLY FRUSTRATED ISING MODEL}
\label{sec:square}

% *************** ?????????????????????????????? ***************************
% DO ALL THIS:
% Introduce the model, talk about the SW algorithm and percolation probs
% and bond/link rations (take all this stuff from the next section), then
% get on to the KBD algorithm, the operations, and then say
% how they do it, how they get clusters (fluke?), and how you can do
% it more generally. Give our results for the angle operation.
% *************** ?????????????????????????????? ***************************

The fully frustrated Ising model (FFIM) was introduced by 
Villain~\cite{villain}
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~\cite{kandel90,kandel91b}.


\subsection{The KBD Algorithm}
\label{subsec:kbd}

In the KBD algorithm there are seven basic freeze/delete operations, which 
are shown in Fig.~\ref{ops_sq}, 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~\ref{tab:probs_sq}.
(SHOW HOW WE GET THIS USING THE MATRIX FORMALISM, AS AN EXAMPLE??)

% \input fig:ops_sq.tex
% \input tab:probs_sq.tex

Note that different orientations of the same basic tolopogy of links, 
such as operations 2 and 3, or operations 4-7, in Fig.~\ref{ops_sq},
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 \times 64$).

The main reason for this is that the square lattice has two energy conserving
operations (numbers 2 and 3 in Fig.~\ref{ops_sq}), 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.


\subsection{Dynamic Critical Exponents}
\label{subsec:dynamics}

KBD measured a dynamic critical exponent of $z \approx 0.55$ 
for their algorithm~\cite{kandel91b}, compared to $z \approx 2$ for the
Metropolis algorithm
(see Refs.~\cite{monte,sokalrev} 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~\cite{paul}.
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 
%
\begin{equation}
\rho_A(t) = {<A(0)A(t)> - <A(0)>^2 \over <A(0)A(0)> - <A(0)>^2}
\label{eq:rho}
\end{equation}
%
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 $\tau_{exp}$, via 
%
\begin{equation}
\rho(t) \sim e^{-t/\tau_{exp}},
\label{eq:tau_exp}
\end{equation}
%
as well as integrated autocorrelation times
%
\begin{equation}
\tau_{int} = {1 \over 2} + \sum_{t=1}^{\infty} \rho_A(t),
\label{eq:tau_int}
\end{equation}
%
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.~\cite{paul}.

In Table~\ref{tab:tau_kbd} 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 \times 10^6$ sweeps. 
Fig.~\ref{tau_kbd} shows a log-log plot of these results. 
Straight line $\chi^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.~\ref{tau_kbd}), so that $z_{int}$ could also be zero.
The same situation also occurs for the autocorrelations of cluster algorithms 
for the ferromagnetic Ising model~\cite{paul}, since it is very 
difficult to differentiate between a logarithm and a small power.

% \input fig:tau_kbd.tex
% \input tab:tau_kbd.tex

Let us also define $\tau_{init}$ to be the autocorrelation time obtained from 
an exponential fit to $\rho(t)$ at only two points, $t=0$ and $t=1$.
This is plotted in Fig.~\ref{tau_init} on a log-linear plot. 
We can see that $\tau_{init}$
grows extremely slowly with $L$ (the dashed line shows the function 
$\log(\log L)$). Thus is in marked contrast to $\tau_{exp}$, which grows as a
substantial power of $L$. Unfortunately this means that it is very 
difficult to get a good asymptotic fit to $\rho(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 $\tau_{exp}$. The fact 
that $\tau_{init}$ increases so slowly means that the $\tau_{int}$ grows
much slower than $\tau_{exp}$, and consequently $z_{int}$ appears to be
smaller than $z_{exp}$.

% \input fig:tau_init.tex

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 $\tau_{int}$ by a factor of 2, however this is not the case. In fact
the autocorrelation function for this algorithm is {\it oscillatory},
as shown in Fig.~\ref{new_kbd}. Thus using the definition for
$\tau_{int}$ (equation~\ref{eq:tau_int}) gives a value close to zero, 
and much less than half of the value for the standard KBD algorithm.

% \input fig:new_kbd.tex

% *************** ?????????????????????????????? ***************************
WHY DOES THIS GIVE NEGATIVE AUTOCORRELATIONS? 
(LOOKS A BIT LIKE WHAT HAPPENS FOR THE SW REPLICA METHOD). 
DO THE FORMULAS BREAK DOWN IF THE AUTOCORRELATION FUNCTION IS NEGATIVE?
CHECK SOKAL'S PAPERS.
% *************** ?????????????????????????????? ***************************



\subsection{A Variant of the KBD Algorithm}
\label{subsec:var_kbd}

% ?????????????
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.
% , the latter being the one used by KBD~\cite{kandel90}. 
% NO - they use 0.5 * p'.
Integrated autocorrelation times for these 
three variations of the KBD algorithm are shown in Fig.~\ref{tau_all}, 
along with results for the Metropolis algorithm, 
and the KBD results from Fig~\ref{tau_kbd}.

% \input fig:tau_all.tex

The results are quite surprising. As expected, $p_{sub} = 0.5$ gives results
which are very similar to the Metropolis algorithm ($z \approx 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 \approx 0.??$, while $p_{sub} = 1.0$ has $z \approx 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.


\section{THE TRIANGULAR LATTICE ISING ANTIFERROMAGNET}
\label{sec:triangle}

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~\cite{wannier,stephenson}.

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~\cite{perc}.
%
% *************** ?????????????????????????????? ***************************
% which is also the ratio of frozen bonds to links for the SW algorithm for the
% Ising ferromagnet at its critical point, since in two dimensions the SW 
% clusters percolate there (REF??).
%
% This doesn't seem to be true for the triangular lattice!
% *************** ?????????????????????????????? ***************************
%
For the square lattice the bond percolation threshold is 0.5~\cite{perc},
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.~\ref{ops_tri}. 
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~\cite{kandel91a}, or by solving
the matrix equation constructed using the method of section~\ref{sec:algorithm},
(DO THIS EXPLICITLY AS AN EXAMPLE??)
gives the probabilities given in Table~\ref{tab:probs_tri}.

% \input fig:ops_tri.tex
% \input tab:ops_tri.tex

% *************** ?????????????????????????????? ***************************
% DO THIS EXPLICITLY!
%
% If there are two satisfied links, we follow KBD and chose one at random as
% the position of the frozen bond. 
% EXPLAIN WHY YOU *HAVE* TO CHOOSE THE BONDS WITH THE SAME PROBABILITY,
% UNLIKE THE SQUARE LATTICE CASE.  IT'S WORTH MAKING THIS POINT FOR THE
% SQUARE LATTICE CASE AS WELL, PERHAPS? 
%
% *************** ?????????????????????????????? ***************************

At the zero temperature critical point, we only perform the energy conserving
freeze/delete operations. This gives a bond/link ratio of $1 \over 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.~\ref{fig:dbltri}a. 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.~\ref{fig:tricover}, 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??).

% *************** ?????????????????????????????? ***************************
%
% is this true? N.B. Have 2 ground state energies here, too!
%
% *************** ?????????????????????????????? ***************************

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

% *************** ?????????????????????????????? ***************************
% Both? Or only this one?
% *************** ?????????????????????????????? ***************************

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.~\ref{fig:dbltri}b, 
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??)

% *************** ?????????????????????????????? ***************************
%
% (except for the plaq + SW version??) Is this true? True in general?
%
% 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 {\it et al.} have also tried to apply their method to this model,
with the same negative result~\cite{danny}.

%
% This idea doesn't satisfy detailed balance! 
% Should I say anything about it?
%
% One idea we tried to improve the algorithm is to bias the direction of the
% bonds, in the hope that providing less connectivity in certain directions 
% might split up the large cluster.
% Thus instead of freezing one of the two satisfied links at random, we choose
% an ordering of the three link directions, say $d_1$, $d_2$, $d_3$. 
% Now, if one of the satisfied links is in direction $d_1$, we delete it and 
% freeze the other one. If on the other hand the unsatisfied link is in 
% direction $d_1$, then
% we choose the link in direction $d_2$ as the one to be deleted with a
% probability $p_{dir}$ which we are free to vary. In this way we are biasing
% the bonds strongly against the direction $d_1$ and weakly against the 
% direction $d_2$ (if $p_{dir}$ is chosen to be greater than 0.5). We might
% therefore hope to split the largest cluster (at least partially) along the 
% direction $d_3$.
% The choice for the ordering is done at random before each iteration 
% (just like the choice of black or white plaquettes from a checkerboard), 
% so as not to bias the configurations.
% 


\section{OTHER 2-D FRUSTRATED MODELS}
\label{sec:others}

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 $\pm J$ spin glass, where the couplings $J_{ij}$ in 
equation~\ref{eq:energy} are chosen to be
ferromagnetic ($+J, J>0$) or antiferromagnetic ($-J$) at random~\cite{glass}.
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.

% *************** ?????????????????????????????? ***************************
% Do full GCA formalism on this - does it give the same thing?
% *************** ?????????????????????????????? ***************************

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 co-ordination 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
co-ordination 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.



\section{THE 3-D FULLY FRUSTRATED ISING MODEL}
\label{sec:cube}

We next attempted to apply the generalized cluster Monte Carlo method to
the three dimensional FFIM~\cite{BMB,grest,diep}.
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~\cite{grest,diep}. 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~\cite{perc}. 
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$~\cite{diep}) 
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.~\ref{ops_cube_Econs} 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 ??)
%
% 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.~\ref{ops_sq}. 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~\ref{sec:square}.

% \input fig:ops_cube_Econs.tex

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.~\ref{ops_cube_others} 
plus one more operation representing any of the energy conserving operations 
of Fig.~\ref{ops_cube_Econs} (or an arbitrary linear combination of these),
we will have a square ($7 \times 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~\ref{tab:numops_cube},
which we have seen in section~\ref{sec:algorithm} is necessary for the 
probabilities to be positive at all values of the coupling. 

% \input fig:ops_cube_others.tex
% \input tab:numops_cube.tex

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.~\ref{ops_cube_others}(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.~\ref{ops_cube_others}(c-f)) 
we chose operations which are simple, and which 
% in most cases
reduce to the basic plaquette
operations of Fig.~\ref{ops_sq} 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.~\ref{ops_cube_Econs} and Fig.~\ref{ops_cube_others},
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~\cite{mathematica} to
solve the matrix equation~\ref{eq:matrix}, 
giving the probabilities in Table~\ref{tab:probs_cube},
all of which can be shown to be positive for $K \ge 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.

% \input tab:probs_cube.tex

A problem with generalized cluster algorithms of this type is that there
is no way of knowing {\it 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.~\ref{ops_cube_Econs}, 
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.~\cite{diep}. 
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~\ref{tab:max_csize_Econs}. 
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.~\ref{ops_cube_others}(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 
% function?
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 PARAMETERS (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?
% *************** ?????????????????????????????? ***************************


\section{CONCLUSIONS}
\label{sec:conclusion}

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
%
\begin{references}
\bibitem{SW}
R.~H.~Swendsen and J.-S.~Wang, Phys. Rev. Lett. {\bf 58}, 86 (1987);
J.-S.~Wang and R.~H.~Swendsen, Physica A {\bf 167}, 565 (1990).
% \bibitem{SWanti}
% J.-S.~Wang, R.H.~Swendsen and R.~Kotecky, 
% Phys. Rev. Lett. {\bf 58}, 86 (1987).
% \bibitem{SWfrust}
% V.S.~Dotsenko, W.~Selke, and A.L.~Talapov,
% Physica A {\bf 170}, 278 (1991).
\bibitem{metrop}
N.~Metropolis {\it et al.}, J. Chem. Phys. {\bf 21}, 1087 (1953).
\bibitem{monte}
K.~Binder ed., {\it Monte Carlo Methods in Statistical Physics}, 
(Springer-Verlag, Berlin, 1986);
K.~Binder and D.W.~Heermann, {\it Monte Carlo Simulation in Statistical 
Physics}, (Springer-Verlag, Berlin, 1988);
H.~Gould and J.~Tobochnik,
{\it An Introduction to Computer Simulation Methods, Vol. 2},
(Addison-Wesley, Reading, Mass., 1988).
\bibitem{sokalrev}
A.~D.~Sokal, in {\em Computer Simulation Studies in Condensed Matter
Physics:  Recent Developments}, eds.\ D.~P.~Landau {\em 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.) {\bf 20}, 55 (1991).
\bibitem{kandel90}
D.~Kandel, R.~Ben-Av and E.~Domany, Phys. Rev. Lett. {\bf 65}, 941 (1990).
\bibitem{kandel91a}
D.~Kandel and E.~Domany, Phys. Rev. A {\bf 43}, 8539 (1991).
\bibitem{villain}
J.~Villain, J. Phys. C {\bf 10}, 1717 (1977).
\bibitem{replica}
R.~H.~Swendsen and J.-S.~Wang, Phys. Rev. Lett. {\bf 57}, 2607 (1986).
\bibitem{repglass}
R.~H.~Swendsen and J.-S.~Wang, Phys. Rev. B {\bf 38}, 4840 (1988).
\bibitem{kandel91b}
D.~Kandel, R.~Ben-Av and E.~Domany, Phys Rev B {\bf 45}, 4700 (1992).
\bibitem{paul}
C.~F.~Baillie and P.~D.~Coddington, Phys. Rev. B {\bf 43}, 10617 (1991).
\bibitem{wannier}
G.~H.~Wannier, Phys. Rev. {\bf 79}, 357 (1950).
\bibitem{stephenson}
J.~Stephenson, J. Math. Phys. {\bf 11}, 420 (1970).
\bibitem{perc}
D.~Stauffer, {\it Introduction to Percolation Theory},
Taylor and Francis, London, 1985).
\bibitem{danny}
D.~Kandel, private communication.
\bibitem{glass}
K.~Binder and A.~P.~Young, Rev. Mod. Phys. {\bf 58}, 801 (1986);
D.~Chowdhury, {\it Spin Glasses and Other Frustrated Systems},
(Princeton University Press, Princeton, NJ, 1986);
M.~Mezard and G.~Parisi, {\it Spin Glass Theory and Beyond},
(World Scientific, Singapore, 1987).
\bibitem{BMB}
D.~Blankschtein, M.~Ma and A.~N.~Berker, Phys. Rev. B {\bf 30}, 1362 (1984).
\bibitem{grest}
G.~S.~Grest, J. Phys. C {\bf 18}, 6239 (1985).
\bibitem{diep}
H.~T.~Diep, P.~Lallemand and O.~Nagai, J. Phys. C {\bf 18}, 1067 (1985).
\bibitem{mathematica}
S.~Wolfram, {\it Mathematica: a system for doing mathematics by computer},
(Addison-Wesley, Redwood City, Calif., 1990).
\end{references}
%
%
% FIGURES --- put these in the appropriate place in the text.
%
% ********* section square **************
%
\figure{Freeze/delete operations for square lattice FFIM.
\label{ops_sq}}
%
\figure{Log-log plot of autocorrelation times of $M^2$ for the standard KBD
algorithm for the square lattice FFIM at zero temperature. 
\label{tau_kbd}}
%
\figure{Log-linear plot of $\tau_{init}$ for the standard KBD algorithm for 
the square lattice FFIM at zero temperature. 
The dashed line shows the function $\log(\log L)$.
\label{tau_init}}
%
\figure{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. 
\label{tau_all}}
%
\figure{Autocorrelation function for the improved KBD algorithm 
for the square lattice FFIM at zero temperature. 
\label{new_kbd}}
% 
%
% ********* section triangle **************
%
\figure{Freeze/delete operations for triangular lattice FFIM.
\label{ops_tri}}
%
\figure{Possible extensions of basic triangular plaquette: (a) double
plaquette, (b) double plaquette plus extra link.
\label{fig:dbltri}}
%

\figure{The triangular lattice cannot be covered using only a double plaquette.
\label{fig:tricover}}
%
%
% ********* section cube **************
%
\figure{All possible energy conserving freeze/delete operations for the 
cubic lattice FFIM.
\label{ops_cube_Econs}}
%
\figure{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.
\label{ops_cube_others}}
%
%
% TABLES --- put these in the appropriate place in the text.
%
\begin{table}
\caption{Probabilities for the freeze/delete operations of Fig.~\ref{ops_sq}
for the square lattice FFIM.
\bigskip}
\label{tab:probs_sq}
\end{table}
%
\begin{table}
\caption{Number of possible occurrences of the freeze/delete operations of 
Figs.~\ref{ops_cube_Econs} and~\ref{ops_cube_others}
for the cubic lattice FFIM.
\bigskip}
\label{tab:numops_cube}
\end{table}
%
\begin{table}
\caption{
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).
\bigskip}
\label{tab:max_csize}
\end{table}
%
\begin{table}
\caption{Autocorrelation times of $M^2$ for the standard KBD algorithm 
applied to the square lattice FFIM at zero temperature.
\bigskip}
\label{tab:tau_kbd}
\end{table}
%
\begin{table}
\caption{Probabilities for the freeze/delete operations of Fig.~\ref{ops_tri}
for the triangular lattice FFIM.
\bigskip}
\label{tab:probs_tri}
\end{table}
%
\begin{table}
\caption{
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.
\bigskip}
\label{tab:max_csize_Econs}
\end{table}
%
\begin{table}
\caption{Probabilities for the freeze/delete operations of 
Figs.~\ref{ops_cube_Econs} and~\ref{ops_cube_others}
for the cubic lattice FFIM at the critical point ($\beta = 0.738$).
These particular probabilities are obtained using only the first of the
possible operations shown in the figures.
\bigskip}
\label{tab:probs_cube}
\end{table}
%

\end{document}


