\documentstyle[med_headings]{article}
\oddsidemargin=0in\textwidth=6.4in\topmargin=-0.80in\textheight=10in
\renewcommand{\textfraction}{0.01}
\pagestyle{empty}\begin{document}\title{}\author{David J.C. MacKay }\date{}
\noindent
\input{psfig.tex}
\input{/home/mackay/tex/newcommands1}
\begin{center}
{\large \bf
Computational Methods in Physics
\\ \medskip
}
\end{center}
%
\newcommand{\bb}{{\bf b}}
\section*{Ising Models---Monte Carlo simulation}
An `Ising model' is an array of spins (\eg, atoms that can take
states $\pm 1$) that are
magnetically coupled to each other. In a ferromagnetic model it
is energetically favourable for neighbours to
be in the same state; in antiferromagnets, it is favourable for
them to be in opposite states.
Let the state $\bx$ of an Ising model with $N$ spins be a vector in which
each component $x_n$ takes values $-1$ or $+1$. If two spins $m$ and
$n$ are neighbours we write $(m,n) \in {\cal N}$. The energy of a state
$\bx$ is
\beq
E(\bx;J,H) = - \left[
% \frac{1}{2}
\sum_{(m,n) \in {\cal N}} J x_m x_n
+ \sum_{n} H x_n \right] ,
\eeq
where $J$ is the coupling between spins,
% $m$ and $n$,
and $H$ is
the applied field. If $J > 0$ then the
model is ferromagnetic, and if $J < 0$ it is antiferromagnetic.
(Note that each pair $(m,n)$ is counted just once in the first sum).
% In Physics we may be interested in the properties of Ising models with
% a large number $N$ of spins having regular geometric neighbourhood
% relationships.
In this problem we will study two-dimensional planar Ising models,
using a simple `Gibbs sampling' or `heat bath Monte Carlo' method.
Starting from some initial state, a spin $n$ is selected at
random, and the probability that it should be $+1$ given the state of
the other spins and the temperature is computed,
\beq
P(+1|h_n)= 1/(1+\exp(- 2 \beta b_n)),
\eeq
where $\beta = 1/k_{\rm B}T$ and $b_n$ is the local field
\beq
b_n = \sum_{m:(m,n) \in {\cal N}} J x_m + H.
\eeq
Spin $n$ is set to $+1$ with that probability, and then a new spin
is selected at random.
After sufficiently many iterations, this procedure
converges to the equilibrium distribution,
$P(\bx)=\frac{1}{Z}\exp(-\beta E(\bx;J,H))$.
\medskip
% \subsection*{Part a}
(A) Write a program that simulates an Ising model with the square
geometry (a) shown below, and with periodic boundary conditions. A
line between two spins indicates that they are neighbours.
%
% To make a bite-sized example, we will set $b$ to 0 throughout,
Set the external field $H=0$
and consider the two cases $J = \pm 1$.
% which are a ferromagnet and antiferromagnet respectively.
Start from a variety of initial
states to check whether the initial state affects the results.
\begin{center}
\setlength{\unitlength}{2pt}
(a) \begin{picture}(70,40)(0,-5)
\newsavebox{\verticalfour}
\savebox{\verticalfour}(0,0)[bl]{
\multiput(0,0)(0,10){4}{\circle{2}} % spins
\multiput(0,5)(0,10){4}{\line(0,-1){3}} % lines down
\multiput(0,-5)(0,10){4}{\line(0,1){3}} % lines up
\multiput(2,0)(0,10){4}{\line(1,0){6}} % lines right
}
\multiput(0,0)(10,0){6}{\usebox{\verticalfour}}
\end{picture}
(b) \begin{picture}(70,45)(0,-5)
\newsavebox{\verticalfourdiag}
\savebox{\verticalfourdiag}(0,0)[bl]{
\multiput(0,0)(0,10){4}{\circle{2}} % spins
\multiput(0,5)(0,10){4}{\line(0,-1){3}} % lines down
\multiput(0,-5)(0,10){4}{\line(0,1){3}} % lines up
\multiput(2,1)(0,10){4}{\line(2,1){6}} % lines rightup
\multiput(2,-1)(0,10){4}{\line(2,-1){6}} % lines rightdown
}
\multiput(0,0)(20,0){3}{\usebox{\verticalfourdiag}}
\multiput(10,5)(20,0){3}{\usebox{\verticalfourdiag}}
\end{picture}
\end{center}
Track the average energy and the standard deviation of the energy
of the system as a function of temperature.
Repeat the simulation for a variety of
values of $N$ and compare (start with $N=9,16,36$).
Also track the mean square magnetization. Do the
results appear to converge to some infinite $N$ limit? Examine both
the cases $J = \pm 1$. Can you extract heat capacities of the Ising
models from your computations? [Hints: $C = k_{\rm B} \beta^2 {\rm
var}(E)$ might be a useful trick. One awkward decision that
has to be made is how soon to start taking
measurements of energy and magnetization
after a new temperature has been established; it
is virtually impossible to detect `equilibrium'. A simple method
is to pick a total number of iterations $I$ (which should obviously be
significantly greater than the
number of spins $N$), and to assume equilibrium had been reached after
$I/3$ iterations; one can assess how well
equilibrium has been reached by running through a sequence of temperatures
first in one direction then in the other, looking to see if properties such
as the mean energy show hysteresis. If they do, then $I$ needs to be bigger.]
% \subsection*{Part b}
(B)
Repeat the above computations for a triangular Ising model.
Do you expect the triangular Ising model with $J = \pm 1$
to show different physical
properties from the square Ising model? Do you see evidence for
differences?
\subsection*{Optional extras}
Switch on the applied field $H$ and find the energy and mean
magnetization as a function of $\beta$. As $\b$ and
$H$ vary do you find evidence for phase transitions?
\medskip
\hfill DJCM \today
\end{document}