
COS 323  Computing for the Physical and Social Sciences

Fall 2013

Assignment 4: Simulating the Ising Model
Due Friday, Dec. 13
In this assignment we will be examining magnetization by simulating the
2D Ising model. In particular we will use simulation to investigate
the phase transition for magnetization.
The 2D Ising Model
In the 2D Ising model there is an array of atomic spins $s_i \in \{1,1\}$ in
a lattice. A state of the model $S$ is the configuration of all spins of the
lattice. For example, if you have a 3x3 lattice, and you number all the
points in this grid from 1 to 9 starting with the upperleft as 1, the point
underneath 1 as 4, and the lowerright as 9, an example configuration of
$s_1,s_2,s_3,s_4,\dots,s_9$ is $\{1,1,1,1,1,1,1,1,1\}$. A spin, $s_i$, has a neighbor $s_j$ if $s_j$ is adjacent to $s_i$ in the lattice. We are also
assuming periodic boundary conditions. By this we mean that if a spin is at the
edge of the lattice, it is neighbors with the other edges of the lattice if we
wrapped the lattice around. In the 3x3 lattice
example $s_1$ is neighbors with $s_2$,$s_4$,$s_3$,$s_7$. Under the assumptions given for this model for this assignment,
all spins have 4 neighbors.
The energy for a configuration, $s$, is given by its Hamiltonian:
$$E(s) = \sum_{(i,j) \in \cal N}J_{i,j}s_is_j  H\sum_{k}s_k$$
where $J_{i,j}$ is the coupling coefficient between neighboring pairs of spin sites in the lattice and
$H$ is a constant for the external magnetic field. The first sum is taken only
over neighboring pairs of spins and each pair is counted only once. For this assignment
we make the following assumptions: $H = 0$ and $J_{i,j} = J > 0$ for all neighboring pairs and 0 otherwise
The likelihood of a configuration is the following:
$$ P(s) = \frac{exp(E(s)/(\kappa T))}{Z(T)}$$
where $T$ is temperature, $\kappa$ is the Boltzmann constant, and $Z(T)$ is the normalization constant where the sum
is taken over all possible configurations, s:
$$Z(T) = \sum_{s}{exp(E(s)/(\kappa T))}$$
The magnetization of the material is the sum of the spins at all the sites
of the lattice. The spins will naturally try to align themselves to be the
same in the ferromagnetic model that we are using ($J > 0$). We are interested in examining the mean magnetization per spin site. The mean magnetization per spin site, $m$, is the mean over all $N$ spin sites in the
lattice of the sum of all of the spins:
$$m = \frac{1}{N} \sum_{i} s_i$$
In 1944 Lars Onsager showed analytically that the 2D Ising model without an external magentic field exhibited a phase transition. A phase transition in this context means that at temperatures below a critical temperature materials will exhibit nonzero spontenous magnetization and above the critical temperature the mean magnetization will go to zero[Gould]. Onsager determined the critical temperature, $T_c$, to be
$$\frac{\kappa T_c}{J}= \frac{2}{\ln (1+\sqrt 2)} \approx 2.269$$
For this assignment
and algorithm as given, we are assuming temperature is measured in units
such that $J/\kappa =1$ [Gould]. In the remainder of this assignment
we will simulate the Ising model using the Metropolis algorithm in
order to get simulated measurements for mean magnetization per spin in
order to investigate the phase transition relative to the critical temperature.
Monte Carlo Simulation: The Metropolis Algorithm
This section summarizes the Metropolis algorithm particularly for the
Ising model [Beichl, Gould] with the assumptions and
clarifications for how we will implement the algorithm. (More general
discussion of the Metropolis algorithm is offered in the brief
survey article by that name in the references.)
Algorithm
 Initialize the configuration by selecting the spin for each site
in the lattice to 1 or 1 at random.
 Choose a grid location uniformly at random and flip its sign.
 Compute the change of energy,$\Delta E$, that would occur as a result of
the flip (Hint: You need only consider neighbors of the chosen spin)
 If $\Delta E < 0$, accept the flip with probability 1.
 If
$\Delta E > 0$ accept the flip with probability
$e^{(\Delta E/\kappa T)}$. (Recall $J$ is constant and $J/\kappa = 1$.)
 Repeat 25 for $n^2$ times for an $n \times n$ lattice.
(The $n^2$ potential flips taken together
constitute a single unit of time.
We would like to have a unit of time in the
simulation such that each spin in the lattice may have
the chance to change in that unit of time and on average each spin is
chosen uniformly.)
 Repeat 26 for the desired number of iterations or time steps of the
simulation.
Your assignment:

Implement the Metropolis algorithm (given above) to simulate the Ising model
in the following matlab function:
function [m] = ising2d[size,temp,iter]
This function will take a lattice size, temperature, and number
of iterations and will simulate the model for the given number
of iterations. It will return the mean magnetization of each step
of the iteration. A step of the iteration will be the unit of time
in which a flip is possible for all spins in the lattice (i.e., all size x size spins).
 To test your implementation simulate the Ising model
for a 100 x 100 lattice for 1000 iterations
 For all temperatures
between [.5,3] in .25 intervals, create a scatter plot of the
mean magnetization per spin against temperature (For each temperature
there should be one point on the scatter plot.)
 For the temperatures T=1,2.269, and 3, create a scatter plot of the
mean magnetization per spin at each iteration of the simulation. (These
scatter plots should have 1000 points each.)
The code to run the trials and create the scatter
plots are to be in a script named ising_trials.m.

Explain
what you observe about the phase transition for magnetization per spin relative to the critical temperature in your scatter plots
You may discuss this entirely from your scatter plots.
How would you determine if you have run enough iteration steps
to feel confident in the validity of your simulation? (This is
for a single simulation temperature, but you can also consider
various simulation runs at different temperatures.)
 Suppose you would like to measure the sample variance of the mean
magnetization per spin of your
simulations. What are the difficulties in doing so for this model?
 Explain how the algorithm as given uses the rejection method.
Specify the specific step(s) in the algorithm above.
 At what temperature in your trials above, is your simulation rejecting more flips on average? (If this is not clear, you can also code for this.)
 Importance sampling could also be used. Explain what
distribution we would need to use importance sampling directly and one challenge with
doing so for this model.
References
 Isabel Beichl and Francis Sullivan. 2000.
The Metropolis Algorithm. Computing in Science and Engg. 2, 1 (January 2000), 6569. DOI=10.1109/5992.814660
 Harvey Gould and Jan Tobochnik.
"STP Textbook Chapter 5: Magnetic Systems." In Thermal and Statistical Physics. 2008.
Submitting
This assignment is due Friday, December 13, 2013 at 11:59 PM.
Please see the course policies
regarding assignment submission, late assignments, and collaboration.
Please submit:
 Wellcommented source files ising2d.m and ising_trials.m
as well as an auxilliary files you may have.
 A writeup in plaintext or PDF format, containing:
 A description of your code
 All 4 scatter plots
 Answers to all bold questions
The Dropbox link to submit your assignment is
here.
Last update
3Dec2013 16:47:20
smr at princeton edu