COS 323 - Computing for the Physical and Social Sciences

Fall 2013
Course home Outline and lecture notes Assignments


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 upper-left as 1, the point underneath 1 as 4, and the lower-right 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

  1. Initialize the configuration by selecting the spin for each site in the lattice to -1 or 1 at random.
  2. Choose a grid location uniformly at random and flip its sign.
  3. 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)
  4. If $\Delta E < 0$, accept the flip with probability 1.
  5. If $\Delta E > 0$ accept the flip with probability $e^{(-\Delta E/\kappa T)}$. (Recall $J$ is constant and $J/\kappa = 1$.)
  6. Repeat 2-5 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.)
  7. Repeat 2-6 for the desired number of iterations or time steps of the simulation.

Your assignment:

  1. 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).
  2. To test your implementation simulate the Ising model for a 100 x 100 lattice for 1000 iterations The code to run the trials and create the scatter plots are to be in a script named ising_trials.m.
  3. 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.)
  4. 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?
  5. Explain how the algorithm as given uses the rejection method. Specify the specific step(s) in the algorithm above.
  6. 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.)
  7. 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

  1. Isabel Beichl and Francis Sullivan. 2000. The Metropolis Algorithm. Computing in Science and Engg. 2, 1 (January 2000), 65-69. DOI=10.1109/5992.814660
  2. 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:

The Dropbox link to submit your assignment is here.

Last update 3-Dec-2013 16:47:20
smr at princeton edu