Project_3_Instructions.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.4"
    },
    "colab": {
      "name": "Project 3 Instructions",
      "provenance": []
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_V7TRkFV8Y0B",
        "colab_type": "text"
      },
      "source": [
        "# Project 3: Markov Chain Monte Carlo (MCMC) and Statistical Mechanics"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1ZA_Sb-V8Y0E",
        "colab_type": "text"
      },
      "source": [
        "## Simulate a biased coin\n",
        "\n",
        "Consider the game we discussed in lecture in which Helena flips a coin.  Helena wins $\\$ 1$ if the coin comes up heads and loses $\\$ 1$ if the coin comes up tails.\n",
        "\n",
        "- Suppose you have a sequence $o_1, o_2, \\dots, o_n$ of real numbers corresponding to, say, measured values of some observable quantity $O$.  The $k^\\mathrm{th}$ **running average** of this sequence is defined as follows:\n",
        "\\begin{align}\n",
        "    \\langle O \\rangle_k = \\frac{o_1 + o_2 + \\cdots + o_k}{k}.\n",
        "\\end{align}\n",
        "Show (type out your proof in LaTeX) that the $k^\\mathrm{th}$ running averages satisfy the following identity:\n",
        "\\begin{align}\n",
        "    \\langle O\\rangle_{k+1} = \\langle O\\rangle_k + \\frac{1}{k+1}(o_{k+1} - \\langle O\\rangle_k).\n",
        "\\end{align}\n",
        "This is a useful identity because it means that if you're generating the sequence of $o_k$'s one after the other in, say, a loop, computing the running average on each iteration doesn't require summing up all of the previous values and dividing by the total number of values every single time -- you can simply update the prior running average by adding a change term that depends on the next value in the sequence and the current running average.\n",
        "- Write a Python function called \n",
        "\n",
        "        weighted_coin \n",
        "\n",
        "  whose inputs are:\n",
        "  - `beta`: a float that gives the probability of the coin landing heads (it's \"bias\" toward heads if you will)\n",
        "  - `n`: a positive integer representing the number of flips of the coin\n",
        "  \n",
        "  and whose output is Helena's average earnings per flip in dollars based on a Metropolis MCMC simulation of `n` flips of the coin as well as a plot that tracks the average earnings per flip during the simulation.  This plot should have the average earnings per flip on the vertical axis and the simulation \"time\" on the horizontal axis.  You may choose to include an additional argument in your function that allows you to only keep track of and plot the running average after having waited a specified number of iterations each time you sample the running average so that the code will run faster.  For example, if you want to flip the coin 10 million times, then you probably don't need to be plotting the running average *every* flip but instead maybe every 10 thousand flips (or some reasonable number of your choosing). \n",
        "- Write a Python function \n",
        "\n",
        "        average_earnings_per_flip \n",
        "\n",
        "  whose input is `beta`, a float that gives the probability of a coin landing heads, and whose output is the average earnings per flip in dollars calculated exactly (following the method a human might use to compute the average earnings per flip on a piece of paper).\n",
        "- Generate a grid of two plots: one that gives the theoretical average earnings per flip as a function of `beta` on the interval $[0, 1]$ and compares it to the average earnings per flip computed via Metropolis MCMC simulation of 1 thousand flips for 20 equally-spaced values of `beta` on the interval $[0, 1]$, and one that's the same except thatt uses 1 million flips for each run of the simulation.\n",
        "\n",
        "\n",
        "## Simulate a weighted die\n",
        "\n",
        "Imagine that you are playing a game with your friend in which you roll a 6-sided die.  The way the game works is that every time the die lands on side 1 or side 2, you win one dollar, and every time the die lands on side 3, 4, 5, or 6, your friend wins one dollar.  You want to win the game at all costs, so naturally you decide to cheat and use a weighted die that's preferentially weighted so that sides 1 and 2 are each three times more likely to land face up than sides 3, 4, 5, and 6.\n",
        "\n",
        "- Write a function\n",
        "\n",
        "        weighted_die\n",
        "\n",
        "  whose input is\n",
        "  - `n`: the number of times the die is rolled\n",
        "  that runs a Metropolis MCMC simulation and outputs the average earnings per roll after rolling the die `n` times. \n",
        "- Deteremine roughly how many steps you need before the simulation's average earnings seems to converge to the nearest $\\$0.01$.  In other words, play around with the simulation to determine how many steps in the simulation you need so that when you run the simulation with that many steps, the result is pretty much always the same value to within one cent.\n",
        "- Figure out the expected value of your earnings analytically and compare this to what your simulation produces as a sanity check."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rOtodzeIHzEL",
        "colab_type": "text"
      },
      "source": [
        "# Your computational task relating to a physical system\n",
        "\n",
        "Write a function\n",
        "\n",
        "    two_dim_ising(L, temp, num_steps)\n",
        "    \n",
        "that simulates the 2-d Ising model on a square, periodic lattice of a specified side length (number of spins) `L` and at a given temperature.  The inputs `temp` and `num_steps` specify the desired temperature of the lattice and the number of MCMC updates you wish to perform during the simulation.  The function should return whatever you think is useful in completing the Jupyter notebook described below."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-oGyuKSE8Y0F",
        "colab_type": "text"
      },
      "source": [
        "# Some details on the classical 2D Ising model\n",
        "\n",
        "The 2D Ising model is a simple model of a ferromagnetic material with a phase transition.\n",
        "The model consists of a 2D lattice $L\\times L$ of spins $s_i\\in\\{-1,+1\\}$. Each spin interacts only with\n",
        "its nearest neighbors.  When the lattice is finite The energy is expressed as\n",
        "\n",
        "\\begin{equation}\n",
        "E(\\{s_i\\}) = - \\sum_{\\langle i,j\\rangle}s_is_j - H\\sum_i s_i,\n",
        "\\end{equation}\n",
        "\n",
        "where $\\{s_i\\}$ is notation for the entire configuration of spins, $H$ is the external magnetic field, and $\\langle i,j\\rangle$ implies a summation over all nearest-neighbour pairs.  A nearest neighbor is a spin either directly to the right or left of a given spin, or directly above or below.\n",
        "We have normalized energy with $J$, spin with $\\hbar/2$, and magnetic field with $J/\\mu$, where $J$ is the \n",
        "exchange energy and $\\mu$ is the atomic magnetic moment.\n",
        "When the system is in contact with a heat bath at temperature $T$ the equilibrium probability density is the \n",
        "Boltzmann distribution\n",
        "\n",
        "\\begin{equation}\n",
        "    p(\\{s_i\\}) = Z^{-1}\\exp\\bigl(-E(\\{s_i\\})/T\\bigr),\n",
        "\\end{equation}\n",
        "\n",
        "where the partition function $Z$ is the sum of exponential factors $\\exp(-E/T)$ over all possible configurations.\n",
        "Below is the critical temperature\n",
        "\n",
        "\\begin{equation}\n",
        "T_c = \\frac{2}{\\log(1+\\sqrt 2)}\\approx 2.2692\n",
        "\\end{equation}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cG3PHL928Y0G",
        "colab_type": "text"
      },
      "source": [
        "## Write a Jupyter notebook that explores and describes the physics of the 2d Ising model\n",
        "\n",
        "\n",
        "Your notebook should apply the function `two_dim_ising` to investigate the physics of the 2d Ising model on an $N = L\\times L$ periodic lattice in the absence of an external magnetic field.  Periodicity of the lattice means that if the lattice configuration is represented as an array $s_{ij}$ of spins\n",
        "\n",
        "\\begin{align}\n",
        "    \\begin{pmatrix}\n",
        "        s_{00} & \\cdots & s_{0, L-1} \\\\\n",
        "        \\vdots & \\ddots & \\vdots \\\\\n",
        "        s_{L-1,0} & \\cdots & s_{L-1, L-1} \\\\\n",
        "    \\end{pmatrix}\n",
        "\\end{align}\n",
        "\n",
        "Then the nearest neigbors of a spin on the edge are identified by \"wrapping the lattice around\" to the opposite edge.  For example, for a spin $s_{i, 0}$ on the left edge of the lattice, the nearest neighbor on the left would be $s_{i, L - 1}$, a spin that is on the right edge of the lattice, and similarly for the bottom and top of the lattice.  In this notebook, you will explore the behaviors of the mean internal energy $U$, \n",
        "magnetization $M$, specific heat $C_H$, and magnetic susceptibility $\\chi_T$ (per lattice site) on the temperature, defined as\n",
        "\n",
        "\\begin{align}\n",
        "U& = \\frac{1}{N}\\langle E\\rangle ,& M& = \\frac{1}{N}\\langle S\\rangle ,\\\\\n",
        "\\chi_T& = \\frac{1}{NT}\\left(\\langle S^2\\rangle  - \\langle S\\rangle^2\\right),&\n",
        "C_H& = \\frac{1}{NT^2}\\left(\\langle E^2\\rangle - \\langle E\\rangle^2\\right),\n",
        "\\end{align}\n",
        "\n",
        "where $N=L^2$ is the total number of sites, and $S = \\sum_i s_i$ is the net magnetization.  One way to implement the Metropolis-Hastings algorithm for the Ising model on a square lattice is the following:\n",
        "\n",
        "1. Pick a random site $i$ on the 2D lattice and compute the energy change $\\Delta E$ due to the change of \n",
        "sign in $s_i$:\n",
        "\n",
        "    \\begin{align}\n",
        "        \\Delta E = 2s_i\\bigl(s_{\\text{top}} + s_{\\text{bottom}} + s_{\\text{left}} + s_{\\text{right}} + H\\bigr),\n",
        "    \\end{align}\n",
        "    \n",
        "where $s_{\\text{top}}$, $s_{\\text{bottom}}$, $s_{\\text{left}}$, and $s_{\\text{right}}$ are the 4\n",
        "nearest neighbours of $s_i$.\n",
        "1. If $\\Delta E\\leq 0$ accept the move. If $\\Delta E>0$ accept the move with probability $A=\\exp(-\\Delta E/T)$.\n",
        "1. Flip the spin $s_i$ if the move has been accepted. \n",
        "1. Repeat steps 1-3 until you generate a large sample of spin configurations.\n",
        "\n",
        "To save memory, use new samples $\\{s_i\\}_{n+1}$ immediately as they arrive to update your estimates according to the rule\n",
        "\n",
        "\\begin{align}\n",
        "\\langle O\\rangle_{n+1} = \\langle O\\rangle_{n} + \\frac{1}{n+1}\\Bigl(O\\bigl(\\{s_i\\}_{n+1}\\bigr) - \\langle O\\rangle_{n}\\Bigr),\n",
        "\\end{align}\n",
        "\n",
        "where $\\langle O\\rangle_n$ is the previous estimate for the mean of some observable $O$ obtained from a \n",
        "sequence of $n$ samples, and $\\langle O\\rangle_{n+1}$ is the new improved estimate.\n",
        "Avoid also updating the energy and magnetization by looping over all sites of the $L\\times L$ lattice. \n",
        "It is necessary to perform this time-consuming operation only once at the very \n",
        "start of the simulation. Afterwards you can keep track of $E$ and $S$ by adding the increments \n",
        "$\\Delta E$ and $\\Delta s$ to the old values.\n",
        "\n",
        "### Plotting time series' of intensive quantities\n",
        "\n",
        "As you run your MCMC simulation, you should find that the \n",
        "for larger system size, the Monte Carlo simulation needs to run longer steps, to reach equilibrium.  In order to monitor the convergence of values of observables during the simulation, it is helpful to plot their values as a function of time step.  In order to make the benchmarking a fair game, the \"updating step per site\" is defined as follows:\n",
        "\n",
        "\\begin{align}\n",
        "    t=\\frac{\\text{total number of updates}}{\\text{number of sites}}\n",
        "\\end{align}\n",
        "\n",
        "**What you need to do:**\n",
        "\n",
        "- Choose system sizes $L = 16, 32$ and a temperature of your choosing not to close to the critical temperature, and run the simulation.  Plot $U$ and $M$ as functions of $t$, the update step per site defined above.  \n",
        "- Roughly how long does it seem to take for the simulation to converge?  \n",
        "- Does the time it takes seem to depend on the system size?\n",
        "\n",
        "### Magnetization curves for different lattice sizes\n",
        "\n",
        "Onsager's exact solution gives the following result for the magnetization as a function of temperature:\n",
        "\n",
        "\\begin{align}\n",
        "M(T)& = \\begin{cases}\\Bigl[1 - \\sinh^{-4}\\bigl(2/T\\bigr)\\Bigr]^{1/8},& T<T_c\\\\\n",
        "0,& T\\geq T_c\n",
        "\\end{cases}\n",
        "\\end{align}\n",
        "\n",
        "obtained in the thermodynamic limit $N\\to\\infty$. \n",
        "\n",
        "**What you need to do:**\n",
        "\n",
        "- Consider system sizes $L = 8, 16, 32, 64$. \n",
        "- For each size, generate an $M(T)$ curve.  Plot these 4 curves, together with the exact curve, *in the same plot.* - We suggest you start the simulations at a high temperature above the critical point $T_c$  and slowly cool the ferromagnet by decreasing the temperature in small steps. After each update of $T$, perform a large number of iterations with the Metropolis-Hastings algorithm without calculating the ensemble averages. Once the system is close to thermal equilibrium start drawing the samples from the equilibrium distribution and calculate the averages along the way.\n",
        "\n",
        "**Notes.**\n",
        "\n",
        "- These are computationally expensive calculations.\n",
        "- There are 4 curves for different system sizes.\n",
        "- Each curve is made of many dots. (you decide the spacing of temperatures)\n",
        "- Each dot is an averaged magnetization per site at a given temperature\n",
        "\n",
        "### Typical spin configuration at different temperatures\n",
        "\n",
        "Once the energy and magnetization have come to equilibrium in your simulation, you can look at the configuration of the system to see what a typical configuration of the system looks like in thermal equilibrium.  One way to do this is to make a plot of the whole array of spins with a black or white square representing spin up or spin down.\n",
        "\n",
        "**What you need to do:**\n",
        "\n",
        "- Choose your system size to be $L = 100$.  \n",
        "- Make sure run your simulation enough time steps, until the system is in equilibrium.\n",
        "- Plot three pictures of spin configuration, at temperature $T=1.8$, $T=2.3$, $T=4.0$.  Your pictures should display the grid of spins with either a black or white square at a given site depending on whether the site is in a spin up state or a spin down state.\n",
        "- Comment on how the patterns differ at different temperatures and what they might physically tell us about how the system behaves at temperatures below, near, and above the critical temperature.\n",
        "\n",
        "### Further exploration if you're interested\n",
        "\n",
        "- Consider the case where the external magnetic field $H$ is nonzero, and investigate the physics of the model.\n",
        "- Explore the behaviors of $C_H$ and $\\chi_T$ in addition to $U$ and $M$ as a function of temperature and lattice size."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "imKOHy3ZSyY0",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}