%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 and L. Han}
\begin{instit}
Physics Department, Syracuse University, Syracuse, NY 13244, USA
\end{instit}
\begin{center}
\bigskip
1 May 1992
\bigskip
\end{center}
% The author will not know the received date when the compuscript
% is first submitted; production will insert this. 
\begin{abstract}
\end{abstract}
% \pacs{PACS numbers 05.50.+q, 75.40.Mg}

\newpage
\narrowtext

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

The cluster update algorithm of Swendsen and Wang~\cite{SW} (SW) provides
a great improvement in the efficiency of Monte Carlo simulations over the
Metropolis algorithm~\cite{metrop} and other local update schemes
by dramatically reducing critical slowing down~\cite{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 is 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}. 
% A bit more waffle here.

In section~\ref{sec:algorithm} we introduce the KBD general cluster algorithm,
and present the method we used for implementing the algorithm for a particular 
spin model. We then apply this to the fully frustrated Ising model on a 
triangular lattice (section~\ref{sec:triangle}) and a cubic lattice 
(section~\ref{sec:FFI3d}).


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

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.

FFIM
link (interaction between neighboring sites)


\subsection{The Metropolis Algorithm}
\label{subsec:met2d}

Ergodic if you sweep over sites chosen at random????
Not ergodic if you loop in order over sites. But in that case, it
gives autocorrs which are smaller than the random algorithm, 
and even gives z < 1???


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

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}.

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 (64x64).

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.


In Ref.~\cite{kandel91b}, KBD give evidence that their algorithm, with a 
Metropolis sweep after every cluster iteration, is ergodic at zero temperature.
It is easy to see that the KBD algorithm without Metropolis sweeps is
not ergodic at zero temperature. It is illuminating to study the algorithm
both with and without Metropolis sweeps, to get some idea of 

Consider any elementary plaquette of a ground state configuration.
Suppose the unsatisfied link for this plaquette is in the vertical 
direction. After a KBD update (using operations 2 or 3 in Fig.~\ref{ops_sq}),
the unsatisfied link must remain in the vertical direction,
so configurations with an unsatisfied link on this plaquette in the horizontal
direction can never be reached (this would require using operations 4-7 in 
Fig.~\ref{ops_sq}), and thus the algorithm is highly non-ergodic.

Because the model does not have rotational symmetry (the antiferromagnetic
couplings are only in one direction), this non-ergodicity can have 
measureable consequences. 
For example, consider a $2 \times 2$ lattice (i.e. a single plaquette). 
It has 8 possible ground states, 4 of which have unsatisfied links in the
direction of the antiferromagnetic bonds, 
% (Fig.~\ref{gs2}(a)), 
and 4 of which have unsatisfied links in the other direction.
% (Fig.~\ref{gs2}(b)). 
The average values of $|M|$ and $M^2$ are both 0.5 for the former set of
configurations,
% for the configurations in Fig.~\ref{gs2}(a), 
but they are 0.5 and 0.25 respectively for the latter set.
% configurations shown in Fig.~\ref{gs2}(b). 
An ergodic algorithm would therefore find the average $M^2$ to be 0.375
for this lattice size,
whereas the KBD algorithm would find 0.25 or 0.5, depending on its starting
configuration, since it is trapped in one or the other set of configurations.

This problem occurs for arbitrarily large lattice sizes. For example, if we
take as an initial configuration the completely ordered state (which is a
ground state of this model), then a KBD update of this configuration
will produce clusters which are horizontal strips (i.e. each cluster 
consists of sites $(i,j)$ with constant $j$). Flipping the spins
in some of these clusters will produce configurations for which the spins
are the same along a horizontal strip. Further updates will only produce 
configurations of the same type, so the algorithm is trapped in a very
small subspace of the possible ground states.

In Table~\ref{tab:ergodic} we give some measured averages of $M^2$
for the KBD algorithm at $T = 0$ with and without 
Metropolis sweeps, for various lattice sizes and initial configurations.
Each of the runs was at least 500,000 sweeps.
The non-ergodicity is clearly seen in the results for smaller lattice 
sizes. However it is interesting to note that the effect decreases as $L$ 
is increased, since the number of ground state configurations and the range
of values for the magnetization of each configuration become much greater.
For lattices sizes greater than about $100^2$ it becomes difficult to 
see any discrepancy between the results from the non-ergodic and ergodic 
algorithms.

We checked that the KBD algorithm with Metropolis sweeps gave results which 
were independent of the starting configuration (within errors) for all lattice
sizes, and also agreed with the exact values which can be calculated
for the $2^2$ and $4^2$ lattices.

Adding the angled operations 4-7 may make the KBD algorithm ergodic,
however the problem then is that the lattice is not guaranteed to split
into at least two large clusters, so the algorithm is much less
effective at reducing critical slowing down.

The KBD algorithm works extremely well for the square lattice FFI, almost
completely eliminating critical slowing down. However it is not at all 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 with 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 general cluster algorithms will work for other
frustrated spin models.


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


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


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

We next attempted to apply the generalized cluster Monte Carlo idea to
the three dimensional FFIM~\cite{BMB,grest,diep}.

Any ground state config is frozen for any local update algorithm at zero 
temperature 



\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.
Even if such an algorithm does exist, the resulting clusters may be too big, 
as in both the cases which we have investigated.
Unfortunately it seems that although the KBD general 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 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{metrop}
N. Metropolis {\it et al.}, J. Chem. Phys. {\bf 21}, 1087 (1953).
\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{villain}
J. Villain, J. Phys. C {\bf 10}, 1717 (1977).
\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{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 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{glass}
reference to 2-d spin glass.
\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{perc}
D. Stauffer,
%
% \bibitem{ising}
% See for example B.~M.~ McCoy and T.~T.~Wu, {\em The Two-Dimensional Ising
% Model,} (Harvard University Press, Cambridge, Mass., 1973).
% \bibitem{binder}
% {\it Monte Carlo Methods in Statistical Physics, Topics in Current 
% Physics 7, Second Edition}, Ed. K. Binder (Springer-Verlag, Berlin, 1986);
% {\it Applications of the Monte Carlo Method in Statistical Physics, Topics in 
% Current Physics 36, Second Edition}, Ed. K. Binder (Springer-Verlag, Berlin,
% 1987).
%
\end{references}
%
%
% FIGURES
%
% ********* section FFI2d **************
%
\figure{Freeze/delete operations for square lattice FFIM.
\label{ops_sq}}
%
\figure{Possible ground states of the FFIM on a $2 \times 2$ lattice. 
The dashed lines indicate unsatisfied links.
\label{gs2}}
%
\figure{Autocorrelation function for the KBD algorithm with and without 
Metropolis sweeps, for the square lattice FFIM at zero temperature. 
\label{rhoM0M1}}
% 
%
% ********* section triangle **************
%
\figure{Freeze/delete operations for triangular lattice FFIM.
\label{ops_tri}}
%
%
% ********* section FFI3d **************
%
\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
%
\begin{table}
\caption{Average of the $M^2$ for different lattice sizes $L$ using the 
KBD algorithm with no Metropolis sweeps,
for four different initial zero temperature configurations C0 (the
ordered state), and C1, C2, C3 (thermalized from 3 different random starts). 
Also shown are results for the KBD algorithm with one Metropolis sweep 
per iteration (KBD+M).
The averages were taken over at least 50,000 iterations, and the 
statistical error in the last digit is given in brackets. 
For the two smallest lattice sizes we also computed the exact result.
\bigskip}
Measured averages of $M^2$
\label{tab:ergodic}
\end{table}
%

\end{document}


