activity-2-1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "K7PbPUVeY-ZD"
   },
   "source": [
    "# Activity 2-1: Simulating the Ising model\n",
    "\n",
    "## The Ising model\n",
    "\n",
    "As a simple model of a ferromagnetic material, we can consider a set of spins arranged on a lattice. The energy of a pair of spins will be lower if they are aligned in the same direction. To make our analysis simple, we'll assume that only **nearest neighbor** spins are coupled together. Spins that aren't nearest neighbors could also interact, but these interactions will be much weaker due to the greater distance between them.\n",
    "\n",
    "Mathematically, we can write the energy of a configuration $\\underline{\\sigma}=\\{\\sigma_1, \\sigma_2, \\ldots, \\sigma_N\\}$ of spins $\\sigma_i\\in\\{-1,1\\}$ in a one-dimensional lattice as\n",
    "\n",
    "$$ \n",
    "E(\\underline{\\sigma}) = -\\sum_{i=1}^{N-1} J \\sigma_i\\,\\sigma_{i+1}\\,.\n",
    "$$\n",
    "\n",
    "The coupling $J$ tells us how strong the interaction is between two neighboring spins.\n",
    "\n",
    "## Monte Carlo simulation\n",
    "\n",
    "If we want to measure properties of an Ising model with $N$, such as the average magnetization of a particular spin, we won't be able to average over all configurations of the system unless $N$ is very small (or unless we are very clever see the [Wikipedia page](https://en.wikipedia.org/wiki/Ising_model) for much more information)!\n",
    "\n",
    "However, we can use a **Monte Carlo simulation** to get information about the Ising model even if we can't average over all configurations exactly. To do this, we'll start in a particular configuration of the system. Then, we'll randomly change the configuration in a way that satisfies **detailed balance**. There are many possible ways to do this, but the following rule (known as the [Metropolis rule](https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm)) is one of the simplest:\n",
    "\n",
    "- Choose a spin at random  \n",
    "- Compute the difference between the energy we would have if the chosen spin was \"flipped\" from $+1$ to $-1$, or vice versa, and the energy of the current configuration, which we call $\\Delta E$    \n",
    "- If $\\Delta E<0$, so that the new energy is less than the current one, then the spin is always flipped  \n",
    "- Otherwise, the spin is flipped with probability $p = e^{-\\beta \\Delta E}$, where $\\beta = 1 / T$ and $T$ is the temperature\n",
    "\n",
    "## 1. Building the simulation step by step\n",
    "\n",
    "Monte Carlo simulation of a large system usually can't be performed in just a few lines of code. We'll make our work simpler by writing down the things that our program needs to do, and then writing a function to execute each of those things.\n",
    "\n",
    "In general, the steps are:\n",
    "\n",
    "1. Prepare an initial configuration of the system  \n",
    "2. Iterate Monte Carlo  \n",
    "    a. Randomly choose a spin to flip  \n",
    "    b. Compute the change in energy if the spin flips  \n",
    "    c. Draw a random number and choose whether or not the spin flips  \n",
    "    d. Update and save the configuration  \n",
    "3. Analyze and visualize the results  \n",
    "\n",
    "### 1.a. Set up the initial configuration\n",
    "\n",
    "Let's call our configuration `spins`. We'll choose a system size of $N=20$. We'll also try starting with a random configuration of spins instead of a fixed order. \n",
    "\n",
    "Fill in the code below to define a function, `initialize_spins`, that returns a list of spins with random orientations (that is, each spin is assigned $+1$ or $-1$ at random with equal probabilities)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.random as rng\n",
    "\n",
    "# Define the system size and initialize random number generator\n",
    "\n",
    "N = 20\n",
    "r = rng.RandomState(0)\n",
    "\n",
    "# Write the function to produce a random initial configuration\n",
    "\n",
    "def initialize_spins():\n",
    "    \"\"\" \n",
    "    Returns a configuration of N spins, where each spin is assigned\n",
    "        +1 with probability 1/2, and\n",
    "        -1 with probability 1/2\n",
    "    \"\"\"\n",
    "    \n",
    "    spins = []\n",
    "    for i in range(N):\n",
    "        # < Your code here! >\n",
    "        # Add a new entry to `spins` that is either +1 or -1 with probability 1/2\n",
    "        # Hint: you could do this using r.rand() combined with if/else statements,\n",
    "        # or you could use the numpy.random.choice function (more on this later)\n",
    "            \n",
    "    # Return the result\n",
    "    \n",
    "    return spins"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.b. Test `initialize_spins`\n",
    "\n",
    "Execute the code block below to test your function. If it works properly, you should see a set of $20$ numbers that are equal to $\\pm 1$, with about an even number of positive and negative spins (though there will be a substantial amount of variation!)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the initial configuration of spins\n",
    "\n",
    "N = 20\n",
    "spins = initialize_spins()\n",
    "\n",
    "# And print the results\n",
    "\n",
    "print('The length of the configuration is {}'.format(len(spins)))\n",
    "print('There are {} + spins and {} - spins'.format(np.sum(np.array(spins)==1), np.sum(np.array(spins)==-1)))\n",
    "print('The configuration is:')\n",
    "print(spins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.c. Choose a random spin\n",
    "\n",
    "Next, during each step of the MC algorithm we need to choose a spin at random, which we will attempt to flip. To do this, we'll use a function `numpy.random.choice`, which selects a random element from a list. By default, the probability of each element is equal. We can also pass this function a number $N$ instead of a list to randomly choose a number between $0$ and $N-1$. For more information, search this function in the `NumPy` docs or use `help(np.random.choice)`.\n",
    "\n",
    "Execute the code below to define a function, `choose_spin`, that will select a spin index uniformly at random from the lattice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def choose_spin():\n",
    "    \"\"\"\n",
    "    Choose spin indices at random. \n",
    "    We assume that the lattice size N has already been defined.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Return a random index between 0 and N-1\n",
    "    \n",
    "    return np.random.choice(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.d. Test `choose_spin`\n",
    "\n",
    "Execute the code block below to test your function. If it works properly, you should see a nearly uniform distribution over the range from $\\left[0, 1, \\ldots, N-1\\right]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "indices = []\n",
    "\n",
    "# Randomly sample a lot of indices\n",
    "\n",
    "for k in range(10**4):\n",
    "    indices.append(choose_spin())\n",
    "\n",
    "# And check their distribution\n",
    "\n",
    "sns.distplot(indices, bins=np.arange(0, N+1, 1), norm_hist=True,\n",
    "             kde=False, hist_kws={'align': 'left'})\n",
    "plt.title(r'Distribution of indices over $10^4$ random samples')\n",
    "plt.xlabel('index')\n",
    "plt.ylabel('frequency')\n",
    "plt.xticks(range(N));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.e. Compute the change energy of a configuration\n",
    "\n",
    "The most important step we need to complete is to compute how the energy would change if we flipped a selected spin $\\sigma_{i}$ to $-\\sigma_{i}$. For example, if $\\sigma_{i}=-1$, how would the energy change if we set $\\sigma_{i}=1$ instead?\n",
    "\n",
    "To do this, we can write a function to compute the energy of each configuration and take the difference in the energy of the two configurations (with and without the spin flipped). Fill in the code below to define a function `compute_E` to compute the energy of a configuration.\n",
    "\n",
    "**Note: you could also compute the difference in energy $\\Delta E$ directly! If you're up for more of a challenge, instead of writing code to compute the energy of a configuration, try writing a function that computes $\\Delta E$.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the function to compute the energy of a configuration\n",
    "\n",
    "def compute_E(spins):\n",
    "    \"\"\" This function computes the energy of a spin configuration assuming J=1. \"\"\"\n",
    "    \n",
    "    energy = 0\n",
    "    \n",
    "    # Fill in code here to add up the energy\n",
    "    \n",
    "    return energy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.f. Test `compute_E`\n",
    "\n",
    "Execute the code block below to test your function. We'll prepare two special configurations and check to see if the computed $\\Delta E$ is the same as what we expect. If you don't get the expected result, return to 1.e. above to revise your function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the test configuration\n",
    "\n",
    "N = 20\n",
    "spins = [1 for i in range(N)] # this is called a list comprehension!\n",
    "spins_flipped = [i for i in spins] # this line creates a copy of `spins`\n",
    "\n",
    "\n",
    "# And test compute_E\n",
    "\n",
    "spins_flipped[-1] = -1\n",
    "dE = compute_E(spins_flipped) - compute_E(spins)\n",
    "if dE == 2:\n",
    "    print('Success! Got delta_E = %.2f, expected delta_E = %.2f' % (dE, 2))\n",
    "else:\n",
    "    print('Test unsuccessful! Got delta_E = %.2f, expected delta_E = %.2f' % (dE, 2))\n",
    "\n",
    "spins[-1] = -1\n",
    "spins_flipped[-2] = -1\n",
    "dE = compute_E(spins_flipped) - compute_E(spins)\n",
    "if dE == 0:\n",
    "    print('Success! Got delta_E = %.2f, expected delta_E = %.2f' % (dE, 0))\n",
    "else:\n",
    "    print('Test unsuccessful! Got delta_E = %.2f, expected delta_E = %.2f' % (dE, 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Code the simulation\n",
    "\n",
    "Now we've defined all the helper functions that we need to run our Monte Carlo simulation. Recall that our steps were\n",
    "\n",
    "1. Prepare an initial configuration of the system  \n",
    "2. Iterate Monte Carlo  \n",
    "    a. Randomly choose a spin to flip  \n",
    "    b. Compute the change in energy if the spin flips  \n",
    "    c. Draw a random number and choose whether or not the spin flips  \n",
    "    d. Update and save the configuration  \n",
    "3. Analyze and visualize the results  \n",
    "\n",
    "To choose whether or not the spin flips, we'll follow the [Metropolis rule](https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm),\n",
    "\n",
    "$$\n",
    "P({\\rm flip}) = \\begin{cases} \n",
    "1 &\\text{if }\\Delta E < 0\\\\\n",
    "e^{-\\Delta E/T} &\\text{if } \\Delta E > 0\n",
    "\\end{cases} \\,.\n",
    "$$\n",
    "\n",
    "\n",
    "Fill in the code below to execute these steps. We will put our Monte Carlo simulation code inside a function `MC_simulate` so that we can reuse it easily. Analysis and visualization will be in part 3 below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.random as rng\n",
    "\n",
    "\n",
    "def MC_simulate(T, n_steps):\n",
    "    \"\"\"\n",
    "    This function runs a Monte Carlo simulation of the 1D Ising model\n",
    "    and returns an array of n_steps configurations of the system.\n",
    "    We assume that the coupling strength J is equal to 1.\n",
    "    The temperature T and number of Monte Carlo steps n_steps are\n",
    "    passed as parameters.\n",
    "    \"\"\"\n",
    "\n",
    "    # Define global variables\n",
    "\n",
    "    N = 20        # system size\n",
    "    configs = []  # container for configurations\n",
    "\n",
    "    # Randomize the starting configuration\n",
    "\n",
    "    spins =  # < Your code here! Use initialize_spins from 1.a to do this. >\n",
    "\n",
    "    # Begin MC simulation\n",
    "\n",
    "    for _step in range(n_steps):\n",
    "\n",
    "        # Record the current configuration\n",
    "        # NOTE: we have to create a copy, or all configurations will be one copy of the last one\n",
    "\n",
    "        configs.append([i for i in spins])\n",
    "        \n",
    "\n",
    "        # Choose a random spin `k` to attempt to flip\n",
    "\n",
    "        # < Your code here! Use choose_spin from 1.c to do this. >\n",
    "        \n",
    "\n",
    "        # Compute the change in energy\n",
    "\n",
    "        # step 1: make a copy of the current spins\n",
    "\n",
    "        # < Your code here! If you're unsure, look back to 1.f above. >\n",
    "        \n",
    "\n",
    "        # step 2: change spin k from +1 to -1 (or vice versa) in the copied configuration\n",
    "\n",
    "        # < Your code here! >\n",
    "        \n",
    "\n",
    "        # step 3: compute the change in energy\n",
    "\n",
    "        # < Your code here! Use compute_E to do this. >\n",
    "        \n",
    "\n",
    "        # Check to see if the spin flips\n",
    "        # Here we use a new if/else statement `elif`, which is short for \"else if\"\n",
    "\n",
    "        if delta_E < 0:\n",
    "            # < Your code here! Energy is lowered, so we flip the spin at index k. >\n",
    "\n",
    "        elif rng.rand() < np.exp(-delta_E/T):\n",
    "            # < Your code here! Energy is higher but check passes, so we flip the spin at index k. >\n",
    "\n",
    "    # Finally, return the stored configurations\n",
    "\n",
    "    return configs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Analyze and visualize the results\n",
    "\n",
    "### 3.a. The Ising model at high temperature\n",
    "\n",
    "First, let's run a simulation at a temperature of $T=5$. We'll take $10^4$ Monte Carlo steps in our simulation. Execute the code below to gather the simulation data and display it. We'll use some functions from `matplotlib` to make an [animation](https://towardsdatascience.com/how-to-create-animated-graphs-in-python-bb619cc2dec1) so that we can visualize how the configuration of the system changes!\n",
    "\n",
    "**Note:** animation can be tricky, and the libraries needed to create animations are not as well standardized as other ones we've been using. **If this code fails**, you may be able to troubleshoot the problem (for example, this might require installing a library like `ffmpeg`). However, I've confirmed that this approach does work on Google Colaboratory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import HTML\n",
    "import matplotlib.animation as animation\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "%matplotlib inline\n",
    "\n",
    "# Gather the simulation data\n",
    "\n",
    "configs = np.array(MC_simulate(5, 10000)).reshape((10000, 1, 20))\n",
    "\n",
    "# Visualize the configurations!\n",
    "\n",
    "xs, ys = np.meshgrid(range(N), range(1))\n",
    "fig = plt.figure(figsize=(5, 2))\n",
    "plt.xlim(0, N-1)\n",
    "plt.ylim(0, 1)\n",
    "plt.xticks(range(0, N, 2))\n",
    "plt.xlabel('index')\n",
    "plt.title('1D Ising model evolution (T = 5)')\n",
    "\n",
    "writer = animation.writers['ffmpeg']\n",
    "writer = writer(fps=20, metadata=dict(artist='Me'))\n",
    "\n",
    "p = plt.pcolormesh(configs[0], cmap=plt.get_cmap('binary'))\n",
    "plt.yticks([])\n",
    "\n",
    "def animate(i):\n",
    "    p.set_array(configs[i].ravel())\n",
    "    return p\n",
    "\n",
    "\n",
    "anim_high = animation.FuncAnimation(\n",
    "    fig, animate, frames=range(0, 10000, 10), interval=10)\n",
    "#anim_high.save('T=5.mp4', writer=writer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.b. Display the animation (high T)\n",
    "\n",
    "The code above will create and save (if you uncomment the bottom line) an animation of the evolution of the system. Execute the code block below to display the animation in the Jupyter notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "HTML(anim_high.to_html5_video())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.c. The Ising model at low temperature\n",
    "\n",
    "Next, let's run a simulation at a temperature of $T=1$. As before we'll take $10^4$ Monte Carlo steps. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import HTML\n",
    "import matplotlib.animation as animation\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "%matplotlib inline\n",
    "\n",
    "# Gather the simulation data\n",
    "\n",
    "configs = np.array(MC_simulate(1, 10000)).reshape((10000, 1, 20))\n",
    "\n",
    "# Visualize the configurations!\n",
    "\n",
    "xs, ys = np.meshgrid(range(N), range(1))\n",
    "fig = plt.figure(figsize=(5, 2))\n",
    "plt.xlim(0, N-1)\n",
    "plt.ylim(0, 1)\n",
    "plt.xticks(range(0, N, 2))\n",
    "plt.xlabel('index')\n",
    "plt.title('1D Ising model evolution (T = 1)')\n",
    "\n",
    "writer = animation.writers['ffmpeg']\n",
    "writer = writer(fps=20, metadata=dict(artist='Me'))\n",
    "\n",
    "p = plt.pcolormesh(configs[0], cmap=plt.get_cmap('binary'))\n",
    "plt.yticks([])\n",
    "\n",
    "def animate(i):\n",
    "    p.set_array(configs[i].ravel())\n",
    "    return p\n",
    "\n",
    "\n",
    "anim_low = animation.FuncAnimation(\n",
    "    fig, animate, frames=range(0, 10000, 10), interval=10)\n",
    "#anim_low.save('T=1.mp4', writer=writer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.d. Display the animation (low T)\n",
    "\n",
    "Execute the code block below to display the evolution of the 1D Ising model at low temperature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HTML(anim_low.to_html5_video())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.e. Assessing the results\n",
    "\n",
    "Compare your simulations at high and low temperatures. What differences do you observe? Can you explain the results?\n",
    "\n",
    "We won't have time to discuss this in the class, but the **2D** Ising model has a [phase transition](https://en.wikipedia.org/wiki/Phase_transition) in between these two temperatures, which occurs at a critical temperature\n",
    "\n",
    "$$ T_c = \\frac{2}{\\log\\left(1+\\sqrt{2}\\right)} \\approx 2.27\\,. $$\n",
    "\n",
    "The simulation at high temperature is above $T_c$, in the disordered phase, and the simulation at low temperature is in the ordered phase below $T_c$. In **1D** there is no phase transition, but you should observe a clear difference in behavior!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "G4FFfTBSY-aJ"
   },
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "name": "homework-6-solution.ipynb",
   "provenance": [],
   "version": "0.3.2"
  },
  "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.7.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}