Workshop Part 1 - Ising Model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 1: The Ising Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Goals of this part\n",
    "1. Learn what the Ising model is, and (roughly) how macroscopic magnetism arises from microscopic components\n",
    "2. Write a simple code to simulate the simplest model of magnetism\n",
    "3. Explore some basic Julia linear algebra features on a single node\n",
    "\n",
    "We'll start with the *classical* Ising model, a simple physical model of magnetism. As we add complexity to the model, we'll be able to represent more interesting physics more accurately. Hopefully, the complexity of our *code* will grow more slowly than the complexity of the physics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction to the Model\n",
    "\n",
    "The Ising model is a powerful but simple model of magnetism in materials. Small constituents, local magnetic moments, interact with each other to give \"macroscopic\" phenomena, like long range correlations we observe as ferromagnetism.\n",
    "\n",
    "A magnetic moment is a local \"unit\" of magnetism. We say that if something (an atom, an electron, an ion) has *non-zero* magnetic moment then it \"feels\" magnetic fields and responds to them. The grains in a hard drive have non-zero magnetic moment, as do all the nuclei in a human being (that's how MRI works).\n",
    "\n",
    "These local magnetic moments will add up (or not) to give the macroscopic moment you'd see in a compass needle. The magnetic moment of a big object, made of many little magnetic pieces, is just the sum of all their magnetic moments (we'll return to this later). If most of the magnetic moments are *positive*, the sum will be a large positive number, and the magnet points north. If most are *negative*, the sum will be a large negative number, and the magnet will point south. All the parts of a human have magnetic moment but they average out to pointing no direction at all.\n",
    "\n",
    "As an analogy, imagine a group of people holding a rope trying to pull a 747. If all of them pull in the same direction, their small individual efforts may combine to enough force to move the jet. However, if everyone is too busy pulling in whatever direction they feel, the jet won't move at all. If the people are on slippery ground and end up flummoxed, prevented from working together, again they cannot generate enough force to move the airplane. We want to find out what combination of community-mindedness and dry dirt will let the people move/not move the airplane. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make this a little concrete. For the workshop, we'll work in *one dimension* on a *chain of spins*. This is both easy for us to visualize and easy to represent on a computer. We have a single \"wire\" crystal, where on each site in the wire, the *spin* (local magnetic moment) can be *up* or *down*.\n",
    "\n",
    "Here's a sample configuration, represented in many equivalent ways:\n",
    "\n",
    "![separatespins](pics/separatespins.svg)\n",
    "![spinstogether](pics/spinstogether.svg)\n",
    "![zerosandones](pics/zerosandones.svg)\n",
    "\n",
    "Another model you might have seen built out of similar constituents is [Conway's Game of Life](https://en.wikipedia.org/wiki/Conway%27s_Game_of_Life). If you're familiar with cellular automata, this stuff might be sounding familiar as well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Physical systems want to minimize their energy. The Ising model is a set of energetic interactions which add up to give the total energy of the system, which we'll write Julia code to minimize. By simulating the model we can learn about what physics is favourable in different parameter regimes. Another way of saying this is:\n",
    "- If the small constituent parts have *this sort* of interaction (a set of parameters), what will we observe *overall* (the macroscopic physics)?\n",
    "- How far can we vary the interactions and keep the same physics (parameter regimes)?\n",
    "- What will happen when we cross from one regime into another?\n",
    "\n",
    "Here is the (classical) energy function, which we call a Hamiltonian:\n",
    "\n",
    "$$ H = -\\sum_{\\langle i, j \\rangle} \\sigma_i \\sigma_j $$\n",
    "\n",
    "There are a couple things we should explain here:\n",
    "  - $\\langle i, j \\rangle$ means \"nearest neighbors\", so look at how site 1 and site 2 interact, but not site 1 and site 8.\n",
    "  - $\\sigma_i$ is the local configuration (the arrow) on site $i$. It can be 0 (down and blue) or 1 (up and red). We can specify the configuration of the entire system $\\{\\sigma\\}$ with a list of all the local configurations.\n",
    "  - Right now there aren't any *parameters* - in a \"real\" simulation, temperature could be one, but for now we just want to encode some physics in a simple matrix\n",
    "  \n",
    "So, this classical Hamiltonian says \"sum up over all the connections on the square grid - neighbors having the same value is good, having different values is bad\".\n",
    "\n",
    "Note that there is *no* preferred direction/value - all being 0 (blue/down) is just as good as all being 1 (red/up)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Now we can write down some code! Right now our Hamiltonian just measures the energy of a given configuration - it is *diagonal* in the spin basis. So:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "L = 10\n",
    "function energy(config)\n",
    "    E = 0.\n",
    "    for i in 1:L^2\n",
    "        # horizontal site\n",
    "        j = mod(i, L) + L*div(i, L)\n",
    "        if j > L^2 # literal corner case :)\n",
    "            j = L^2 - L\n",
    "        end\n",
    "        E += -1 + 2*xor(config[i], config[j])\n",
    "        # vertical site\n",
    "        j = i + L\n",
    "        if j > L^2\n",
    "            j = mod(i + L, L^2) + 1\n",
    "        end\n",
    "        E += -1 + 2*xor(config[i], config[j])\n",
    "    end\n",
    "    return E\n",
    "end\n",
    "\n",
    "# make a matrix H\n",
    "H = zeros(2^L, 2^L)\n",
    "basis = Vector{Bool}[]\n",
    "\n",
    "# translate the integer into a bit representation\n",
    "for element in 1:2^L\n",
    "    bit_rep = falses(L^2)\n",
    "    for site in 1:L^2\n",
    "       bit_rep[site] = (element >> (site - 1)) & 1\n",
    "    end\n",
    "    push!(basis, bit_rep)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could have also used the Julia `digits` function to achieve the same outcome:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Bool,1}:\n",
       " false\n",
       " false\n",
       "  true\n",
       " false\n",
       "  true"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Vector{Bool}(digits(20, 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The point is, we're going to use whether the bit is `0` or `1` on each \"site\" (site `i` is indexed as bit `i` in the integer) to encode the \"state\" of up/down on each site, then make a matrix describing the energy of each state. We'll deal with small systems for now, so whatever method is clearest to you is completely fine to use! First, let's check we did this correctly (each element of `basis` should be unique)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@assert length(unique(basis)) == 2^L \"Basis elements aren't unique!\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for (index, configuration) in enumerate(basis)\n",
    "    H[index, index] = energy(configuration)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.733284062"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@elapsed eig(H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.011597495"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@elapsed eig(Diagonal(H))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since this model is completely diagonal in the *simulation basis* (of lists of up/downs) it's pretty boring. In quantum mechanics, the *eigensystem* of the Hamiltonian is extremely important. The *eigenvectors* are the quantum states and the *eigenvalues* are the corresponding energies of the states. The lowest energy state is called the *groundstate* and it's the one we \"see\" at very low temperatures. The next lowest is the *first excited state*.\n",
    "\n",
    "\n",
    "Time to add some complications to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0-rc1",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}