ising_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this lab we will look at modeling an Ising system with Monte Carlo.\n",
    "If you are unfamiliar with Python, a set of tutorials can be found [here](https://www.learnpython.org/en/Welcome).\n",
    "From the tutorials linked above listed under Python basics, you are to be familiar with:\n",
    "- Variables & Types\n",
    "- Lists\n",
    "- Basic Operators\n",
    "- Conditions\n",
    "- Loops\n",
    "- Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ising Model \n",
    "\n",
    "The Ising model is a lattice of $k$ interacting sites. Each site has a spin $\\sigma_k$ which can have a value of +1 or -1. The spins are allowed to interact with their nearest-neighbor spins. The Hamiltonian for the system is given by:\n",
    "\n",
    "$$\n",
    "H(\\sigma)=-\\sum_{ a }^{k} \\left(\\sum_{b}^{\\langle a\\rangle} J \\sigma_{a} \\sigma_{b}\\right)- h \\sigma_{a}\n",
    "$$\n",
    "\n",
    "where $\\langle a\\rangle$ indicate a sum over only the nearest neighbors of site $a$. The interaction between spins is captured in $J$. If $J$ is positive, the spins will want to align with their neighbors or system will behave ferromagnetically. If $J$ is negative, the spins will want to antialign with their neighbors or system will behave antiferromagnetically. A magnetic field aligned with the spins is represented by $h$.\n",
    "\n",
    "We will start by introducing some code for this model. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To begin, we need to import the Python modules we will need. If you're not familiar with these, here is an overview: \n",
    "- `numpy` allows us to perform the necessary math\n",
    "- `matplotlib` will allows us to visualize our results\n",
    "- `pandas` will help up tabulate data cleanly\n",
    "- `pyblock` handles reblocking, which will be discussed later in the notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cell below will ensure you have all necessary packages installed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install numpy matplotlib pyblock pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pyblock\n",
    "import pandas as pd\n",
    "%matplotlib inline\n",
    "from IPython.display import set_matplotlib_formats\n",
    "set_matplotlib_formats('png', 'pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll also import some helper utilities that we wrote to make things easier for you"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ising_utils as iu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ising Model Parameters\n",
    "We will store the information about a $M$ by $N$ 2D lattice of spins using a dictionary like below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ip = {\"J\": None, \"h\": None, \"M\": None, \"N\": None}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Initializing a lattice\n",
    "The function below takes a dictionary of Ising model parameters `ip` and initializes a lattice to be all spin up.\n",
    "\n",
    "The lattice state will be stored as a $M$ by $N$ array of 1 and -1 which will represent spin up and down respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initialize_lattice_state(ip):\n",
    "    \"\"\"\n",
    "    Initializing the lattice of size MxN to all ones\n",
    "    \"\"\"\n",
    "    lattice_state = np.ones((ip[\"M\"], ip[\"N\"]))\n",
    "    return lattice_state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code below initializes an example $3\\times3$ lattice with $J = 1, h = 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initializing a model of 3x3 spins, with J = 1 and h = 0\n",
    "ip = {\"J\": 1, \"h\": 0 , \"M\": 3, \"N\": 3}\n",
    "example_lattice = initialize_lattice_state(ip) \n",
    "# visualize the lattice\n",
    "iu.plot_lattice(example_lattice)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1: Initialize an example lattice\n",
    "\n",
    "Use the function above to initialize an example $5\\times8$ lattice with $J = 2, h = 0.01$.\n",
    "\n",
    "1. Create a dictionary of Ising model parameters named `ip`\n",
    "2. Call the function `initialize_lattice_state` above with `ip` as an argument. Store the result in a variable named `lattice`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# place code to initialize here\n",
    "ip = # YOUR CODE HERE\n",
    "lattice = # YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Don't modify this, it will check your answer above\n",
    "iu.check_exercise_1(ip, lattice)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Coding functions for the model\n",
    "Next, you will be coding your first functions for the model. \n",
    "\n",
    "The functions we will need to create are: \n",
    "- the ability to flip a specific spin within the lattice\n",
    "- the ability to go through each lattice site to calculate the energy of the lattice.\n",
    "- a function to calculate the energy per site.\n",
    "\n",
    "\n",
    "\n",
    "```python\n",
    " def function_name(ARGUMENTS HERE):\n",
    "        CODE HERE\n",
    " ```\n",
    " \n",
    " \n",
    " We will first start with the function to flip (i.e. at a specific site in the lattice, $ i $, $ j $, the spin will flip and that lattice site will be the negative of itself.)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def flip_spin(lattice_state, i, j):\n",
    "    \"\"\"\n",
    "    Flip spin (i, j)\n",
    "    i.e. -1 ---> 1\n",
    "          1 ---> -1\n",
    "    \"\"\"\n",
    "    # Make a copy of the lattice state to modify\n",
    "    new_lattice_state = np.copy(lattice_state)\n",
    "    \n",
    "    # YOUR CODE HERE\n",
    "    return new_lattice_state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Excercise 2\n",
    "As you did above initaite a new 3x3 lattice, with J = 1 and h = 0. \n",
    "Place the lattice in a variable called `example_lattice`.\n",
    "Use the `flip_spin` function you wrote above twice to flip the spin at (0,1) and then at (2,2).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ip = # YOUR CODE HERE\n",
    "example_lattice = # YOUR CODE HERE\n",
    "# YOUR CODE HERE\n",
    "iu.plot_lattice(example_lattice)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Don't modify this, it will check your answer above\n",
    "iu.check_exercise_2(ip, example_lattice)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating the contributions to a lattice sites\n",
    "\n",
    "The next function we want to write is the function to calculate the energy contributions of a specific site.\n",
    "The total energy will be the sum of this function over all sites and correcting for double counting.\n",
    "\n",
    "Our model will implement periodic boundary conditions (PBC) to approximate a large infinite system with a small unit cell.\n",
    "Breifly stated, a single lattice state will also feel the affect of having a repeated lattice state to the top, bottom, left, and right.\n",
    "We will build our understanding of PBC using the diagram below. \n",
    "\n",
    "<img src=\"pbc.png\" width=\"60%\">\n",
    "\n",
    "Here, we have described a unit cell as a 3 x 3 lattice, which we visualize as the center block of this image (the one containing the red dot).\n",
    "Each dot represents a lattice site which either has a +/-1 spin.\n",
    "Let's explore the interactions felt by the center unit cell above, specifically the red dot.\n",
    "Where the red dot in the individual unit cell will interact with other sites within the unit cell such as the yellow sites.\n",
    "Additionally, to model the bulk system, it will also feel the effect of sites outside of the reference cell, indicated by the blue dots.\n",
    "The situations where a site is feeling contributions from other unit cells will happen on all edges. \n",
    "As such, they are colloquially reffered to as edge cases.\n",
    "Conveniently, instead of storing multiple versions of our lattice, we can instead realize that the blue points are also found in the unit cell and we can thus stick to referencing a single unit cell.\n",
    "This will then require some clever coding on our part.\n",
    "\n",
    "\n",
    "# Exercise: calculate the contributions of a site\n",
    "\n",
    "We will first focus on writing the energy contributions of the single site.\n",
    "The function below will return those contributions.\n",
    "This will require knowledge of the spins on all neighbors of the central site $i$, $j$.\n",
    "\n",
    "There will be boundary cases to consider:\n",
    "   - If you are along the top edge of the lattice, $j == 0$.\n",
    "   - If you are along the bottom edge of the lattice, $j == N-1$.\n",
    "   - If you are along the left edge of the lattice, $i == 0$.\n",
    "   - If you are along the right edge of the lattice, $i == M-1$.\n",
    "\n",
    "Recall that python indexing starts from 0 which is why $M-1$ and $N-1$ describe the right and bottom edge respectively rather than $M$ and $N$.\n",
    "\n",
    "Once we have the values of the neighboring spins, the interactions of a specific lattice site with the external field and the neighboring sites is expressed as:\n",
    "\n",
    "$$E = -h\\sigma_{ij} - J\n",
    "\\sigma_{ij}(\\sigma_{i+1,j}+\\sigma_{i-1,j}+\\sigma_{i,j+1}+\\sigma_{i,j-1})$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_energy_contrib_of_site(lattice_state, ip, i, j):\n",
    "    \"\"\"\n",
    "    Calculate energy of spin (i, j)\n",
    "\n",
    "    Periodic boundary conditions implemented\n",
    "    \"\"\"\n",
    "    spin_here = lattice_state[i, j]  # value of spin here\n",
    "\n",
    "    # value of spin above, below, left, and right of spin (i, j)\n",
    "    # for each, if on boundary, we wrap around to the other side\n",
    "    # of the lattice for periodic boundary conditions\n",
    "    \n",
    "    # create variables spin_left spin_right spin_above and spin_below using PBC\n",
    "    \n",
    "    # YOUR CODE HERE\n",
    "    \n",
    "    return (-ip[\"h\"] * spin_here) - (ip[\"J\"] * spin_here * (\n",
    "        spin_above + spin_below + spin_left + spin_right))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Excercise 3\n",
    "The code in the cell below will test the `calculate_energy_contrib_of_site` function you wrote above. Pay attention to the output to see if you made any errors in `calculate_energy_contrib_of_site`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Don't modify this, it will check your answer above\n",
    "iu.check_exercise_3(initialize_lattice_state,flip_spin, calculate_energy_contrib_of_site)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Additional functions for observables of the model\n",
    "\n",
    "However, looping over all sites and calculating the energy contributions of a site suffers from double counting of the sites.\n",
    "You can think about this as the site $\\sigma_{i,j}$ will show up in this expression but also for the the iteration assessing $\\sigma_{i+1,j}$. Thus, we reduce the contribution of J by a factor of 2 when calculating the total energy. Thus the final expression for the energy of a site is: \n",
    "\n",
    "$$E = - \\frac{-2h\\sigma_{ij} - J\n",
    "\\sigma_{ij}(\\sigma_{i+1,j}+\\sigma_{i-1,j}+\\sigma_{i,j+1}+\\sigma_{i,j-1})}{2}$$\n",
    "\n",
    "This lets us loop over all sites and calculate the total energy as a sum over the site energies.\n",
    "\n",
    "We have provided a set of functions to calculate other obervables.\n",
    "These calculate the instantaneous value of the observable, and we will be interested in the ensemble average of these values.\n",
    "\n",
    "- The total spin $s$\n",
    "- The total spin per particle $\\frac{s}{M\\times N}$\n",
    "- The total energy $E$\n",
    "- The total energy per particle $\\frac{E}{M \\times N}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_total_spin(lattice_state):\n",
    "    \"\"\"\n",
    "    Calculate total spin of the lattice\n",
    "    \"\"\"\n",
    "    return np.sum(lattice_state)\n",
    "\n",
    "\n",
    "def calculate_total_spin_per_spin(lattice_state, ip):\n",
    "    \"\"\"\n",
    "    Calculate total spin of the lattice\n",
    "    \"\"\"\n",
    "    return calculate_total_spin(lattice_state) / (ip[\"M\"] * ip[\"N\"])\n",
    "\n",
    "\n",
    "def calculate_total_energy(lattice_state, ip):\n",
    "    \"\"\"\n",
    "    Calculate total energy of the lattice\n",
    "    \"\"\"\n",
    "    E = 0.0\n",
    "    for i in range(ip[\"M\"]):\n",
    "        for j in range(ip[\"N\"]):\n",
    "            E += calculate_energy_contrib_of_site(lattice_state, ip, i, j)\n",
    "    # This sum double counts spin-spin interactions so add the interaction \n",
    "    # with the field again and divide by 2\n",
    "    return (E - ip['h']*calculate_total_spin(lattice_state))/2\n",
    "\n",
    "\n",
    "def calculate_total_energy_per_spin(lattice_state, ip):\n",
    "    \"\"\"\n",
    "    Calculate energy of lattice normalized by the number of spins\n",
    "    \"\"\"\n",
    "    return calculate_total_energy(lattice_state, ip) / (ip[\"M\"] * ip[\"N\"])\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monte Carlo Simulations\n",
    "\n",
    "Now that we have a model for the ising lattice, we are ready to calculate the observables of interest. \n",
    "An observable can be calculated as a weighted sum over states \n",
    "$$\n",
    "\\langle O\\rangle=\\sum_{\\alpha} O(\\alpha) P(\\alpha)\n",
    "$$\n",
    "For example, for the expectation value of the energy of an $M \\times N$ lattice, the expectation value would be\n",
    "\n",
    "$$\n",
    "\\langle E\\rangle=\\sum_{\\alpha} E(\\alpha) P(\\alpha)\n",
    "$$\n",
    "\n",
    "where the probability distribution $P(\\alpha)$ is\n",
    "\n",
    "$$\n",
    "P(\\alpha) = \\frac{e^{-\\beta E_{\\alpha}}}{\\sum_{i} e^{-\\beta E_{i}}}\n",
    "$$\n",
    "\n",
    "where $E(\\alpha)$ is the energy of a fixed state $\\alpha$, and  $P(\\alpha)$ is the probability of being in that fixed state.\n",
    "It is worth noting that this approach is true for any other property as well, not just the energy.\n",
    "However, the number of fixed states grows as $2^{(N*M)}$ where $N*M$ is the total number of lattice points.\n",
    "Exactly calculating the contributions for every possible lattice state quickly becomes impractical as the lattice size grows.\n",
    "Additionally we don't know the probability of being in each state since we cannot compute the denominator.\n",
    "\n",
    "One way to deal with this is Monte Carlo sampling, which can be used to estimate the probability distribution by using a random distribution of $n$ points in the configuation space (i.e. the space of all possible lattices): \n",
    "$$\n",
    "\\langle E\\rangle = \\frac{\\sum_{n} E(n)e^{-\\beta E_{n}}}{\\sum_{n} e^{-\\beta E_{n}}}\n",
    "$$\n",
    "However, by simply choosng $n$ random lattices from the configuration space, there is no guarantee that we have sampled enough points that are representative of $P(\\alpha)$.\n",
    "\n",
    "To address this, we sample lattices drawn from the probability distribution $P(\\alpha)$.\n",
    "This may seem weird since we just said that we do not know the distribution $P(\\alpha)$, however using the Metropolis algorithm we can sample from an unknown distribution.\n",
    "\n",
    "To understand this consider starting at a certain lattice configuration $X$. The probability of being at state $X$ is:\n",
    "$$\n",
    "P(X) = \\frac{e^{-\\beta E_{X}}}{\\sum_{i} e^{-\\beta E_{i}}}\n",
    "$$\n",
    "If we generate a new configuration $Y$, the probability of being at lattice $Y$ is given by P($Y$):\n",
    "$$\n",
    "P(Y) = \\frac{e^{-\\beta E_{Y}}}{\\sum_{i} e^{-\\beta E_{i}}}\n",
    "$$\n",
    "The ratio of these two probabilities is:\n",
    "$$\n",
    "P(Y \\leftarrow X) = e^{-\\beta (E_{Y} - E_{X})}\n",
    "$$\n",
    "Notice that the denominators now drop out! Therefore if we generate lattice configurations and accept new configurations according to this ratio then the generated configurations will be distributed according to $P(\\alpha)$. One issue is that this ratio can be greater than 1 so it is no longer a probability distribution function. \n",
    "We therefore take the acceptance probability of a new move to be:\n",
    "$$\n",
    "P_{\\mathrm{accept}} = min\\left(1,\\frac{P(Y)}{P(X)}\\right)\n",
    "$$\n",
    "where $min$ means we take the minimum of the two values.\n",
    "\n",
    "This also indicates the need for an equilibration period.\n",
    "If the first sample $X$ is not from a region of high probability, samples must be drawn from $P(\\alpha)$ until some steady state is reached.\n",
    "\n",
    "The algorithm is now:\n",
    "- Generate an initial configuration\n",
    "- equilibrate\n",
    "- for N samples\n",
    "    - generate a trial configuration\n",
    "    - calculate the metropolis acceptance probability $e^{-\\beta (E_{Y} - E_{X})}$\n",
    "    - draw a uniform random number from [0,1]\n",
    "    - if the random number <= P_accept\n",
    "        - update the initial configuration to be the current configuration\n",
    "    - else\n",
    "        - stay at the current configuration (i.e. don't update the initial configuration)\n",
    " \n",
    "Notice that this means if a trial configuration is lower in energy it is always accepted. If a trial configuration is higher in energy than it is accepted with probability $P_{accept}$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulation parameters and simulation results will be stored in dictionaries. The simulation parameters are the number of equilibration sweeps (`num_equil_sweeps`), the number of production sweeps(`num_sample_sweeps`), and kT (`kT`).\n",
    "\n",
    "During the simulation we want to record the observables that we would like ensemble averages of: the energy per site `energy_list` and the total spin per site `total_spin_list`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# simulation parameters\n",
    "sp = {\"num_equil_sweeps\": None, \"num_sweeps\": None, \"kT\": None}\n",
    "# simulation results:\n",
    "sr = {\"energy_list\": [], \"total_spin_list\": []}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4 Metropolis test\n",
    "Above the probability to accept a new configuration was given by $P_{\\mathrm{accept}}$.\n",
    "This function will return a boolean value (True or False) depending on if the new configuration is accepted.\n",
    "This function will receive the energy of the current configuration (`E_old`) and the trial configuration (`E_new`) and the simulation parameters (`sp`).\n",
    "The code to generate a random number sampled from a uniform distribution between 0 and 1 in python using `numpy` is `np.random.random()` and the exponential function is `np.exp`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def metropolis_test(sp, E_new, E_old):\n",
    "    # YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Don't modify this, it will check your answer above\n",
    "iu.check_exercise_4(metropolis_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When generating new configurations, we could do a number of things.\n",
    "For example we could pick a spin at random and flip it and test if we accept it. \n",
    "Or we could sweep through the lattice and flip each spin and test if we accept it.\n",
    "Here we choose the latter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sweep(lattice_state, ip, sp):\n",
    "    for site_i in range(ip[\"M\"]):\n",
    "        for site_j in range(ip[\"N\"]):\n",
    "            E_old = calculate_energy_contrib_of_site(lattice_state, ip, site_i, site_j)\n",
    "            # flip spin i and j\n",
    "            new_lattice_state = flip_spin(lattice_state, site_i, site_j)\n",
    "            # calculate updated energy\n",
    "            E_new = calculate_energy_contrib_of_site(new_lattice_state, ip, site_i, site_j)\n",
    "            # Monte Carlo step\n",
    "            if metropolis_test(sp, E_new, E_old):\n",
    "                # accept move\n",
    "                E_old = E_new\n",
    "                lattice_state = np.copy(new_lattice_state)\n",
    "            else:\n",
    "                # reject move\n",
    "                # do nothing just continue\n",
    "                continue\n",
    "    return lattice_state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These functions calculate some of the statistics for you. \n",
    "The mean and standard error of the mean of the energy per site and spin per site."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_average_energy(sr):\n",
    "    return np.mean(sr[\"energy_list\"])\n",
    "\n",
    "\n",
    "def get_stderr_energy(sr, sp):\n",
    "    return np.std(sr[\"energy_list\"]) / np.sqrt(sp[\"num_sample_sweeps\"])\n",
    "\n",
    "\n",
    "def get_average_spin(sr):\n",
    "    return np.mean(sr[\"total_spin_list\"])\n",
    "\n",
    "\n",
    "def get_stderr_spin(sr, sp):\n",
    "    return np.std(sr[\"total_spin_list\"]) / np.sqrt(sp[\"num_sample_sweeps\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, these expression are serially correlated. For example, the energy for one sweep is not statistically independent from previous or following sweeps. To remove this correlation, one can average over short blocks of data, a process known as reblocking. This removes serial correlation because each block is now statistically independent. The resulting averages are equivalent when averaging over all the data or blocks of data, but the standard error will differ.\n",
    "\n",
    "The reblocking implementation of mean and standard error function are included in the ising_utils as:\n",
    "\n",
    "`iu.get_reblocked_avg_stderr_spin`\n",
    "\n",
    "`iu.get_reblocked_avg_stderr_energy`\n",
    "\n",
    "You will see these functions being called later in the simulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next function will run a simulation.\n",
    "The arguments required are `ip` for the Ising model parameters and `sp` for the simulation parameters.\n",
    "\n",
    "This function will return a dictionary, `sr` containing the simulation results.\n",
    "The two statistics to collect per simulation step are the energy of the lattice (`sr['energy_list']`) and the total sping of the lattice (`sr['total_spin_list']`)\n",
    "\n",
    "The simulation takes the following steps:\n",
    "1. initialize the lattice state\n",
    "2. A number of equilibration steps are done in which statistics are not collected.\n",
    "3. A number of production steps are done.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulation(ip, sp):\n",
    "    sr = {\"energy_list\": [], \"total_spin_list\": []}\n",
    "    lattice_state = initialize_lattice_state(ip)\n",
    "    for step in range(sp[\"num_equil_sweeps\"]):\n",
    "        lattice_state = sweep(lattice_state, ip, sp)\n",
    "    for step in range(sp[\"num_sample_sweeps\"]):\n",
    "        lattice_state = sweep(lattice_state, ip, sp)\n",
    "        # calculate total energy\n",
    "        e = calculate_total_energy_per_spin(lattice_state, ip)\n",
    "        sr[\"energy_list\"].append(e)\n",
    "\n",
    "        # calculate total s2\n",
    "        s = calculate_total_spin_per_spin(lattice_state, ip)\n",
    "        sr[\"total_spin_list\"].append(s)\n",
    "    return sr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Applications"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will study the ising model under two sets of parameters. You will be tasked with making some observations about the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will be comparing our simulation to that of the analytical result. The analytical result is defined as:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analytical(x, J):\n",
    "    analytical_solution = []\n",
    "    for i in x:\n",
    "        if i < 2.269:\n",
    "            analytical_solution.append(\n",
    "                (1 - np.sinh((2 * J) / i) ** (-4)) ** (1.0 / 8.0)\n",
    "            )\n",
    "        else:\n",
    "            analytical_solution.append(0)\n",
    "\n",
    "    return np.array(analytical_solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application 1\n",
    "We want to simulate a system ($M$ = 30, $N$ = 30) with $J = 1$ and $h = 0$.\n",
    "We want to find the transition temperature, so we will run over a list of kT values around the transition temperature. \n",
    "We will look at the kT range between 1.2 and 3.0.  \n",
    "The ising model parameters are defined first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# bad sampling bad lattice for quick running will remove\n",
    "ip = {\"N\": 10, \"M\": 10, \"J\": 1, \"h\": 0.0}\n",
    "num_equil_sweeps = 100\n",
    "num_sample_sweeps = 1500\n",
    "# good sampling good lattice\n",
    "# ip = {\"N\": 30, \"M\": 30, \"J\": 1, \"h\": 0.0}\n",
    "# num_equil_sweeps = 1000\n",
    "# num_sample_sweeps = 15000\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we calculate the analytical solution to this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "analytical_kT_list = np.arange(1.2, 3.0, 0.001)\n",
    "s_analytical = analytical(analytical_kT_list, ip[\"J\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we set up list of kT values to simulate, the data structures for storing observables:\n",
    "- `avg_s_list` : average spin per site \n",
    "- `stderr_s_list` : standard error of the mean of s per site\n",
    "- our simulation parameters are defined for each kT, but the `num_equil_sweeps` and `num_sample_sweeps` will be the same for each kT so we define them here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "simulation_kT_list = np.arange(1.2, 3.0, 0.1)\n",
    "\n",
    "avg_s_list = []\n",
    "stderr_s_list = []\n",
    "\n",
    "for kT in simulation_kT_list:\n",
    "    sp = {\"num_equil_sweeps\": num_equil_sweeps, \"num_sample_sweeps\": num_sample_sweeps, \"kT\": kT}\n",
    "\n",
    "    sr = run_simulation(ip, sp)\n",
    "    mean, std_err = iu.get_reblocked_avg_stderr_spin(sr, sp)\n",
    "    print(\"kT={:>4.2f}    mean {:>4.2f}    std_err {:>4.2f}\".format(kT, mean, std_err))\n",
    "\n",
    "    avg_s_list.append(mean)\n",
    "    stderr_s_list.append(std_err)\n",
    "    fig, ax = plt.subplots(tight_layout=True)\n",
    "    fig.set_size_inches(8, 5)\n",
    "    n_bins = 25\n",
    "    bins = np.linspace(-1, 1, n_bins)\n",
    "    ax.hist(sr[\"total_spin_list\"], bins=bins)\n",
    "    ax.set_title(\n",
    "        \"kT={:>4.2f} m={:>02d}, n={:>02d}, J={:>4.2f}, h={:>4.2f} samp={}\".format(\n",
    "            sp[\"kT\"], ip[\"M\"], ip[\"N\"], ip[\"J\"], ip[\"h\"], sp[\"num_sample_sweeps\"]\n",
    "        )\n",
    "    )\n",
    "    ax.set_xlim(-1.1, 1.1)\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(\n",
    "        \"metropolis_s_hist_kT_{:4.2f}_n_{:>02d}_m_{:>02d}_J_{:>4.2f}_h_{:>4.2f}_samp_{}.png\".format(\n",
    "            sp[\"kT\"], ip[\"M\"], ip[\"N\"], ip[\"J\"], ip[\"h\"], sp[\"num_sample_sweeps\"]\n",
    "        ),\n",
    "        dpi=200,\n",
    "    )\n",
    "    # plt.show()\n",
    "#     s = pd.Series(sr[\"total_spin_list\"])\n",
    "#     print(\"avg s autocorr = {}\".format(s.autocorr()))\n",
    "plt.figure(figsize=(8, 5))\n",
    "plt.errorbar(simulation_kT_list, avg_s_list, stderr_s_list, label=\"avg S\",color='darkorchid')\n",
    "plt.plot(analytical_kT_list, s_analytical, label=\"analytical\",color='dodgerblue')\n",
    "plt.plot(analytical_kT_list, -s_analytical,color='dodgerblue')\n",
    "plt.title(\n",
    "    \"m={:>02d}, n={:>02d}, J={:>4.2f}, h={:>4.2f} samp={}\".format(\n",
    "        ip[\"M\"], ip[\"N\"], ip[\"J\"], ip[\"h\"], sp[\"num_sample_sweeps\"]\n",
    "    )\n",
    ")\n",
    "\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.savefig(\n",
    "    \"metropolis_avg_s_m_{:>02d}_n_{:>02d}_J_{:>4.2f}_h_{:>4.2f}_samp_{}.png\".format(\n",
    "        ip[\"M\"], ip[\"N\"], ip[\"J\"], ip[\"h\"], sp[\"num_sample_sweeps\"]\n",
    "    ),\n",
    "    dpi=200,\n",
    ")\n",
    "# plt.show()\n",
    "\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finite field (h=0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application 2\n",
    "Consider the previous Ising system ($M$ = 30, $N$ = 30) with $J = 1$ and now an applied field $h = 0.01$.\n",
    "This should bias the system to align with the applied field\n",
    "\n",
    "We will look at the kT range between 1.2 and 3.0.  \n",
    "The ising model parameters are defined first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# bad sampling bad lattice for quick running will remove\n",
    "ip = {\"N\": 10, \"M\": 10, \"J\": 1, \"h\": 0.01}\n",
    "num_equil_sweeps = 100\n",
    "num_sample_sweeps = 1500\n",
    "# good sampling good lattice\n",
    "# ip = {\"N\": 30, \"M\": 30, \"J\": 1, \"h\": 0.01}\n",
    "# num_equil_sweeps = 1000\n",
    "# num_sample_sweeps = 15000\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the simulation as we did above. Note we still plot the analytical results, but those are for zero applied field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "simulation_kT_list = np.arange(1.2, 3.0, 0.1)\n",
    "\n",
    "avg_s_list = []\n",
    "stderr_s_list = []\n",
    "\n",
    "for kT in simulation_kT_list:\n",
    "    sp = {\"num_equil_sweeps\": num_equil_sweeps, \"num_sample_sweeps\": num_sample_sweeps, \"kT\": kT}\n",
    "\n",
    "    sr = run_simulation(ip, sp)\n",
    "    mean, std_err = iu.get_reblocked_avg_stderr_spin(sr, sp)\n",
    "    print(\"kT={:>4.2f}    mean {:>4.2f}    std_err {:>4.2f}\".format(kT, mean, std_err))\n",
    "\n",
    "    avg_s_list.append(mean)\n",
    "    stderr_s_list.append(std_err)\n",
    "    fig, ax = plt.subplots(tight_layout=True)\n",
    "    fig.set_size_inches(8, 5)\n",
    "    n_bins = 25\n",
    "    bins = np.linspace(-1, 1, n_bins)\n",
    "    ax.hist(sr[\"total_spin_list\"], bins=bins)\n",
    "    ax.set_title(\n",
    "        \"kT={:>4.2f} m={:>02d}, n={:>02d}, J={:>4.2f}, h={:>4.2f} samp={}\".format(\n",
    "            sp[\"kT\"], ip[\"M\"], ip[\"N\"], ip[\"J\"], ip[\"h\"], sp[\"num_sample_sweeps\"]\n",
    "        )\n",
    "    )\n",
    "    ax.set_xlim(-1.1, 1.1)\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(\n",
    "        \"metropolis_s_hist_kT_{:4.2f}_n_{:>02d}_m_{:>02d}_J_{:>4.2f}_h_{:>4.2f}_samp_{}.png\".format(\n",
    "            sp[\"kT\"], ip[\"M\"], ip[\"N\"], ip[\"J\"], ip[\"h\"], sp[\"num_sample_sweeps\"]\n",
    "        ),\n",
    "        dpi=200,\n",
    "    )\n",
    "    # plt.show()\n",
    "#     s = pd.Series(sr[\"total_spin_list\"])\n",
    "#     print(\"avg s autocorr = {}\".format(s.autocorr()))\n",
    "plt.figure(figsize=(8, 5))\n",
    "plt.errorbar(simulation_kT_list, avg_s_list, stderr_s_list, label=\"avg S\",color='darkorchid')\n",
    "plt.plot(analytical_kT_list, s_analytical, label=\"analytical\",color='dodgerblue')\n",
    "plt.plot(analytical_kT_list, -s_analytical,color='dodgerblue')\n",
    "plt.title(\n",
    "    \"m={:>02d}, n={:>02d}, J={:>4.2f}, h={:>4.2f} samp={}\".format(\n",
    "        ip[\"M\"], ip[\"N\"], ip[\"J\"], ip[\"h\"], sp[\"num_sample_sweeps\"]\n",
    "    )\n",
    ")\n",
    "\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.savefig(\n",
    "    \"metropolis_avg_s_m_{:>02d}_n_{:>02d}_J_{:>4.2f}_h_{:>4.2f}_samp_{}.png\".format(\n",
    "        ip[\"M\"], ip[\"N\"], ip[\"J\"], ip[\"h\"], sp[\"num_sample_sweeps\"]\n",
    "    ),\n",
    "    dpi=200,\n",
    ")\n",
    "# plt.show()\n",
    "\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}