%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 that works
extremely well for the fully frustrated Ising model on a square lattice. 
Here we study some properties of this algorithm and some simple extensions
to it.  We introduce a practical methodology for applying the generalized 
cluster algorithm to a given spin model, and investigate applications of
the algorithm to some other frustrated Ising models in two dimensions.
We find that such algorithms work poorly or not at all for these models.
\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 computational efficiency of Monte Carlo 
simulations over the
Metropolis algorithm~\cite{metrop,monte} and other standard 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, the systems for which 
critical slowing down can be the most extreme~\cite{glass}. The reason for
this failure is that at the critical point, the SW algorithm produces a cluster
that encompasses a large fraction (if not all) of the sites in
the lattice, so that updating the spins in this cluster is virtually the
same as a trivial global spin update.

Kandel, Ben-Av and Domany (KBD) introduced a generalized cluster
algorithm~\cite{kandel90,kandel91} that works very well for the two 
dimensional fully frustrated Ising model~\cite{kandel90,kandel92}. 
The only other cluster algorithm that 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}.
It is therefore of great interest to determine whether the generalized
cluster algorithm can be applied to other frustrated spin models.
We have investigated some variants of the KBD algorithm for the 2-$d$
fully frustrated Ising model (FFIM), and attempted to apply the method
to some other frustrated two dimensional Ising models.

We introduce a practical method for formulating a generalized cluster 
algorithm for a given spin model in section~\ref{sec:algorithm}. 
We then apply this methodology to 
the FFIM on a square lattice in section~\ref{sec:square}, 
the triangular lattice Ising anti-ferromagnet in 
section~\ref{sec:triangle}, and 
some other 2-$d$ frustrated Ising models in section~\ref{sec:others},
and analyze the performance of generalized cluster algorithms for each of
these models.

\newpage

\section{THE GENERALIZED CLUSTER ALGORITHM}
\label{sec:algorithm}

First we 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 > 0$.
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, bonds are introduced 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 failure
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 on average, 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 some other group of sites which can give a tiling of the 
lattice, that is, they can be replicated over the lattice in such a way
that each link is part of one (and only one) subregion. 
We will refer to such a subregion as a {\it tile}.
For each tile $l$ we stochastically assign one of $n$ possible operations with
a probability $P_{i}^{l}(u)$, that depends on the spin configuration $u$ 
and the operation $i$. These operations involve placing bonds between
certain sites in the region $l$. 
We refer to the placement (or non-placement) of bonds as {\it freezing} 
(or {\it deleting}) the connections between two neighboring sites.
The frozen bonds produce 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$.
KBD show that detailed balance is satisfied if the probabilities 
also satisfy~\cite{kandel91}
%
\begin{equation}
E(u) - {1 \over \beta} \log P_{i}(u) + C_i = 
\left\{ \begin{array}{ll}
	0      \quad & \mbox{for an allowed operation $i$} \\
        \infty \quad & \mbox{otherwise}
	\end{array}
\right.                                                       \label{eq:detbal}
\end{equation}
%
Here $E(u)$ is the energy for the configuration $u$, 
$C_i$ is a constant (independent of the configuration), 
$\beta = 1 / kT$ is the inverse temperature, and we have
supressed the subscript $l$ labeling the particular tile.
For our purposes we will define an operation as ``allowed'' if it does
not freeze any unsatisfied links of the configuration $u$
(for a more general definition, see Ref.~\cite{kandel91}). 
If we define $a_i(u)$ to be 1 if operation $i$ is allowed for configuration 
$u$, and 0 otherwise, 
then equation~\ref{eq:detbal} can be written as
%
\begin{equation}
P_{i}(u) = a_{i}(u) \, e^{\beta(E(u) + C_i)}.               \label{eq:probs}
\end{equation}
%
Imposing the normalization condition of equation~\ref{eq:norm}, we have
\begin{equation}
\sum_{i} a_i(u) \, c_i = e^{-\beta E(u)}  \quad \forall \, u\, , 
\qquad \quad c_i = e^{\beta C_i}.                           \label{eq:normprob}
\end{equation}
This can be conveniently written as a matrix equation
%
\begin{equation}
{\bf A} \, \underline{\bf c} = \underline{\bf b},           \label{eq:matrix}
\end{equation}
%
where the $c_i$ are the elements of the vector \underline{\bf c}, 
the Boltzmann factors $b_j = e^{-\beta E(u)}$ are the elements of the vector
\underline{\bf b}, and the $a_i(u)$ make up the elements $A_{ij}$ of the matrix {\bf A}, 
where $j$ labels the different possible configurations $u$.

Formulating a generalized cluster algorithm involves identifying some
candidate operations, constructing the matrix {\bf A}, and then solving 
equation~\ref{eq:matrix} for \underline{\bf c}, from which the probabilities
for each operation can be calculated via equation~\ref{eq:probs}. 
Thus to guarantee a solution, the number of operations should be greater than
or equal to the number of possible configurations. If the number of operations
is greater than the number of configurations, there is a possibility of
having free parameters in the solution. 
Also, if an operation $i$ only occurs for a single configuration, there are
no constraints on the choice of $C_i$, which also allows for a free parameter
in choosing the probabilities.
(IS THIS TRUE IN GENERAL, OR ONLY TRUE FOR ENERGY CONSERVING OPERATORS,
AS KDB CLAIM???).
This may allow some optimization of the algorithm,
such as tuning the size of the largest cluster.

There is no guarantee that a solution will exist for which the probabilities 
are all between 0 and 1 for an arbitrarily chosen set of operations.
% However there is a simple method of choosing the operations to 
% ensure (???????) that this constraint is satisfied.
% (WHAT ABOUT NON-ZERO TEMP??) 
% Let us consider what happens as the temperature approaches zero 
% (i.e. $\beta \rightarrow \infty$).
% From equation~\ref{eq:probs} we find that in order to avoid having either
% infinite probabilities or all probabilities being zero, we must choose
% $C_i = -E_{max}^i$, where $E_{max}^i$ is the maximum
% energy of all configurations $u$ for which operation $i$ is allowed
% (i.e. $a_i(u) \neq 0$). Thus at zero temperature $P_i(u)$ is 0
% if $E(u) < E_{max}^i$, and 1 if $E(u) = E_{max}^i$. To satisfy 
% equation~\ref{eq:norm}, we must therefore choose only operations for which
% there is only a single configuration of energy $E_{max}^i$ for which
% the operation $i$ is allowed.

Another constraint on the algorithm is that the number of operations
which conserve energy must be greater than or equal to the number of
possible ground state configurations, otherwise the equations will be
over-determined at zero temperature. 
%
Finally, the matrix {\bf A} must have non-zero determinant in order for
a solution to exist. (BUT DOESN'T DET=0 MEAN NO SOLUTION OR MULTIPLE
SOLUTIONS?????)


% BUT WHAT HAPPENS IN THE CASE OF CONFIGS OF HIGHER ENERGY, WHICH DON'T
% OCCUR AT T=0? THEY STILL OCCUR WITH FINITE PROBABILITY AS T -> 0, SO 
% THIS ARGUMENT STILL GOES THROUGH.


\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 lacking the extra complication 
of disorder that is present in systems such as spin glasses.
On a square lattice, the FFIM has ferromagnetic couplings on all links
except for a line of anti-ferromagnetic couplings on every second column 
of links. This means that every plaquette is frustrated, having 3
ferromagnetic links and one anti-ferromagnetic link.
 
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 bonds all the sites in the lattice together into
a single cluster at the critical point, 
so the SW algorithm does not work for this model.
However the generalized cluster algorithm of KBD, that deletes bonds using
information on the state of plaquettes rather than links, works extremely 
well in this case~\cite{kandel90,kandel92}.


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

We will derive the KBD algorithm in a slightly different way, so as to
make the generalization to other frustrated models more obvious, and
to highlight the simple matrix methodology outlined in 
section~\ref{sec:algorithm}. 
 
In the KBD algorithm, the tile used is a plaquette. 
There are seven freeze/delete operations, as
shown in Fig.~\ref{ops_sq}, where the bold lines indicate frozen bonds.
Operations 2-7 conserve the energy of a configuration.
% Such energy conserving operations must be included in any list of operations,
% so that zero temperature updates can occur.
Note that operations 2-3 and operations 4-7 are just different orientations 
of the same basic topology of links.
To simplify matters, we will classify the operations into the 3 basic 
topologies (A, B, and C in Fig.~\ref{ops_sq}) which we will refer to as 
{\it operators}. Each of the operators has some number
of possible orientations: 1 for operator A, 2 for operator B, and 4 for
operator C. 
However, the important feature of these operators for constructing a 
generalized cluster algorithm is the number of orientations that are possible 
for any given configuration of spins in a plaquette. 
For example, consider a plaquette having the lowest energy, $-2J$. This has one
unsatisfied link. Since we can only freeze satisfied links, this means there
is only one allowed orientation of operator B (the
one with the deleted bonds corresponding to the unsatisfied link and the
link parallel to it), and two allowed orientations of operator C (deleting
the unsatisfied link and one of the links at right angles to it on either side).

% \input fig:ops_sq.tex

Likewise, the important
property of a plaquette configuration is not the spins, but rather which of
the links are satisfied, since that alone determines which freeze/delete
operations can be applied. Again, there are only a few basic topologies of
satisfied and unsatisfied links -- for the square lattice FFIM there are
only two (either 1 or 3 unsatisfied links), in general there will be at
least as many as there are possible energy states of the tile. All
other configurations are just different orientations of these basic link
topologies. From now on we will use the term ``configuration'' to mean a 
configuration of satisfied and unsatisfied links, rather than spins.
% , and ``link topology'' to refer to these basic link topologies.

To simplify matters, we will consider different orientations of the 
configurations and the freeze/delete operators to be equivalent. 
In order to use the matrix formalism outlined in section~\ref{sec:algorithm}, 
we need to know the number of different link topologies and the
number of allowed orientations $n_{ij}$ of each operator $i$ for any link 
topology $j$. If we adopt this approach, the 
elements of the matrix {\bf A} will be $n_{ij}$.
The matrix for the KBD algorithm is shown in Table~\ref{tab:matrix_sq}.

Solving the matrix equation~\ref{eq:matrix} 
(which is easily done using Mathematica~\cite{mathematica})
gives the probabilities $P_{i}(j)$ given in Table~\ref{tab:probs_sq}. 
Since for this model we have 3 operators and only 2 link topologies,
the solution to the matrix equation has a free parameter $p$. 
KBD use this freedom in the choice of probabilities to prove 
(using a subtle geometrical argument) that for a particular choice of the
probabilities, the lattice is always split into at least two large clusters
at the zero temperature critical point~\cite{kandel92}.
For the probabilities chosen by KBD, the average maximum cluster size is 
neither too large nor too small (the average being 0.6432(5) for a 
$64^2$ lattice), so the KBD cluster algorithm works extremely well.

% \input tab:matrix_sq.tex
% \input tab:probs_sq.tex

Note that the probabilities in Table~\ref{tab:probs_sq} refer to the 
probability of choosing operator $i$. Each operator has a number of
allowed orientations, which would then be assigned with equal probability,
$1 / n_{ij}$. For example, if we chose operator C, then we would randomly
choose with probability $1 \over 2$ to delete either of the 2 bonds 
perpendicular to the unsatisfied link.

When we set up the matrix equation, we could have treated all possible 
orientations of the configurations and the operators as different, rather
than grouping them together as we have done. This would give 4 
lowest energy configurations (rather than just 1), corresponding to the 
4 possible positions of the single unsatisfied link in the plaquette.
It would also give the 6 energy conserving operations of 
Fig.~\ref{ops_sq} (rather than the 2 basic operators). 
So instead of having a single free parameter, from having 2 operators
and only 1 configuration, we now have 2 free parameters, from having
6 operations and only 4 configurations. If we solve the extended
matrix equation, we find that the free parameter enters (as one might expect) 
in the probability of choosing which of the bonds perpendicular to the
unsatisfied link are to be deleted for operator C.  
We now have an extra probability $p_{sub}$ of choosing the operations from
the subgroup (4-5) of the angled operations, 
and $1-p_{sub}$ of choosing the operations from the subgroup (6-7), 
where previously we had $p_{sub} = {1 \over 2}$.

Both these approaches -- using the basic topologies or using all possible
orientations --  are correct, i.e. satisfy detailed balance.
The first is simpler and allows a simpler implementation
of more complex problems where there are a large number of possible
topologies and orientations; 
the other may provide some extra freedom in the choice of probabilities,
that can allow some tuning to improve the performance of the algorithm, 
as we shall see in section~\ref{subsec:var_kbd}.

Although the KBD algorithm satisfies detailed balance, it is not ergodic
at zero temperature -- that is, starting from a particular ground state
configuration, it cannot generate all other ground state configurations.
In fact, it can readily be seen for small lattice sizes that the KBD algorithm
cycles between subsets of possible ground state configurations. In order
to ensure ergodicity, each KBD update must be followed by a Metropolis
sweep. Note also that for the Metropolis update to be ergodic for this model,
the site to be updated must be chosen in random, not serial, order.


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

KBD measured the dynamic critical exponent for the exponential autocorrelation
time to be $z \approx 0.55$ for their algorithm~\cite{kandel92},
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 $z$ from
such data~\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 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}
%
that 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 $7.5 \times 10^6$ iterations
(more for smaller lattice sizes). 
Fig.~\ref{tau_kbd} shows a log-log plot of these results. 
Straight line $\chi^2$ fits give 
$z_{int} = 0.28(1)$ and
$z_{exp} = 0.66(5)$.
The integrated autocorrelations also fit fairly well to a logarithm,
so that $z_{int}$ could also be zero.
The same situation 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 define $\tau_{0}$ 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_0} on a log-linear plot. 
We can see that $\tau_{0}$
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 makes it 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_{0}$ increases so slowly means that the $\tau_{int}$ grows
much slower than $\tau_{exp}$, and consequently $z_{int}$ is much 
smaller than $z_{exp}$.

% \input fig:tau_0.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, so there is no change in
the configuration apart from a trivial global spin flip. Only half of the
iterations make a non-trivial update of the configuration by 
updating one of the clusters and not the other. 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 {\it not} the case. In fact
the autocorrelation function for this algorithm alternates in sign,
as shown in Fig.~\ref{new_kbd}. 
The autocorrelation function for a Markov process can generally be expressed
asymptotically as $\rho(t) = a^t$, where $|a| < 1$~\cite{timeseries}.
If $a > 0$, then redefining $a = e^{-1/\tau_{exp}}$ gives the standard 
exponential form for $\rho(t)$.
However it is also possible to have $a < 0$, in which 
case the autocorrelation function still falls off exponentially, but 
is proportional to $(-1)^t$. This indicates that the original data forms an
alternating time series, with successive measurements falling predominantly
on one side and then the other of the mean value~\cite{timeseries}.
%
Since $\rho(t)$ alternates in sign, 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
with an exponential autocorrelation function.
(SO HOW DOES THIS WORK????)

% \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.
% *************** ?????????????????????????????? ***************************

Another interesting result concerns the autocorrelation function $\rho(t)$ for
the KBD algorithm {\it without} the Metropolis sweep needed for ergodicity,
which is shown in Fig.~\ref{tau_ergodic}, along with $\rho(t)$ for the
algorithm including Metropolis sweeps. The latter falls off asymptotically
as an exponential, as expected. However $\rho(t)$ for the non-ergodic 
algorithm asymptotes to a logarithm. This implies that the autocorrelation
time is infinite, which is consistent for a non-ergodic algorithm,
although we do not know why $\rho(t)$ should have a perfectly logarithmic form.

WHY DOES THE AUTOCORR FUNCTION FOR JUST KBD FALL OFF LIKE A LOG? 


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

In order to guarantee that the lattice was broken up into at least 2 large
clusters at zero temperature, KBD used the freedom in choosing the 
probabilities so that only operations 2 and 3 had non-zero probability
at ${T \! = \!0}$. 
However in order to investigate the effect of tuning the extra free 
parameter $p_{sub}$, we used an alternate version of the KBD algorithm
that uses only the angle operations (4-7) at zero temperature, 
rather than the operations (2,3) used by KBD. 
We measured the cluster sizes and autocorrelations for this algorithm
for different values of $p_{sub}$: 1.0, 0.8, and 0.5. As noted by KBD, for 
this choice of probabilities the lattice is no longer split into two
large clusters, rather we get one very large cluster and a number of very
small clusters. 
Tuning this extra free parameter does have an effect
on the largest cluster size, although it is quite small:
the average maximum cluster size is 
0.91762(1) of the lattice for $p_{sub} = 0.5$,
0.89978(6)  for $p_{sub} = 0.8$, and
0.87371(10) for $p_{sub} = 1.0$. 
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 a value
$z = 1.97(2)$ which is very similar to the Metropolis algorithm,
since the largest cluster encompasses most of the lattice and consequently
the cluster update does little apart from a global spin flip.
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 substantially reduce the autocorrelations -- 
$p_{sub} = 0.8$ has $z = 0.62(2)$, while $p_{sub} = 1.0$ has $z = 0.48(1)$. 
It is very surprising that in going from $p_{sub} = 0.5$ to 1.0 the biggest 
cluster size is only reduced by about $5\%$, yet $z$ is reduced from 
around 2.0 to around 0.5.
This result is very promising, since it implies that critical slowing down
can be substantially reduced even when the biggest cluster size is quite 
large (87\% 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 generalized cluster 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 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 works well only because
it has a tunable parameter $p_{sub}$ that 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 would 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. Hence, 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 that 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
(triangular) plaquette, 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. 
Since the ground state configuration has a single unsatisfied link, and
the energy conserving operations (2-4 of Fig.~\ref{ops_tri}) freeze only a
single bond, there are 3 possible orientations of each. Thus, unlike the
square lattice case, there are no free parameters when we solve the matrix
equation~\ref{eq:matrix}, even when all possible orientations are used.
The probabilities for each operation are given in Table~\ref{tab:probs_tri}.

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

At the zero temperature critical point, we only perform the energy conserving
operations (2-4). This gives a bond/link ratio of $1 \over 3$, 
which is below the percolation threshold for the triangular lattice, 
although in this case the bonds are not placed randomly, 
but rather one per plaquette.
However we expect to at least be able to create multiple clusters,
as with the KBD algorithm with angled operations for the square lattice FFIM. 

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}$) that 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 for the triangular lattice 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
introducing a bias in the choice of the operations (2-4) does indeed reduce 
the largest cluster size,
but of course it also gives incorrect results due to the violation of 
detailed balance.

Since the ground state for this model is paramagnetic, we measured the
paramagnetic susceptibility ${\chi = V <m^2>}$, 
where $m$ is the average magnetization per site,
and $V$ is the lattice volume. At the zero temperature critical point, this
quantity is expected to remain of $O(1)$ as the lattice size is increased.
We measured the autocorrelations in $\chi$ and in the spin correlation
function $<\sigma_0 \sigma_R>$ at a distance $R=2$
for both the plaquette algorithm and the standard Metropolis algorithm.
As expected, we found that the
plaquette algorithm gave very little improvement over standard Metropolis,
as can be seen in Fig.~\ref{tri_corrs}, which shows the autocorrelation
function for the paramagnetic susceptibility at zero temperature 
on a $128^2$ lattice.
We have plotted the results in such a way that we compare two standard 
Metropolis sweeps with one iteration of the plaquette algorithm, which is a 
plaquette update followed by a Metropolis update to ensure ergodicity. 

To try to improve on this algorithm, we investigated making the tile
something larger than a plaquette. We first tried using the
the double plaquette shown in Fig.~\ref{dbl_tri}a. 
(DO WE REAALLY WANT TO LIST THE POSSIBLE OPERATIONS FOR THIS TOPOLOGY???
I THINK NOT).
%
% This has 5 possible link topologies, corresponding to 4 possible energy
% states, as shown in Fig.~\ref{dbl_conf}.
% There is a subtletly here, in that the ground state configuration of the
% whole lattice can induce 2 different energies for the double plaquette.
% IS THIS A PROBLEM???
% There are also 5 different operators, as shown in Fig.~\ref{dbl_ops}.
% This offers the possibility of reducing the bond/link ratio by sometimes 
% freezing only the link common to both the triangular plaquettes. 
% As with the single plaquette case, the number of energy conserving operations
% is the same as the number of possible ground state configurations, so there
% is again no free parameter.
%
There is a major problem with this approach, however, in that the 
double plaquette cannot tile the lattice -- that is, it is not possible 
to cover all the links in the lattice such that each link is uniquely 
assigned to a double triangular plaquette. Since the double plaquette
has only 5 links and the number of links in the lattice is a multiple of 3, 
there is an extra link required. Fig.~\ref{dbl_tri}b shows a set of links
which can tile a triangular lattice.

We might try to approach this problem by treating the double plaquettes and
the extra links separately, by applying the generalized cluster
algorithm to them as separate entities.
For the single links, this is equivalent to applying the SW algorithm, 
i.e. deleting the link if it is unsatisfied, 
and freezing it with probability $e^{-2K}$ if it is satisfied.
However this method ignores the interplay between the links and the double
plaquettes, and can lead to non-conservation of energy at zero temperature.
As shown in Fig.~\ref{bogus}, for some configurations it is possible to 
choose to freeze only the diagonal links for two double plaquettes
bordering a single unsatisfied link. This results in all
bonds being deleted in the triangular plaquette bordered by the link and 
the two double plaquettes, which allows the possibility of the update 
flipping all spins in this triangular plaquette independently, and thus
changing the energy. This is not a proof that the energy is not conserved,
since it is still possible that the geometry may force these spins to end
up in the same cluster.
However when we implemented this algorithm we found that there were
instances for which the energy was not conserved. (IS THIS TRUE???)

This failure can be rectified by just applying the generalized cluster 
algorithm to the basic object in Fig.~\ref{dbltri}b. The possible link
configurations and freeze/delete operations are shown in Fig.~\ref{dblink_ops}.
Here there are 8 configurations of the tile which can occur
in the ground state (a ground state configuration has two satisfied links
in each triangular plaquette). 
Unlike the algorithm using the triangular plaquette, these configurations
of the tile can have different energies.

An operation conserves energy if it freezes at least one satisfied bond
for each triangular plaquette.
For this object, there are ?? possible energy conserving operations which
freeze 3 or fewer bonds, of which only 8 satisfy the constraint described
in section~\ref{sec:algorithm}. This means we can construct an algorithm 
which is valid at zero temperature, but it has no free parameters, since the
number of ground state configurations is equal to the number of energy
conserving operations. 
By choosing the energy non-conserving operations to be those in 
Fig.~\ref{dblink_ops} and solving the matrix equation, 
we obtain the probabilities of Table~\ref{tab:dblink_ops}, which give
a valid algorithm for any temperature.

At zero temperature, this algorithm gives a reduced average maximum cluster
size of ?????(?) for the $64^2$ lattice.
The autocorrelations are ??????????

% *************** ?????????????????????????????? ***************************
%
% and different numbers of frustrated and satisfied links, for this object,
% depending on which links it covers. In thi particular case it is not 
% possible to find a freeze/delete operation that will conserve
% the energy for all possible ground state configurations.
% (except for the plaq + SW version??) Is this true? True in general?
%
% *************** ?????????????????????????????? ***************************

Thus although it is possible to construct a simple generalized cluster 
algorithm for the triangular lattice FFIM that 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}.


\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$) 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 links, we use
the KBD algorithm, while for plaquettes with an even number of such links,
we use the SW algorithm on each of the links (i.e. just freeze the satisfied
bonds at zero temperature).
This is necessary at the zero temperature critical point, since satisfied
bonds cannot be deleted without allowing non-conservation of energy for 
plaquettes with an even number of ferromagnetic links.
For the spin glass there is not the regular geometry that allows the KBD
algorithm to split the lattice into two clusters, and at zero temperature
we again have the problem that the lattice freezes 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 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. 

The three shapes shown in Fig.~\ref{honeytile} are the simplest non-trivial
objects which can tile the honeycomb lattice.
Using tile A gives an algorithm which is equivalent to SW (just as using the
general cluster formalism on a tile which is two links at right angles 
gives the SW algorithm for the square lattice FFIM).
Using tile B ???????????
Using tile C is really complicated???????????.
DO THIS??

% 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{CONCLUSIONS}
\label{sec:conclusion}

THIS NEEDS IMPROVEMENT.

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. We have seen that we can construct
a cluster algorithm for a spin model based on a larger group of spins
as long as it can tile of the lattice. 
However for a particular spin model,
there is no guarantee that a valid algorithm with sensible probabilities 
for each possible freeze/delete operation can be constructed. Even
if such an algorithm does exist, the resulting clusters may be too big, 
as is the case for the triangular lattice FFIM and the Ising spin glass,
so that the cluster algorithm performs no better than traditional 
algorithms such as Metropolis. 
We have also found that this problem also occurs with a generalized cluster
algorithm for the 3-d FFIM~\cite{cubic}.
Unfortunately it appears that although the 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.

Here we have looked only at the case where the bonds are either frozen or
deleted. This produces clusters that can be updated independently of
one another. However the generalized cluster algorithm allows for another
possibility, in which there is an interaction energy between the clusters.
This could allow the creation of smaller clusters that may be updated with
a probability dependent on the energies of the other clusters, in a fashion
similar to the algorithms of Niedermayer~\cite{niedermayer} and 
D'Onorio de Meo {\it et al.}~\cite{demeo}.
However these kinds of interacting cluster
algorithms have not proven to be effective for non-frustrated systems,
mainly because the probability of updating a cluster is approximately 
inversely proportional to the exponential of the cluster size. 
Thus updating large clusters occurs with a very small probability, so 
these algorithms tend to be not much better than standard local Monte Carlo 
methods, at a much greater computational cost.

\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{Wfrust}
% 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{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{kandel90}
D.~Kandel, R.~Ben-Av and E.~Domany, Phys. Rev. Lett. {\bf 65}, 941 (1990).
\bibitem{kandel91}
D.~Kandel and E.~Domany, Phys. Rev. A {\bf 43}, 8539 (1991).
\bibitem{kandel92}
D.~Kandel, R.~Ben-Av and E.~Domany, Phys Rev B {\bf 45}, 4700 (1992).
\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{villain}
J.~Villain, J. Phys. C {\bf 10}, 1717 (1977).
\bibitem{mathematica}
S.~Wolfram, {\it Mathematica: a system for doing mathematics by computer},
(Addison-Wesley, Redwood City, Calif., 1990).
\bibitem{paul}
C.~F.~Baillie and P.~D.~Coddington, Phys. Rev. B {\bf 43}, 10617 (1991).
\bibitem{timeseries}
M.B.~Priestley, {\it Spectral Analysis and Time Series}, 
(Academic Press, London, 1981).
\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{cubic}
P.D.~Coddington and L.~Han, ``A generalized cluster algorithm for the 3-d
fully frustrated Ising model'', in preparation.
\bibitem{niedermayer}
F.~Niedermayer, Phys. Rev. Lett. {\bf 61}, 2026 (1988).
\bibitem{demeo}
M. D'Onorio de Meo, D.W. Heermann and K. Binder,
J. Stat. Phys. {\bf 60}, 585 (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_{0}$ for the standard KBD algorithm for 
the square lattice FFIM at zero temperature. 
The dashed line shows the function $\log(\log L)$.
\label{tau_0}}
%
\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}}
%
\figure{Autocorrelation function for the KBD algorithm with and without
Metropolis sweeps for the square lattice FFIM at zero temperature. 
\label{tau_ergodic}}
% 
%
% ********* 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{dbltri}}
%
\figure{The triangular lattice cannot be covered using only a double plaquette.
\label{tricover}}
%
%
% TABLES --- put these in the appropriate place in the text.
%
\begin{table}
\caption{Matrix containing the number of possible orientations of the
freeze/delete operators of Fig.~\ref{ops_sq}
for each energy state of the the square lattice FFIM.
}
\label{tab:matrix_sq}
\end{table}
%
\begin{table}
\caption{Probabilities for the freeze/delete operations of Fig.~\ref{ops_sq}
for the square lattice FFIM.
}
\label{tab:probs_sq}
\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).
}
\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.
}
\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.
}
\label{tab:probs_tri}
\end{table}


\end{document}


