Project_4.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Project 4",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kw62klmxOLuu",
        "colab_type": "text"
      },
      "source": [
        "**Name**: Brandon Loptman\n",
        "\n",
        "**Total Score**: 20/20\n",
        "\n",
        "Code Tests:\n",
        "\n",
        "Completed Tests: 3/3 \n",
        "\n",
        "Tests Explained: 2/2\n",
        "\n",
        "Application to Physical Systems Modeling:\n",
        "8/8 \n",
        "\n",
        "General:\n",
        "\n",
        "Run without problems: 1/1\n",
        "\n",
        "Correct structure: 1/1\n",
        "\n",
        "Logical flow: 1/1\n",
        "\n",
        "Good graphical representaions: 3/3 \n",
        "\n",
        "Graphs have axes, labels, ticks, titles, etc.: 1/1 \n",
        "\n",
        "**Comments**:\n",
        "\n",
        "**Description**:\n",
        "\n",
        "Code Tests\n",
        "\n",
        "(3 pts) If code tests are suggested, all such tests are completed.  If code tests aren’t suggested, student has articulated their own code tests and executed those tests to show that code is working properly in at least some small number of obvious cases.  The more robust the code tests are, the better.\n",
        "\n",
        "(2 pts) All code tests are explained in complete sentences in text cells (along with any relevant mathematics).  There is a reasonable amount of explanation of how one knows that the code tests have been passed.\n",
        "\n",
        "Application of Code to Physical Systems Modeling\n",
        "\n",
        "(8 pts) Notebook contains a mixture of explanation text in complete sentences (with math where appropriate) and graphs/other visual representations illustrating how the computer model exhibits the relevant phenomena predicted by the mathematical model of the system.  All prompts have been answered.\n",
        "\n",
        "General\n",
        "\n",
        "(1 pts) The whole notebook can be re-started and re-run without halting due to any coding errors.\n",
        "\n",
        "(1 pts) The notebook has the correct general structure (code tests placed in the beginning and application of code to analysis of relevant physical systems after).\n",
        "\n",
        "(1 pts) Flow of notebook is logical, clear, and concise.\n",
        "\n",
        "(3 pts) Graphical/visual representations are used to concisely illustrate results of code tests and simulations.\n",
        "\n",
        "(1 pts) All graphs with axes have labels on axes with units and tick marks where appropriate.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mRi4aS5w-BrU",
        "colab_type": "text"
      },
      "source": [
        "# Project 4 (Path Integral Monte Carlo)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hFno37Bu9bFw",
        "colab_type": "text"
      },
      "source": [
        "##Abstract"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "a--_NfCw80gb",
        "colab_type": "text"
      },
      "source": [
        "Our project was about numerically approximating ensemble averages of a system using Markov Chain-Monte Carlo (MCMC) simulations. In particular, we applied the Feynmann path integral formulation of quantum mechanics to the quantum harmonic oscillator to determine the ensemble average energy $\\langle \\hat{H} \\rangle$ using the Metropolis algorithm. For our MCMC simulations we used a paths consisting of five points over a range of $\\beta \\in [.01,10]$. The values of $\\langle \\hat{H} \\rangle$ obtained from our MCMC simulations agreed very strongly with the analytical result, differing on average by only 2.89%."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "R8OVQ0_M9f4X",
        "colab_type": "text"
      },
      "source": [
        "## Physical System and Model"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ulu38Qp4Jz_o",
        "colab_type": "text"
      },
      "source": [
        "Equilibrium statistical mechanics is based upon the idea of a partition function which contains all of the essential information about the system under consideration. The partition function $Z$ is the sum over all states of the Boltzmann factor for each state. Recall, that the Boltzmann factor for a state $s$ is given by $e^{-E_{s}/k_{b}T}$ where $E_{s}$ is the energy of state $s$, $k_{b}$ is the Boltzmann's constant, and $T$ is the temperature in Kelvin. This is often written as $e^{E_{s}\\beta}$, where $\\beta = \\tfrac{1}{k_{b}T}$. Thus, the partition function is:\n",
        "\n",
        "$$\n",
        "  Z = \\sum_{\\mbox{all states}} e^{-E_{s}\\beta}\n",
        "$$\n",
        "\n",
        "The above sum is over all states of the system and thus depends on size of the system and the number of degrees of freedom for each particle. For systems consisting of only a few interacting particles the partition function can be written down exactly, thus the properties of the system can be calculated in a closed form. \n",
        "\n",
        "A classic example of a system for which the partition function can be computed exactly is the quantum harmonic oscillator. Recall, that the energy of the quantum harmonic oscillator is given by\n",
        "\n",
        "$$\n",
        "  E_{n} = (\\tfrac{1}{2} + n) \\hbar \\omega\n",
        "$$\n",
        "\n",
        "where $\\omega$ is the angular frequency of the oscillator and \n",
        "$n \\in (0,1,2,\\dots)$. Plugging this into the equation for the parition function we get\n",
        "\n",
        "$$\n",
        "  Z = \\sum_{n=0}^{\\infty} e^{-E_{n}\\beta} = e^{-\\tfrac{1}{2} \\beta \\hbar \\omega} \\sum_{n=0}^{\\infty} e^{-n \\beta\\hbar\\omega}\n",
        "$$\n",
        "\n",
        "The sum over n above is an infinite geometric series and can be evaluated immediately. The result is\n",
        "\n",
        "$$\n",
        "  \\sum_{n=0}^{\\infty} e^{-n \\beta \\hbar \\omega} = \\frac{1}{1 - e^{\\beta\\hbar\\omega}}\n",
        "$$\n",
        "\n",
        "Thus, the partition function takes the form \n",
        "\n",
        "$$\n",
        "  Z_{Oscillator} = -\\frac{1}{2}\\frac{\\beta \\hbar \\omega}{1 - e^{\\beta\\hbar\\omega}}\n",
        "$$\n",
        "\n",
        "We can also recall from statistical mechanics that the average energy of a system $\\langle \\hat{H} \\rangle$ is given by\n",
        "\n",
        "$$\n",
        "  \\langle \\hat{H} \\rangle = - \\frac{\\partial}{\\partial \\beta} \\ln{Z}\n",
        "$$\n",
        "\n",
        "Thus, the average energy of the quantum harmonic oscillator is given by \n",
        "\n",
        "$$\n",
        "  \\langle \\hat{H} \\rangle_{Oscillator} = \\hbar\\omega (\\tfrac{1}{2} + \\tfrac{1}{e^{\\beta\\hbar\\omega}-1})\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "c6JpxKGCBnft",
        "colab_type": "text"
      },
      "source": [
        "In general, the partition function cannot be evaluated exactly. Not only this, but if the number of states is large (as it is for real systems) the partition function becomes impossible to compute on a computer. Fortunately, there are methods by which we can still compute properties of systems without explicitly knowing its partition function.\n",
        "\n",
        "One such way is using Feynmann's path integral formulation of quantum mechanics in conjunction with MCMC. This is called path integral Monte-Carlo. The derivation is laborious (and done by the professor in his notes) so I will omit it here. The relevant results are:\n",
        "\n",
        "$$\n",
        "  \\langle \\hat H\\rangle_N = \\frac{\\int_{\\mathbb R^N}d^Nq\\,\\Phi(\\vec q) \\,e^{-U(\\vec q)}}{\\int_{\\mathbb R^N}d^Nq\\,e^{-U(\\vec q)}}, \\qquad \\Phi(\\vec q) = \\frac{1}{N}\\sum_{i=1}^N\\left[V(q_i) + \\frac{1}{2}q_i V'(q_i)\\right], \\qquad U(\\vec q) = \\sum_{i=1}^N\\left[\\frac{m}{2(\\beta/N)\\hbar^2}\\left(q_{i+1}-q_i\\right)^2 + (\\beta/N) V(q_i)\\right]\n",
        "$$\n",
        "\n",
        "where $\\vec q = (q_1, q_2, \\dots, q_N)$. To apply the this to MCMC, one needs to sample points from this space according to the un-normalized probability distribution $e^{-U(\\vec q)}$ using the Metropolis-Hastings algorithm (this is discussed more in the next section).  \n",
        "\n",
        "In order to do this in the case of the quantum harmonic oscillator we need to recall the potential of the quantum harmonic oscillator in one dimension is given by\n",
        "$$\n",
        "  V = \\tfrac{1}{2}m x^{2} \\omega^{2}\n",
        "$$\n",
        "\n",
        "where $m$ is the mass and $x$ is the displacement from equalibrium of the oscillator."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_eBKPKai9iiB",
        "colab_type": "text"
      },
      "source": [
        "## Algorithms"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0ebAHPbCJuwt",
        "colab_type": "text"
      },
      "source": [
        "The main algorithm we are using to numerically determine $\\langle \\hat{H} \\rangle$ is the Metropolis algorithm. The Metropolis algorithm is a simple algorithm for producing samples from distributions that may otherwise be difficult to sample from.  \n",
        "\n",
        "Suppose we want to sample from a distribution $\\pi$. The Metropolis algorithm works by simulating a Markov Chain, whose stationary distribution is $\\pi$. This means that, in the long run, the samples from the Markov chain look like the samples from $\\pi$.   \n",
        "\n",
        "The algorithm runs as follows:\n",
        "1. Initiate the sequence of samples in a state $s_{1}$ of your own choosing.\n",
        "2. If $s$ is the current state, propose a new state $s'$ according to a conditional probability distribution of our own choosing, defined such that $r(s'|s)$ is the probability of proposing $s'$ given the current state is $s$.\n",
        "3. Accept the proposed state, namely add it as the next state in the sample sequence, with probability\n",
        "\n",
        "$$\n",
        "  A(s|s') = \\mbox{min}(1,\\tfrac{p_{s'}}{p_{s}})\n",
        "$$\n",
        "\n",
        "If the proposed state is not accepted, add the current state $s$ to your sequence instead as the next state.\n",
        "\n",
        "4. Repeat from step 2 for a large number of samples.\n",
        "\n",
        "For our implementation of the Metropolis algorithm for the path integral we closely follow the procedure of Å. Larsen and F. Ravndal in the paper provided by the professor. It is as follows:\n",
        "\n",
        "1. Start from some arbitrary path $q$.\n",
        "2. For each point $q_{i}$ in the path, we propose a new point $q_{i}'$ by generating a random value uniformly on $(q_{i} - \\Delta,q_{i} + \\Delta)$, where $\\Delta = 1$.\n",
        "3. Accept or reject the proposed state according to probability distribution $e^{-U(q)}$. This step is repeated $\\bar{n} = 5$ times for each $q_{i}$.\n",
        "4. Return to step 2 and repeat for each $q_{i}$ in $q$.\n",
        "\n",
        "Going through all points in the path constitute on Monte Carlo iteration. The reason we repeat step 3 $\\bar{n}$ times is because Å. Larsen and F. Ravndal claim this improves convergence to equilibrium."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "s4Zh2F_xbNur",
        "colab_type": "text"
      },
      "source": [
        "For the problem of the quantum harmonic oscillator, as we have already seen, the partition function can be found analytically. Thus, we are able to compute whatever ensemble averages we want to exactly. For this reason this elaborate algorithmic approach is not necessary in this specific case, though it does serve as a good example of how to apply the path integral Monte Carlo method to more complex problems where the partition function is not known/cannot be written in a closed form easily."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "j3ernjrr9lsP",
        "colab_type": "text"
      },
      "source": [
        "## Code"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7A-BvNscNSkM",
        "colab_type": "text"
      },
      "source": [
        "The following cell contains the main functions and modules used to obtain our results."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "43T5tnBtiI7p",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy import constants\n",
        "plt.style.use(\"ggplot\")\n",
        "\n",
        "h_bar = 1 #constant of nature, set to 1 for convenience \n",
        "\n",
        "def V(q,m=1,omega=1):\n",
        "  \"\"\"\n",
        "  \"\"\"\n",
        "  return 0.5 * m * (omega**2) * (q**2)\n",
        "\n",
        "def dV(q,m=1,omega=1):\n",
        "  \"\"\"\n",
        "  \"\"\"\n",
        "  return m * (omega**2) * q\n",
        "\n",
        "def Phi(q,m=1,omega=1):\n",
        "  \"\"\"\n",
        "  \"\"\"\n",
        "  N = np.size(q)\n",
        "  term = V(q,m,omega) + (0.5 * q * dV(q,m,omega))\n",
        "  return np.sum(term)/N\n",
        "\n",
        "def U(q,beta=1,m=1,omega=1):\n",
        "  \"\"\"\n",
        "  \"\"\"\n",
        "  N = np.size(q)\n",
        "  term = np.zeros(N)\n",
        "\n",
        "  for i in range(N-1):\n",
        "    term[i] = (m*(q[i+1]-q[i])**2)/(2*(beta/N)*(h_bar**2)) + (beta/N)*V(q[i])  \n",
        "\n",
        "  term[N-1] = (m*(q[0]-q[N-1])**2)/(2*(beta/N)*(h_bar**2)) + (beta/N)*V(q[i])\n",
        "\n",
        "  return np.sum(term)\n",
        "\n",
        "def prob(U):\n",
        "  \"\"\"\n",
        "  \"\"\"\n",
        "  return np.exp(-U)\n",
        "\n",
        "def exact_H(beta,omega=1):\n",
        "  \"\"\"\n",
        "  \"\"\"\n",
        "  return (h_bar * omega) * (0.5 + (np.exp(beta*h_bar*omega)-1)**-1)\n",
        "\n",
        "def inital_config(N):\n",
        "  \"\"\"\n",
        "  \"\"\"\n",
        "  return np.random.uniform(-1,1,N)\n",
        "\n",
        "def path_integral_mcmc(q,m,beta,omega,num_steps,delta=1,n_bar=5):\n",
        "  \"\"\"\n",
        "  \"\"\"\n",
        "  N = np.size(q)\n",
        "  q_proposed = np.zeros(N)\n",
        "\n",
        "  p_current = prob(U(q,beta,m,omega))\n",
        "  p_proposed = 0\n",
        "\n",
        "  H_avg = np.zeros(num_steps)\n",
        "  H_avg[0] = Phi(q,m,omega)\n",
        "\n",
        "  for i in range(num_steps-1):\n",
        "    for j in range(N):\n",
        "      for k in range(n_bar):\n",
        "        q_proposed = np.copy(q)\n",
        "        q_proposed[j] = np.random.uniform(q[j]-delta,q[j]+delta)\n",
        "\n",
        "        p_proposed = prob(U(q_proposed,beta,m,omega))\n",
        "\n",
        "        if np.random.rand() < min(1,p_proposed/p_current):\n",
        "          q[j] = q_proposed[j]\n",
        "          p_current = prob(U(q,beta,m,omega))\n",
        "\n",
        "    H_avg[i+1] = H_avg[i] + (Phi(q,m=1,omega=1)-H_avg[i])/(i+2)\n",
        "\n",
        "  return q, H_avg"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "i91Dt5m39suK",
        "colab_type": "text"
      },
      "source": [
        "## Results"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XfOi1ESQOiom",
        "colab_type": "text"
      },
      "source": [
        "In the following cells we run code to obtain our results."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "iclISDqKPDsK",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#Running MCMC to approximate <H> for various values of beta.\n",
        "#Takes awhile to run...\n",
        "betas_1 = np.array([.1,.3,.5,.7,.9,2,4,6,8,10])\n",
        "H_avg = np.zeros(np.size(betas_1))\n",
        "\n",
        "m = 1\n",
        "omega = 1\n",
        "\n",
        "N = 5\n",
        "delta = 1\n",
        "n_bar = 5\n",
        "num_steps = np.array([80000,80000,80000,80000,80000,50000,50000,50000,50000,50000])\n",
        "\n",
        "for i in range(np.size(betas_1)):\n",
        "  q = inital_config(N)\n",
        "  q, H = path_integral_mcmc(q,m,betas_1[i],omega,num_steps[i])\n",
        "  H_avg[i] = H[-1]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1GxO1soHPSnW",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "#Getting exact values of <H> for various values of beta.\n",
        "betas_2 = np.linspace(.01,10,100)\n",
        "H_avg_exact = exact_H(betas_2,omega)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "W9fqTIaXPzz6",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 314
        },
        "outputId": "43e59f4f-c42d-4e23-aee6-c236317642f4"
      },
      "source": [
        "#Plotting our results.\n",
        "plt.figure(figsize=(8,4)) \n",
        "\n",
        "plt.scatter(betas_1,H_avg,label=\"MCMC\")\n",
        "plt.plot(betas_2,H_avg_exact,linestyle=\"--\",label=\"Exact\",color=\"blue\")\n",
        "\n",
        "plt.title(r\"$\\langle \\hat{H} \\rangle$ for various values of $\\beta$\",fontsize=18)\n",
        "plt.xlabel(r\"$\\beta$ (1/$k_{b}$T)\")\n",
        "plt.ylabel(r\"$\\langle \\hat{H} \\rangle$ (Joules)\")\n",
        "\n",
        "plt.ylim(0,12)\n",
        "\n",
        "plt.legend()\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfkAAAEpCAYAAACQiFvLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0\ndHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deVxU5f4H8M+ZGYZh3xFBEBX3BTHR\ntHLFrVTIW1pumXktxbyZG3lNy36WlaZ5w2xTc7tXqxukphaRmtqigvuagqiI7DsMzMz5/YGMzWWR\nYQbnDHzer5ev4jxn+Z6H0e8853nO8wiiKIogIiKiRkdm6QCIiIioYTDJExERNVJM8kRERI0UkzwR\nEVEjxSRPRETUSDHJExERNVJM8kRERI0UkzxZjUuXLuGNN97ApUuXLB0KEZFVEDgZDlkDURTx2GOP\n4fTp0+jevTsOHjwIQRAsHRYRkaSxJU9W4YsvvkBycjJOnjyJq1ev4ssvv7R0SEREkseWPEleZmYm\nOnbsiO3bt2PIkCH44YcfMGnSJFy4cAHu7u6WDo+ISLKY5ImIiBopPq4nIiJqpJjkyeKysrKgUCjw\n0UcfmXSe9evXQyaT4fbt20Ydl5SUhIiICHh5eUEQBEyZMsWkOKTqwIEDEAQBmzZtsnQoFtNY66C+\nn+FDhw5hyJAhcHZ2hkqlQvfu3fHtt982bLD0QDHJk8Xt3r0bWq0Wo0ePNtien58PmUyG4cOH13hs\nhw4d4OLiAlEU9cfv2rXLqOtPmTIFBw8exMKFC7Flyxa8+OKLxt8EkQXV5zO8b98+DBo0CKmpqVi6\ndCnee+89FBQUYOzYsbh27doDiJoeBIWlAyCKiYlBSEgIAgICDLYnJCRAFEU8/PDD1R6Xm5uLy5cv\no1+/fhAEAb6+vujZsydiYmIwffr0Ol1brVbjl19+waxZszBv3jyT70XK+vXrh5KSEtjY2Fg6FDKj\n+nyGdTodZsyYgaCgIBw/fhx2dnYAgMDAQISHh+PXX39F69atGzJsekDYkieLKikpwQ8//IDw8PAq\nZQkJCQCA3r17V3vsH3/8AVEU0aNHD/22iIgIxMfHo7CwsE7Xv3PnDkRRNPsofa1Wi+LiYrOes74q\nY5HJZFCpVJDL5ZYOicyoPp/hI0eOIDk5GbNnz9YneAD6L4B/3UbWjUmeLCouLg7FxcXVJvkTJ04A\nAHr16lXtsb///jsAGCT58PBwqNVq7Nu3777XnjJlClq2bAkAePPNNyEIAgRBwIEDBwBUvLoXGRkJ\nf39/KJVK+Pv7IzIyEllZWQbn2bRpEwRBQFxcHN566y20adMGKpUKO3furPa6e/fuhSAIWLt2bbXl\nffr0gZeXF8rLywEABQUFWLx4MXr37g1PT0/Y2toiKCgIUVFRVb5I1BZLTf3Rdb3PN954A4IgIDk5\nuUrMgYGBGDBggMG20tJSvPHGG2jfvj3s7e3h6uqKrl27Yv78+dXe94Oon+oYe19qtRpvv/02Onfu\nDJVKBVdXV4waNQqJiYlmuf9Kdfm93O8zXJOjR48CAAYPHmywPT4+HoDh3ymybnxcTxYVExODgIAA\ndO/evUpZQkIC/Pz8kJOTg5ycnCrlhw4dAmD4D1Lnzp0RFBSEmJgYPPXUU7Ve+8UXX0T37t0xZ84c\nPPnkkxgzZgwAoGPHjsjLy0Pfvn3x559/YurUqejRowcSExPx8ccfIz4+Hn/88QecnJwMzjdv3jyU\nl5fj73//O5ydndG+fftqrzt06FD4+Phg8+bNmD17tkHZlStX8Ntvv2H27Nn6VtWtW7fw+eef429/\n+xvGjx8PhUKBgwcP4r333kNiYiL2799f5RrVxaJWq6vsV5/7rKvIyEhs2LABkydPxquvvgqNRoMr\nV67oE0lNHkT91Fd5eTmGDx+Oo0ePYtKkSZg1axby8vLw2Wef4ZFHHsGhQ4fQs2dPk+4fqPvvpbbP\ncG0SExPh6OiIdu3a6bfdvHkTGzZsQGhoKAIDA+tfSSQtIpGFaLVa0dvbW5w1a1aVsoKCAlEmk4kA\nav1jb28varVag2Pnzp0rurm5ieXl5feNISkpSQQgLl261GD7okWLRABidHS0wfaPPvpIBCAuXrxY\nv23jxo0iALFdu3ZiUVFRne593rx5IgDx3LlzBtsXL14sAhBPnDih36ZWq8WysrIq56jc9/fff69T\nLD///LMIQNy4cWO97nPp0qUiADEpKalKLC1bthT79+9vsM3NzU0cMWJEjXVQm4aqn+rqwJj7+uCD\nD0QA4r59+wz2y8vLE/39/Q32NeX+jfm91PQZrk379u3FRx99VBRFUbx06ZK4efNmsVWrVqKNjY14\n+PDhesVM0sTH9WQxv/76K9LT06t9VH/y5EnodDq88MIL+Oqrr6r8WbduHQCge/fukMkMP8bh4eHI\nycnRt/Tr49tvv4WXl1eVAXwvvvgivLy8qn3NaMaMGbC3t6/T+Z977jkAwObNm/XbRFHE1q1b0aVL\nF4OnE0qlUt9q1Wg0yMnJQWZmJsLCwgDc67aoTyz1uc+6cnFxwblz53D27Fmjj23o+qmvrVu3okOH\nDnjooYeQmZmp/1NWVoYhQ4bg8OHDKCkpAWDa/Tfk76WoqAhXrlxBjx49cP36dbRv3x6TJ09GUlIS\nZs6ciT59+tT73CQ9TPJkMTqdDgCqHQhW2R//zDPP4Kmnnqryp3Lkb3V9h5Xn02q19Y4tKSkJ7du3\nh0Jh2KOlUCjQrl27al8x+uujz/upTFTbtm3T18OhQ4eQnJyMyZMnV9l/3bp16NatG2xtbeHu7g4v\nLy99X3F1XRl1jaU+91lXa9asQU5ODrp27Yo2bdpg2rRpiI2N1d9vbRq6furrwoULuHjxIry8vKr8\n2bBhA7RaLTIzM02+/4b8vZw+fRo6nQ4PPfQQnJ2dsWvXLkRHR2Po0KH48MMP8fzzz9f73CQ9TPJk\nMX379oWnpydiY2OrlFWOrO/WrVu1x546dQpA9Uk+JiYGLi4uVQZMNbS6tuIrTZ48GTdv3tT30W7e\nvBlyuRwTJ0402O+DDz5AZGQkmjdvjk8++QR79uzBjz/+qB9AV13SMDaWuqht1T+NRlNlW3h4OJKT\nk7FlyxYMGjQIP/30EyIiIjBgwACUlZXd93oNWT/1vS9RFNG1a1f8+OOPNf7x8vIyy/03lJMnTwKo\n+Lvj5uaGkSNHYubMmdi/fz/CwsKwdetW/aBGsn4ceEcWI5fLMXLkSMTGxmLNmjUGZSdOnICPjw+8\nvb2rPfb06dMAqk/ysbGxePzxx016H7x169a4dOkSNBqNQWtKo9Hg8uXLZnmHePz48Zg/fz42b96M\nRx55BF9//TWGDBmC5s2bG+y3ZcsWBAYGYu/evQZdE3V5g+B+jLnPyle0srOzDQZmlZaW4vbt2wgK\nCqpyfnd3d0ycOBETJ06EKIqIiorCe++9h9jYWDz99NO1xvag6seY+2rbti0yMjIwaNCgKt1ENZ27\nPvffkJ+/kydPws7OrtrBeTqdDt7e3pxLoRFhS54sqrK1U5m0AaC4uBgXL15EcHBwjcedOnUKtra2\n6Ny5s8H2Cxcu4PLly9X28xsjIiICGRkZ+Pzzzw22f/bZZ8jIyMCTTz5p0vkBwMvLCyNGjMB///tf\nbNu2Dfn5+fq+6L+Sy+UQBAHiX9aS0mg0WLFihckxGHOflV0AcXFxBvuuXr26SmtZq9UiNzfXYJsg\nCAgJCQFQkVDv50HVjzH3NXnyZKSlpeGDDz6o9lx37twBYPr9N+TnLzExETqdDklJSQbbjx49igMH\nDmDSpEn1PjdJD1vyZFFDhw6FnZ0dYmNj9Y/mT506Ba1WW2OSV6vVuHjxIrp3716lzzI2NhZKpRIj\nRowwKa4FCxbgq6++QmRkJBISEhASEoLExER88cUXaN++PRYsWGDS+Ss999xz+O677zB37ly4uLgg\nIiKiyj5PPfUUXnvtNYwYMQJjxoxBfn4+tm/fbpbWljH3GRYWhvbt22PJkiXIyspCq1atcPjwYfz2\n22/w9PQ0OG9BQQGaN2+O0aNHIyQkBN7e3khKSsLHH38MNzc3jBo1qk7xPYj6Mea+/vGPf+DHH3/E\n/PnzER8fj0GDBsHZ2RkpKSn46aefoFKp8PPPP5t8/w31+dNqtTh79izUajUee+wxREZGolmzZjhz\n5gw2bNiA4OBgvP766/U6N0mUBUf2E4miKIqjR48We/bsqf+58jWhrVu3Vrt/QkKCCECcPn16lbKH\nH35YHDZsWJ2vXdvrR+np6eKMGTNEPz8/UaFQiH5+fuLMmTPFjIwMg/0qX1v7+eef63zdSmq1WnR3\ndxcBiNOmTat2H41GI7799ttimzZtRKVSKQYEBIjz588Xz58/XyX22mKp7vUxY+5TFCtetxo2bJho\nZ2cnuri4iE8//bR48+bNKq+aqdVqMSoqSgwNDRXd3d1FpVIptmzZUnz++efFy5cvW6x+aqqDut6X\nKIpieXm5+OGHH4o9e/YU7e3tRXt7ezEoKEgcP368uH//frPdf11/L8a8Qnf27FkRgPjss8+K3bt3\nF5VKpejq6ip2795dfOedd8S8vLw6xUbWg+vJk8Vt2LAB06ZNw40bN+Dn51fv86SlpcHX1xfR0dGY\nMWOGGSMkahy2bduGiRMn4siRI+jbt6+lw6EHgH3yZHGjRo2CIAj47rvvTDpP5fGm9scTNVYnT56E\nIAjo2rWrpUOhB4QteZIEtVoNuVxepY/dGBqNBlqtFra2tmaMjKjxCAsLw9WrV6sMuqPGSxID79at\nW4eEhAS4uLhg1apVACpeizlx4gQUCgWaNWuGmTNnwsHBwcKRUkMxR2JWKBQmfUkgauxOnTrFGe2a\nGEm05M+fPw+VSoXo6Gh9kj916hS6dOkCuVyOrVu3AkCVSTCIiIioZpLok+/UqRMcHR0NtgUHB+un\nJ23Xrl2d3qslIiKie6zi2WZ8fHytI0Hj4uL0E1mYY4IQIiKixkDySf6///0v5HI5HnvssRr3CQsL\n0684BQCpqalmu76np6d+wYlly5yxZYs9rlxJM9v5m4K/1iHVD+vQdKxD82A9ms7cdejr61tjmaST\n/IEDB3DixAksWbKk1kUkHpQxY4rRvbvlFpYgIiIyhmST/MmTJxEbG4s333xTMq9EdemiQZcuVVfb\nIiIikiJJJPk1a9bg/PnzKCgowEsvvYSxY8fi22+/hUajwVtvvQWgYvWn6dOnWzTO5GQ5bt+Wo08f\ntuaJiEj6JJHkX3nllSrbBg0aZIFIardtmz02bHDE1au3LR0KEVGTIYoiSktLodPpJNF1a6o7d+5A\nrVYbdYwoipDJZFCpVEbVgSSSvDWx/KwCRERNS2lpKWxsbBrNZFcKhUL/irgxNBoNSktLYWdnV+dj\nJPGevLVoBF8giYisjk6nazQJ3hQKhQI6nc6oY5jkjcSWPBHRg9UYHtGbi7F1wSRvBH7OiIjImvD5\nhxGeeqoEoaEcWU9E1NT4+/ujQ4cO+p/Dw8Mxa9Yss5z77NmzuHPnDgYPHmyW8/0Vk7wR2rbVoG1b\nvidPRNTUqFQq/Pjjjw1y7nPnzuH06dMNkuT5uN4If/6pwI8/SmNiHiIisqz8/Hw89thj+PPPPwEA\nM2fOxLZt2wAAUVFRGDFiBAYOHIiVK1fqjzl58iSeeOIJhIWF4YknnkB+fj5WrlyJ7777DkOGDEFs\nbKxZY2RL3gjffGOH6GhHpKTwPXkiIkt56imPKttGjizBlCnFKCkRMGmSe5Xyp58uxrhxJcjOlmH6\ndDeDsq+/zrrvNUtLSzFkyBD9z7NmzUJ4eDiWL1+OOXPmYNq0acjLy8OECRMAAAsXLoSbmxu0Wi3G\njRuH8+fPIygoCDNmzMCnn36Krl27oqCgAHZ2dpg3bx5Onz6N5cuXG1sV98UkbySOricianpqelzf\nr18/7N69G4sWLTIo37VrF7Zt2watVos7d+7gypUrEAQB3t7eCAkJgUajgZOTU4PHzSRvBI6uJyKy\nvNpa3nZ2Yq3l7u66OrXc60qn0+HKlSuws7NDXl4efH19kZKSgk8++QR79uyBq6srXnnlFZSWlprt\nmsZgn7yR2JInIqJKn376Kdq2bYvo6Gi8+uqrKC8v1z+Gd3Z2RkZGBn7++WcAQJs2bZCeno7ExEQA\nQGFhITQaDRwdHVFYWNgg8bElbwS25ImImqb/7ZMfOHAgxo4di3//+9/Ys2cPHB0d0bt3b3z44YeY\nN28eunTpgn79+sHX1xehoaEAAKVSiY8//hiLFi1CSUkJVCoVduzYgb59+yI6OhpDhgzR9/WbiyCK\nja9tmpqaarZzeXp6IjMzEwCQkiLHrVtchc5Yf61Dqh/WoelYh+ZhiXosLi6Gvb39A71mQ1IoFNBo\n6vc6dnV14evrW/O16nWVJiogQIuAAK2lwyAiIqoT9skb4cIFBb791o798kREZBWY5I2wd68Ks2a5\n3X9HIiIiCWCSrwe25ImIyBowyRuhcnQ9kzwREVkDJvk60mWkQTxxBACg/Xw1dBlpFo6IiIiodkzy\ndaDLSIO4egmEW8kAAPGPQxBXL2GiJyJqIvz8/PDyyy/rf9ZoNOjatSsmT56s3xYfH48RI0ZgwIAB\nGDp0KN58800AwKpVq+Dn54ekpCT9vp999hn8/Pxw6tQpAEBRUREWLFiAvn37Yvjw4XjqqaeQkJBg\nctx8ha4uYrcBGWkY1yIGA7yOQC5ogYy0iu3T5lo6OiIiamD29va4ePEiSkpKYGdnh0OHDsHHx0df\nfvHiRSxevBibN29GUFAQtFottm7dqi/v2LEjYmNj8corrwAAdu/ejfbt2+vL582bh4CAABw+fBgy\nmQwpKSm4fPmyyXGzJV8HYm42AKCZKgvdXC7e65u/u52IiKRDl5EG3eeroF35T+g+X2W2p66DBg3C\nTz/9BACIiYlBRESEvmzdunWYPXs2goKCAAByuRzPPfecvnzYsGH44YcfAADJyclwdnaGu7u7/ufE\nxEQsWLAAMllFWg4ICEBYWJjJMTPJ14HgWvGLOJvXHltTxkArygy2ExGRNFR2r4q/HwQunYH4+0Gz\nda+Gh4cjNjYWpaWluHDhAkJCQvRlly5dQteuXWs81snJCb6+vrh48SK+/fZbjBo1Sl92+fJldO7c\nGXK53OQY/xeTfF2ETwC8fHAgoy8WnV0EjU4OePlUbCciIum4271qoLJ71USdOnXCzZs3ERsbi0GD\nBhl9/OjRoxEbG4t9+/ZhxIgRJsdTF0zydSDz8oEwZxkE/0AAgNDrMQhzlkHm5VP7gURE9EDV1I1q\nru7VoUOHYtmyZQaP6gGgXbt2OHPmTK3HDhkyBN988w38/PwM1pJv164dzp8/D63W/NOmM8nXkczL\nB0LoYwAAYcorTPBERBJUUzequbpXx40bh1dffRUdO3Y02D5jxgz861//wtWrVwFUrDO/efNmg33s\n7Ozwz3/+Uz/4rlJgYCC6deuGlStXonLNuBs3biAuLs7keDm63gicDIeISOLCJwDXLhk+sjdj96qv\nry9eeOGFKts7deqEN954A5GRkSgpKYEgCNUOnAsPD692FbqVK1di2bJleOSRR6BSqeDu7o7Fixeb\nHK8klppdt24dEhIS4OLiglWrVgEACgsLsXr1amRkZMDLywtz5syBo6Njnc7XUEvNRkc74u23nfHn\nn6mwszPbJRo9LvFpOtah6ViH5mENS83q7vbBi7nZFS348AmSevra5JaaHTBgAIYPH47o6Gj9tpiY\nGHTt2hURERGIiYlBTEwMJk6caMEogfHjizB8eAlsbS0aBhER1ULm5cM5TO6SRJ98p06dqrTSjx07\nhv79+wMA+vfvj2PHjlkiNANubiLatNFCJolaIyIiqp1k01VeXh7c3CqWdXV1dUVeXp6FIwJOnrTB\n+vUOKC+3dCRERE2HBHqVJcPYupDE4/r7EQQBQuWot2rExcXpRyGuWLECnp6eZru2QqHQn+/UKRne\nekuBV16xQx2HBxAM65Dqh3VoOtaheViiHgVBgE6ng42NzQO9bkNSKIxPv+Xl5XB0dISHh0fdr2P0\nVR4QFxcX5OTkwM3NDTk5OXB2dq5x37CwMINRjOYcFPLXQSYlJQ4AXJCVlYXSUn6zrCsOeDId69B0\nrEPzsEQ9iqKI0tJSFBcX19rgsxa2trZQq9VGHSOKImQyGVQqVZX6l/zAu+r07NkTBw8eREREBA4e\nPIjQ0FBLh6THJ0dERA+OIAiwa0SvND3IL0qSSPJr1qzB+fPnUVBQgJdeegljx45FREQEVq9ejfj4\neP0rdFLBJE9ERNZAEkn+f2f/qbRkyZIHHEntGsFTIiIiakIkO7peiiZOLEZCQhocHNiUJyIi6ZNE\nS95aODiITPBERGQ12JI3wvHjNli50gklJZaOhIiI6P6Y5I2QmKjE6tVOUKvZOU9ERNLHJG8EDrwj\nIiJrwiRfD3yFjoiIrAGTvBG4njwREVkTJnkj8HE9ERFZEyZ5I0ycWITLl2/DzY1NeSIikj6+J28E\npRJQKpngiYjIOrAlb4Rjx5R4801nFBXxuT0REUkfk7wRzp1T4NNPHVFayiRPRETSxyRvBLm84r/l\n5ZaNg4iIqC6Y5I3g6FjRH19YyJY8ERFJH5O8ERwcdACAwkJWGxERSR+zlRGcnEQIgojiYrbkiYhI\n+vgKnRF69y5DSsptyPjViIiIrACTvBGY3ImIyJowbRmhqEjA/PkuOHjQ1tKhEBER3ReTvBFkMhHb\ntzvgzBkbS4dCRER0X0zyRlCpALlcREEBB94REZH0MckbQRAqRthzWlsiIrIGTPJGcnDQoaCA1UZE\nRNLHbGUkb28dR9kTEZFV4Ct0Rtq9O9PSIRAREdUJ26RERESNFJO8kTZscMBrr7lYOgwiIqL7kvzj\n+t27dyM+Ph6CIMDf3x8zZ86EUqm0WDznzilw4IAK77yTZ7EYiIiI6kLSLfns7Gzs3bsXK1aswKpV\nq6DT6XD06FGLxuToKHKpWSIisgqSTvIAoNPpUFZWBq1Wi7KyMri5uVk0HkfHivfkdTqLhkFERHRf\nkn5c7+7ujlGjRmHGjBlQKpUIDg5GcHCwRWNydNRBFAUUFwtwdBQtGgsREVFtJJ3kCwsLcezYMURH\nR8Pe3h4ffPABDh06hH79+hnsFxcXh7i4OADAihUr4OnpabYYFAqFwflatZKhVSsRDg4eMONlGrX/\nrUMyHuvQdKxD82A9mu5B1qEgiqJkm6O//vorTp48iRkzZgAADh48iCtXrmDatGm1Hpeammq2GDw9\nPZGZyXfjTcE6NB3r0HSsQ/NgPZrO3HXo6+tbY5mk++Q9PT1x5coVqNVqiKKIM2fOwM/Pz9JhERER\nWQVJJ/m2bdvi4YcfxsKFCzFv3jyIooiwsDCLxnT5sgLjxnng1CkuN0tERNIm6T55ABg7dizGjh1r\n6TD0ysqAw4dtcfu2HMHB5ZYOh4iIqEb1TvJ5eXk4deoUkpOTUVxcDHt7ewQGBqJbt25wdXU1Z4yS\nUjminu/KExGR1Bmd5G/evIkdO3bg3LlzaN26Nfz8/ODq6oqSkhIcOnQImzZtQufOnTFu3Di0aNGi\nIWK2KCZ5IiKyFkYn+XXr1mH06NGYPXs2bGyq9kuXl5fj+PHj+Pjjj7F8+XKzBCkljo4Vs+BwTXki\nIpI6o5P822+/XWu5jY0N+vTpgz59+tQ7KCmztQU6dSqHszOnvCMiImkzaeDd2bNn4e3tDW9vb+Tk\n5GDbtm2QyWQYP358o+2XFwTgxx8zLB0GERHRfZn0zPmLL76ATFZxis2bN0Or1UIQBHzyySdmCY6I\niIjqz6SWfHZ2Njw9PaHVanHq1CmsW7cOCoUCL774ornik6TZs13h6qrDsmX5lg6FiIioRiYleTs7\nO+Tm5uLGjRto0aIFVCoVNBoNNBqNueKTpOvXFbhzR7KzARMREQEwMckPHz4cr732GjQaDaZMmQIA\nuHjxYqOfetbJSYfcXI6uJyIiaTMpyUdERKBXr16QyWTw8fEBULE87EsvvWSW4KTKwUHEjRt8T56I\niKTN5Glt/3f1m9pWw2ksnJx0KCxkS56IiKTN5CR/+vRpHDlyBHl5eYiKisLVq1dRUlKCLl26mCM+\nSWrfXoOsrDJLh0FERFQrk5qje/fuxWeffYbmzZvjwoULAAClUon//Oc/ZglOqv7+9yJs3Jhj6TCI\niIhqZVKS//777/H6668jIiJC/768n58fUlNTzRIcERER1Z9JSb6kpASenp4G2zQaDRQKya9ga5Ld\nu1Xo29cbd+6wX56IiKTLpCzVsWNHxMTEGGzbu3cvOnfubFJQUqfRCLh+XcFFaoiISNJMylJTp07F\nH3/8gcjISJSWluIf//gHfv31Vzz33HPmik+S7q1Ex9foiIhIukx6ru7m5oZ33nkHV69eRUZGBjw8\nPBAUFKTvn2+suKY8ERFZA6OT/NmzZ6vd7uTkhLKyMpw/fx4AGvUrdJUt+aKixv1lhoiIrJvRSf7j\njz++7z6CIOCjjz6qV0DWwMNDhyFDSuHmxjXliYhIuoxO8tHR0Q0Rh1Vp3lyHTZuyLR0GERFRrfi8\nmYiIqJEyaeDdjh07aiwbN26cKaeWvNBQbzzzTAnmzi2wdChERETVMinJZ2VlGfycm5uL8+fPo1ev\nXiYFZQ2Ki2XIzuaDECIiki6TkvzMmTOrbDt58iQOHz5symmtgqOjju/JExGRpJm9KdqtWzccO3bM\n3KeVHCcnEUVFTPJERCRdJrXk79y5Y/CzWq3G4cOHq8xn3xg5OIic1paIiCTNpCQ/e/Zsg5+VSiVa\ntWqFyMhIk4KyBmFhpRBFS0dBRERUswYbXW8uRUVFWL9+PW7cuAFBEDBjxgy0a9euwa97Py+/XGjp\nEIiIiGpl8pqwt2/fxpEjR5CdnQ13d3c88sgjaN68uTliAwBs3LgR3bt3x9y5c6HRaKBWq812blOJ\nIiCwW56IiCTKpE7l48ePIyoqCrdu3YKjoyNSU1MRFRWF48ePmyW44uJiXLhwAYMGDQIAKBQKODg4\nmOXcpnrnHSd06uRj6TCIiIhqZFJL/t///jfmz59vsBjNuXPnsGHDBvTs2dPk4NLT0+Hs7Ix169bh\n+vXraN26NaZMmQKVSmWwX9Pq1y0AACAASURBVFxcHOLi4gAAK1asMOvAP4VCUe35XFzkyM+Xwd3d\nE4180T2T1VSHVHesQ9OxDs2D9Wi6B1mHJiX57OxsdOzY0WBbhw4dqkySU19arRZJSUmYOnUq2rZt\ni40bNyImJgbPPPOMwX5hYWEICwvT/5yZmWmW6wOAp6dnteeTyRwAuOD69Sw4OXEEXm1qqkOqO9ah\n6ViH5sF6NJ2569DX17fGMpPaoIGBgdi1a5fBtt27dyMwMNCU0+p5eHjAw8MDbdu2BQA8/PDDSEpK\nMsu5TcU15YmISOpMaslPmzYN7777Lvbu3QsPDw9kZWVBqVRi4cKFZgnO1dUVHh4eSE1Nha+vL86c\nOYMWLVqY5dymupfkZQC45CwREUmPSUnez88Pq1evxpUrV/Sj64OCgqBQmDxoX2/q1KlYu3YtNBoN\nvL29q51K1xLati3HCy8UwsGBCZ6IiKRJEMXGN6VLamqq2c7F/ifTsQ5Nxzo0HevQPFiPpnuQffL1\nanLHx8fXWi4IApycnBAYGNhoR2GKIlBaCsjlgFJp6WiIiIiqqleS/+WXX+67T3FxMVJTUzFhwgQM\nHz68PpeRtJQUOfr2bYbVq3MwdmyJpcMhIiKqol5JfunSpXXa78aNG3j77bcbZZKvfG2OK9EREZFU\nGf0K3ffffw+NRlPrPuXl5fj+++/h7++PRx99tN7BSVnlgDuuREdERFJldEs+NzcXL7/8MkJCQtCp\nUyf4+vpCpVKhtLQUqampOH/+PBITE9G/f38AwIQJE8wetBTY2gJKJdeUJyIi6TI6yY8fPx4jR47E\ngQMHEB8fj5SUFBQVFcHR0REBAQEICQnBs88+Cycnp4aIV1IcHHRsyRMRkWTVq0/e2dkZo0ePxujR\no80dj1WZObMIQUHllg6DiIioWuabtaYJmjmTa8oTEZF08VmzCfLzBWRksAqJiEia2JI3QWSkG7Kz\nZdizh7M/ERGR9LAZagIHBxEFBRxdT0RE0sQkbwInJx2KiliFREQkTSY/ri8pKUFiYiJSUlIQHByM\nDh06QBCaRuuWLXkiIpKyeiX57OxsHD9+HMeOHUNycjK6du2KgIAA7Ny5E6mpqQgODkZoaCiCg4Oh\nbMSrtzg5iSgqkkGnA2Rs0BMRkcTUK8kvWbIEPXv2xOjRo9GpUyfI5XIAQEREBIqKipCQkIBffvkF\nx44dk8z67w1hwIBSuLjooNUyyRMRkfRwPfn74NrJpmMdmo51aDrWoXmwHk33INeTr1f789NPP8UP\nP/yAy5cvQ61W1zswa1dcLODyZQVKStgvT0RE0nPfx/WFhYXYuXMnpk6dqt9mY2ODX3/9FV999RUK\nCwvh4+ODli1bolWrVggMDERgYCBcXFwaNHApOH3aBn/7mye2b89C//5N98sOERFJU61JPjExETt2\n7MD48eMNtj///PP6/8/NzcWcOXPQq1cvXL16FT/88AMyMzOxY8eOholYQlq3rlhy99o1Oe4uukdE\nRCQZ923Ji6IIWS2jylxdXSGTyTBs2DC4u7sDAPLz880XoYR5eeng6KjD1aucOJCIiKSn1uwUEhKC\noKAg7Ny5E126dKnzSZ2dnU0OzBoIQkVr/to1JnkiIpKe+w68c3JywgsvvGCwLSMjo8ECsjZM8kRE\nJFX1yk6zZs2CSqVCixYt4O/vj7KyMly7dg0ODg6wtbU1d4ySNnVqEf72txJLh0FERFRFvZL8hg0b\nkJycjOTkZCQlJcHHxwerVq0CALRo0QJt2rTBSy+9ZNZApeqhh8otHQIREVG16pXkHRwc0LlzZ3Tu\n3Fm/TaPRICUlBUlJSbh+/brZApQ6tRo4csQWrVpp0KqV1tLhEBER6dVrMpwvv/wS58+fx18ny1Mo\nFGjdujUGDx6MZ599FkePHsUPP/xgtkClqrxcwKRJHtizx87SoRARERmoV0vey8sLX3/9NW7duoXg\n4GD07NkTAQEBOH36NI4dO4br16+jS5cuGDRokLnjlRxHRxHNmmk5+I6IiCSnXpnp8ccfx+OPP46i\noiIkJibiyJEj+M9//oNu3bohPDwcnTp1qvXdemPpdDpERUXB3d0dUVFRZjuvuVSMsJdbOgwiIiID\nJjU/HRwc8Oijj+LRRx81VzzV+v777+Hn54eSEmmOYm/dWoN9+1SWDoOIiMiA5BdIzcrKQkJCAgYP\nHmzpUGrUurUGWVly5OZyoRoiIpIOyXckb9q0CRMnTqy1FR8XF4e4uDgAwIoVK+Dp6Wm26ysUivue\nb+pUIDy8HIGBHlBIvkYfvLrUIdWOdWg61qF5sB5N9yDrUNIp6cSJE3BxcUHr1q1x7ty5GvcLCwtD\nWFiY/mdzrtNbl3V/lUqgeXMgN9dsl21UuP606ViHpmMdmgfr0XQPcj15SSf5S5cu4fjx40hMTERZ\nWRlKSkqwdu1azJ4929KhVbFjhx28vXUYOJBLzhIRkTRIOsmPHz9ev8ztuXPnsGvXLkkmeABYu9YJ\nXbqUM8kTEZFkSH7gnbXgQjVERCQ1VpOV/ncaXalp3VqDX39VQqcDzDhFABERUb0xHZlJmzYalJTI\nkJbGKiUiImlgRjKT1q01AIDkZKt5OEJERI0cM5KZhIaW4ezZ23BzE++/MxER0QPAJG8mtraArS0T\nPBERSQcf15vRli32iI52tHQYREREAJjkzerIEVts325v6TCIiIgAMMmbjS4jDa3yDiPlugyl69dA\nl5Fm6ZCIiKiJY5I3A11GGsTVS9Aq7xfoRBmuH7wKcfUSJnoiIrIoJnlziN0GZKShreM1AMDZvA5A\nRlrFdiIiIgthkjcDMTcbANDZ+TKaq9Jwu7SZwXYiIiJL4Ct0ZiC4ukMEIBd0ODIgHAqZVr+diIjI\nUtiSN4fwCYCXDwDoE7zo6VOxnYiIyEKY5M1A5uUDYc4yCL37Q9u2G/529it8oPoMsruJn4iIyBL4\nuN5MZF4+wLS5UAKQ/eaBuCMyLECGpcMiIqImjC35BjB4sBrnztng9m1WLxERWQ6zUAMYPLgUABAf\nr7JwJERE1JQxyTeA9u018PPT4KefbC0dChERNWHsk28AggBERhZyVToiIrIoJvkG8txzxZYOgYiI\nmjg+rm9AGRkyJCTYWDoMIiJqotiSb0CvvuqKa9cUOHw4HYJg6WiIiKipYUu+AQ0eXIrkZAWuXpVb\nOhQiImqCmOQb0JAhagiCiK++srd0KERE1AQxyTcgPz8tRo4sxaZNDsjN5fN6IiJ6sJjkG9js2QUo\nLhbwyy98Z56IiB4sDrxrYJ06afD773fg66uzdChERNTEMMk/AL6+Ougy0lCw42s4lt6uWGc+fAJX\nqSMiogYl6SSfmZmJ6Oho5ObmQhAEhIWF4fHHH7d0WEbTZaThX1PPYOPZV3F4YDjs5Grg2iXo5ixj\noiciogYj6T55uVyOSZMmYfXq1Vi+fDn279+PmzdvWjos48VuQy/lIWSUeWJ7ypMV2zLSgNhtlo2L\niIgaNUkneTc3N7Ru3RoAYGdnBz8/P2RnZ1s4KuOJudno5X4KD7ufwPprk6HW2ui3ExERNRRJJ/m/\nSk9PR1JSEoKCgiwditEEV3cAwMtBX+CO2hubr4812E5ERNQQBFEUJb9UWmlpKZYuXYoxY8agd+/e\nVcrj4uIQFxcHAFixYgXKysrMdm2FQgGNRmPSOTRpqch94x/QpN3C1BMf4GxeBxwZG4nm//c+FD6+\nZopUusxRh00d69B0rEPzYD2aztx1qFQqayyTfJLXaDR49913ERwcjJEjR9bpmNTUVLNd39PTE5mZ\nmSafR3e3Dz7rdhnK7L3QfPLIJjPozlx12JSxDk3HOjQP1qPpzF2Hvr41NxYlPbpeFEWsX78efn5+\ndU7wUiXz8gGmzYXX3Z91OuDnn20xcKDaonEREVHjJek++UuXLuHQoUM4e/Ys5s+fj/nz5yMhIcHS\nYZnFzp12mDjRA7GxKkuHQkREjZSkW/IdOnTAzp07LR1Gg/jb30qwbZsDFi50RUhIBlrY3QJit0HM\nzeZkOUREZBaSTvKNmY0NsG5dDoYO9cL4cS74T68FaF58DgAgApwsh4iITCbpx/WNnb+/Fps3ZyHz\njojn9i+GVvzLr4OT5RARkYnYkrew0NBybB/zPoqv3YJcMFzEhpPlEBGRKdiSl4DuHfLQx+MEAGDn\njVG4WtgSACfLISIi07AlLwXhE4Brl5CXWogVl2ZBrbPF233WISI8zNKRERGRFWNLXgJkXj4Q5iyD\n66MP4btJ76B98wy8/MsCzH27A4qKBEuHR0REVooteYmonCynJYBvlgAffFCAtWsdceqUDfbvz4CN\njaUjJCIia8MkL0EKBbBgQQEefVSNpCQFbGwAbXoasrbHwkOTwvfoiYioTvi4XsL69i3DhAnF0GWk\n4YdX96DvB1FYtasPCo/8AXH1kor58ImIiGrAJG8NYreho/gHwrx/wYd/TscjB77DqiPhSN/6naUj\nIyIiCWOStwJibjb87W8jOmQRYvpMwUOup7H2zxcw9osXIe01BImIyJLYJ28FBFd3VObyHm5n8UXP\nuUgq8kdaqyEQhAiUlgJjx3oiLKwUTz5ZAn9/rUXjJSIiaWBL3hqETwD+Z5Bdq8By9H3lYQBARoYc\nMpmId991xsMPN0N4uCc++sgRt2/z10tE1JSxJW8FZF4+0M1ZVuMqdf7+WsTEZOHGDTm+/dYOe/ao\n8M47znjkETWaN9fhzBkbXL2qQO/eFT8TEVHTwCRvJSrfo6+Nv78Ws2cXYvbsQqSlyeDtXZHQ//tf\nO3z6qSMAoGVLDXr1KkNISBkmTiyGXN7goRMRkYUwyTdSPj73WuyLF+djzJgS/PabEr//rkR8vC0O\nHLDF5MnFAID33nNCaqocHTqUo0MHDdq3L4ePjw6CiZPt6e6upJddVACdgxPf7SciesCY5JsAuRzo\n2rUcXbuW4+9/L4IoAjk5Mn0ST0+X4eBBW3z1lb3+mN691fjvf7MAAJs328PWVkRAgBb+/lo0b669\n7xMAXUYaxNVLgIw0lFduvHYJujnLmOiJiB4QJvkmSBAAd/d7Lf2VK/MA5CE7W4ZLlxS4eFEBB4d7\n7+a9954TcnLuZXWFQsSzzxZjxYo8AMC77zrB1VWHZs108PLSolkzHbz3fw2H/52s527L/n7dDkRE\nZB5M8qTn7q5Dnz5l6NOnzGD7iRN3kJoqx40bcty4oUBKihwdOmgAAGo18OmnDigtNRzJP/2hfljc\n7AcUaeww5tcv4KbMhZtNHtxuyuGe5YRBg0oRGlqOkhIgIUEJFxcdnJzEu390nKv/LnZ5EJEpmOTp\nvmxtgVattGjVSgugrErZn3+mIS9PQEaGHHfuyJCRIUfL0xeBm0CZTokWdreRXe6KiwVtkZPvidyz\nDnBz0yE0tBy3bskxdqxnlWu++24uJk4sxqVLCsya5QZHRx0cHETY2YmwtxcxZUoRQkLKceOGHN99\nZwc7OxEq1b0/oaFl8PLSITdXwM2bcqhUgFIpQqkUYWtb8WVCIfFPP7s8zINflMyD9Wg6S9ShxP+Z\nI2sgCICrqwhXVw3atq3Ypnu0H8TVcXDLSMMXPe8+nr+7pC48fKC9O1+Pr68OX32Vibw8GQoKBBQU\nVPy3e/eKLxNyOdCihQaFhTLk5sqQmiqgqEjAqFElAICrVxV4+23nKjFt356F/v3VOHLEFtOnu1cp\nj4nJRGhoGb76yg5RUa6wsRFhYyNCqazojti2LRtBQRrExqrw6aeOUCgAG5uKLwYKhYjVq3Ph5aXD\nvn0q7N6tglxesV0ur4j59dfzYW8vIj7eFr//roRMhrtlImQyYNasQsjlwOHDSly8aAO5XIQgVOxj\nYyPimWdKgNhtOH7ZGzeKu0MmiBCgg+y2Dqr3EjHs/REAgFOnbJCeLoNMBv0fOzsRvXpV1N/58wrk\n51eMv5DJAKDiS1LnzhVPYq5cUaC4WLhbXtFF4+Ag3v1CByQlyVFWVlFe+cfeXgdf34runhs35NBo\noB/fIQgVx3t6VpSnpsogigJwdzqniuNFuLpW/Hznzr0nQJXnsLcX4egoQhSB7GyZQRkAqFQV96DT\nAXl5QrXlKhWg0wH5yekQ160AMu/oy22vXIVq3hKI7j4oKbl3YOU5lEqxYlEoLaBWVx19qlRWfA60\nWqCs7K/l4t3yit+jVguUl1c5HEplxe9CqwU0mqrlNjb3yrXVzGtlY1MRa3XlglCxwFVlua6aN2Zr\nK688Hqgor5xRU5eRBvHDN4CMOxBlFRfV/nkF+McbVZJUdcdXV17dvVf+Hag8vjqmlN/7e1B7uShW\nX3emlOsy0oAPl0KWefuBfmkXRLHxTYyamppqtnN5enoiMzPTbOdrSiq/tSqKCqBpoG+togiUlgIl\nJTKUlAClpQJKSwUEBGjh5CQiLU2GkyeVKC0VoFZX/KOsVld8SWjWTIdTp2ywe7cKarWA8nIBGk3F\nPlFR+WjeXIf9+1XYssVeX6bRVPx306ZseHnpsGmTPT77zPHuP+gCdLqKf7wOH06Hi4uIFSuc8Mkn\njtBoAJ3uXkK4fj0VCgXw2msu2LzZweCeVCoRV6/ehnblP/GPHU/j29THDco97PJx+s9CAMALL7hh\n3z47g/KAAA1+/TUdADBunAcOH7Y1KO/YsRxxcRkAgJEjPZGYqDQoDw1VIyamYtDlgAFeuHLFsO9k\n4MBSbN2aDQDo2bMZbt82HIU5cmQJPvkkBwDQqZMP8vIMu3KeeaYIq1bl3Y21ObRaw0T6wguFWLYs\nHyUlQFCQL/7X7NkFWLiwAJmZMgQHV/08LVqUj8jIQly/Lkffvs2qlP9f5xV4bmwGLvRZiKFDvauU\nr1mTg6efLsEffyjx5JNVnzJ9/nk2RowoRXy8LSZN8qhS/u9/Z6JfvzLs2qXCSy9V/YK5a1cGevQo\nx3/+Y4e5c92qlMfHp6N9ew0+/9wBS5e6VCn/44878PPTYs0aR7z/ftUvuOfO3Yarq4jly52wbp1T\nlfLKz96iRS748svqP3sAMHu2K775xt6g3F2Zg5NhQwAAfz/xPvbfGWhQbs7P3sCBXrh8uebPXmio\nN1JTDdupTzxRgk8/rfmzN25cMT74IPdurPX/7GVlydCtm/Gfvbc6vYvnAr8y2Cb07g+ZieOUfH2r\nxlqJLXlqMJXv9rs34BclQQDs7AA7u+on+fHx0WH48NIajw8OLkdwcDXNrbuGDSvFsGE1Hz9lSjGm\nTCmusTwqqgBRUQUA7n2712rvtTaWLMnHggX5EEVB37KqbAEIru5Y3HENXmn7GXSiDDoI0IlyyIN7\nAJgEoOL1yJdfLtSfW6cToFTe+97++ut5yM2V3T1nxZeQvw6qXLIkH7m5wt34KspdXe/V5Rtv5CM/\nX9C3yERRgLf3vSbQsmV5KC6uPL7iT4sW98rfeitP39qtPEfr1veacG+/nac/rnKfjh0rym1sgOXL\nc6u0Brt1q/h9OTiIWLYsr0p5aGjFUwxXVx2W9P83kH7boLyn2ymIuSo0a6bD66/n4X9Vnr9FCw0W\nL65a3q5dRXmbNhr885/5BvcGQP8UpGNHDV57Lb/K8b6+FeVdu5YjKqpqeeVTkJ49y7BwYdVyJ6eK\n8r59DcsrY7C1rfif/v3VcHSs2oarbMkOGVKKZs3u/a5EEQZdWE88Uar/XYlH4oDMO7CT3fu7EO67\nD13b5EB4JEy/zcXl3mdn7Nhi9OmjrvbeAGDixCKEhRn+3Wre/F48U6YU6Z/kVAoMvFf+4otFKCgw\nTNJt2977bL38cgFKSw3LO3e+93f91VcLqnx2QkIqyhUKYN68qnVf+YTMzk6stdzFRVelXDwSh+7C\nuSrHiLnZVbaZE1vy98GWvOlYh/Xz1z55vbtdHuwLrRvd56sg/n6wynZztJ6aEtaj6RqyDmtryXNy\ncyKJkt1N6ELv/rDp0gNC7/5M8MaqZt0HePlUbKe6Yz2azkJ1yJb8fbAVajrWoelYh/X3IMaGNAWs\nR9M1VB3W1pJnkr8P/uNqOtah6ViHpmMdmgfr0XTmrkOrHnh38uRJbNy4ETqdDoMHD0ZERISlQyIi\nIrIKku6T1+l0+OKLL7Bo0SKsXr0aR44cwc2bNy0dFhERkVWQdJL/888/4ePjg2bNmkGhUKBv3744\nduyYpcMiIiKyCpJO8tnZ2fDwuDfZhIeHB7KzG/adQiIiosZC8n3ydREXF4e4uDgAwIoVK2odhFAf\n5j5fU8Q6NB3r0HSsQ/NgPZruQdWhpFvy7u7uyMrK0v+clZUFd/eq00SGhYVhxYoVWLFihdljiIqK\nMvs5mxrWoelYh6ZjHZoH69F0D7IOJZ3k27Rpg9u3byM9PR0ajQZHjx5Fz549LR0WERGRVZD043q5\nXI6pU6di+fLl0Ol0GDhwIPz9/S0dFhERkVWQdJIHgB49eqBHjx4Wu35YWNj9d6JasQ5Nxzo0HevQ\nPFiPpnuQddgoZ7wjIiIiiffJExERUf1J/nG9pXA6XdNlZmYiOjoaubm5EAQBYWFhePzxxy0dllXS\n6XSIioqCu7s7RzfXQ1FREdavX48bN25AEATMmDED7dq1s3RYVmX37t2Ij4+HIAjw9/fHzJkzoVQq\nLR2W5K1btw4JCQlwcXHBqlWrAACFhYVYvXo1MjIy4OXlhTlz5sDR0bFBrs+WfDU4na55yOVyTJo0\nCatXr8by5cuxf/9+1mM9ff/99/Dz87N0GFZr48aN6N69O9asWYP333+fdWmk7Oxs7N27FytWrMCq\nVaug0+lw9OhRS4dlFQYMGIBFixYZbIuJiUHXrl2xdu1adO3aFTExMQ12fSb5anA6XfNwc3ND69at\nAQB2dnbw8/PjjIX1kJWVhYSEBAwePNjSoVil4uJiXLhwAYMGDQIAKBQKODg4WDgq66PT6VBWVgat\nVouysjK4ublZOiSr0KlTpyqt9GPHjqF///4AgP79+zdofuHj+mpUN53ulStXLBiR9UtPT0dSUhKC\ngoIsHYrV2bRpEyZOnIiSkhJLh2KV0tPT4ezsjHXr1uH69eto3bo1pkyZApVKZenQrIa7uztGjRqF\nGTNmQKlUIjg4GMHBwZYOy2rl5eXpvyS5uroiLy+vwa7Fljw1uNLSUqxatQpTpkyBvb29pcOxKidO\nnICLi4v+iQgZT6vVIikpCUOHDsV7770HW1vbBn082hgVFhbi2LFjiI6OxieffILS0lIcOnTI0mE1\nCoIgQBCEBjs/k3w16jqdLt2fRqPBqlWr8Nhjj6F3796WDsfqXLp0CcePH0dkZCTWrFmDs2fPYu3a\ntZYOy6p4eHjAw8MDbdu2BQA8/PDDSEpKsnBU1uXMmTPw9vaGs7MzFAoFevfujcuXL1s6LKvl4uKC\nnJwcAEBOTg6cnZ0b7FpM8tXgdLrmIYoi1q9fDz8/P4wcOdLS4Vil8ePHY/369YiOjsYrr7yCLl26\nYPbs2ZYOy6q4urrCw8MDqampACoSVosWLSwclXXx9PTElStXoFarIYoizpw5w8GLJujZsycOHjwI\nADh48CBCQ0Mb7FqcDKcGCQkJ+PLLL/XT6Y4ZM8bSIVmdixcvYsmSJQgICNA/jnr22WctOoOhNTt3\n7hx27drFV+jqITk5GevXr4dGo4G3tzdmzpzZYK8sNVY7d+7E0aNHIZfLERgYiJdeegk2NjaWDkvy\n1qxZg/Pnz6OgoAAuLi4YO3YsQkNDsXr1amRmZjb4K3RM8kRERI0UH9cTERE1UkzyREREjRSTPBER\nUSPFJE9ERNRIMckTERE1UkzyREREjRSTPBERUSPFJE9EVWzfvh179uyp9/EzZsww29Sxr732Gm7c\nuGGWcxE1NVyFjqiJSE9Px+eff44rV65AqVTimWeewcCBA6vsl5+fj4MHD+Jf//qXftu+fftw4MAB\npKSk4JFHHkFkZKTBMdnZ2Vi0aBHWr1+PwsJC5OTk1Hna00mTJun/v6ysDAqFAjJZRftj+vTpGDVq\nFHbs2IF58+bV57aJmjQmeaImonKhoKioKJw8eRIbN26sNskfOHAAISEhUCqV+m1ubm4YM2YMTp06\nhbKysirHJCYmonv37gCAlJQU+Pj4GBxfmy1btuj/PzIyEi+++CK6deum31ZWVobPPvsMubm5cHV1\nrfP9EhEf1xM1CdevX0dhYSFGjhypbyXXtPJVYmIiOnXqZLCtd+/e6NWrF5ycnGo8JiQkBEBFkvf3\n9wcAqNVqfPjhh1i5ciVKS0vrFbtSqUTr1q1x8uTJeh1P1JQxyRM1AZcuXUL79u2h0+lw7do1bN68\nGUOGDKl235SUFPj6+tb53BqNBhcuXNC3vlNSUhAQEID09HS8/vrr8PX1xdy5c6FSqeodv5+fH65f\nv17v44maKj6uJ2oCkpOT0aZNG7z55pu4cOECWrVqhV69elW7b3FxMezs7Op87gsXLqBly5b6YyqT\n8ZtvvokpU6YYLKO5ZcsW9OnTB0FBQUbFb2dnp19/m4jqji15oiagMskvXboUa9euhaOjI7Zt21bt\nvg4ODigpKanzuf/6qF4URdy4cQPHjh3DkCFDqqyTfePGjXqt5V5SUgIHBwejjyNq6pjkiRo5nU6H\nW7duoVWrVpDJZPDx8UH79u1r3L9ly5a4fft2nc//1ySfnp4OAHj99dexe/duXL161WDfnJwcfPnl\nl1i4cCH2799f52vcunULLVu2rPP+RFSBSZ6okUtNTYVarUZiYiJ0Oh2Sk5MRHx+P/v37V7t/SEgI\nzp8/b7BNq9WirKwMOp0OOp0OZWVl0Gq1SE9PR3l5ub51fv36dbRs2RIBAQGYPn06Vq5cqX/Mnp+f\nj8LCQkyYMAFvvfUWfvrppzrFX1ZWhmvXrhmMuCeiumGfPFEjl5SUhBYtWmDLli1Yt24dmjVrhuef\nfx7t2rWrdv9+/fphwYIFKCsr078G98033+Drr7/W7/PLL7/gqaeegrOzs74VD9wbdAcAvXr1QkpK\nCt5//3288cYbSElJFFKDNQAAAQtJREFUwWOPPQZHR0cAgK2tbZ3iP3HiBDp37gx3d/d63T9RUyaI\noihaOggiajhbtmyBo6MjnnzyyTofs337dri4uOCJJ56odb933nkHw4YNQ48ePe57zj179qCgoADP\nPPMMfvvtN6SlpSEiIuK+xy1atAgvvfSS/ssDEdUdkzxRI/fWW29h2LBhNY6mN0VsbCxGjBhRp4lv\nPv74Y9ja2iI3NxcqlQrTp0+HQsGHiUQNiX/DiBq55OTkOk8xa6zw8PA67ztjxowGiYGIasaWPBER\nUSPF0fVERESNFJM8ERFRI8UkT0RE1EgxyRMRETVSTPJERESNFJM8ERFRI8UkT0RE1EgxyRMRETVS\n/w9R6gIiam4ORgAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 576x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dc_oqOdy5JnP",
        "colab_type": "text"
      },
      "source": [
        "From the above code cell, we can qualitatively verify that the values of $\\langle \\hat{H} \\rangle$ attained from our path integral Monte Carlo simulation agree very well with the exact values of $\\langle \\hat{H} \\rangle$ found analytically. A simple way we can quantitatively show the accuracy of our results is by calculating the percent error between our approximate values and the exact values of $\\langle \\hat{H} \\rangle$. This is done in the following cell."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rQwM4-KbB1g-",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 102
        },
        "outputId": "ee8f5919-280c-49a9-86ed-617baabea950"
      },
      "source": [
        "#Calculating percent errors.\n",
        "H_avg_exact_2 = exact_H(betas_1,omega)\n",
        "\n",
        "percent_diff = (np.abs(H_avg_exact_2 - H_avg)/H_avg_exact_2)*100\n",
        "avg = np.sum(percent_diff)/np.size(percent_diff)\n",
        "\n",
        "print(\"Percent difference for each beta value: \\n\", percent_diff,\"\\n\")\n",
        "print(\"Average percent differance: \",avg)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Percent difference for each point: \n",
            " [4.95034444 5.01568641 1.1105442  5.4308563  3.32181376 3.45835661\n",
            " 0.18165948 2.28718827 1.92065033 1.21397751] \n",
            "\n",
            "Average percent differance:  2.8891077314405402\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hrmuBzT_FLzt",
        "colab_type": "text"
      },
      "source": [
        "Even though we are using a relatively small number of points to form our paths ($N = 5$) the results of our simulation are still very accurate, having an average percent error of 2.89% and a maximum percent error of 5.43%. Because the difference between the approximate and exact values of $\\langle \\hat{H} \\rangle$ is very small (that is, the approximate values are very nearly equal to the exact values) we can be fairly confident that our results are trustworthy and accurate.\n",
        "\n",
        "If we were to increase the value of $N$ our results would most likely become even closer to the true values of $\\langle \\hat{H} \\rangle$."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "73uCcP4V9udl",
        "colab_type": "text"
      },
      "source": [
        "## References"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jOYTcyOqmsBl",
        "colab_type": "text"
      },
      "source": [
        "**Textbooks:**  \n",
        "* [*A Guide to Monte-Carlo Simulations in Statistical Physics*](https://el.us.edu.pl/ekonofizyka/images/6/6b/A_guide_to_monte_carlo_simulations_in_statistical_physics.pdf) by David P. Lanadu and Kurt Binder\n",
        "\n",
        "**Papers:**\n",
        "* [*The harmonic oscillator at finite temperature using path integrals*](https://courses.physics.ucsd.edu/2017/Spring/physics142/Lectures/Lecture14/Larsen.pdf) by Å. Larsen and F. Ravndal\n",
        "\n",
        "**Lecture Notes:**\n",
        "* [Professor Samani's notes on MCMC and Feynmann Path Integral](https://drive.google.com/drive/u/1/folders/1nergRg0IK9ks8kbN0rmyAjaaxSTvAQY6)\n",
        "\n",
        "* [Hitoshi Murayama's (UC Berkeley) note on the Feynmann Path Integral](http://hitoshi.berkeley.edu/221A/pathintegral.pdf)"
      ]
    }
  ]
}