
\section{Computer Simulation}

\begin{slide}

\subsection{\bf The Metropolis Update}

Note that the Metropolis algorithm does not specify {\em how} the changes 
to the configuration should be made --- it just says that any {\em proposed}
change to the system should be {\em accepted} with a certain probability
that depends on the change in energy.

How the changes are made depends on the variables and the model being 
studied. The only constraints on the update procedure are:
%
\begin{itemize}
\item  it should be ergodic 
\item  it should not be biased in such a way as to violate detailed balance
\end{itemize}
%
Another issue is efficiency --- the procedure should sample the 
configuration space as effectively as possible. There is often some 
freedom in tuning the algorithm to improve efficiency and performance.

\end{slide}


\begin{slide}


\subsection{\bf Metropolis Algorithm for the Ising Model}

For the Ising model, the obvious change in the configuration is to try to 
update (or ``flip'') a spin, i.e. flip the sign (or direction) of the spin
variable. If we try to change many spins at once, $\Delta E$ will be 
large, so the probability to accept the change, $e^{-\Delta E/kT}$, will be 
small. So update a single spin at a time.

For a single spin flip, $\Delta E$ depends only on the spin values at the 
site and its nearest neighbors, i.e. the update is {\em local}.

% Figure ??

If $\Delta E \le 0$, we make the proposed change.

If $\Delta E > 0$, we make the change with probability $e^{-\Delta E/kT}$.

\end{slide}


\begin{slide}

\subsection{\bf Programming the Metropolis Algorithm}

Update a single site at a time.  Choose sites at random to ensure no bias
(in most cases, this is not really necessary, and a simple contiguous loop 
over lattice sites is sufficient).

A loop over all sites is referred to as one Metropolis {\em sweep} 
(or iteration).

The probabilistic part of the algorithm is done using a random number 
generator (hence the name -- Monte Carlo -- in reference to games of
chance). Generate a random floating point number $r$ in [0,1). 
A process with probability $p$ is performed if $r < p$.

The Metropolis algorithm is very simple to program.

\end{slide}


\begin{slide}

\subsection{\bf Code for Metropolis Algorithm (in C)}

substitute this for figure in portrait mode

\end{slide}


% Ising simulation configuration picture ??


\begin{slide}

\subsection{\bf Measurements}

We generate configurations of the Ising model using a Monte Carlo algorithm, 
such as Metropolis.  Then what do we do with them?  
We want to numerically measure quantities of interest, such as:
 
\begin{tabular}{ll}
\smallskip
energy & $E = - \frac{J}{V}  {\displaystyle 
\sum_{\langle i,j\rangle} } S_i S_j $ \qquad  $V$ = volume (\# sites) \\
\smallskip
magnetization & $M = \frac{1}{V} {\displaystyle \sum_i} S_i$ \\
\smallskip
correlation function & $\Gamma (n) = \frac{1}{V} {\displaystyle \sum_i }
S_i S_{i+n} - \left( \frac{1}{V} {\displaystyle \sum_i } S_i \right)^2$ \\
\medskip
correlation length & $\Gamma (n) \sim e^{-n/\xi}$ \\
\medskip
specific heat & $C/V =          % check factors of V and T ??
\frac{1}{V} \frac{\partial E}{\partial T} = 
\left\langle \left( E - \langle E\rangle \right)^2 \right\rangle = 
\left\langle E^2\right\rangle - \langle E\rangle^2$ \\
\medskip
susceptibility & $\chi/V =      % check factors of V and T ??
\frac{1}{V} \frac{\partial M}{\partial T} =  
\left\langle \left( M - \langle M\rangle \right)^2 \right\rangle = 
\left\langle M^2\right\rangle - \langle M\rangle^2$ \\
\end{tabular}
 
Quantities such as these are measured for each configuration, and the
averages and statistical errors calculated.

\end{slide}


\begin{slide}

\subsection{\bf Phase Transitions}

Phase transitions are described by an {\em order parameter} which 
differentiates the two phases. Usually it is zero in the disordered 
phase, non-zero in the ordered phase, e.g. the magnetization for a 
magnetic system.

Phase transitions are classified in part by their {\em order}\/:

\begin{itemize}
\item First order --- discontinuity in the order parameter or energy, 
i.e., first derivative of the partition function.

\item Second order --- divergence in the susceptibility or specific heat, 
i.e., second derivative of the partition function.

\item Third order, etc.
\end{itemize}

First order transitions are classified by quantities such as the 
{\em latent heat}, the discontinuous change in the energy at the 
critical point $T_c$.

\end{slide}


\begin{slide}

\subsection{\bf Critical Exponents}

Second order transitions are classified by their {\em critical exponents}, 
which measure how quantities diverge at the critical point. \\ 
N.B.\ These are asymptotic $(T\rightarrow T_c)$ results.
% 
\begin{eqnarray*}
M & \sim & \vert T-T_c \vert^\beta 
\hspace{2em} \mbox{\Smallsize magnetization}
\\
\chi & \sim & \vert T-T_c \vert^\gamma 
\hspace{2em} \mbox{\Smallsize susceptibility} 
\\
C & \sim & \vert T-T_c \vert^\alpha 
\hspace{2em} \mbox{\Smallsize specific heat} 
\\
\xi & \sim & \vert T-T_c \vert^\nu 
\hspace{2em} \mbox{\Smallsize correlation length} 
\end{eqnarray*}
 
In many cases, these exponents ($\alpha, \beta, \gamma, \nu$) 
are {\em universal}, i.e., they do 
not depend on details of the model, but only on gross features such as the 
dimension of the space and the symmetries of the energy function.  
This explains the success of very simple spin models like the Ising model in 
providing a quantitative description of real magnets with complex interactions.

\end{slide}


\section{Problems, Subtleties, and Tricks of the Trade}


\begin{slide}

\subsection{Problems and Subtleties in Monte Carlo Simulations}

Although Monte Carlo simulation (using the Metropolis algorithm, for example)
appears very straightforward, there are in fact many problems and 
subtleties that can trap the unwary and produce unreliable results.

The following are some examples of common problems that arise in Monte Carlo
simulations, along with examples of techniques to attempt to overcome them.

\end{slide}


\begin{slide}

\subsection{\bf Thermalization}

Before taking measurements, we must be sure that the Markov process has 
{\em thermalized}, or reached {\em equilibrium}, i.e. $P(A,t) = p(A)$ 
and the arbitrary initial configuration has been ``forgotten.''

Approach to equilibrium is exponential $\sim e^{t/\tau}$.
The {\em autocorrelation time} $\tau$ can be very large. 
May thermalize for 100,000 sweeps at a certain $T$, 
but if $\tau$ is of the same order, then the measurements are suspect.

There are 3 main techniques for checking for thermalization:

{\bf Binning:} 
Check to see whether measurements are converging to a constant 
average value.  Can do this by {\em binning}\@ the data, i.e.
splitting measurements up into a number of large contiguous 
segments ({\em bins}).  Averages are taken over each bin, 
and data from initial deviant bins is discarded. 
Binning can also be used to calculate errors (see later).

\end{slide}


\begin{slide}

\subsection{\bf Thermalization (cont.)}

{\bf Autocorrelations:}  Measure the autocorrelation time
$\tau$ over the whole simulation, and then throw away data 
for at least the first $10 \tau$ sweeps.

{\bf Hot and cold starts:} 
Start from hot (random) configuration and also a cold (ordered) configuration.
When the two have converged to the same results, then the system 
is equilibrated. 
(But what about first order phase transitions at $T_c$? 
Usually will not get convergence ---
in fact this is a good way to look for first order transitions.)

% Figure ??

\end{slide}


\begin{slide}

\subsection{\bf Ergodicity}

Can be very difficult to prove that the method of choosing new
configurations is ergodic (especially at $T \! = \! 0$). 
Can get trapped in subsets of configuration space, such as periodic loops.

May also have updates that are {\em quasi-ergodic}\@ --- there may be a 
high energy barrier between subspaces, 
so it is theoretically possible to reach all possible states, 
but \Underline{highly improbable} in a finite time simulation.

{\bf Check:} Compare results from multiple runs with different initial 
configurations.

\end{slide}


\begin{slide}

\subsection{\bf Fitting}

Once measurements are taken, generally need to fit the data to find 
parameters such as critical exponents, using 
a multi-parameter fit to data which usually has large errors.  
This is very difficult! 

For example
$ C = A + B \left\vert T-T_c \right\vert^{-\alpha} $
requires a four parameter fit, and often have so much freedom that a wide 
range of values are possible.  Usually need to assume $A\approx 0$, and 
get $T_c$ from other measurements.

Above formula is only true ``asymptotically close'' to $T_c$. 
How close is that?  Not always clear.  
May see ``apparent'' or {\em pseudo-critical} exponents, 
which change or disappear as the size of the system is increased and 
we move closer to $T_c$.

There may also be $\log\left\vert T-T_c \right\vert$ corrections 
to the asymptotic formula, which are hard to distinguish from a small 
power-law exponent. In this case may need input from theory.

\end{slide}


\begin{slide}

\subsection{\bf Finite Size Effects}

When studying any macroscopic system with a very large number of degrees of 
freedom, invariably make an approximation and simulate a smaller and/or 
discretized model system.  This introduces systematic errors called 
{\it finite size effects}.  
Have to understand these, and be able to extrapolate to an infinite system, 
usually by doing a number of simulations at different system sizes.

For spin models, we have a finite $d$-dimensional lattice of $V=L^d$ sites. 
But only get a true phase transition (i.e., divergence) when 
$L\rightarrow\infty$.  
For a finite system, get rounded peaks rather than divergences. 
The peaks narrow and increase in height as $L$ is increased, and the
location of the peak shifts slightly. 

Many problems require an empirical extrapolation to an infinite system.
But for phase transitions in statistical mechanics, some elegant and
useful theoretical results exist.

\end{slide}



\begin{slide}

\subsection{\bf Illustration of Finite Size Effects}

Substitute for scanned figure ??

\end{slide}



\begin{slide}

\subsection{\bf Finite Size Scaling}

In some instances, e.g. second order phase transitions, the form 
of the asymptotic (large $L$) behavior is known for finite size systems.
%
\begin{eqnarray*}
C_{\max} (L) & \sim & L^{\alpha/\nu} 
\hspace{2.0em} \phantom{^-}
\mbox{\Smallsize $C_{\max}(L)$ is height of specific heat peak}
\\
\left\vert T_c (L) - T_c\right\vert & \sim & L^{1/\nu} 
\hspace{2.0em} \phantom{^-}
\mbox{\Smallsize $T_c (L)$ is where $C_{\max} (L)$ occurs, 
$T_c = T_c (\infty)$} 
\\
M\left( T_c \right) & \sim & L^{-\beta/\nu} 
\hspace{2.0em}
\mbox{\Smallsize Magnetization at $T_c$ is zero only for $L \rightarrow \infty$}
\\
\chi\left( T_c \right) & \sim & L^{\alpha/\nu} 
\hspace{2.0em} \phantom{^-}
\mbox{\Smallsize Susceptibility diverges only for $L \rightarrow \infty$}
\end{eqnarray*}
 
Example of simple (non-rigorous) ``proof'' (for specific heat):
%
\begin{eqnarray*}
\xi_{\max} (L) & \sim & L \hspace{2em} 
\mbox{\Smallsize correlation length cannot exceed length of lattice!} 
\\
\xi_{\max} (L) & \sim & \left\vert T_c (L) - T_c \right\vert^{-\nu} 
\hspace{2em} \mbox{\Smallsize from definition of $\nu$} 
\\
\mbox{Therefore} \qquad &  & 
\\
L^{1/\nu} & \sim & \left\vert T_c (L) - T_c \right\vert  
\\
C_{\max} (L) & \sim & \left\vert T_c (L) - T_c \right\vert^{-\alpha} 
\\
& \sim & L^{\alpha/\nu} 
\end{eqnarray*}

\end{slide}


\begin{slide}

\subsection{\bf Finite Size Scaling (cont.)}

These relations are extremely useful, and even allow us to extract exponents 
and $T_c (\infty )$.  But only if we are at ``large enough'' $L$.  
If $\xi$ is very large (but not infinite), may get incorrect
``pseudo-critical'' exponents.

This is often the case at a ``weakly'' first order transition 
(very small latent heat).  
So it is often very difficult to distinguish between a first and 
second order phase transition, or even to say for sure if there is a 
phase transition at all --- is $\xi = \infty$, or just $\xi \gg L$?
May need extremely large system sizes to get correct results.

In general, don't know the form of the finite size scaling, so try to fit 
data to a {\em scaling function}, or go to large enough systems so 
that results $\rightarrow$ constant or can do a simple linear extrapolation.

\end{slide}


\begin{slide}

\subsection{\bf Boundary Conditions}

Generally use {\em periodic boundary conditions}, where the lattice wraps 
around on itself to form a torus.  This has been shown to give the 
smallest finite size effects.

But then must be careful that in measuring the correlation function, 
can only take measurements to $L/2$ for a lattice of size $L$, since
a distance $n$ is equivalent to a distance $n + L/2$ on a lattice that
wraps around on itself.

% Figure ??

\end{slide}


\begin{slide}

\subsection{\bf Tuning the Acceptance}

For the Ising model, there is only one possible change to the state\@: $S_i 
\rightarrow -S_i$.  More generally (e.g., continuous spin models) we are free 
to make the change in the variable as large or small as we like.

A small change means $\Delta E$ is small and so the {\em acceptance} 
(proportion of proposed changes that are accepted) is close to 1.  
But then we move slowly through configuration space.  
Large changes usually have large $\Delta E$ and thus low acceptance,
which again means that configuration space is not sampled efficiently. 

General practice is to do a tradeoff and tune the size of the change 
so that the acceptance $\approx 0.5$.

\end{slide}


\begin{slide}

\subsection{\bf Random Number Generators}

Generating ``good'' random numbers can be a major problem in 
Monte Carlo simulations.  
Many basic random number generators (RNGs) supplied with computers (especially 
parallel computers) are not sufficient.  Particularly true in the past, but 
unfortunately still true far too often today!

There are many statistical tests for uniform distribution, correlations, etc. 
Monte Carlo measurement of the energy of the 2-$d$ Ising model is a good 
test --- the result is known exactly for $L^2$ lattice.  
Some RNGs pass statistical tests, but fail the Ising model test!

Even some ``good'' RNGs having adequate randomness properties are not 
good enough for Monte Carlo simulation because their period of repetition 
is too small.

\end{slide}


\begin{slide}

\subsection{\bf Random Number Generators (cont.)}

Multiplicative linear congruential generators (MLCG)
using 32-bit integers have a period of at most 
$2^{31} \approx 10^{10}$. This many random numbers
can be generated in seconds on a modern workstation.

Period must be much greater than the number of random numbers used in 
the simulation, or else the results can be incorrect.
Modest Ising model simulation, $1024^2$ lattice, $10^7$ sweeps, 
uses $\sim 10^{13}$ random numbers.

Using 64-bit words, or combining two 32-bit MLCGs, can give periods $\sim 
10^{18} \sim$ \# nanoseconds per year. These are adequate for largest
simulations done today (Gigaflop-year).
Other generators (lagged Fibonacci, etc.) have even longer periods.
These or combined 64-bit MLCGs will be required for Teraflop-year 
simulations.

Random number generators will be the subject of a later lecture.

\end{slide}


\begin{slide}

\subsection{\bf Critical Slowing Down}

A problem with Monte Carlo (MC) simulation is that the error decreases only 
as $1/\sqrt{N}$, so 
need to do 100 times more work to get one more significant figure accuracy.
Unfortunately, things are even worse.  Error is only $1/\sqrt{N}$ for $N$ 
{\em independent} measurements.  But MC configurations are usually 
correlated --- $C_{i+1}$ depends on $C_i$.  

Measure the {\em autocorrelation function} 
(similar to correlation function except in computer time, not real space) 
for an operator (measurement) $A$:
%
\begin{eqnarray*}
\rho(t) & = & \frac{\langle A(t) A(0)\rangle - \langle A(0)\rangle^2}{\langle 
A(0)^2 \rangle - \langle A(0)\rangle^2} \\
& \sim & e^{-t/\tau_{\exp}} 
\end{eqnarray*}
%
$\tau_{\exp}$ is the (exponential) {\em autocorrelation time}, 
and is equivalent to 
the $\tau$ which measures equilibration time of the Markov process.  
Need to thermalize for time $t \gg \tau_{\exp}$.

\end{slide}


\begin{slide}

\subsection{\bf Critical Slowing Down (cont.)}

The {\em integrated autocorrelation time}

\begin{displaymath}
\tau_{\rm int} = \frac{1}{2} + \sum^\infty_{t=1} \rho (t)
\end{displaymath}

is a measure of how correlations affect statistical errors. 
It is normalized so that $\tau_{\rm int} 
\approx \tau_{\rm exp}$ if $\rho\approx 
e^{-t/\tau_{\exp}}$. It can be shown that for MC measurements
%
\begin{eqnarray*}
\mbox{error in mean} & = & \sqrt{\frac{2\tau_{\rm int}}{N}} \cdot \sigma \\
\sigma & = & \sqrt{\rm variance} 
\end{eqnarray*}

A major problem for MC simulations is that 
$\tau_{\rm int}$ \Underline{diverges} at $T_c$!!  
This is known as {\em critical slowing down}.

\end{slide}


\begin{slide}

\subsection{\bf Dynamic Critical Exponents}

To decorrelate a configuration, need to make changes on a scale $>$ 
spatial correlation length.  For \Underline{local} MC 
algorithms (such as Metropolis) the effect of a change is only seen by nearest 
neighbors. Since the update is a stochastic process, the changes should 
propagate through the system like a random walk, so the number 
of iterations required to propagate a distance $\xi$ goes like $\xi^2$. 
At the critical point, $\xi \sim L$, so we expect $\tau \sim L^2$.  
This is a major problem for MC simulations, since we need large $L$.

More generally, $\tau (L) \sim L^z$ where $z$ is called the {\em dynamic 
critical exponent}.
For a local algorithm $z \geq 1$ and usually $z \approx 2$.
One of the active areas of research in MC simulation is 
trying to find new algorithms with $z<2$.  
These are usually non-local, multi-scale algorithms, 
c.f. multigrid methods for solving PDEs, 
where the same type of problem occurs.

\end{slide}


\section{Other Monte Carlo Algorithms}


\begin{slide}

\subsection{\bf Acceleration Algorithms}

A number of different methods have been used to attempt to ``accelerate'' the 
dynamics of Monte Carlo updates (i.e., reduce $z$).
\begin{itemize}
\item Over-relaxation
\item Multigrid
\item Fourier acceleration
\item Cluster algorithms
\end{itemize}

The first two are commonly used to improve numerical solutions of PDEs.
In some cases, over-relaxation can reduce $z$ to $\approx 1$.  
Multigrid Monte Carlo and Fourier acceleration work well in certain cases.
Cluster algorithms work very well, but currently are only 
applicable to a small class of spin models.

\end{slide}


\begin{slide}

\subsection{\bf Cluster Algorithms}

The idea behind cluster MC algorithms is very simple.  Local MC algorithms 
perform poorly because they update only one spin at a time.  Cluster 
algorithms use a clever way of finding large clusters of sites which can be 
updated at once.  (Remember that trying to update a group of sites chosen at 
\Underline{random} gives a large $\Delta E$ and thus a small 
acceptance).

Cluster algorithms were invented by Swendsen and Wang in 1987, 
for Ising and Potts spin models. 
They have since been generalized to other models, and other algorithms.

\end{slide}


\begin{slide}

\subsection{\bf The Potts Model}

The Potts model is an obvious generalization of the Ising model to an 
arbitrary number, $Q$, of discrete states.  Proposed by Domb to his grad 
student, Potts.
%
\begin{eqnarray*}
E & = & -J \sum_{\langle i,j\rangle} \delta \left( s_i, s_j \right) \\
 & & \\
\delta & = & 1 \hspace{1em} \mbox{if}\ s_i = s_j \quad \Rightarrow E = -J \\
& = & 0 \hspace{1em} \mbox{if}\ s_i \not= s_j \quad \Rightarrow E = 0 
\end{eqnarray*}

As for Ising model, if $J > 0$, the energy is lower if neighboring spins 
are the same (aligned).

2-state Potts $\Leftrightarrow$ Ising.

$Q>2$ model has different critical exponents, interesting phase structure.  
$Q \ge 5$ has first order transition in two dimensions.

\end{slide}


% ADD picture??


\begin{slide}

\subsection{\bf Percolation}

Swendsen and Wang derived their algorithm from a simple relation which 
transforms a Potts spin model into a bond percolation model.
What is percolation?

\begin{minipage}[h]{4.5in}
Suppose we have a two-dimensional lattice of sites, and we put a bond between 
every pair of neighboring sites with probability $p$.  Sites are then bonded 
together into connected clusters.  
\end{minipage} \ \
%
\begin{minipage}[h]{3.5in}
\begin{center}
\begin{picture}(151,151)
%%BoundingBox: 0 0 151 151
\psfig{figure=percolation.eps}
\end{picture}
\end{center}
\end{minipage} 

\vspace{1em}

At a critical value of $p$, will get {\em percolation}, 
i.e., one cluster will span the lattice, meaning that there is a connected
path across the lattice through a single {\em spanning cluster}.  
For an infinite two-dimensional lattice, $p_c = \frac{1}{2}$.

\end{slide}


% ADD picture??


\begin{slide}

\subsection{\bf The Swendsen-Wang Algorithm}

Swendsen-Wang (SW) found that by placing bonds between neighboring sites 
with the same spin value, the Potts model can be written in terms of 
interacting {\em clusters} of connected spins.  

For the particular choice $p = 1 - e^{-\beta}$
($\beta=\frac{J}{kT}$)
the interaction energy is zero, so the clusters are independent.  This means 
we can update all the clusters (by changing the current spin values to a 
random new spin value) completely independently.

Since the clusters can be very large, the update is highly non-local, 
and $z$ can be very small (even zero).

\end{slide}


\begin{slide}

\subsection{\bf The Wolff Algorithm}

Wolff introduced a variant of the SW algorithm, in which a site is chosen at 
random, and a single cluster is grown around the site.  This tends to favor 
larger clusters, and thus smaller autocorrelation times.  Wolff also 
generalized the cluster algorithms to continuous spin models.

Cluster algorithms work amazingly well for some models, e.g., 2-$d$ 
Ising model, where $z \approx 0$.  
Even for a large $512^2$ lattice, $\tau\approx 10$, which 
is 0(1000) times smaller than $\tau$ for the Metropolis algorithm.

\end{slide}


% sw.ps goes here ??

% wolff.ps goes here ??


\begin{slide}
\subsection{\bf Correlations}

Correlations fig goes here. % ??

\end{slide}



\begin{slide}

\subsection{\bf Connected Component Labeling}

The cluster algorithms are very different to Metropolis.  
The main computational 
task is to take local information (bond connections) and work out global 
information (clusters of sites).  

This is an example of {\em connected component labeling}, which is
a standard graph problem ---
given a graph of vertices and edges, the problem is to identify those
vertices which form a connected set, where a vertex is in the set if it
has an edge connecting it to another vertex in the set.
%
For cluster algorithms, the vertices are the lattice sites, the 
edges are the bonds, and the connected components are the clusters.

Component labeling is commonly used in image processing, to join neighboring 
pixels into connected regions which are the shapes in the image.  
Efficient sequential algorithms exist for this procedure, however
implementing an efficient parallel algorithm is difficult.

\end{slide}


