QMC_Uniform_reference.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reference implementation of the QMC algorithm on a uniform Cartesian grid (1D version) #\n",
    "\n",
    "### Abstract ###\n",
    "\n",
    "This notebook presents a reference implementation of a quantum Monte-Carlo algorithm designed for sampling the ground state of the molecular Hamiltonian represented on a uniform Cartesian grid. Details of the theory as well as the aplications to model systems in 1D and 3D can be found in the accompaning paper. Here, we provide a simplified version of the code for 1D Hamiltonians and demonstrate its perfomance for several test cases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation of the quantum Monte-Carlo method ##\n",
    "\n",
    "**Random_walk** class provided below implements a grid propagator defined by the following expression for the transition probability:\n",
    "\n",
    "$$ p_d = \\frac{1}{2\\pi} \\int_{ -\\pi}^{\\pi} \\cos(k d) e^{-\\frac{2\\tau}{\\delta^{2}}\\sin^2 (k/2)} dk, $$\n",
    "\n",
    "where d is the length of the step on the 1-D integer lattice ($d = 0, \\pm 1, \\pm 2, \\pm 3,...$), $\\tau$ is the timestep of the diffusion random walk, and $\\delta$ is a uniform grid spacing. The integral is calculated numerically using a general 1-D quadrature subroutine *quad* from *NumPy* for a series of d such that $p_d < 10^{-8}$ (the default value of thresh). The absolute integration accuracy is set to thresh/100 and the relative accuracy is $10^{-3}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Functor class generating grid displacements\n",
    "\n",
    "class Random_walk:\n",
    "        \n",
    "    \n",
    "    def __init__(self, delta, tau, thresh = 1e-8):\n",
    "        \n",
    "        # The names of the input parameters should be self-explanatory\n",
    "        \n",
    "        import numpy as np\n",
    "        from scipy.integrate import quad\n",
    "        \n",
    "        assert delta > 0. and tau > 0.\n",
    "        self.delta = delta\n",
    "        self.tau = tau\n",
    "                \n",
    "        # Calculate the cumulative probabilites\n",
    "        disp = [0]\n",
    "        pd = []\n",
    "        \n",
    "        def kernel(k, d, t):\n",
    "            return 0.5 * np.exp(-2. * t * np.sin(k/2.)**2 / d**2) / np.pi\n",
    "        \n",
    "        terminate = False\n",
    "        while not terminate:\n",
    "            d = disp[-1]\n",
    "            r = quad(kernel, -np.pi, np.pi, (self.delta, self.tau), weight='cos', wvar=d, epsabs=thresh/100.,epsrel=1e-3)\n",
    "            if r[1] > max(thresh/100., r[0] * 1e-3):\n",
    "                print(\"WARNING! Insufficient integration accuracy for displacement {0}!\\n\".format(r[1]))\n",
    "                \n",
    "            pd.append(r[0])\n",
    "            if abs(r[0]) < thresh:\n",
    "                terminate = True\n",
    "            else:\n",
    "                disp.append(d + 1)\n",
    "                \n",
    "        self.pd = np.array(pd[-1:0:-1] + pd)\n",
    "        self.cumm_pd = np.cumsum(self.pd)\n",
    "        self.disp = np.array([-disp[j] for j in range(len(disp) - 1, 0, -1)] + disp)\n",
    "        \n",
    "    def __call__(self, num_samples):\n",
    "        \n",
    "        import numpy as np\n",
    "        \n",
    "        r = np.random.choice(self.disp, size=num_samples, p=self.pd)\n",
    "        \n",
    "        return r\n",
    "        \n",
    "         \n",
    "    def __call_old__(self, num_samples):\n",
    "        \n",
    "        import numpy as np\n",
    "        \n",
    "        r = np.random.rand(num_samples)\n",
    "        #print len(r)\n",
    "        compare_mask = self.cumm_pd[np.newaxis,:] > r[:,np.newaxis]\n",
    "        #print compare_mask.shape\n",
    "        return np.array(self.disp[np.argmax(compare_mask, axis=1)])\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It can be shown that $p_d$ resembles Gaussian distribution in the limit $\\delta \\to 0$ ($\\delta << \\sqrt(\\tau)$). However, the two may differ significantly if $\\delta$ is comparable to $\\sqrt(\\tau)$ as shown by the following example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1520ec0e90>"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmYFOW59/HvzYAgigQFcwRkBgweAw4MOIBbxLggioK7ICoucYJHYpSYI4pxRjzJG5eoIXJixiUxcaJGo0iixqgBlxiVAVEENQICGeEoskoQZeR+/6iapml6umuWmvX3ua6+upanqu6q7ul76nmqnjJ3R0REJJM2jR2AiIg0fUoWIiKSlZKFiIhkpWQhIiJZKVmIiEhWShYiIpKVkoWIiGSlZCEiIlkpWYiISFZtGzuA+tK1a1fPy8tr7DBERJqVefPmferu3bKVazHJIi8vj/Ly8sYOQ0SkWTGzFVHKqRpKRESyUrIQEZGslCxERCSrFtNmIdLSbdu2jYqKCrZu3drYoUgz1KFDB3r27Em7du1qtbyShUgzUVFRQadOncjLy8PMGjscaUbcnbVr11JRUUHv3r1rtQ5VQ4nEqKwM8vKgTZvgvays9uvaunUr++yzjxKF1JiZsc8++9TprFRnFiIxKSuDoiLYsiUYX7EiGAcYP75261SikNqq63dHZxYiMZk6dUeiqLJlSzBdpLlRshCJycqVNZveHJgZP/jBDxLjt912GyUlJQ0aQ3l5OVdccUWtlj366KOb7c27M2fOZPHixY22fSULkZj06lWz6fWtPttLqrRv357HH3+cTz/9tO4rq4XKykoKCwuZPn16o2y7MdUmWdRnzEoWIjH58Y+hY8edp3XsGEyPW1V7yYoV4L6jvaSuCaNt27YUFRVxxx137DLvwgsv5LHHHkuM77nnngDMmTOH4cOHc/bZZ3PggQcyZcoUysrKGDp0KPn5+SxduhSANWvWcMYZZzBkyBCGDBnC3//+dwBKSkooKipixIgRXHDBBcyZM4eTTz4ZgM2bN3PRRReRn5/PgAED+OMf/wjAZZddRmFhIf3796e4uDjrfuXl5XHNNdcwdOhQhg4dypIlSxL7NHnyZL797W9zzTXXsG7dOk499VQGDBjAoYceyttvv52I8fzzz+eYY46hb9++3HPPPYn4jj32WAYPHkx+fj5PPvlkYps33XQTBx10EMcffzzjxo3jtttuA+Cee+5hyJAhDBw4kDPOOIMtW7bw6quvMmvWLH74wx9SUFDA0qVLWbBgAYceeigDBgzgtNNOY/369UBw9nTdddcxfPhwfv7zn0f9aLNz9xbxOuSQQ1ykqXnwQffcXHez4P3BB2u/rsWLF0cum5vrHqSJnV+5ubXfvrv7Hnvs4Rs3bvTc3FzfsGGD33rrrV5cXOzu7hMmTPBHH310p7Lu7rNnz/bOnTv7qlWrfOvWrd69e3e/4YYb3N39zjvv9O9///vu7j5u3Dh/+eWX3d19xYoVftBBB7m7e3FxsQ8ePNi3bNmSWN+oUaPc3f2///u/E8u7u69bt87d3deuXevu7pWVlT58+HB/66233N19+PDhPnfu3DTHK9f/53/+x93dH3jggcT6J0yY4KNGjfLKykp3d580aZKXlJS4u/sLL7zgAwcOTMQ4YMAA37Jli69Zs8Z79uzpH330kW/bts03btzo7u5r1qzxAw44wLdv3+5z5871gQMH+pYtW3zTpk3+jW98w2+99VZ3d//0008TcU2dOtWnT5+e9vjm5+f7nDlz3N39Rz/6UeI4DB8+3C+77LK0n1+67xBQ7hF+Y3U1lEiMxo+v/ZVPdRFne8lee+3FBRdcwPTp09l9990jLTNkyBD2228/AA444ABGjBgBQH5+PrNnzwbg+eef36maZdOmTXz22WcAjB49Ou22nn/+eR5++OHEeJcuXQD4wx/+QGlpKZWVlaxevZrFixczYMCAjDGOGzcu8X7VVVclpp911lnk5OQA8MorryTOXo455hjWrl3Lxo0bARgzZgy77747u+++O9/+9rd54403GDVqFNdddx0vvfQSbdq04aOPPuLjjz/mlVdeSZQHOOWUUxLbe+edd7j++uvZsGEDmzdv5oQTTtgl1o0bN7JhwwaGDx8OwIQJEzjrrLMS888555yM+1obShYiLVCvXkHVU7rp9eHKK69k8ODBXHTRRYlpbdu2Zfv27UBQY/Hll18m5rVv3z4x3KZNm8R4mzZtEvXq27dv5x//+EfapLDHHnukjcPdd7kk9MMPP+S2225j7ty5dOnShQsvvDDS/QXJ60keTt528I94+uVS4zAzysrKWLNmDfPmzaNdu3bk5eWxdevWtOupcuGFFzJz5kwGDhzIb37zG+bMmZM19lTVHa+6iLXNwsxGmtn7ZrbEzKakmT/RzBaa2QIze8XM+oXT88zs83D6AjO7O844RVqauNtL9t57b84++2zuu+++xLS8vDzmzZsHwJNPPsm2bdtqtM4RI0Zw1113JcYXLFhQ42XWr1/Ppk2b2GOPPejcuTMff/wxzzzzTKTtP/LII4n3ww47LG2Zo446irKw4WfOnDl07dqVvfbaCwj2eevWraxdu5Y5c+YwZMgQNm7cyL777ku7du2YPXs2K8IMfuSRR/KnP/2JrVu3snnzZp566qnENj777DP2228/tm3bltgWQKdOnRJnWp07d6ZLly68/PLLAPzud79LnGXEJbYzCzPLAWYAxwMVwFwzm+Xuyc35v3f3u8Pyo4HbgZHhvKXuXhBXfCItWVXV19SpQdVTr15BoqjPKrEf/OAHO/1QX3rppYwZM4ahQ4dy7LHH1vi/2+nTp3P55ZczYMAAKisrOeqoo7j77sz/J15//fVcfvnlHHzwweTk5FBcXMzpp5/OoEGD6N+/P3369OGII46ItP0vvviCYcOGsX37dh566KG0ZUpKSrjooosYMGAAHTt25IEHHkjMGzp0KKNGjWLlypX86Ec/onv37owfP55TTjmFwsJCCgoKOOigg4CgWm706NEMHDiQ3NxcCgsL6dy5MxA0fA8bNozc3Fzy8/MTCWLs2LFceumlTJ8+nccee4wHHniAiRMnsmXLFvr06cOvf/3rSPtZW5bpdKhOKzY7DChx9xPC8WsB3P3/VVN+HHCBu59oZnnAn9394KjbKyws9OZ6/bRIFO+++y7f/OY3GzuMFqnq4Wldu3at1fIlJSXsueeeXH311ZGX2bx5M3vuuSdbtmzhqKOOorS0lMGDB9dq+1Gl+w6Z2Tx3L8y2bJxtFj2AfyWNVwDDUguZ2eXAZGA34JikWb3N7E1gE3C9u78cY6wiIg2qqKiIxYsXs3XrViZMmBB7oqirOJNFuo5IdjmNcfcZwAwzOxe4HpgArAZ6uftaMzsEmGlm/d19004bMCsCigB6NdSdTiLS4ixfvrxOy9fmLvbf//73ddpmQ4uzgbsC2D9pvCewKkP5h4FTAdz9C3dfGw7PA5YCB6Yu4O6l7l7o7oXdumV93riIiNRSnMliLtDXzHqb2W7AWGBWcgEz65s0Ogr4IJzeLWwgx8z6AH2BZTHGKiIiGcRWDeXulWY2CXgWyAHud/dFZjaN4I7BWcAkMzsO2AasJ6iCAjgKmGZmlcBXwER3XxdXrCIiklmsN+W5+9PA0ynTbkga/n41y/0R+GOcsYmISHTqSFBEIsvJyaGgoID+/fszcOBAbr/99sRd23XpOjyKOXPm8OqrrybG7777bn7729/Web15eXmN1otuqvrapziouw+Rlq6kJHjVg9133z1xZ/Unn3zCueeey8aNG7nxxhspLCyksDDr5foZVVZW0rZt+p+lOXPmsOeee3L44YcDMHHixDptqylqyvukMwuRlu7GG2NZ7b777ktpaSl33XUX7r5T1+EvvvgiBQUFFBQUMGjQoMRdyLfccgv5+fkMHDiQKVOCHoBSu9RO11X58uXLufvuu7njjjsoKCjg5ZdfpqSkhNtuu41Vq1YltlVQUEBOTg4rVqyotsvztWvXMmLECAYNGsR3v/vdavtpuu+++zjwwAM5+uijufTSS5k0aRIAf/rTnxg2bBiDBg3iuOOO4+OPPwZIxFPl4IMPZvny5fz73/9m1KhRDBw4kIMPPjjRrciUKVPo168fAwYMSNzMl7yOdF2VQ9B31BVXXMHhhx9Onz59duoWPk46sxCRWuvTpw/bt2/nk08+2Wn6bbfdxowZMzjiiCPYvHkzHTp04JlnnmHmzJm8/vrrdOzYkXXrdlyzsmHDBl588UUAzj33XK666iqOPPJIVq5cyQknnMC7777LxIkTd7pL+oUXXgCge/fuibOdGTNm8OKLL5Kbm1vtem688UaOPPJIbrjhBp566ilKS0t32a9Vq1Zx0003MX/+fDp16sQxxxzDwIEDgaBfp9deew0z49577+WWW27hZz/7WbXH6C9/+Qvdu3dP9P+0ceNG1q1bxxNPPMF7772HmbFhw4Zdljv99NO59NJLgaBbk/vuu4/vfe97AKxevZpXXnmF9957j9GjR3PmmWdG+LTqRslCpCUqKdn5jKKqR9Ti4nqrkqqS7j/zI444gsmTJzN+/HhOP/10evbsyfPPP89FF11Ex7CHw7333jtRPrlL7UxdlWfy97//nXvvvTfRuV5163nppZd4/PHHARg1alSiW/Nkb7zxBsOHD0/EeNZZZ/HPf/4TgIqKCs455xxWr17Nl19+Se/evTPGlZ+fz9VXX80111zDySefzLe+9S0qKyvp0KED3/nOdxg1alTijCxZpq7KTz31VNq0aUO/fv0SZzZxUzWUSENo4OdUU1Ky45lHsGO4nuNYtmwZOTk57LvvvjtNnzJlCvfeey+ff/45hx56KO+9917a7sSrJHc6WNVV+YIFC1iwYAEfffQRnTp1yhjH6tWrueSSS3jkkUcST+jLtJ7q4qiSqc+8733ve0yaNImFCxfyq1/9KtH9eXIX7UBi+oEHHsi8efPIz8/n2muvZdq0abRt25Y33niDM844g5kzZzJy5MhdtnPhhRdy1113sXDhQoqLi3fqZj25y/e4+vdLpWQh0hBiajdoTGvWrGHixIlMmjRplx/fpUuXkp+fzzXXXENhYSHvvfceI0aM4P7770/UvSdXQyWrrqvy5C66k23bto2zzz6bm2++mQMPPDDrepK7GX/mmWcSjyNNNnToUF588UXWr19PZWVl4oFHEFQj9ejRA2CnXmfz8vKYP38+APPnz+fDDz8Egiqtjh07ct5553H11Vczf/58Nm/ezMaNGznppJO4884703bHXl1X5Y1F1VAiLV2EZ1BH9fnnn1NQUMC2bdto27Yt559/PpMnT96l3J133sns2bPJycmhX79+nHjiibRv354FCxZQWFjIbrvtxkknncRPfvKTXZatrqvyU045hTPPPJMnn3ySX/ziF4nyr776KnPnzqW4uDjxvO2nn3662vUUFxczbtw4Bg8ezPDhw9P2K9ejRw+uu+46hg0bRvfu3enXr1+iC/GSkhLOOussevTowaGHHppICmeccQa//e1vKSgoYMiQIYnEtXDhQn74wx/Spk0b2rVrxy9/+Us+++wzxowZk3gQUrpnmlfXVXljia2L8oamLsqlyUltN6hSy3YDdVHesKq6EK+srOS0007j4osv5rTTTmvssOqkLl2UqxpKJC4N1G4g8SgpKaGgoICDDz6Y3r17c+qppzZ2SI1K1VAiImkk3zMhOrMQaRj11G7QUqqNpeHV9bujZCHSEOqh6qlDhw6sXbtWCUNqzN1Zu3YtHTp0qPU6VA0l0kz07NmTiooK1qxZ09ihSDPUoUMHevbsWevllSxEmol27dplvVtYJC6qhhIRkayULEREJCslCxERySrWZGFmI83sfTNbYmZT0syfaGYLzWyBmb1iZv2S5l0bLve+mZ2QuqyIiDSc2JKFmeUAM4ATgX7AuORkEPq9u+e7ewFwC3B7uGw/YCzQHxgJ/G+4PhERaQRxnlkMBZa4+zJ3/xJ4GBiTXMDdNyWN7gFUXUA+BnjY3b9w9w+BJeH6RESkEcR56WwP4F9J4xXAsNRCZnY5MBnYDTgmadnXUpbtEU+YIiKSTZxnFumeLrLLrafuPsPdDwCuAa6vybJmVmRm5WZWrhuVRETiE2eyqAD2TxrvCazKUP5hoKpbx0jLunupuxe6e2G3bt3qGK6IiFQnzmQxF+hrZr3NbDeCButZyQXMrG/S6Cjgg3B4FjDWzNqbWW+gL/BGjLGKiEgGsbVZuHulmU0CngVygPvdfZGZTQPK3X0WMMnMjgO2AeuBCeGyi8zsD8BioBK43N2/iitWERHJTE/KExFpxfSkPBERqTdKFiIikpWShYiIZKVkISIiWSlZiIhIVkoWIiKSVdZkYWZfa4hARESk6YpyZjHPzB4ysxGxRyMiIk1SlGTRF/gtcKmZfWBm08zsgJjjEhGRJiRrsnD37e7+jLufBVwKXAIsMLMXzEzPmBARaQWy9g0VtlmMBy4g6L/pKuAJ4BDgEaB3nAGKiEjji9KR4Fzg98DZ7r4iafprZnZPPGGJiEhTEqXNYqq7FycnCjM7HcDdfxJbZCIi0mRESRZT0kybWt+BiIhI01VtNZSZnQCMBHqY2e1Js/YCtscdmIiINB2Z2iw+Ad4BtgKLkqZ/RvqzDRERaaGqTRbu/ibwppk96O5fNGBMIiLSxGSqhnrI3ccRXPW0y+P03H1wrJGJiEiTkaka6ofh+5m1XbmZjQR+TvAM7nvd/acp8ycD3yF4zvYa4OKqq67M7CtgYVh0pbuPrm0cIiJSN5mqoSrC96W1WbGZ5QAzgOOBCmCumc1y98VJxd4ECt19i5ldBtwCnBPO+9zdC2qzbRERqV+ZqqHWA8nVTxaOG+DuvneWdQ8Flrj7snB9DwNjgESycPfZSeVfA86rUfQiItIgMlVDda3junsA/0oarwCGZSh/CfBM0ngHMysnqKL6qbvPTF3AzIqAIoBevXrVMVwREalOpmqor6qGzWwAcCTBmcUr7r6wuuWSWLrVpi1odh5QCAxPmtzL3VeZWR/gb2a2MLVKzN1LgVKAwsLCtOsWEZG6i/Lwo6nAQwRnCj2Bh8zs2gjrrgD2TxrvCaxKs/7jCO4IH518ia67rwrflwFzgEERtikiIjGI0t3HecAQd5/q7lMJ2iIuiLDcXKCvmfU2s92AscCs5AJmNgj4FUGi+CRpehczax8OdwWOIKmtQ0REGlaUXmdXpJRrCyzLtpC7V5rZJOBZgktn73f3RWY2DSh391nArcCewKNmBjsukf0m8Csz206Q0H6achWViIg0oExXQ91B0MawBVhkZs+G4yOAV6Ks3N2fBp5OmXZD0vBx1Sz3KpAfZRsiIhK/TNVQ7xD0CfUUUAL8g+Dy1mnA32KPTKSpKSlpXdsVSWLuLeMiosLCQi8vL2/sMKQlM4PG+HtprO1Kq2Bm89y9MFu5KI9VPQD4MdAP6FA13d0PrFOEIiLSbES5Guo3wK8J7ps4EfgD8HCMMYk0HSUlwX/2Ft42VDUcd9VQY21XpBpZq6HCU5RDwpvi8sNpL7v7txokwohUDSWxUzWUtED1Vg0FfGHBda1LzWwi8BGwb10DFBGR5iNKsriK4F6IKwjaLjoDF8cZlEiTVFzcurYrkiTy1VBm1pGgt9nP4w2pdlQNJSJSc1GroaL0DTXYzN4E/gl8YGbzwm46RFqFsjLIy4M2bYL3srKWvV2RdKJUQ/0auLLq2RNmdjTBFVID4wtLpGkoK4OiItiyJRhfsSIYBxg/vuVtV6Q6Ua6GetXdD0+Z9nd3PyLWyGpI1VASh7y84Ic6VW4uLF/e8rYrrU+dr4YKn2EB8LqZzSDoptwJHns6u7rlRFqSlStrNr25b1ekOpmqoWakjA9IGtZF39Iq9OqV/j/8uB/M2FjbFalOpiflNamb7kQaw49/vHPbAUDHjsH0lrhdkepEuRqqk5ndYmavha+bzaxTQwQn0tjGj4fS0qCtwCx4Ly2Nv5G5sbYrUp0oDdyPElw2+0A46Xzgm+5+Zsyx1YgauEVEaq4+u/vo6+5nJY3/yMwW1D40ERFpbqL0OrvVzA6rGjGzQ4Gt8YUkIiJNTZRkcRlwr5ktMbMlwD3AxCgrN7ORZvZ+uOyUNPMnm9liM3vbzF4ws9ykeRPM7IPwNSHqDomISP3LWA1lZjlAH3fvb2Z7E7RxrI2y4nDZGcDxQAUw18xmufvipGJvAoXuvsXMLgNuAc4Jt1UMFBJcpjsvXHZ9TXdQRETqLuOZhbt/BVwZDq+LmihCQ4El7r7M3b8keGDSmJT1z3b3qosDXwN6hsMnAM+F21wPPAeMrMG2RUSkHkWphnrWzK40s/3MbK+qV4TlegD/ShqvCKdV5xLgmZosa2ZFZlZuZuVr1qyJEJKIiNRGlKuhvhu+/yBpmgPZ7iW1NNPSXqdrZucRVDkNr8my7l4KlEJw6WyWeEREpJayJgt337+W664AkpftCaxKLWRmxwFTgeHu/kXSskenLDunlnGIiEgdRbmDu72ZXWFmfzCzR8xskpm1j7DuuUBfM+ttZrsBY4FZKeseBPwKGO3unyTNehYYYWZdzKwLMCKcJiIijSBKNdQDwBcEl8wCjAunjc20kLtXmtkkgh/5HOB+d19kZtOAcnefBdxK8MjWR4PHfLPS3Ue7+zozu4kg4QBMc/d1Ndw3ERGpJ1G6+3jb3QekTHvL3ZvUw4/U3YeISM3V22NVgQVmNiRpxYcA/6hLcCIi0rxEqYYaDLxmZh+G472BReFzud3dB8cWnYiINAlRksWY7EVERKQli3Lp7NKGCERERJquKG0WIiLSyilZiIhIVkoWIiKSVZQ7uMeY2btmttHMNpnZZ2a2qSGCExGRpiHK1VA/A05z94VxByMiIk1TlGqoj5UoRERatyhnFnPNrAyYSdBHFABh304iItIKREkW+wDbgdFJ05yUHmRFRKTlinJT3vkNEYiIiDRdUa6G6m5mj5rZ6vD1iJl1b4jgRESkaYjSwP1r4K9AXvh6LpwmIiKtRJRk8XV3v8fdvwhf9wJfjzswERFpOqIki3VmNtZ2OAfQU+tERFqRKMniYuAC4FNgDXA+cEmcQYmISNOSNVm4+3J3P8nd93H3ru5+srt/mG05ADMbaWbvm9kSM5uSZv5RZjbfzCrN7MyUeV+Z2YLwpct0RUQaUbWXzprZD9z9Z2Z2B8F9FTtx98mZVmxmOcAM4HigguDmvlnuvjip2ErgQuDqNKv43N0Lsu+CiIjELdN9FlUPPXqnluseCixx92UAZvYwwVP3EsnC3ZeH87bXchsiItIAqk0W7j4zHFzv7o8nzzOz0yOsuwfwr6TxCmBYDWLrYGblQCXw06R4kuMoAooAevXqVYNVi4hITURp4L4+zbSpEZazNNN2qc7KoJe7FwLnAnea2QG7rMy91N0L3b2wW7duNVi1iIjURKY2ixOAkUAPM7s9adZeBH1FZVMB7J803hNYFTUwd18Vvi8zsznAIHZUjYmISAPKdGbxCUF7xVZgUdLrr8CJEdY9F+hrZr3NbDdgLBE7HzSzLmbWPhzuChxBUluHiIg0rExtFm8Cb5rZg+7+RXXlMixfaWaTgGeBHOB+d19kZtOAcnefZWZDgCeALsApZnaju/cHvgn8Kmz4bkPQZqFkISLSSMw9fTOCmT3k7uPM7E3SXzo7OO7gaqKwsNDLy8sbOwwRkWbFzOaF7cMZZbp09ofh+5kZyoiISCtQbZuFu1eEg6uAZe5e1bj8n8CKuAMTEZGmI8qlsy8Du5vZfsCLwGXA/bFGJSIiTUqUZNHG3bcAZwB3ufspwIB4wxIRkaYkUrIIr1o6F/hzOC0nvpBERKSpiZIsJgM3Ak+5+ztm1oegakpERFqJTFdDAeDufwP+ljS+DPivOIMSEZGmJWuyMLNvEJxd5CWXd/cR8YUlIiJNSdZkATwG3Ac8CHwVbzgiItIURUkW2939F7FHIiIiTVaUBu4nzazIzLqZ2V5Vr9gjExGRJiPKmcV3wvcfJU1zQE8bEhFpJaJcDbV/tjIiItKyRTmzwMwOAvoBHaqmufvv4wpKRESaliiXzl4PjAAOIng2xQnAK4CShYhIKxGlgfsc4NvAanc/HxhIxDMSERFpGaIki8/d/Sug0sw6Af8H9Ik3LBERaUqiJIs3zexrBN2SlwNvAPOjrNzMRprZ+2a2xMympJl/lJnNN7NKMzszZd4EM/sgfE2Isj0REYlHxuokMzOgxN03ADPM7FlgL3fPmizMLAeYARwPVABzzWxWyrO0VwIXAlenLLs3UAwUElymOy9cdn3kPRMRkXqT8czCgwd0/zlpfEmURBEaCixx92Xu/iXwMDAmZf3L3f1tYHvKsicAz7n7ujBBPAeMjLhdERGpZ1Gqod4ws8G1WHcP4F9J4xXhtLiXFRGRelZtNZSZtXX3SuBI4FIzWwr8GzCCk45sCcTSTPOIcUVa1syKgCKAXr10Q7mISFwytVm8AQwGTq3luiuA5Lu/ewKrarDs0SnLzkkt5O6lQClAYWFh1EQkIiI1lClZGIC7L63luucCfc2sN/ARMJbg0axRPAv8xMy6hOMjgGtrGYeIiNRRpmTRzcwmVzfT3W/PtGJ3rzSzSQQ//DnA/e6+yMymAeXuPit8tvcTQBfgFDO70d37u/s6M7uJIOEATHP3dTXZMRERqT+ZkkUOsCfp2w8icfengadTpt2QNDyXoIop3bL3E9zbISIijSxTsljt7tMaLBKRuJWUBK/WpDXus8Qi06WztT6jEGmSbryxsSNoeK1xnyUWmZLFsQ0WhYiINGnVJgs1KEuLUFICZsELdgy35KqZ1rjPEjsLevRo/goLC728vLyxw5CmzAxayPc9sta4z1IjZjbP3QuzlYvS3YdIs1ZWBnl5wXBeXjDe0rXGfZZ46SFG0qKVlUFREWzZAiUUs2JFMA4wfnzjxhaX1rjPEj9VQ0mLlpcHK1bsOj03F5Yvb+hoGkZr3GepPVVDiQArV9ZsekvQGvdZ4qdkIS1adZ0Rt+ROilvjPkv8lCykRfvxj6Fjx52ndewYTG+pWuM+S/yULKR7vtfDAAANtElEQVRFGz8eSkuD+nqz4L20tGU39LbGfZb4qYFbRKQVUwO3iIjUGyULERHJSslCRESyUrIQEZGslCxERCSrWJOFmY00s/fNbImZTUkzv72ZPRLOf93M8sLpeWb2uZktCF93xxmniIhkFltHgmaWA8wAjgcqgLlmNsvdFycVuwRY7+7fMLOxwM3AOeG8pe5eEFd8IiISXZxnFkOBJe6+zN2/BB4GxqSUGQM8EA4/BhxrZnqcq4hIExNnsugB/CtpvCKclraMu1cCG4F9wnm9zexNM3vRzL4VY5wiIpJFnM+zSHeGkHq7eHVlVgO93H2tmR0CzDSz/u6+aaeFzYqAIoBe6iVNRCQ2cZ5ZVAD7J433BFZVV8bM2gKdgXXu/oW7rwVw93nAUuDA1A24e6m7F7p7Ybdu3WLYBRERgXiTxVygr5n1NrPdgLHArJQys4AJ4fCZwN/c3c2sW9hAjpn1AfoCy2KMVUREMoitGsrdK81sEvAskAPc7+6LzGwaUO7us4D7gN+Z2RJgHUFCATgKmGZmlcBXwER3XxdXrCIikpl6nRURacXU66yIiNQbJQsREclKyUJERLJSshARkayULEREJCslCxERyUrJQkREslKyEBGRrJQsREQkKyULaRbKyiAvD9q0Cd7Lyho7opZNx1tSKVlIk1dWBkVFsGIF3OAlrFgRjOsHLB463pKO+oaSJi8vL/jhAnAMCx+LkpsLy5c3Wlgtlo5366K+oaTFWLmyZtOlbnS8JR0lC2nybt+rBMfw8MGKVcO371XSuIG1UDreko6ShTR53WaUsEdHT1SHGM4eHZ1uM0oaN7AWSsdb0lGykCZv/HgoLQ3qzCF4Ly0Npkv90/GWdNTALc1LSUnwkoah493iqYFbmpx6uXZfP1wNqxbHW/dotEyxJgszG2lm75vZEjObkmZ+ezN7JJz/upnlJc27Npz+vpmdEGecEj9du9866HNuwdw9lheQAywF+gC7AW8B/VLK/Bdwdzg8FngkHO4Xlm8P9A7Xk5Npe4cccojXxoMPuufmupdQ7Lm5wXhLXraxtp2b6w7ByyExnJtbs9ilaavr56y/qYaN290dKPcov+lRCtXmBRwGPJs0fi1wbUqZZ4HDwuG2wKeApZZNLlfdqzbJ4sEH3Tt23PmL3bFjtAPeHJdtzG2bpf8RMYsWtzQPdfmc9TfVsHFXaQrJ4kzg3qTx84G7Usq8A/RMGl8KdAXuAs5Lmn4fcGam7dUmWdTlv6DmuGxjbvuOzsU7Fk563dG5OFrg0izU5XPW31TDxl0larKIs83C0kzziGWiLIuZFZlZuZmVr1mzpsYBXrSiJO3NRxetKGmRyzbmtnXtfutQl89Zf1MNG3eNRckotXnRDKqhmuN/BM31vyD3HfWrDrWuX5Wmr7afs/6mmvaZRZzJoi2wjKCBuqqBu39KmcvZuYH7D+Fwf3Zu4F5GDA3czbGusbnWr+6kuLiGC0izVMPPWX9TrbTNIoiBk4B/ErRFTA2nTQNGh8MdgEeBJcAbQJ+kZaeGy70PnJhtW7oaqnlsWyQT/U013auhdAe3iEgrpju4RUSk3ihZiIhIVkoWIiKSlZKFiIhkpWQhIiJZtZirocxsDbCiDqvoSnBTYFOjuGpGcdWM4qqZlhhXrrt3y1aoxSSLujKz8iiXjzU0xVUziqtmFFfNtOa4VA0lIiJZKVmIiEhWShY7lDZ2ANVQXDWjuGpGcdVMq41LbRYiIpKVzixERCSrVpUszOwsM1tkZtvNrDBl3rVmtsTM3jezE6pZvreZvW5mH5jZI2a2WwwxPmJmC8LXcjNbUE255Wa2MCwXew+KZlZiZh8lxXZSNeVGhsdwiZlNaYC4bjWz98zsbTN7wsy+Vk25Bjle2fbfzNqHn/GS8LuUF1csSdvc38xmm9m74ff/+2nKHG1mG5M+3xvijivcbsbPxQLTw+P1tpkNboCY/jPpOCwws01mdmVKmQY5XmZ2v5l9YmbvJE3b28yeC3+HnjOzLtUsOyEs84GZTahzMFG6pm0pL+CbwH8Cc4DCpOn92Pn5GUtJ8/wM4A/A2HD4buCymOP9GXBDNfOWA10b8NiVAFdnKZMTHrs+7HiGSb+Y4xoBtA2HbwZubqzjFWX/gf9i52e4PNIAn91+wOBwuBPBYwNS4zoa+HNDfZ+ifi4Ejzl4huChaIcCrzdwfDnA/xHci9Dgxws4ChgMvJM07RZgSjg8Jd13Htib4DlAewNdwuEudYmlVZ1ZuPu77v5+mlljgIfd/Qt3/5Dg+RpDkwuYmQHHAI+Fkx4ATo0r1nB7ZwMPxbWNGAwFlrj7Mnf/EniY4NjGxt3/6u6V4ehrQM84t5dFlP0fQ/DdgeC7dGz4WcfG3Ve7+/xw+DPgXaBHnNusR2OA33rgNeBrZrZfA27/WGCpu9flht9ac/eXgHUpk5O/Q9X9Dp0APOfu69x9PfAcMLIusbSqZJFBD+BfSeMV7PrHtA+wIemHKV2Z+vQt4GN3/6Ca+Q781czmmVlRjHEkmxRWBdxfzalvlOMYp4sJ/gtNpyGOV5T9T5QJv0sbCb5bDSKs9hoEvJ5m9mFm9paZPWNm/RsopGyfS2N/p8ZS/T9sjXG8AL7u7qsh+EcA2DdNmXo/bm3rsnBTZGbPA/+RZtZUd3+yusXSTEu9TCxKmUgixjiOzGcVR7j7KjPbF3jOzN4L/wuptUxxAb8EbiLY55sIqsguTl1FmmXrfLldlONlZlOBSqCsmtXU+/FKF2qaabF9j2rKzPYE/ghc6e6bUmbPJ6hq2Ry2R80E+jZAWNk+l8Y8XrsBo4Fr08xurOMVVb0ftxaXLNz9uFosVgHsnzTeE1iVUuZTglPgtuF/hOnK1EuMZtYWOB04JMM6VoXvn5jZEwRVIHX68Yt67MzsHuDPaWZFOY71HlfYeHcycKyHFbZp1lHvxyuNKPtfVaYi/Jw7s2s1Q70zs3YEiaLM3R9PnZ+cPNz9aTP7XzPr6u6x9oMU4XOJ5TsV0YnAfHf/OHVGYx2v0Mdmtp+7rw6r5D5JU6aCoF2lSk+CttpaUzVUYBYwNrxSpTfBfwhvJBcIf4RmA2eGkyYA1Z2p1NVxwHvuXpFuppntYWadqoYJGnnfSVe2vqTUE59WzfbmAn0tuGpsN4JT+FkxxzUSuIbgue5bqinTUMcryv7PIvjuQPBd+lt1Ca6+hG0i9wHvuvvt1ZT5j6q2EzMbSvDbsDbmuKJ8LrOAC8Krog4FNlZVwTSAas/uG+N4JUn+DlX3O/QsMMLMuoRVxiPCabUXd2t+U3oR/MhVAF8AHwPPJs2bSnAly/vAiUnTnwa6h8N9CJLIEuBRoH1Mcf4GmJgyrTvwdFIcb4WvRQTVMXEfu98BC4G3wy/rfqlxheMnEVxts7SB4lpCUDe7IHzdnRpXQx6vdPsPTCNIZgAdwu/OkvC71KcBjtGRBFUQbycdp5OAiVXfM2BSeGzeIrhQ4PAGiCvt55ISlwEzwuO5kKSrGGOOrSPBj3/npGkNfrwIktVqYFv423UJQRvXC8AH4fveYdlC4N6kZS8Ov2dLgIvqGovu4BYRkaxUDSUiIlkpWYiISFZKFiIikpWShYiIZKVkISIiWSlZSItmZl+FvYIuCrtmmGxmbcJ5hWY2vZbrXW5mXes32vpnZnlmdm5jxyHNX4u7g1skxefuXgAQdinxe4K7povdvRyIvXv3RpYHnEuw3yK1pjMLaTXc/ROgiKBDRAufSfBnADMbnvRsgjfNrFM4/yULnpOx2MzurjorSWZmM8OO8BYld4ZnwbMt5odnNC+E0/YIO2KcG25nTDj9wnA9fzKzD81sUngW9KaZvWZme4flDjCzv4Tbe9nMDgqn/8aC5z68ambLzKyqp4GfAt8K9+uqOI+vtHANcTekXno11gvYnGbaeuDrJD2TAPgTQad2AHsSnHUfDWwluNM4h6Cb5zPDMssJn8PAjjtodyformIfoBvBneW9U8r8BDgvHP4awd3eewAXEtxp2ylcdiM77ha+g6DzPwju2O0bDg8j6C4Egrv+HyX4B7AfQVfp0EjPqdCr5b1UDSWtUboeOf8O3G5mZcDj7l4Rdv3zhrsvAzCzhwi6zngsZdkrzOy0cHh/gr7FugEvefB8FNy9qrPAEcBoM7s6HO8A9AqHZ3vwvInPzGwjQQKDoJuLAWGvsYcDj9qOR2C0T4pjprtvBxab2dcjHguRSJQspFUxsz7AVwQ9dX6zarq7/9TMniLoM+k1M6vq6Ta1P5ydxs3saIKOHw9z9y1mNocgAViaZQmnn+EpD+Eys2EEfZZV2Z40vp3gb7UNwTNVCqrZveTlY32gkrQ+arOQVsPMuhE8Dvcud0/90T/A3Re6+80Ejd4HhbOGhr3ItgHOAV5JWW1nYH2YKA4iePQnwD+A4WEvxlS1ORD0/Pm9pB5LB0WN34NusT80s7PCZc3MBmZZ7DOCqi2ROlGykJZu96pLZ4Hngb8CN6Ypd6WZvWNmbwGfs+OJe/8gaCR+B/gQeCJlub8Abc3sbYKHQr0G4O5rCBrTHw/X+UhY/iagHfC2mb0TjtfEeOCScJ2LyP7Y2reByrCRXQ3cUmvqdVakGmEV09XufnJjxyLS2HRmISIiWenMQkREstKZhYiIZKVkISIiWSlZiIhIVkoWIiKSlZKFiIhkpWQhIiJZ/X/V73I4bwLYXwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "delta = 0.1\n",
    "tau = 0.015 \n",
    "rw = Random_walk(delta=delta, tau=tau)\n",
    "disp = rw.disp\n",
    "pd = rw.pd\n",
    "#Discretized Gaussian\n",
    "pd_gau = delta / np.sqrt(2. * np.pi * tau)  * np.exp(-delta**2 * disp**2 / (2. * tau))\n",
    "\n",
    "plt.plot(disp, pd, 'bo', label='Numerical propagator')\n",
    "plt.plot(disp, pd_gau,'r+', label='Discretized gaussian')\n",
    "\n",
    "plt.xlabel(\"Displacement\")\n",
    "plt.ylabel(\"Transition probablity\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The instances of the **Random_walk** class can be used to generate random walks on integer lattices by invoking the '__call__' method. In order to test the validity of our implementation we generate a series of *n_samples* and demonstrate that their distribution agrees with the analytic result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of sampled displacements 100000\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1520f9bb10>"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEKCAYAAADjDHn2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt8FdW5//HPQwIEUNAqKodAEnrwAgQCBAQvIIKCWsVWPSKxora/FCsePbS+xHKOAh5Orba1tbUi9nipRqDaatHq8VYQbVEJFUWwKFeJIHJRLnIxl+f3x0zSTbKTPbnsJITv+/WaV/asWWvmyeyd/WRmzawxd0dERKQmrZo6ABERaf6ULEREJCElCxERSUjJQkREElKyEBGRhJQsREQkISULERFJSMlCREQSUrIQEZGEUps6gIZy7LHHemZmZlOHISJySFm6dOk2d++cqF6LSRaZmZkUFhY2dRgiIocUM9sQpZ5OQ4mISEJKFiIikpCShYiIJNRi+ixEWpri4mKKiorYv39/U4ciLUBaWhrp6em0bt26Tu2VLESaqaKiIo488kgyMzMxs6YORw5h7s727dspKioiKyurTuvQaSiRJCpYXkDmLzJpNb0Vmb/IpGB5QeS2+/fv55hjjlGikHozM4455ph6HaUqWYgkScHyAvKfzWdD0QH84QVsKNpP/rP5tUoYShTSUOr7WVKyEEmSqa9OZW/xXnjtv+DjM+C129hbvJepr05t6tBEak3JQiRJNtz8AUxzKPw+eErwc5oH5YeILVu2MH78eHr06MHAgQMZOnQoTz/9dJPFs3DhQv72t7/Vul1mZibbtm2rsc4jjzzCpEmTAJg1axa/+93v6hzH/PnzufPOOwG4+uqreeqpp2oV7//8z/8cNH/aaafVqn0yKFmIJEnX/xwGfQog9cugIPVLyH6c9P8alpTt1ad/JB535+KLL2bYsGGsXbuWpUuXMnfuXIqKihoo4vhKSkqqXVbXZFFbEydO5KqrrqpTHCUlJVx00UVMmTKlztuvnCwa43dORMlCJEl+8q2bSGm/F0rTIHUflKaR0m4vd37zpgbfVkX/yM4NOM6GnRtq3T9S2V/+8hfatGnDxIkTK8oyMjK44YYbACgtLeXmm29m0KBB9O3blwceeAAIvkjPOussLr30Uk4++WTy8vJwdwCWLl3K8OHDGThwIKNHj2bz5s0AnHXWWfzoRz9i+PDh/PKXv+TZZ5/l1FNPpX///owaNYotW7awfv16Zs2axT333ENOTg6vv/46W7du5ZJLLmHQoEEMGjSIv/71rwBs376dc889l/79+/O9732vYvuVPfzww5x44okMHz68oi3AtGnT+OlPfwrAvffeS69evejbty/jxo2LG8fVV1/N5MmTGTFiBLfccstBRykAr7zyCmeeeSYnnngizz33HECVOt/4xjdYuHAhU6ZMYd++feTk5JCXlwfAEUccAQQJ/Oabb6ZPnz5kZ2czb968hPu8wbh7i5gGDhzoIs1N7sgNfsTpjzoT+/kRpz/quSM3RG67cuXKyHUz7slwplFlyrgnow5RB375y1/6TTfdVO3yBx54wO+44w53d9+/f78PHDjQ165d6wsWLPCOHTv6xo0bvbS01IcMGeKvv/66f/XVVz506FD/7LPP3N197ty5fs0117i7+/Dhw/26666rWPeOHTu8rKzM3d0ffPBBnzx5sru733777X733XdX1Lviiiv89ddfd3f3DRs2+Mknn+zu7jfccINPnz7d3d2fe+45B3zr1q0Hxb9p0ybv1q2bf/bZZ37gwAE/7bTT/Prrr6+ynS5duvj+/fvd3f3zzz+PG8eECRP8ggsu8JKSEnd3f/jhhyvWNWHCBB89erSXlpb6hx9+6F27dvV9+/YdVMfd/YILLvAFCxa4u3uHDh0OirV8/qmnnvJRo0Z5SUmJf/rpp96tWzfftGlTtfu8snifKaDQI3zH6j4LkSRa8kp34KpwSp6Pd35cq/K6uP7663njjTdo06YNS5Ys4aWXXuK9996rOB+/c+dOPvroI9q0acPgwYNJT08HICcnh/Xr13PUUUfx/vvvc8455wDBkUmXLl0q1n/55ZdXvC4qKuLyyy9n8+bNfPXVV9XeG/DKK6+wcuXKivldu3axe/duFi1axB//+EcALrjgAo4++ugqbd966y3OOussOnfuXLH9Dz/8sEq9vn37kpeXx8UXX8zFF19c7f657LLLSElJibvs3/7t32jVqhU9e/akR48e/OMf/6h2PTV54403uOKKK0hJSeH4449n+PDhLFmyhI4dO8bd52eccUadthOPkoVIC9C9U3c27Kw6eGj3Tt3rvM7evXvzhz/8oWL+vvvuY9u2beTm5gLBWYlf/epXjB49+qB2CxcupG3bthXzKSkplJSU4O707t2bxYsXx91ehw4dKl7fcMMNTJ48mYsuuoiFCxcybdq0uG3KyspYvHgx7dq1q7IsyqWiUer8+c9/ZtGiRcyfP5877riDFStWJIw/0XbMjNTUVMrKyirKotwD4TWcWoq3zxtSUvsszGyMma0ys9VmVqW3x8wmmtlyM1tmZm+YWa+wPNPM9oXly8xsVjLjFDnUzRw5k/at2x9U1r51e2aOnFnndZ599tns37+f+++/v6Js7969Fa9Hjx7N/fffT3FxMQAffvghX375ZbXrO+mkk9i6dWtFsiguLq72i3fnzp107doVgEcffbSi/Mgjj2T37t0V8+eeey6//vWvK+aXLVsGwLBhwygoCPprXnjhBT7//PMq2zj11FNZuHAh27dvp7i4mCeffLJKnbKyMjZu3MiIESO46667+OKLL9izZ0+VOBJ58sknKSsrY82aNaxdu5aTTjqJzMxMli1bVrGNt99+u6J+69atK/ZrrGHDhjFv3jxKS0vZunUrixYtYvDgwZHjqI+kJQszSwHuA84DegFXlCeDGE+4e7a75wB3AT+PWbbG3XPCaSIiUq287DxmXzibjE4ZGEZGpwxmXzibvOy8Oq/TzHjmmWd47bXXyMrKYvDgwUyYMIGf/OQnAHz3u9+lV69eDBgwgD59+vC9732vxv9m27Rpw1NPPcUtt9xCv379yMnJqfYqn2nTpnHZZZdx5plncuyxx1aUX3jhhTz99NMVHcv33nsvhYWF9O3bl169ejFrVvB/5e23386iRYsYMGAAL730Et27Vz3C6tKlC9OmTWPo0KGMGjWKAQMGVKlTWlrKlVdeSXZ2Nv379+c//uM/OOqoo6rEkchJJ53E8OHDOe+885g1axZpaWmcfvrpZGVlkZ2dzQ9/+MODtp+fn19x+ivWN7/5Tfr27Uu/fv04++yzueuuuzjhhBMSbr8hWE2HNfVasdlQYJq7jw7nbwVw9x9XU/8K4Cp3P8/MMoHn3L1P1O3l5ua6Hn4kLckHH3zAKaec0tRhSAsS7zNlZkvdPTdR22SehuoKbIyZLwrLDmJm15vZGoIji3+PWZRlZu+Y2WtmdmYS4xQRkQSSmSzi9RxVOYxx9/vc/evALcB/hsWbge7u3h+YDDxhZh2rbMAs38wKzaxw69atDRi6iIjESmayKAK6xcynA5tqqD8XuBjA3Q+4+/bw9VJgDXBi5QbuPtvdc909t/zyNxERaXjJTBZLgJ5mlmVmbYBxwPzYCmbWM2b2AuCjsLxz2EGOmfUAegJrkxiriIjUIGn3Wbh7iZlNAl4EUoCH3H2Fmc0guGNwPjDJzEYBxcDnwISw+TBghpmVAKXARHffkaxYRUSkZkm9Kc/dnweer1R2W8zrG6tp9wfgD/GWiYhI49NAgiJSrZSUFHJyciqm8mG3G8KyZct4/vnnE1dspr744gt+85vfNHUYjUbJQqQF2bwZhg+HTz9tmPW1a9eOZcuWVUz1GXa7svomC3c/aLiMxlaXZFFaWpqkaJJPyUKkBbnjDnjjDZgxI3nb2LlzJyeddBKrVq0C4IorruDBBx8E4LrrriM3N5fevXtz++23V7RZsmQJp512Gv369WPw4MHs3LmT2267jXnz5pGTk1Mx1Ha5Rx55hLFjxzJmzBhOOukkpk+fDsD69es55ZRT+P73v8+AAQPYuHEjc+bMITs7mz59+nDLLbdUrOOII47gBz/4AQMGDGDkyJGUX17/4IMPMmjQIPr168cll1xSMYTJmjVrGDJkCIMGDeK2226rGBZ8z549jBw5kgEDBpCdnc2f/vQnAKZMmcKaNWvIycnh5ptvrnH48BEjRjB+/Hiys7Mb/P1oNFGGpj0UJg1RLi1NbYYoT0tzh6pTWlr9YmjVqpX369evYpo7d667u7/00ks+ZMgQnzNnjo8ePbqi/vbt293dvaSkxIcPH+7vvvuuHzhwwLOysvztt992d/edO3d6cXFxlSG6Yz388MN+wgkn+LZt23zv3r3eu3dvX7Jkia9bt87NzBcvXuzu7p988knFMOPFxcU+YsQIf/rpp93dHfDHH3/c3d2nT59esa1t27ZVbGfq1Kl+7733unswRPgTTzzh7u73339/xbDgxcXFvnPnTnd337p1q3/961/3srIyX7dunffu3btiXTUNH96+fXtfu3Ztnd+HhlKfIcp1ZCHSAqxdC+PHQ/twLMH27SEvD9atq996K5+GKh9G/JxzziE7O5vrr7+e3/72txX1f//73zNgwAD69+/PihUrWLlyJatWraJLly4MGjQIgI4dO5KamvjamnPOOYdjjjmGdu3a8a1vfYs33ngDCB7ANGTIECA4YikfZjw1NZW8vDwWLVoEQKtWrSrivfLKKyvav//++5x55plkZ2dTUFBQMZjh4sWLueyyywAYP358RRzuzo9+9CP69u3LqFGj+OSTT9iyZUuVeKsbPhxg8ODB1Q6zfqjQEOUiLUCXLtCxI+zfD2lpwc+OHSFZY8yVlZXxwQcf0K5dO3bs2EF6ejrr1q3jpz/9KUuWLOHoo4/m6quvZv/+/bh7pKHAK4s3rDccPBS412Jsu/L2V199Nc888wz9+vXjkUceYeHChTW2KygoYOvWrSxdupTWrVuTmZkZdzjxmmKpafjyQ4WOLERaiC1bYOJEePPN4GdDdXLHc88993DKKacwZ84crr32WoqLi9m1axcdOnSgU6dObNmyhRdeeAGAk08+mU2bNlX8l717925KSkoSDvP98ssvs2PHDvbt28czzzzD6aefXqXOqaeeymuvvca2bdsoLS1lzpw5DB8+HAgSWvmDmZ544omKBwHt3r2bLl26UFxcXDGMOcCQIUMqnt8xd+7civKdO3dy3HHH0bp1axYsWMCGDcFzQyrH35TDhzcGHVmItBDhg+EAuO++hlln+bOgy40ZM4Zrr72W3/72t7z99tsceeSRDBs2jP/+7/9m+vTp9O/fn969e9OjR4+KL/c2bdowb948brjhBvbt20e7du145ZVXGDFiBHfeeSc5OTnceuutBz0pD+CMM87g29/+NqtXr2b8+PHk5uayfv36g+p06dKFH//4x4wYMQJ35/zzz2fs2LFA8N/8ihUrGDhwIJ06darocL7jjjs49dRTycjIIDs7u+IL/xe/+AVXXnklP/vZz7jgggvo1KkTAHl5eVx44YXk5uaSk5PDySefDMAxxxzD6aefTp8+fTjvvPO46667WLx4Mf369cPMKoYPr+tT8ZqbpA1R3tg0RLm0NIfzEOWPPPIIhYWFBz3YqLaOOOII9uzZE7n+3r17adeuHWbG3LlzmTNnTsWVTy1FfYYo15GFiAiwdOlSJk2ahLtz1FFH8dBDDzV1SM2KjixEmqnD+chCkqO5PvxIROqppfwzJ02vvp8lJQuRZiotLY3t27crYUi9uTvbt28nLS2tzutQn4VIM5Wenk5RURF6CqQ0hLS0NNLT0+vcXslCpJlq3br1IX/Xr7QcOg0lIiIJKVmIiEhCShYiIpJQUpOFmY0xs1VmttrMqjw1xcwmmtlyM1tmZm+YWa+YZbeG7VaZ2ehkxikiIjVLWrIwsxTgPuA8oBdwRWwyCD3h7tnungPcBfw8bNsLGAf0BsYAvwnXJyIiTSCZRxaDgdXuvtbdvwLmAmNjK7j7rpjZDkD5BeVjgbnufsDd1wGrw/WJiEgTSOals12BjTHzRcCplSuZ2fXAZKANcHZM2zcrte2anDBFRCSRZB5ZxHvaSZVbUd39Pnf/OnAL8J+1aWtm+WZWaGaFunFJRCR5kpksioBuMfPpwKYa6s8FLq5NW3ef7e657p7buXPneoYrIiLVSWayWAL0NLMsM2tD0GE9P7aCmfWMmb0A+Ch8PR8YZ2ZtzSwL6Am8ncRYRUSkBknrs3D3EjObBLwIpAAPufsKM5sBFLr7fGCSmY0CioHPgQlh2xVm9ntgJVACXO/upcmKVUREaqbnWYiIHMb0PAsREWkwShYiIpKQkoWIiCSkZCEiIgkpWYiISEJKFiIikpCShYiIJKRkISIiCSlZiIhIQkoWIiKSUMJkYWYdzKxV+PpEM7vIzFonPzQREWkuohxZLALSzKwr8CpwDfBIMoMSEZHmJUqyMHffC3wL+JW7f5PgmdoiInKYiJQszGwokAf8OSxL5uNYRUSkmYmSLG4CbgWeDp8z0QNYkNywRESkOUl4hODurwGvmVmHcH4t8O/JDkxERJqPKFdDDTWzlcAH4Xw/M/tN0iMTEZFmI8ppqF8Ao4HtAO7+LjAsmUGJiEjzEummPHffWKko0vOwzWyMma0ys9VmNiXO8slmttLM3jOzV80sI2ZZqZktC6f5UbYnIiLJEeWqpo1mdhrgZtaGoL/ig0SNzCwFuA84BygClpjZfHdfGVPtHSDX3fea2XXAXcDl4bJ97p5Ti99FRESSJMqRxUTgeqArwZd+TjifyGBgtbuvdfevgLnA2NgK7r4gvIcD4E0gPWrgIiLSeKJcDbWN4B6L2uoKxJ6+KgJOraH+d4AXYubTzKwQKAHudPdn6hCDiIg0gGqThZn9CvDqlrt7ostnLV6zarZ1JZALDI8p7u7um8L7Ov5iZsvdfU2ldvlAPkD37t0ThCMiInVV05FFYT3XXQR0i5lPBzZVrmRmo4CpwHB3P1Be7u6bwp9rzWwh0B84KFm4+2xgNkBubm61iU1EROqn2mTh7o/GzptZx6DYd0dc9xKgp5llAZ8A44DxldbZH3gAGOPun8WUHw3sdfcDZnYscDpB57eIiDSBhH0WZpYLPAwcGczaF8C17r60pnbuXmJmk4AXgRTgoXC4kBlAobvPB+4GjgCeNDOAj939IuAU4AEzKyPohL+z0lVUIiLSiMy95rM3ZvYecL27vx7OnwH8xt37NkJ8keXm5nphYX3PnImIHF7MbKm75yaqF+XS2d3liQLA3d8Aop6KEhGRFiDKTXlvm9kDwByCq5kuBxaa2QAAd/97EuMTEZFmIEqyKL+L+vZK5acRJI+zGzQiERFpdqLclDeiMQIREZHmK8rVUEcBVwGZsfUj3JQnIiItRJTTUM8TjNu0HChLbjgiItIcRUkWae4+OemRiIhIsxXl0tnHzOz/mVkXM/ta+ZT0yEREpNmIkiy+IrjTejGwNJx095scNgqWF5A+fRCW+RrdZgyiYHlBo21782YYPhw+/bTRNikSV5RkMRn4V3fPdPescOqR7MBEmoOC5QXkP5vPJ89dAx+fQdGz15D/bH6jJIyC5QWceNljLHq9lJ6XPtaoSUqksijDfcwHxsU8pKhZ0nAfkgytWu/HS9KqlFvqfsqKq5Y3lDZtSyj+qmqXYus2JXx1IEpXo0g0DTncRymwzMweMLN7y6f6hyjS/PmNPaBPAaR+GRSkfgnZj+M3ZiV1u8f9aGjc7R4/dWhStytSnSj/ojwTTiKHnYz0NmxouwtK0yB1X/Cz7S4y0tsmdbubWApxtvtJzYM9iyRNlDu4H01UR6SlmjlyJlc9eARlA2dB7mwozKfV3nRmjpyZ1O1279SdDV8eBzHbZc8JdO+kJ0JK04hyB3dP4MdAL6DiJK06ueVwkJedB3MLmPrq3Xy882O6593NzJEzg/IkmjlyJvl7r2JvcdhV+I1JtG/dnpkjZyd1uyLViXIa6mGCQQTvAUYA1xD/+doiLVJedl7Sk0O8bQJMfXVqkKQ6dW+UJCVSnShXQy1194Fmttzds8Oy1939zEaJMCJdDSUiUntRr4aKcmSx38xaAR+Fj0n9BDiuvgGKiMihI8qlszcB7YF/BwYC3wYmRFm5mY0xs1VmttrMpsRZPtnMVprZe2b2qpllxCybYGYfhVOk7YmISHJEuRpqSfhyj5l9BzjC3XclamdmKcB9wDlAEbDEzOa7+8qYau8Aue6+18yuA+4CLg/HnrodyCV4wNLSsO3ntfnlRESkYSQ8sjCzJ8yso5l1AFYCq8zs5gjrHgysdve17v4VMBcYG1vB3RfE3Bn+JpAevh4NvOzuO8IE8TIwJtqvJCIiDS3Kaahe4ZHExQTPtuhOcCoqka7Axpj5orCsOt8BXqhjWxERSaIoHdytzaw1QbL4tbsXm1nNl1AF4l1eG7edmV1JcMppeG3amlk+kA/QvbtuVhIRSZYoRxYPAOuBDsCisBM6YZ8FwdFAt5j5dGBT5UpmNgqYClzk7gdq09bdZ7t7rrvndu7cOUJIIiJSFwmThbvf6+5d3f18D27K+Jjg5rxElgA9zSzLzNoA44D5sRXMrD9BMrrI3T+LWfQicK6ZHW1mRwPnhmUiItIEaj3WcZgwSiLUKwnvy3gRSAEecvcVZjYDKHT3+QQPVToCeNLMAD5294vcfYeZ3UGQcABmuPuO2sYqIiINI+Ed3IcK3cEtIlJ79X6ehZldFv5M7sD9IiLS7NXUZ3Fr+PMPjRGIiIg0XzX1WWw3swVAVvho1YO4+0XJC0tERJqTmpLFBcAA4DHgZ40TjoiINEfVJotwiI43zew0d99qZkcGxb6n8cITEZHmIMpNeceb2TvA+8BKM1tqZn2SHJeIiDQjUZLFbGCyu2e4e3fgB2GZiIgcJqIkiw7uvqB8xt0XEgz9ISIih4kod3CvNbP/IujoBrgSWJe8kEREpLmJcmRxLdAZ+GM4HQtck8ygRESkeYnypLzPCR6pKiIih6koRxYiInKYU7IQEZGElCxERCShhH0WZtYZ+H9AZmx9d782eWGJiEhzEuXS2T8BrwOvAKXJDUdERJqjKMmivbvfkvRIRESk2YrSZ/GcmZ2f9EhERKTZipIsbiRIGPvNbHc47YqycjMbY2arzGy1mU2Js3yYmf3dzErM7NJKy0rNbFk4VXmehoiINJ4oN+UdWZcVm1kKcB9wDlAELDGz+e6+Mqbax8DVwA/jrGKfu+fUZdsiItKwovRZYGYXAcPC2YXu/lyEZoOB1e6+NlzHXGAsUJEs3H19uKysFjGLiEgjS3gayszuJDgVtTKcbgzLEukKbIyZLwrLokozs0Ize9PMLq4mtvywTuHWrVtrsWoREamNKEcW5wM57l4GYGaPAu8AVfogKrE4ZV6L2Lq7+yYz6wH8xcyWu/uag1bmPpvw2Rq5ubm1WbeIiNRC1Du4j4p53SlimyKgW8x8OrApYlvcfVP4cy2wEOgfta2IiDSsKEcWPwbeMbMFBEcLw4BbI7RbAvQ0syzgE2AcMD5KUGZ2NLDX3Q+Y2bHA6cBdUdqKiEjDi3I11BwzWwgMIkgWt7j7pxHalZjZJOBFIAV4yN1XmNkMoNDd55vZIOBp4GjgQjOb7u69gVOAB8KO71bAnZWuohIRkUZk7vFP9ZvZye7+DzMbEG+5u/89qZHVUm5urhcWFjZ1GCIihxQzW+ruuYnq1XRkMRnIB34WZ5kDZ9cxNhEROcRUmyzcPT98eZ67749dZmZpSY1KRESalShXQ/0tYpmIiLRQ1R5ZmNkJBDfRtTOz/vzzvomOQPtGiE1ERJqJmvosRhOM25RO0G9Rnix2AT9KblgiItKc1NRn8SjwqJld4u5/aMSYRESkmYnSZzHQzCru4Dazo83sv5MYk4iINDNRksV57v5F+Yy7f04wXpSIiBwmoiSLFDNrWz5jZu2AtjXUFxGRFibK2FCPA6+a2cMEN+NdCzya1KhERKRZiTI21F1mthwYSXBF1B3u/mLSIxMRkWYj0pPy3P0F4IUkxyIiIs1UlCflDTGzJWa2x8y+MrNSM9vVGMGJiEjzEKWD+9fAFcBHQDvgu8CvkhmUiIg0L1FPQ602sxR3LwUeNjONDSUichiJkiz2mlkbYJmZ3QVsBjokNywREWlOopyG+nZYbxLwJcFztS9JZlAiItK81HhkYWYpwEx3vxLYD0xvlKhERKRZqfHIIuyj6Byehqo1MxtjZqvMbLWZTYmzfJiZ/d3MSszs0krLJpjZR+E0oS7bFxGRhhGlz2I98Fczm09wGgoAd/95TY3Co5L7gHOAImCJmc1395Ux1T4mGAb9h5Xafg24HcgluGt8adj28wjxiohIA4uSLDaFUyvgyFqsezCw2t3XApjZXGAsUJEs3H19uKysUtvRwMvuviNc/jIwBphTi+2LiEgDqelJeY+5+7eBL9z9l3VYd1dgY8x8EXBqPdp2jRNjPpAP0L179zqEKCIiUdTUZzHQzDKAa8NnWHwtdoqwbotT5hHjitTW3We7e66753bu3DniqkVEpLZqOg01C/g/oAewlIO/wD0sr0kRwWW25dIJTmdFUQScVantwohtRUSkgVV7ZOHu97r7KcBD7t7D3bNipkSJAmAJ0NPMssKrqcYB8yPG9SJwbnhEczRwblgmIiJNIOFNee5+XV1W7O4lBDfyvQh8APze3VeY2QwzuwjAzAaZWRFwGfCAma0I2+4A7iBIOEuAGeWd3SIi0vjMPWo3QvOWm5vrhYWFTR2GiMghxcyWuntuonpRhvsQEZHDnJKFiIgkpGQhIiIJKVmIiEhCShbS4hUsLyB9+iAs8zW6zRhEwfKCpg6pUWzeDMOHw6efNnUk0hIoWUiLVrC8gPxn8/nkuWvg4zMoevYa8p/Nb/EJo2B5ASde9hiLXi+l56WPtfjfV5JPl85Ki9aq9X68JK1KuaXup6y4anlL0KZtCcVfVR2coXWbEr46EOlJynIY0aWzIoDf2AP6FEBqOLp+6peQ/TgTI+DIAAAPd0lEQVR+Y1bTBpZEx/1oaNzf+fipQ5s2MDmkKVlIi5aR3gba7oLSNEjdF/xsu4uM9LZNHVrSbGJp3N/5E1/a1KHJIUzJQlq0mSNn0mrfv8DAWfDdITBwFq32dmXmyJlNHVrSdO/UHb487qDfmT3HB+UidaQTmNKi5WXnwdwCpr56Nx/v/JjueXczc+TMoLyFmjlyJvl7r2Jv8d6g4BuTaN+6PTNHzm7awOSQpg5ukRaoYHkBU1+dGiTITt1bfIKUuovawa1kISJyGNPVUCIi0mCULEREJCElCxERSUjJQkREEkpqsjCzMWa2ysxWm9mUOMvbmtm8cPlbZpYZlmea2T4zWxZOs5IZp4iI1Cxp91mYWQpwH3AOUAQsMbP57r4yptp3gM/d/V/NbBzwE+DycNkad89JVnwiIhJdMo8sBgOr3X2tu38FzAXGVqozFng0fP0UMNLMLIkxiYhIHSQzWXQFNsbMF4Vlceu4ewmwEzgmXJZlZu+Y2WtmdmYS4xQRkQSSOdxHvCOEyncAVldnM9Dd3beb2UDgGTPr7e67Dmpslg/kA3TvrnFvRESSJZlHFkVAt5j5dGBTdXXMLBXoBOxw9wPuvh3A3ZcCa4ATK2/A3We7e66753bu3DkJv4KIiEByk8USoKeZZZlZG2AcML9SnfnAhPD1pcBf3N3NrHPYQY6Z9QB6AmuTGKuIiNQgaaeh3L3EzCYBLwIpwEPuvsLMZgCF7j4f+F/gMTNbDewgSCgAw4AZZlYClAIT3X1HsmIVEZGaaSBBEZHDmAYSFBGRBqNkISIiCSlZiIhIQkoWIiKSkJKFiIgkpGQhIiIJKVmIiEhCShYiIpKQkoWIiCSkZCGHhILlBaRPH4Rlvka3GYMoWF7Q1CG1aNrfUlkyhygXaRAFywvIfzafvc/dDR+fQdGz15Cfmg9AXnZeE0fX8mh/SzwaG0qavVat9+MlaVXKLXU/ZcVVy6V+tL8PLxobSloMv7EH9CmA1C+DgtQvIftx/Maspg2shdL+lniULKTZy0hvA213QWkapO4LfrbdRUZ626YOrUXS/pZ4lCyk2Zs5ciat9v0LDJwF3x0CA2fRam9XZo6c2dShtUja3xKPkoU0e3nZefxu7h4y8u7GTlhORt7d/G7uHnW2Jon2t8SjDm4RkcOYOril2dG1+4cHvc8tU1KThZmNMbNVZrbazKbEWd7WzOaFy98ys8yYZbeG5avMbHQy45TkK792/5PnrvnntfvP5uuLpIXR+9yCuXtSJiAFWAP0ANoA7wK9KtX5PjArfD0OmBe+7hXWbwtkhetJqWl7AwcO9LratMl92DD3zZsPjbaPv/e4d52W62Qs9PTpuf74e483Wvu6trXUfQ5eZbLUfbWKXZq3+r7PTfHZrG/bhmjfVN9B7u5AoUf5To9SqS4TMBR4MWb+VuDWSnVeBIaGr1OBbYBVrhtbr7qpPsniuuvcW7UKfjb3to+/97i3n9neyb3PsRIn9z5vP7N95A9nfdrXpy0/6OL0edxJ3RN8gaTucbIfc35wQqS45dBQn/e5qT6bTfk3Va6pvoPcm0eyuBT4bcz8t4FfV6rzPpAeM78GOBb4NXBlTPn/ApfWtL26JIu0tKr/AUFQ3lzb1vc/t/q0r0/bjHsynIG/Cf6YUvdW/FFl3JMRKW45NNTnfW6qz2ZT/k011fdIrKjJIpl9FhanrPKlV9XVidIWM8s3s0IzK9y6dWutA1y7FsaPh/btg/n27SEvD9ata75t63t3bX3a16etrt0/PNTnfW6qz2ZT/k011fdIXSQzWRQB3WLm04FN1dUxs1SgE7AjYlvcfba757p7bufOnWsdYJcu0LEj7N8PaWnBz44d4YQTmm/b+t5dW5/29Wmra/cPD/V5n5vqs9mUf1NN9T1SF8lMFkuAnmaWZWZtCDqw51eqMx+YEL6+FPhLeFg0HxgXXi2VBfQE3k5GkFu2wMSJ8Oabwc9PP23ebev7H3p92td323nZeay/aT1lt5ex/qb1ShQtVF3f56b6bDbl3xQ03XdQrUU5V1XXCTgf+JCgL2JqWDYDuCh8nQY8CawmSAY9YtpODdutAs5LtK36dHAfah5/73HPuCfDbZp5xj0Zdbpyo67t67ttkZo01WezKf+mmhoR+yx0B7eIyGFMd3CLiEiDUbIQEZGElCxERCQhJQsREUlIyUJERBJqMVdDmdlWYEM9VnEswdhUzY3iqh3FVTuKq3ZaYlwZ7p7wruYWkyzqy8wKo1w+1tgUV+0ortpRXLVzOMel01AiIpKQkoWIiCSkZPFPs5s6gGoortpRXLWjuGrnsI1LfRYiIpKQjixERCShwypZmNllZrbCzMrMLLfSslvNbLWZrTKz0dW0zzKzt8zsIzObFw693tAxzjOzZeG03syWVVNvvZktD+slfQRFM5tmZp/ExHZ+NfXGhPtwtZlNaYS47jazf5jZe2b2tJkdVU29RtlfiX7/cNj9eeHyt8wsM1mxxGyzm5ktMLMPws//jXHqnGVmO2Pe39uSHVe43RrfFwvcG+6v98xsQCPEdFLMflhmZrvM7KZKdRplf5nZQ2b2mZm9H1P2NTN7OfweetnMjq6m7YSwzkdmNiFenVqJMjRtS5mAU4CTgIVAbkx5L+BdoC2QRTA0ekqc9r8HxoWvZwHXJTnenwG3VbNsPXBsI+67acAPE9RJCfddD6BNuE97JTmuc4HU8PVPgJ801f6K8vsD3wdmha/HAfMa4b3rAgwIXx9J8NiAynGdBTzXWJ+nqO8LwWMOXiB4euYQ4K1Gji8F+JTgXoRG31/AMGAA8H5M2V3AlPD1lHifeeBrwNrw59Hh66PrE8thdWTh7h+4+6o4i8YCc939gLuvI3i+xuDYCmZmwNnAU2HRo8DFyYo13N6/AXOStY0kGAysdve17v4VMJdg3yaNu7/k7iXh7JsET1VsKlF+/7EEnx0IPksjw/c6adx9s7v/PXy9G/gA6JrMbTagscDvPPAmcJSZdWnE7Y8E1rh7fW74rTN3X0Tw9NBYsZ+h6r6HRgMvu/sOd/8ceBkYU59YDqtkUYOuwMaY+SKq/jEdA3wR88UUr05DOhPY4u4fVbPcgZfMbKmZ5ScxjliTwlMBD1Vz6BtlPybTtQT/hcbTGPsryu9fUSf8LO0k+Gw1ivC0V3/grTiLh5rZu2b2gpn1bqSQEr0vTf2ZGkf1/7A1xf4CON7dN0PwjwBwXJw6Db7fUuvTuDkys1eAeE+hneruf6quWZyyypeJRakTScQYr6Dmo4rT3X2TmR0HvGxm/wj/C6mzmuIC7gfuIPid7yA4RXZt5VXEaVvvy+2i7C8zmwqUAAXVrKbB91e8UOOUJe1zVFtmdgTwB+Amd99VafHfCU617An7o54heJxxsiV6X5pyf7UBLgJujbO4qfZXVA2+31pcsnD3UXVoVgR0i5lPBzZVqrON4BA4NfyPMF6dBonRzFKBbwEDa1jHpvDnZ2b2NMEpkHp9+UXdd2b2IPBcnEVR9mODxxV23n0DGOnhCds462jw/RVHlN+/vE5R+D53oupphgZnZq0JEkWBu/+x8vLY5OHuz5vZb8zsWHdP6jhIEd6XpHymIjoP+Lu7b6m8oKn2V2iLmXVx983hKbnP4tQpIuhXKZdO0FdbZzoNFZgPjAuvVMki+A/h7dgK4ZfQAuDSsGgCUN2RSn2NAv7h7kXxFppZBzM7svw1QSfv+/HqNpRK54m/Wc32lgA9LbhqrA3BIfz8JMc1BriF4Lnue6up01j7K8rvP5/gswPBZ+kv1SW4hhL2ifwv8IG7/7yaOieU952Y2WCC74btSY4ryvsyH7gqvCpqCLCz/BRMI6j26L4p9leM2M9Qdd9DLwLnmtnR4Snjc8Oyukt2b35zmgi+5IqAA8AW4MWYZVMJrmRZBZwXU/488C/h6x4ESWQ18CTQNklxPgJMrFT2L8DzMXG8G04rCE7HJHvfPQYsB94LP6xdKscVzp9PcLXNmkaKazXBudll4TSrclyNub/i/f7ADIJkBpAWfnZWh5+lHo2wj84gOAXxXsx+Oh+YWP45AyaF++ZdggsFTmuEuOK+L5XiMuC+cH8uJ+YqxiTH1p7gy79TTFmj7y+CZLUZKA6/u75D0Mf1KvBR+PNrYd1c4Lcxba8NP2ergWvqG4vu4BYRkYR0GkpERBJSshARkYSULEREJCElCxERSUjJQkREElKykBbNzErDUUFXhEMzTDazVuGyXDO7t47rXW9mxzZstA3PzDLNbHxTxyGHvhZ3B7dIJfvcPQcgHFLiCYK7pm9390Ig6cO7N7FMYDzB7y1SZzqykMOGu38G5BMMiGjhMwmeAzCz4THPJnjHzI4Mly+y4DkZK81sVvlRSSwzeyYcCG9F7GB4Fjzb4u/hEc2rYVmHcCDGJeF2xoblV4fredbM1pnZpPAo6B0ze9PMvhbW+7qZ/V+4vdfN7OSw/BELnvvwNzNba2blIw3cCZwZ/l7/kcz9Ky1cY9wNqUlTU03AnjhlnwPHE/NMAuBZgkHtAI4gOOo+C9hPcKdxCsEwz5eGddYTPoeBf95B245guIpjgM4Ed5ZnVarzP8CV4eujCO727gBcTXCn7ZFh2538827hewgG/4Pgjt2e4etTCYYLgeCu/ycJ/gHsRTBUOjTRcyo0tbxJp6HkcBRvRM6/Aj83swLgj+5eFA7987a7rwUwszkEQ2c8Vantv5vZN8PX3QjGFusMLPLg+Si4e/lggecCF5nZD8P5NKB7+HqBB8+b2G1mOwkSGATDXPQNR409DXjS/vkIjLYxcTzj7mXASjM7PuK+EIlEyUIOK2bWAyglGKnzlPJyd7/TzP5MMGbSm2ZWPtJt5fFwDpo3s7MIBn4c6u57zWwhQQKwOG0Jyy/xSg/hMrNTCcYsK1cWM19G8LfaiuCZKjnV/Hqx7ZP6QCU5/KjPQg4bZtaZ4HG4v3b3yl/6X3f35e7+E4JO75PDRYPDUWRbAZcDb1RabSfg8zBRnEzw6E+AxcDwcBRjyvscCEb+vCFmxNL+UeP3YFjsdWZ2WdjWzKxfgma7CU5tidSLkoW0dO3KL50FXgFeAqbHqXeTmb1vZu8C+/jnE/cWE3QSvw+sA56u1O7/gFQze4/goVBvArj7VoLO9D+G65wX1r8DaA28Z2bvh/O1kQd8J1znChI/tvY9oCTsZFcHt9SZRp0VqUZ4iumH7v6Npo5FpKnpyEJERBLSkYWIiCSkIwsREUlIyUJERBJSshARkYSULEREJCElCxERSUjJQkREEvr/R85Q/p32prUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Check whether the class produces correct displacement statistics\n",
    "n_samples = 100000\n",
    "displacements = rw(n_samples)\n",
    "print(\"Number of sampled displacements {0}\".format(displacements.shape[0]))\n",
    "assert displacements.shape[0] == n_samples\n",
    "#print(displacements)\n",
    "\n",
    "hist = {}\n",
    "for i in displacements:\n",
    "    if i in hist:\n",
    "        hist[i] += 1\n",
    "    else:\n",
    "        hist[i] = 1\n",
    "\n",
    "# Normalize and compare with the exact distribution\n",
    "plt.plot(hist.keys(), np.array(hist.values())/float(n_samples), 'go', label='Generated distribution')\n",
    "plt.plot(disp, pd, 'b*', label='Exact propagator')\n",
    "\n",
    "plt.xlabel(\"Displacement\")\n",
    "plt.ylabel(\"fraction of samples\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Walker_ensemble** class implements the core of the grid-based quantum Monte-Carlo algorithm. It uses instances of **Random_walk** class to generate diffusion random walks on spatial grids. The imaginary time propagation is executed by invoking a *one_time_step* method which in turn calls *_diffuse_and_branch, _update_phases* and *_annihilate_walkers* methods. The energy is then 'measured' by the growth estimator. The class also provides a functionality to implement constraints which can for example be used to perform fixed-node simulations and impose spatial symmetry. Present code is restricted to 1-D systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Walker_ensemble:\n",
    "    \n",
    "    def __init__(self, random_walk, potential, n_walkers, n_alpha, n_beta, constraint=None, pop_dyn=False, S=0.0):\n",
    "        \n",
    "        import numpy as np\n",
    "        from sympy.combinatorics.permutations import Permutation\n",
    "        import itertools as it\n",
    "        \n",
    "        L = 5. # Initial walker population will be (uniformly) distributed in the interval [-L, L]\n",
    "        \n",
    "        # Are we performing a regular run or recording population dynamics?\n",
    "        self.pop_dyn = pop_dyn\n",
    "        \n",
    "        # System specification\n",
    "        self.disp_generator = random_walk # An instance of random walk class\n",
    "        self.V = potential # A callable python object capable of operating on a walker (=set of particle coordinates)\n",
    "        self.S = S # initial guess for the growth estimator\n",
    "        self.apply_constraint = False\n",
    "        self.constraint = None\n",
    "        \n",
    "        if constraint != None:\n",
    "            self.constraint = constraint # A callable python object capable of operating on a walker and returning True/False\n",
    "            self.apply_constraint = True\n",
    "        \n",
    "        # Walker positions and weights; \n",
    "        # The former are expressed by the coordinates on the integer lattice.\n",
    "        # The latter should ALWAYS be \\pm 1!\n",
    "        self.n_alpha = n_alpha # number of alpha electrons\n",
    "        self.n_beta = n_beta   # number of beta electrons\n",
    "        self.lookup_a = {} # all permutations of alpha electrons with corresponding parities\n",
    "        self.lookup_b = {} # same for beta\n",
    "        \n",
    "        pa = it.permutations(np.arange(n_alpha))\n",
    "        pb = it.permutations(np.arange(n_beta))\n",
    "        \n",
    "        for i in pa:\n",
    "            parity = 1 if Permutation(i).is_even else -1\n",
    "            self.lookup_a[i] = parity\n",
    "        for i in pb:\n",
    "            parity = 1 if Permutation(i).is_even else -1\n",
    "            self.lookup_b[i] = parity\n",
    "        \n",
    "        self.walker_weights = np.ones(n_walkers, dtype=int)\n",
    "        self.walker_coords = np.zeros((n_walkers, n_alpha + n_beta), dtype=int)\n",
    "        \n",
    "        # Timing info for reference purposes\n",
    "        self.time = {'total': 0.0, 'branching': 0.0, 'phase_update': 0.0, 'annihilation' : 0.0, 'n_steps' : 0}\n",
    "                \n",
    "        # Initialize walker ensemble \n",
    "        # and make sure that the walkers conform to the Pauli principle and don't violate constraints\n",
    "        d = self.disp_generator.delta\n",
    "        grid_init = np.arange(-L/d, L/d + 1, dtype=int)\n",
    "        kwargs = {'size' : n_alpha, 'replace' : False} # replace = False makes sure that the Pauli principle is not violated\n",
    "        alpha_coord = np.apply_along_axis(np.random.choice, 1, np.broadcast_to(grid_init, (n_walkers, grid_init.shape[0])), **kwargs)\n",
    "        alpha_coord = np.apply_along_axis(np.sort, 1, alpha_coord)\n",
    "        kwargs = {'size' : n_beta, 'replace' : False}\n",
    "        beta_coord = np.apply_along_axis(np.random.choice, 1, np.broadcast_to(grid_init, (n_walkers, grid_init.shape[0])), **kwargs) \n",
    "        beta_coord = np.apply_along_axis(np.sort, 1, beta_coord)\n",
    "        \n",
    "        # Check if the coordinates of alpha and beta particles indeed satisfy the Pauli principle\n",
    "        all_distinct = lambda ar : len(ar) == len(np.unique(ar))\n",
    "        a = np.apply_along_axis(all_distinct, 1, alpha_coord)\n",
    "        b = np.apply_along_axis(all_distinct, 1, beta_coord)\n",
    "        \n",
    "        if self.apply_constraint:\n",
    "            constr_a = np.apply_along_axis(self.constraint, 1, alpha_coord)\n",
    "            constr_b = np.apply_along_axis(self.constraint, 1, beta_coord)\n",
    "            a = np.logical_and(a, constr_a)\n",
    "            b = np.logical_and(b, constr_b)\n",
    "            \n",
    "        \n",
    "        alpha_ok = True if n_alpha == 0 else np.all(a)\n",
    "        beta_ok = True if n_beta == 0 else np.all(b)\n",
    "        \n",
    "        assert alpha_ok and beta_ok\n",
    "        \n",
    "        # Combine the coordinates of alpha and beta electrons (the resulting list may contain duplicates)\n",
    "        self.walker_coords = np.hstack((alpha_coord, beta_coord))\n",
    "        \n",
    "    def one_time_step(self):\n",
    "        \n",
    "        # High level grid QMC driver \n",
    "        \n",
    "        import numpy as np\n",
    "        import time\n",
    "        \n",
    "        start = time.clock()\n",
    "        \n",
    "        n_walkers_old = self.get_number_of_walkers()\n",
    "                \n",
    "        self._diffuse_and_branch()\n",
    "        \n",
    "        self.time['branching'] += time.clock() - start\n",
    "        \n",
    "        t1 = time.clock()\n",
    "        \n",
    "        self._update_phases()\n",
    "                \n",
    "        self.time['phase_update'] += time.clock() - t1\n",
    "        \n",
    "        t2 = time.clock()\n",
    "        \n",
    "        self._annihilate_walkers()\n",
    "        \n",
    "        self.time['annihilation'] += time.clock() - t2\n",
    "        self.time['total'] += time.clock() - start\n",
    "        self.time['n_steps'] += 1\n",
    "        \n",
    "        n_walkers_new = self.get_number_of_walkers()\n",
    "        #print(\"Number of walkers at the end of the step is {0}\".format(n_walkers_new))\n",
    "        \n",
    "        if not self.pop_dyn:\n",
    "            self.S = self.S + 1./ self.disp_generator.tau * np.log(n_walkers_old/float(n_walkers_new))\n",
    "    \n",
    "    def _diffuse_and_branch(self):\n",
    "        \n",
    "        import numpy as np\n",
    "        \n",
    "        # Perform diffusion displacement\n",
    "        \n",
    "        random_walk =  self.disp_generator\n",
    "        d = self.disp_generator.delta # grid spacing\n",
    "        t = self.disp_generator.tau   # imaginary time step\n",
    "        \n",
    "        # the code snippet below will work correctly only if all weights are \\pm 1\n",
    "        n_walkers = self.get_number_of_walkers()\n",
    "        assert n_walkers == self.walker_weights.shape[0]\n",
    "        walker_size = self.n_alpha + self.n_beta\n",
    "        disp = random_walk(n_walkers*walker_size)\n",
    "        disp = disp.reshape(n_walkers, walker_size)\n",
    "        \n",
    "        # Some displacements give rise to the particles with\n",
    "        # identical spin labels sitting on top of each other;\n",
    "        # Those violate the Pauli principle and should be removed\n",
    "        \n",
    "        new_walker_coords = self.walker_coords + disp\n",
    "        new_coords_a, new_coords_b = np.hsplit(new_walker_coords, [self.n_alpha])\n",
    "        all_distinct = lambda ar : len(ar) == len(np.unique(ar))\n",
    "        a = np.apply_along_axis(all_distinct, 1, new_coords_a)\n",
    "        b = np.apply_along_axis(all_distinct, 1, new_coords_b)\n",
    "        #print(a)\n",
    "        #print(b)\n",
    "        alive = np.logical_and(a,b) # positions of alive walkers\n",
    "        \n",
    "        # Apply constraints if any\n",
    "        \n",
    "        if self.apply_constraint:\n",
    "            constr_a = np.apply_along_axis(self.constraint, 1, new_coords_a)\n",
    "            constr_b = np.apply_along_axis(self.constraint, 1, new_coords_b)\n",
    "            c = np.logical_and(constr_a, constr_b)\n",
    "            alive = np.logical_and(c, alive)\n",
    "        \n",
    "       \n",
    "        self.walker_coords = np.array(self.walker_coords[alive,:], dtype=int)\n",
    "        self.walker_weights = np.array(self.walker_weights[alive], dtype=int)\n",
    "        new_walker_coords = np.array(new_walker_coords[alive,:], dtype=int)\n",
    "        \n",
    "        # After this step one may still have all-electron walkers sitting atop of each other\n",
    "        # but the algorithm is blind to that\n",
    "        # Perform branching\n",
    "        # NOTE: only alive walkers are participating in branching\n",
    "               \n",
    "        m = np.exp(-0.5 * t * ( np.apply_along_axis(self.V, 1, d*self.walker_coords) + np.apply_along_axis(self.V, 1, d*new_walker_coords) - 2. * self.S))\n",
    "        r = np.random.rand(m.shape[0])\n",
    "        \n",
    "        # Assign integer weights and remove dead walkers\n",
    "        \n",
    "        new_walker_weights = np.asarray(m + r, dtype=int) * self.walker_weights\n",
    "        alive = (new_walker_weights != 0)\n",
    "        self.walker_coords = np.array(new_walker_coords[alive,:], dtype=int)\n",
    "        self.walker_weights = np.array(new_walker_weights[alive], dtype=int)\n",
    "        \n",
    "        # Some walkers can acquire abs. weights > 1\n",
    "        # Those clusters need to be separated\n",
    "        # That will be done after the annihilation step\n",
    "    \n",
    "    def _update_phases(self):\n",
    "        \n",
    "        import numpy as np\n",
    "               \n",
    "        # Update walker signs based on whether the coordinates got interchanged\n",
    "        # Split alpha and beta subsets first\n",
    "        \n",
    "        alpha_coord, beta_coord = np.hsplit(self.walker_coords, [self.n_alpha])\n",
    "        p_alpha = np.argsort(alpha_coord, axis=1)\n",
    "        p_beta = np.argsort(beta_coord, axis=1)\n",
    "        \n",
    "        parity = lambda p: self.lookup_a[tuple(p)]\n",
    "        parity_a = np.apply_along_axis(parity, 1, p_alpha)\n",
    "        parity = lambda p: self.lookup_b[tuple(p)]\n",
    "        parity_b = np.apply_along_axis(parity, 1, p_beta)\n",
    "        \n",
    "        self.walker_weights = self.walker_weights * parity_a * parity_b\n",
    "        alpha_coord = np.sort(alpha_coord, axis=1)\n",
    "        beta_coord = np.sort(beta_coord, axis=1)\n",
    "        \n",
    "        self.walker_coords = np.hstack((alpha_coord, beta_coord))\n",
    "        \n",
    " \n",
    "    def _annihilate_walkers(self, catch_collisions=False):\n",
    "                \n",
    "        # Find and collect duplicates (perform annihilation)\n",
    "        # A \"quick and dirty\" way to accomplish that is to hash\n",
    "        # the walker coordinates & perform sorting;\n",
    "                \n",
    "        converter = lambda a: hash(a.tostring())\n",
    "        walker_coords_str = np.apply_along_axis(converter, 1, self.walker_coords)\n",
    "        \n",
    "        if catch_collisions:\n",
    "            \n",
    "            # Check if we have any collisions, i.e. all the hashes are distinct for distinct walkers\n",
    "            \n",
    "            unique_hashes = set(walker_coords_str)\n",
    "            walker_index = np.arange(walker_coords_str.shape[0], dtype=int)\n",
    "        \n",
    "            if len(unique_hashes) != len(walker_index):\n",
    "            \n",
    "                # Check if there is a collision\n",
    "                for h in unique_hashes:\n",
    "                    walkers_with_hash_h = np.array(self.walker_coords[(h == walker_coords_str), :], dtype=int)\n",
    "                    if len(walkers_with_hash_h) == 1:\n",
    "                        continue\n",
    "                    else:\n",
    "                        assert np.all(walkers_with_hash_h[0][:] == walkers_with_hash_h) # broadcasted!\n",
    "        \n",
    "            \n",
    "        unique_walkers_str_, idx = np.unique(walker_coords_str, return_index = True)\n",
    "        combined_weights = []\n",
    "        \n",
    "        for i in idx:\n",
    "            idx_duplicates = (walker_coords_str[i] == walker_coords_str)\n",
    "            combined_weights.append(np.sum(self.walker_weights[idx_duplicates]))\n",
    "            \n",
    "        self.walker_coords = np.array(self.walker_coords[idx, :], dtype=int)\n",
    "        self.walker_weights = np.array(combined_weights, dtype=int)\n",
    "        \n",
    "        self._expand()\n",
    "    \n",
    "    def _expand(self):\n",
    "        \n",
    "        # Simplified version of expand function\n",
    "        # implemented without using fancy NumPy features\n",
    "                \n",
    "        n = self.walker_weights.shape[0]\n",
    "        split_weights = []\n",
    "        split_coords = []\n",
    "        \n",
    "        for i in range(n):\n",
    "            num_of_instances = abs(self.walker_weights[i])\n",
    "            sign = 1 if self.walker_weights[i] > 0 else -1\n",
    "            for j in range(num_of_instances):\n",
    "                split_weights.append(sign)\n",
    "                split_coords.append(self.walker_coords[i, :])\n",
    "        \n",
    "        self.walker_coords = np.vstack(tuple(split_coords))\n",
    "        self.walker_weights = np.array(split_weights, dtype=int)\n",
    "        \n",
    "    def get_number_of_walkers(self):\n",
    "        import numpy as np\n",
    "        return np.sum(np.abs(self.walker_weights))\n",
    "    \n",
    "    def get_energy(self):\n",
    "        return self.S\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applications ##"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fermions in harmonic trap ###\n",
    "\n",
    "Non-interacting fermions in harmonic trap provide a convenient test case for the algorithm of the previous section for the following reasons:\n",
    "\n",
    "1. Numerically exact energies for finite (non-zero) grid spacings can be readily calculated as sums of one-particle eigenvalues.\n",
    "2. The fermions are confined by soft non-divergent potential so that the branching factors in QMC are always finite.\n",
    "\n",
    "In the following we will assume that the potential energy oprator is $V = \\frac{1}{2} \\sum_{i=1}^{N_e} x_i^{2}$ and consider the cases where $N_e = 1, 2, 3$ and $4$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "def V_harm(x):\n",
    "    return 0.5 * np.sum(x**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Case $N_e = 1$ ###\n",
    "\n",
    "As a first step, we determine the reference energies of one-particle states via direct diagonalization of the Hamiltonian obtained by approximating the kinetic energy operator using a three point finite difference  formula on a **uniform** grid with sapcing $\\delta = 0.1$ a.u. spanning the interval $[-L, L]$:\n",
    "\n",
    "$$ \\hat{T}\\psi(x_i) = \\frac{2\\psi(x_i) - \\psi(x_{i + 1}) - \\psi(x_{i - 1})}{2\\delta^{2}}, $$\n",
    "\n",
    "where $x_{i\\pm1} = x_i \\pm \\delta$. Scince our goal is to approximate the solutions of the stationary Shrodinger equation for $L \\to +\\infty$ we perform a series of diagonalizations for increasing $L$ requiring that the energies of a selected subset of low-lying states have **at least 6 digits of accuracy**. The calculations are performed assuming the boundary condition $\\psi(L) = \\psi(-L) = 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energies of harmonic oscillator states on a uniform grid with spacing 0.1\n",
      "[ 0.4996873   1.49843574  2.49593063  3.49216962  4.48715031]\n"
     ]
    }
   ],
   "source": [
    "from numpy import linalg as LA\n",
    "\n",
    "def construct_H(L, delta):\n",
    "    # Input:\n",
    "    # L     - grid span\n",
    "    # delta - grid spacing \n",
    "    # Output:\n",
    "    # $H_{\\delta}$ - matrix of the discretized Hamiltonian\n",
    "    x = np.linspace(-L, L, int(2. * L / delta) + 1) # grid points\n",
    "    V = np.diag(0.5 * x**2)                        # potential energy\n",
    "    T = -1./ (2. * delta**2) * ( np.eye(len(x), k = 1) - 2. * np.eye(len(x)) + np.eye(len(x), k =-1) ) # kinetic energy\n",
    "    return T + V\n",
    "\n",
    "delta = 0.1\n",
    "L = 4. # Initial grid span\n",
    "eps = 1e-6\n",
    "converged = False\n",
    "e5 = np.zeros(5)\n",
    "\n",
    "while not converged:\n",
    "    H = construct_H(L, delta)\n",
    "    e = LA.eigvalsh(H)\n",
    "    e5_new = np.array(e[0:5])\n",
    "    if max(abs(e5_new - e5)) < eps:\n",
    "        converged = True\n",
    "    else:\n",
    "        L *= 2.\n",
    "    \n",
    "    e5 = np.copy(e5_new)\n",
    "\n",
    "print(\"Energies of harmonic oscillator states on a uniform grid with spacing {0}\".format(delta))\n",
    "print(e5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, the ground state of the 1-D harmonic oscillator can be sampled using quantum Monte-Carlo scheme:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of walkers at the end of equilibration run is 2004\n",
      "Starting production run...\n",
      "Average energy at the end of the production run:      0.500943\n"
     ]
    }
   ],
   "source": [
    "rw0101 = Random_walk(delta=0.1, tau=0.1) # define a new random walk on a grid with spacing delta = 0.1 and tau = 0.1\n",
    "D1 = Walker_ensemble(random_walk=rw0101, potential=V_harm, n_walkers=2000, n_alpha=1, n_beta=0, S=0.5)\n",
    "\n",
    "N_equil = 500\n",
    "N_prod = 2048\n",
    "E = np.zeros(N_prod)\n",
    "N = np.zeros(N_prod)\n",
    "\n",
    "#Perform N_equil equilibration steps\n",
    "for i in range(N_equil):\n",
    "    D1.one_time_step()\n",
    "\n",
    "print(\"Number of walkers at the end of equilibration run is {0}\".format(D1.get_number_of_walkers()))\n",
    "print(\"Starting production run...\")\n",
    "\n",
    "for i in range(N_prod):\n",
    "    D1.one_time_step()\n",
    "    E[i] = D1.get_energy()\n",
    "    N[i] = D1.get_number_of_walkers()\n",
    "    \n",
    "print(\"Average energy at the end of the production run: {0:13.6f}\".format(np.mean(E)))\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to calculate statistical errors of the energy estimators we used blocking algorithm as described in [H.Flyvbjerg et al. *J. Chem. Phys.* 91, 461 (1989)](https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=2ahUKEwiBjaaf68ndAhXH6YMKHSuyCdgQFjAAegQIBBAB&url=https%3A%2F%2Faip.scitation.org%2Fdoi%2F10.1063%2F1.457480&usg=AOvVaw0DEeYW79PHKKZfPb61xsKZ). The notebook relies on the python implementation of the bloking analysis due to J.Spencer which is freely available from the Github page [pyblock](https://github.com/jsspencer/pyblock)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "pyblock_is_installed = True\n",
    "try:\n",
    "    import pyblock as pb\n",
    "except ImportError as e:\n",
    "    print(\"PyBlock was not found! The error estimates will not be available.\")\n",
    "    pyblock_is_installed = False\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal block: [5]\n",
      "BlockTuple(block=5, ndata=64, mean=array(0.5009431527081165), cov=array(0.00014861497030582347), std_err=array(0.0015238467478813254), std_err_err=array(0.00013575505810765735))\n",
      "< E > =   0.500943 +-   0.001524\n"
     ]
    }
   ],
   "source": [
    "# If PyBlock is found we can proceed to calculate statistical errors\n",
    "if pyblock_is_installed:\n",
    "    reblock_data = pb.blocking.reblock(E)\n",
    "    opt = pb.blocking.find_optimal_block(len(E), reblock_data)\n",
    "    print(\"Optimal block: {}\".format(opt))\n",
    "    print(reblock_data[opt[0]])\n",
    "    print(\"< E > = {0:10.6f} +- {1:10.6f}\".format(float(reblock_data[opt[0]][2]), float(reblock_data[opt[0]][4])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Error analysis demonstrates that the energy estimate agrees with the exact diagonlaization result within the statistical error margins."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Case $N_e = 2$ ###\n",
    "\n",
    "Admissible values of the total spin for 2 electrons in a harmonic trap are 0 (singlet) and 1 (triplet). Singlet state corresponds to the doubly occupied ground state of one particle problem and is nodeless. Triplet state exhibits a node at $x_1 = x_2$. The reference energies of the states are calculated below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(singlet) =      0.999375 a.u.\n",
      "E(triplet) =      1.998123 a.u.\n"
     ]
    }
   ],
   "source": [
    "# Singlet state occupation numbers\n",
    "singlet_occ = np.array([2, 0, 0, 0, 0])\n",
    "# Triplet state occupation numbers\n",
    "triplet_occ = np.array([1, 1, 0, 0, 0])\n",
    "\n",
    "print(\"E(singlet) = {0:13.6f} a.u.\".format(np.dot(singlet_occ, e5)))\n",
    "print(\"E(triplet) = {0:13.6f} a.u.\".format(np.dot(triplet_occ, e5)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of walkers at the end of equilibration run is 2002 (singlet)\n",
      "Number of walkers at the end of equilibration run is 2004 (triplet)\n",
      "Starting production runs...\n",
      "Average energy at the end of the production run (S = 0):      1.000598\n",
      "Average energy at the end of the production run (S = 1):      1.998883\n"
     ]
    }
   ],
   "source": [
    "S2 = Walker_ensemble(random_walk=rw0101, potential=V_harm, n_walkers=2000, n_alpha=1, n_beta=1, S = 1.0)\n",
    "T2 = Walker_ensemble(random_walk=rw0101, potential=V_harm, n_walkers=2000, n_alpha=2, n_beta=0, S = 2.0)\n",
    "\n",
    "E_S = np.zeros(N_prod)\n",
    "N_S = np.zeros(N_prod)\n",
    "\n",
    "E_T = np.zeros(N_prod)\n",
    "N_T = np.zeros(N_prod)\n",
    "\n",
    "for i in range(N_equil):\n",
    "    S2.one_time_step()\n",
    "    T2.one_time_step()\n",
    "\n",
    "print(\"Number of walkers at the end of equilibration run is {0} (singlet)\".format(S2.get_number_of_walkers()))\n",
    "print(\"Number of walkers at the end of equilibration run is {0} (triplet)\".format(T2.get_number_of_walkers()))\n",
    "print(\"Starting production runs...\")\n",
    "\n",
    "for i in range(N_prod):\n",
    "    S2.one_time_step()\n",
    "    T2.one_time_step()\n",
    "\n",
    "    E_S[i] = S2.get_energy()\n",
    "    N_S[i] = S2.get_number_of_walkers()\n",
    "    E_T[i] = T2.get_energy()\n",
    "    N_T[i] = T2.get_number_of_walkers()\n",
    "    \n",
    "print(\"Average energy at the end of the production run (S = 0): {0:13.6f}\".format(np.mean(E_S)))\n",
    "print(\"Average energy at the end of the production run (S = 1): {0:13.6f}\".format(np.mean(E_T)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal block: [5]\n",
      "Singlet: < E > =   1.000598 +-   0.001908\n",
      "Optimal block: [5]\n",
      "Triplet: < E > =   1.998883 +-   0.002471\n"
     ]
    }
   ],
   "source": [
    "# A helper function\n",
    "def reblock(a):\n",
    "    if pyblock_is_installed:\n",
    "        reblock_data = pb.blocking.reblock(a)\n",
    "        opt = pb.blocking.find_optimal_block(len(a), reblock_data)\n",
    "        print(\"Optimal block: {}\".format(opt))\n",
    "        #print(reblock_data[opt[0]])\n",
    "        return \"< E > = {0:10.6f} +- {1:10.6f}\".format(float(reblock_data[opt[0]][2]), float(reblock_data[opt[0]][4]))\n",
    "    else:\n",
    "        print(\"PyBlock is not installed!\")\n",
    "        return None\n",
    "\n",
    "print(\"Singlet: \" + reblock(E_S))\n",
    "print(\"Triplet: \" + reblock(E_T))\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Case $N_e = 3$ ###\n",
    "\n",
    "Similar to $N_e = 2$ the system can have two values of total spin: $S = 1/2$ and $S = 3/2$. However, in contrast to the previous case none of the states is nodeless. The reference energies of the states are computed below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E(doublet) =      2.497810 a.u.\n",
      "E(quartet) =      4.494054 a.u.\n"
     ]
    }
   ],
   "source": [
    "# Doublet state occupation numbers\n",
    "doublet_occ = np.array([2, 1, 0, 0, 0])\n",
    "# Quartet state occupation numbers\n",
    "quartet_occ = np.array([1, 1, 1, 0, 0])\n",
    "\n",
    "print(\"E(doublet) = {0:13.6f} a.u.\".format(np.dot(doublet_occ, e5)))\n",
    "print(\"E(quartet) = {0:13.6f} a.u.\".format(np.dot(quartet_occ, e5)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quantum Monte-Carlo calculations performed for S=1/2 and S=3/2 produce energy estimates that are too low compared to the exact results if the number of random walkers is around 2000. This example demonstrates that **similar to FCIQMC present method guarantees convergence to the target state with correct nodal structure only provided the total number of walkers exceeds the critical value ($N_c$).** The $N_c$ values can be estimated from the calculations with the fixed energy onset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of walkers at the end of equilibration run is 2065 (doublet)\n",
      "Number of walkers at the end of equilibration run is 2222 (quartet)\n",
      "Starting production runs...\n",
      "Average energy at the end of the production run (S = 1/2):      2.255693\n",
      "Average energy at the end of the production run (S = 3/2):      3.378452\n"
     ]
    }
   ],
   "source": [
    "D3 = Walker_ensemble(random_walk=rw0101, potential=V_harm, n_walkers=2000, n_alpha=2, n_beta=1, S = 2.5)\n",
    "Q3 = Walker_ensemble(random_walk=rw0101, potential=V_harm, n_walkers=2000, n_alpha=3, n_beta=0, S = 4.5)\n",
    "\n",
    "E_D = np.zeros(N_prod)\n",
    "N_D = np.zeros(N_prod)\n",
    "\n",
    "E_Q = np.zeros(N_prod)\n",
    "N_Q = np.zeros(N_prod)\n",
    "\n",
    "for i in range(N_equil):\n",
    "    D3.one_time_step()\n",
    "    Q3.one_time_step()\n",
    "\n",
    "print(\"Number of walkers at the end of equilibration run is {0} (doublet)\".format(D3.get_number_of_walkers()))\n",
    "print(\"Number of walkers at the end of equilibration run is {0} (quartet)\".format(Q3.get_number_of_walkers()))\n",
    "print(\"Starting production runs...\")\n",
    "\n",
    "for i in range(N_prod):\n",
    "    D3.one_time_step()\n",
    "    Q3.one_time_step()\n",
    "\n",
    "    E_D[i] = D3.get_energy()\n",
    "    N_D[i] = D3.get_number_of_walkers()\n",
    "    E_Q[i] = Q3.get_energy()\n",
    "    N_Q[i] = Q3.get_number_of_walkers()\n",
    "    \n",
    "print(\"Average energy at the end of the production run (S = 1/2): {0:13.6f}\".format(np.mean(E_D)))\n",
    "print(\"Average energy at the end of the production run (S = 3/2): {0:13.6f}\".format(np.mean(E_Q)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting population dynamics run for the doublet state...\n",
      "Starting population dynamics run for the quartet state...\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa8AAAEYCAYAAADrpHnMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXecFdX1wL+HIiKiVFEpYsECdtdeg4pYEhKjBvNLJGqCiSZq1CRoiiZGjYmxxViwG41giUJUVMQSK0qxUFRQUZoU6ShIOb8/zh3e7O57u29337x6vp/PfGbmzp2ZM7tz3pl77rnniqriOI7jOKVEs0IL4DiO4zgNxY2X4ziOU3K48XIcx3FKDjdejuM4TsnhxstxHMcpOdx4OY7jOCVHRRovETlCRGY14fxbReT3uZQpw32aJGdSiMhkETkix9c8S0Suz+U184mIXCsiPy20HPnGdalpuC7VJmtdUtWCLsAM4CtgBTAPuBvYNOF7HgHMyrLuj4BXCvS3yVrOUl6AjYCZQNdY2ZnA+8Dy8F48CbTN8X1/DowDVgP3ZKhzCXAlcAAwGlgELAAeBraK1dsqPMNGBfw7ui7lQM5SXgqoS/cDc4FlwIfAj9PUyakuFUvL65uquimwN7Av8LsCy+PklwHA+6o6G0BEDsde8lNVtS2wC/BQAvedA/wZuKuOOscBTwHtgaFAT2Ab7Ifg7qiSqs7FfiC+lYCcDcF1qbIplC5dBfRU1c0wHfiziOxTo05OdalYjBcA4Q8+CtgVQES2FpGRIrJIRKaLyE+iuiJymYg8IiLDRWS5iEwQkT1ix1VEdojt3yMif053XxEZIiIfhetMEZHvhPJdgFuBA0VkhYgsSXctEflJkG9RkHfrGnL8VESmichiEfmniEgGOVqHay8WkSnYj0907Fci8miN+v+I3AMi8qKIXC4ir4bneFZEOsXqPiwin4vIUhH5n4j0qfG3uVlERoXnfFVEthSR64Ms74vIXrH6M0TkqLDdXEQuif39xotIdzGuE5H54Z7visiu6Z4bOBZ4Kba/L/C6qk4EUNVFqnqvqi7PcH6jUNX/qOrjwBfpjotIe2DHIMsoVX1YVZep6pfATcDBNU55ETg+lzI2Ftcl16VAvnRpsqqujnbDsn3sOXOuS0VlvESkO2adJ4aiB4FZwNbAScCVInJk7JQBWJOzA/Bv4HERadmIW38EHApsDvwRuF9EtlLVqcBPsT/4pqraLo3MfbGvjlOw5u6nwLAa1U7AXqI9Qr1jMshxKfYP3z7UGRQ7dj/QX0Tahfu2AL4H/CtW5/vA6cAWmPvgotixUUCvcGwC8ECNe5+CfaV3wtxor4d6nYBHgGszyHwBcCr2f9sMOAP4EugHHIa9sO2CrGmNBLAb8EFsfyxwjIj8UUQOFpFWGc4DIPxYLMmwvFvXufVwDDBGVdelOXYYMLlG2VTsf1xwXJdclwJ506Vw7pdYq2ku1sqKyL0uFYGPdgbmo1+Cvaw3A62B7sA6Yr5Z7MW+J2xfBrwRO9Ys/MEODfsK7BA7fg/w57B9BHX4v4G3gQFh+0fU8NPXuNadwF9jxzYF1mBN6EiOQ2LHHwKGZLjvx0D/2P7guJyY0vwkbJ8ATIkdexH4XWz/bODpDPdpF+TaPPY8t8eO/wKYGtvfDVhS4392VNj+IPpb1bhHX8z3fQDQrJ53YFr8uUPZscB/w3uxAlP45gm9g38mTZ8X9mP2wzTlu2P++kNrlB8NfOy65LpUwbrUHDgEM94tk9SlYml5fVtV26nqNqp6tqp+hX0hLtLqzdtPga6x/ZnRhqquJ/Vl2SBE5DQReTv6wsBcLZ3qOy+wdZArkmMF9lUUl/Pz2PaXmFJmutbM2P6nNY7fC/wgbP+A6l+KGe8T3BF/Ce6IZZjCQPVnnBfb/irNfiaZu2Nf29VQ1ecxd8A/gXkiMlRENstwjcVA2xrnj1LVb2ItgQHYD9+PM5yfc0SkGaZAT9co3wH74TtPVV+ucVpb7AeikLgupa7lukR+dUlV16nqK0A34GeQnC4Vi/FKxxygg4jE/xE9gNmx/e7RRvgDdQvngb1wm8TqbpnuJiKyDXA7FnnWUc2dMQmIfOn1pd2fg3U6RtdrA3SsIWe2zCX2TNjzxnkc2D34u0+gtrsiE9/HXtqjMHdOz0jcRshYk5nEfNtxVPVGVd0H6IO5PH6V4RrvhuPprrFeVccAzxP6b2oiFm69IsNS0x2RLfsCM1R1Qew+2wDPAZeras0fO7DO8Hcaeb8kcV1yXcq3LrWIPUsiulS0xktVZwKvAVeJyMYisjsW8hl/yfYRkRODz/p8zL/8Rjj2NvD98KXUHzg8w63aYEq1AEBETqf6P3Ye0E1ENspw/r+B00Vkz+BPvhIYq6ozGvbEgLlBLhaR9iLSDXM5bEBVV2E+838Db6rqZ1lety32t/kC+xG6shGyZeIO4HIR6RU6lncXkY4isq+I7B/6TVYCqzDXVTqeIvb/EZEBIjIw/B1ERPYLx99Id7Kq/lStHyXd0ifdOeE+LURkY8zV0Ty8Zy3C4eOJ+exFpCum9P9U1VszXPJw7EuyqHBdcl1KUpdEZItwj03DO3IM1nf3fKiSiC4VrfEKnIp92cwBHgMuVdXRseMjsM7LxcAPgRNVdU04dh7wTazp+X/Yl1YtVHUK8HesU3Ue5pN+NVbleawz8XMRWZjm/DHA74FHsa+97YGBDX9UwDq4PwU+AZ6ltisDzN2xW4ZjmbgvXHc2MIUML24juRb7oXgWG+NxJ9bPshn2Fb443PsL4JoM1/gvsLOkIssWAz/B/PfLsA72v6lqtl/H2fI7zI0zBHMdfUUqtDwK6434MbAdcGn8azQ6KCJbAb3J8J4VAa5LtXFdyg2KuQhnhftdA5yvqiPC8UR0SULnWMkhIpdhncg/qK9uOSEiPbBoni1VdVmh5ckVIjIY6K2q5xeBLF2w1sbWmqWCiMjfgY9U9eZEhUsA1yXXpQRlSUyXWtR10CkuQl/EBcCwclI2AFUdWmgZYmwOXJCtsgGo6oUJyuPkGNelvJGYLrnxKhFC5/U8zG3Qv8DilDWq+iEWmuyUIa5L+SNJXSpZt6HjOI5TuRR7wIbjOI7j1KIs3YadOnXSnj17FloMpwQYP378QlXtXGg5ihXXJScbCqFHZWm8evbsybhx4wothlMCiEjNzAtODNclJxsKoUfuNnQcx3FKDjdejuM4TsnhxstxHMcpOdx4OY7jOCWHGy/HcRyn5HDj5TiO45QciRovEWknIo+IyPsiMlVEDhSRDiIyWkSmhXX7UFdE5EYRmS4i74rI3rHrDAr1p4nIoMx3dBzHcSqBpFteN2DTZ+8M7AFMxaafGKOqvYAxYR9squpeYRkM3AIgIh2AS4H9gf2wNPrtE5bbKRc+/BDuuQeWL6+3quM4dTBiBLz4YqGl2EBixktsmurDsDlpUNWvVXUJNgvpvaHavcC3w/YA4D413gDahXldjgFGq+oiVV0MjMaTaTrZ8tJLcPrpsHRpoSVxnNLmkkvgmkzTiOWfJFte22Ezqt4tIhNF5I6QzbmLqs4FCOstQv2u2DTYEbNCWabyaojIYBEZJyLjFixYUPOwU6msXGnrNm0KK4fj5JPhw+HVV+uvlw2qsGQJTJ0K+++fm2vmgCSNVwtgb+AWVd0Lm756SB31JU2Z1lFevUB1qKpWqWpV586eqs4JrAgTtLrxciqFd9+FgQPhkEPgr3+Fjz+GV16BP/4x8zlr14IIXH99qmzlSrtGs2bQvr0ZsYMOSl7+LEnSeM0CZqnq2LD/CGbM5gV3YDTd8/xY/e6x87thU5ZnKnec+hk/3tYbbVRYORwnX7z9dmr7N7+BCy+EQw+Fyy6zj7hFi2qfs+++tr744lTZf/9bu/V2wAE5F7exJGa8VPVzYKaI7BSKjgSmACOBKGJwEDAibI8ETgtRhwcAS4Nb8Rmgn4i0D4Ea/UKZ42Rm0SL7knz88UJL4jj55YknrLUUEdeBL7+Es8+2VlScyODtumuqbObM6nVaty4qD0bSWeV/ATwgIhsBHwOnYwbzIRE5E/gMODnUfQo4DpgOfBnqoqqLRORy4K1Q70+qmubTwXFifPJJanv33Qsnh+Pkk3fegYcfhv32gzffTJX/5CfW+rrvPusPu+IKeOst6NULXnstVW/cOJg/H7bYAoYOtbLJk6FPH/jb3/L7LPWQqPFS1beBqjSHjkxTV4FzMlznLuCu3ErnlDVr16a2JV23qeOUIW+Fb/zLL4edd7Y+qtmz4aKLYMcdYbvt4Lnn4IMP4NRTre7mm8Mee5jhA2u59egB06fbfu/e8PXX0LJl/p+nDjzDhlMeTJsGd9yR2r/lltR2pJSOU058/TVceWX1MYxTp5p776ijzAA9/zy8/LIZLkit33svdc7SpXDttRbUAXDmmXD00bZ94YW2LjLDBSBa0/dZBlRVValPoFdhbLeduQqXLYO2bau3tvr1g2fSd5OKyHhVTecdcHBdKmqGDbPWU9eucNtt5gKsqoLtt4eJE9Ofo2qRg/Fxj61bW1Rus2a1vRRZ2odC6FFZzqTsVCBRH9dll8Gmm6bKy/DjzHEA65sCcwuecAJ06mStsO23z3yOSPVWVJcu1uKKAjxefRXGjIE//CE5uXOEGy+n9IkbqGuvLZwcjpMPbroJfvGL1P7GG8OqVbBwoe2ffXbd53/3u9ZS69YNpkwxT0XEQQfZcuKJ0KK4zYP3eTmlT9RJXZOBA/Mrh+Mkzbx51Q3X8OE2CDni7behb9+6r3HzzdZamzmzuuGK06cP7LRT+mNFghsvp3S58EJzg2RKWXPPPXkVx3ES5+qrU9tPPw2nnAJbbZUqy2ZYSLNmsPXWuZctzxR3u9Bx6qKmi7BVK1i9uvq+45QLX39t47TAUjdtsknq2Icf2vEKGhbixsspLWbNskGT775b+9icOdCxo21XeQChU0aoWl/VF1/As89WN1xgkYYVhhsvp3T4wQ/ggQfSH3vvPejQwZR83Li6I64cp5RYv97m0XriCduPxmBVON7n5ZQOmQwXVM/JVlVlY1kcp9R59llo3twGGzvVcOPlOI5TrDz9tK2vuMLW8ZydFY4bL6c0+cUvrIMaoGfPgoriOIlRc6yVv+sb8D4vp/gZORIGDEjt7767TZrXrJl1YG++eeFkc5wkmDEDtt3Wtnv2hEGDLFO8swE3Xk7xEzdcAKNGpdLZdOiQf3kcJ2ni47UOOsjSnjnVcLehU/zEXSc33VQWAywbioj8UkQmi8gkEXlQRDYWkW1FZKyITBOR4WHePESkVdifHo73jF3n4lD+gYgcU6jnceohnim+vnRPFYobL6f4Oe641PYXXxROjgIhIl2Bc4EqVd0VaA4MBK4GrlPVXsBi4MxwypnAYlXdAbgu1ENEeofz+gD9gZtFpHk+n8Wph9mz4Ze/tO3DDrMJIQ86qLAyFSluvJziZ+XK1PYOOxROjsLSAmgtIi2ATYC5QF/gkXD8XuDbYXtA2CccP1JEJJQPU9XVqvoJNmu5d6QUC2vXWrLc66+3/QsusBmQKyhrRkNw4+UUP3HjVV/S0TJEVWcD1wCfYUZrKTAeWKKq0ZTRs4CuYbsrMDOcuzbU7xgvT3POBkRksIiME5FxCxYsyP0DOen54IPU9hZbwDHu1a0LN15O8bNiBRx/vE1LvuWWhZYm74hIe6zVtC2wNdAGODZN1WhumHSf6lpHefUC1aGqWqWqVZ07d26c0E7DefZZW/fta+O5Nt64sPIUOW68nOJn5Upo166SUz4dBXyiqgtUdQ3wH+AgoF1wIwJ0A+aE7VlAd4BwfHNgUbw8zTlOIVGFu+6CPfe0ySBr5i50auHGyyl+Vq6sPjty5fEZcICIbBL6ro4EpgAvACeFOoOAEWF7ZNgnHH9eVTWUDwzRiNsCvYA38/QMTl384x8waRK0bl1oSUoGH+flFD/Ll0ObNoWWomCo6lgReQSYAKwFJgJDgSeBYSLy51B2ZzjlTuBfIjIda3ENDNeZLCIPYYZvLXCOqq7L68M4tbn7bjjvPNs+4YTCylJCuPFyipsVK+Crr6wDu4JR1UuBS2sUf0yaaEFVXQWcnOE6VwBX5FxAp3FMngxnnGHbp50Gv/lNYeUpIdxt6BQ3c+faOj5brOOUA6pw5JGp/YsusgzyTlYkarxEZIaIvCcib4vIuFDWQURGh6wAo0MkFWLcGEb/vysie8euMyjUnyYigzLdzylD3Hg55cj8+TBsGMybB6efDjNnwm67FVqqkiIfLa9vqOqeqhpNbTsEGBOyAowJ+2Chv73CMhi4BczYYe6S/TEXyaWRwXMqgM8/t3UFhsg7Zcp778E228D3v2+pzy66yAYnOw2iEG7D+Oj/mlkB7lPjDSwMeCvgGGC0qi5S1cXAaCy1jVMJeMvLKTeefRZWrYJbb7UIw969Cy1RSZJ0wIYCz4qIArep6lCgi6rOBVDVuSIS9cRnGv2fdVYArMVGjx49cv0cTj559FE46SRYutSMV8uWnj3eKU3eeMPSPfXpA1OmwIkn2gSTG28Mgwd76qcmkLTxOlhV5wQDNVpE3q+jbpOzAmDhw1RVVdU67pQQJ4WhS/36wc47Q5cuqSlQHKeU+OlP4Z13UvvDhqW23XA1iUR/EVR1TljPBx7D+qzmBXcgYT0/VM80+t+zAlQqY8dax3aFh8k7CfPaa/DrX8PVV1t/VGN56SV46qnU/uzZZriGDIHf/7563V13bfx9HCBB4yUibUSkbbQN9AMmUX30f82sAKeFqMMDgKXBvfgM0E9E2odAjX6hzCk35syxqKs4o0bBkiWFkcfJH+vWWR9QfB6rfHHaafC3v5mR2X13c+81hiOOsBycGhw/Y8fa+tvfhj/9CUaMSNW99dYmiewk2/LqArwiIu9gKWieVNWngb8AR4vINODosA/wFDbocjpwO3A2gKouAi4H3grLn0KZU068+y507Qrp+is//jj/8jj55cEH4Wc/s9ZPPlmzpvYH0623wsKF8O9/w+rV2V1HYz0VH35o6zfftP7aPfaw/fiYLvcmNJnE+rxU9WNgjzTlX2C52WqWK3BOhmvdBdyVaxmdImKPWq+KUymomvECeOwxuPzyVH/QW2+ZEfnVrxo2g/aqVbBsmfWVNm8O7cPomvXrrWzOHMvc3rw5fP21GaxTT7X0TGPGmEw33WSuv1/9qvq1ly+Hl1+2IIyFC2GffVJRsQAffWSh8A89ZIl2o+zwbdrAgAHWAvOhH03G00M5heWDD9Ircv/+FpUF1pfglCcTJtiPf8SUKfZ/b9ECevaEo4+2qNPPP08ZuPp44QU49lhrNW2yCXz5paUZ+8c/LHP7G29A9+5myCK++U3YbDP40Y/gzDNTraf777dW0sKFcOGFJsuZZ1pEbMRjj5m8Eccfb3J/8om5C+MMH27vfNu2DfkrOelQ1bJb9tlnH3VKgHXrVEF1551tHS0nnaS6cKHqz36mOnFioiIA47QI3tliXRLXpT59Uv/3YcOqvwe77qraurXqAQeotmqlumxZ/debP7/6NaJlzBjV3Xaz7R13rH5s++1T5y9frtq8efprPPZY+nJQHTIkffmaNcn97YqIQuiRxx87+WfNGvuK3m47238/jKB44AF45hl4+GHo2BFuvtncLk55smgRTJ1qrruXX66dZWLSJEtUe8011ooaOTJ1bOFCa0VpjVExV15pLsfRo61V97OfWfkll6QiCT/80Fx5S5da62jChNT5m24KN95oMtV0F37nO7WfIZL5L3+BnXaCHXdMHbvoouotMie35Nta5mPxlleRc/nl6b9SX38976LgLa/C6dLw4fZ/f+012586tfY78dBD1kLfemvVE09MnXvSSXb81VerX3OvvVSPOqp62a67pq533XX2/k2dWrdsK1aoLl6c/j0F1dNPt/WIEamyc89VnTdPdfx41Y8/Vl2/vul/oxKhEHrkLS8n/2QaS9OlS37lcArL009bIMW++9r+TjtZsMakSXDMMdYy79vXAiyOOAImTkydG/VJ9e2b2la17V12qX6fQw6x9UEHwfnnw+9+Z4Pf66JNG5u9W9Vafa+9ZueDtfhuvhkeecT6yi680II9LrvM+sf23hu23dYHISeMt2md/JMp1VPPnnkVwykwkyZBVVXKtSZihgXMsKmmDECPHuZOjqIFly618tWrLUBizBh4/XWbdbtm5Gr//hZNeOaZjZNzo43gwAPh1Vdt/q3evU2u737Xjl9zTeOu6zQJN15O/kk3QHPRIv9SrTRmzbJBwZmIvw/dullf6aefWhj6nDnWP7ViBUyfDvvtZ62zLl1gUI1Zk771LTN2m23WdJn79Gn6NZyc4MbLKTyffZYah+NUBmvWWPh7tlOBRK7F3//eWkFr1sB119l4rosusnmxhg+3vJg1gyREcmO4nKLCjZeTX+Kpd156Ce65x8bcOJXFBx+YW3D77bOrv99+cMEFcO21FpUKFv235ZbWN/Xzn5sha8hAZqek8YANJ7/E3S6HHWad305l8fnnqVmDoyCIbDj//NT200+boWrWzPqyngnpTvfeO/25TtnhLS8nf2hsTM6LLxZMDKfAPP98ajsa65cN3btbDszXX7fpcuIcdpi5n7vWmurPKVPceDn5Yfny6mmADj+8cLI4heXmm209fnzDg3R22y3VaquJu58rCncbOulZssQiu3LFxRfDtGm2/c1v5u66TmmxYIGFnF92mbv4nCbhxstJz/77p8ZdzZljX8giloG7Mfzzn6nts85qsnhOiXL77bYeMKCwcjgljxsvJz1R1oLVq+Gqq1LlEyZYSPJBB9XOK5eJ3/62+v7xx+dGRqf0ePllG0TsOSudJlKv8RKRk2MzIv9ORP4jIt7eLyUWL7asFmPGNPzcWbNsXqOIAw+EgQOt03z+/LrPXbfOjN+VVzb8vmVIxevSV1/Z/Fx1DUx2nCzJpuX1e1VdLiKHAMcA9wK3JCuWk1MmTDADVl96nC+/hHPPtf6uiJqtpjj1Tah38smpifgieveu+5zyprJ1afRo+OILOOWUQkvilAHZGK91YX08cIuqjgA2Sk4kJ+dccYWtP/207tbXsGE2Yd8f/pAqGz7c1oMHpz9n3rz05evW2SR9cVq2rB4mXXlUti599pmto2wZjtMEsjFes0XkNuAU4CkRaZXleU6x8MILqe0LLshcb6PwO/qPf9Q+ttVWtWeFBZuivSarV1vfRk1at670zPGVrUuzZtkHTOfOhZbEKQOyUZxTgGeA/qq6BOgA/KruU5yi5d13Mx9r3br6fjyUuUsXOOOM2uc0q/EKrV1rrsJvfKN23V//Ons5y5PK1qWZM20Qcc13xnEaQZ2DlEWkGfCmqu4alanqXGBu0oI5CTJmDBx5ZO3yL7+svn/eeakM3VtuWT1v3IQJZtwWLKh+zuzZ1ferquxLe/hwywJeobguATNm+LQ3Ts6o8xNIVdcD74hIjzzJ4+SatWttfcopcP31th0PyIj4y1/gtNOql7Vpk+pc79GjejaEHXYwg1Yz4nDOnOr7Y8fCU09B27YVPeWJ6xJuvJyckk16qK2AySLyJrAyKlTVbyUmlZM72rSx9Xe/a2Hu559v00b07w8jR1ofxPr1lgGjJlVV1kIbMCDlQhw40FpubdvarLE1AzZqGkZ3EcWpXF1avdo+bLbdttCSOGVCNsbrj4lL4SRHlBHjoIOqu+2eftpCl6uqaueKW7fOlpYtbf/7308de/DB1HbHjhb6HPHEE6nUTxMneq652lSuLkWpxrzl5eSIeo2Xqr4kItsAvVT1ORHZBGievGhOTuna1QxSnEyZLpo1y67FtNFGFskYTc0ez1nYvbsZN2cDFa1LM2bY2lteTo7IJsPGT4BHgNtCUVfg8WxvICLNRWSiiDwR9rcVkbEiMk1EhovIRqG8VdifHo73jF3j4lD+gYgck/3jVTiqNqvsxRdbf1PNGWbTseOO2V//lVdsHWUJj+OGqxZN1aWS5pNPbO0tLydHZNMhcQ5wMLAMQFWnAVs04B7nAVNj+1cD16lqL2AxEKV9OBNYrKo7ANeFeohIb2Ag0AfoD9wsIpXxtVqTNWvgnXfMEGWT6mnZMgvY6NQpVTZxYuZxNuPGpQxSNrz2mq1/8Qv4+OPsz6tcGq1LItJORB4RkfdFZKqIHCgiHURkdPgQHC0i7UNdEZEbwwffu/EUVCIyKNSfJiKDEnjG9MyYYW5on+nYyRHZGK/VqrohlbiItACyysgqIt2wbAJ3hH0B+mJfn2Dpcb4dtgeEfcLxI0P9AcAwVV2tqp8A04H9srl/2XHGGamEpulmIF61ylpbzz0H06enginixmrPPWtHCPboAaNG2XxbDRlAGs9RF5/OPd4v5sRptC4BNwBPq+rOwB7YB+EQYEz4EBwT9gGOBXqFZTAhBZWIdAAuBfbHdOjSyOAlzowZ9p41r8zvTif3ZGO8XhKRS4DWInI08DDw3yyvfz3wa2B92O8ILFHVEL/NLMx1QljPBAjHl4b6G8rTnLMBERksIuNEZNyCmmOPyoFFi+D++1P7q1ZZsMWjj9r+W2+lMlgcfTT06mWDQiF94EQ0NQVY2p7+/XMj5wMPWESik45G6ZKIbAYcBtwJoKpfh0HO8Q++mh+C96nxBtBORLbC8imOVtVFqroYGI15M5Jn9myf5djJKdkYryHAAuA94CzgKVWtI1urISInAPNVdXy8OE1VredYXeekClSHqmqVqlZ1Lsf0M1GfQcTs2TYT7Ukn2f5+oTEaN9xzw/jXrbaqfb0f/zg3ctW8zg475Oa65UmjdAnYLpx3d+g/vkNE2gBdwkDnaMBz5ILM9MFXuA/BBQtsaIXj5IhsjNdeqnq7qp6sqiep6u0iks1UuAcD3xKRGcAwzF14PfYVGEUOdAOiUa2zgO6wwZ2yObAoXp7mnOIj2zmuGsKyZdbKijN2bP3nTZxo60z5BEeNsnXUemsMt90GN9yQ2vdosrporC61APbGkvnuhY0RG1JH/eL7EFywwHMaOjklG+N1u4hsGAgkIqcCv6vvJFW9WFW7qWpPLODieVX9P+AFIDQXGASMCNsjwz7h+POqqqF8YIipRK7HAAAfkUlEQVRG3Bbz47+Zhdz5Z8oUCxmPjEJjeOklC8iYGftA3nzzus/JZDAffdTC2TOd378/rFwJJ57YOFnBnjeecb5Dh8Zfq/xplC5hH3CzVDX6YnkEM2bzgjuQsJ4fq5/ug68wH4Lr1pnb242Xk0OyMV4nAfeKyC4h1PdsoF8T7vkb4AIRmY71ad0Zyu8EOobyCwhflqo6GXgImAI8DZyjqutqXbUY6NPH1g2ZfHHlSut/uuEGC7Q44ggr32svW48bV73+m2ns9uLF6a/96af2g1FXWqZNNsle1kxsvLH1mz33nHfI102jdElVPwdmishOoehITB/iH3w1PwRPC1GHBwBLg1vxGaCfiLQPgRr9QlmyLFhgH1juNnRyiarWuwA7YsryDNA6m3MKueyzzz5aEExFbVm50srWrlX9059Uv/gi/TkdOlQ/L1ratLHj7dqlyu6+28oefLB63eOPT3+NaHEyAozTPL6bjdUlYE9gHPAuNjasPfbxNwaYFtYdQl0B/gl8hPWvVcWucwYWsTsdOL2+++ZEl958097DESOafi2nKMm3Hqlq5gwbIvIe1f3hHbBsAGNFBFX1ubzjLFtWff/RR+GHP4Qnn7TJHT/8EP71r9rnLVqU/no7hY/sKFfgZZfBj35k2wMHWlqm1avhwgvtHk7RkgtdUtW3gao0h2pNDxB+TM7JcJ27gDTjLBIkcoF365bX2zrlTV0pF07ImxTFykcfwW9/C3fcUf90HjUHDUfplaZPt/X996c3XpmoObVIzcwX55wDkyeb8arJj39sMgO0apX9PZ2kqGxdit5lN15ODqmrz2t5PUv5c8UVNg/VyJH1133iCVv//e+2jvqhvv46ff2ITOmY5s+Hhx9O7acbO1WzD2HBAvjrX2HoUJNnjz1SU687haSydWn+fPuYi2d6cZwmUlfLazx1h9dul4hExcLatSlXYDY5AaOMF4MHW2soMl7xYIrnnoOjjqp+3tKldk6fPjb5Y+/e1lqaODE1l1amoItOnWDIEDjhBBuI3KkT/CpMzHv88ZkT7zr5prJ1af58ezd9ehwnh2T8VVbVyh2w89prcPDBqf2vvqpdZ599rFX13nupshYtzL3Yti3ccouNr/rrX1PHR42qbrzWr7cUTuvWwbnn2gLwv//B4Yen6p2QweskAldd1fDnc/JKResS+BgvJxGy+hQKobX7ichh0ZK0YAVl6NDq+8OGmaGIktaqwoQJMGmSGZ5hw6w86l9avtyyW5x1VvXrLFtm17n2Wtt//nlb33ln9XqHHFJ9/5ZbmvY8TtFQcboE1vLyMHknx2QzJcqPgf9hob1/DOvLkhWrwDz3XPX9p5+29aGHmjsx7go87jg49VTbPuCA9Nd76CFbR0EUF15ora7IGN54Y/X6cffKTjt50EWZUJG6BN7ychIhm5bXecC+wKeq+g1gLyzPWvlSM9IvzuTJqYHEAM8+m9qOsqnHFfW22+Dkk2tfZ8YM+GOYWLdmSwvgwANt7f0E5UTl6RJ4XkMnEbL5ZVylqqvAJoxU1feBneo5p3wZNCh98MQPfpAyWrNnm1tx/frqqZPi3HZbajtdPsCRI+H002tn2HBKmcrTpTVrzFPhLS8nx2QRRscsEWmHjeofLSKLKebEuE1leT2Ry++8k748nrm9ZctUeqdMxAM52rWrfbxTp/RzdjmlTGXpEsDChbZ24+XkmHpbXqr6HVVdoqqXAb/HchB+u+6zSpj337f1OedY5vZXX4Xvfc8iC+vKEbjRRnVf98UXbd2mDXTsmCqPZiN2yp6K0yWwTDBQ/Z13nBxQb8tLRP4EvAy8pqovJS9SgYnmzTrrLNgtJAA/6CBbv/9+Km3T3XfbV2U0rqq+yRwPPzw103Hr1qnyqG/LKXsqTpcgld6sfX4mbHYqh2zchjOAU4EbRWQ5pnz/U9URdZ5VqkQDk9O58nbc0dyGu+xirkGAiy4yg1RXqyzCowYrnRlUki5Bynil0yfHaQLZuA3vUtUzgG8A9wMnh3V5EhmvzTZLf3z33VOGKyIbwxUnGth8++0NO88paSpOl8CNl5MY2bgN7wB6A/OwL8WTgAkJy1U4IuNVXyLeprDrrsnMuOwUNRWnS+DGy0mMbELlO2LTNywBFgELVXVtolIVkuXLzXD5pIpO7qksXQI3Xk5i1NvyUtXvAIjILsAxwAsi0lxVy3N+g6lTbXZjx8kxFadLYMarTZvarnbHaSLZuA1PAA4FDsNmb30ec3mUJ6NGFVoCp0ypOF0CM17e6nISIJtow2OxfGw3qGp5D6h0nGSpPF1avNiNl5MI2bgN004nXpasWWPrmtngHScHVJQuRXjLy0kIz/oaZ+lSW/fuXVg5HKdccOPlJIQbrziR8dp888LK4TjlwtKlrk9OImQ0XiIyJqyvzp84BSZKItqpU2HlcMqKitSliOXLbWZxx8kxdfV5bSUihwPfEpFhQLU0EqpafoMrF4SplTwDtpNbKk+XItx4OQlRl/H6AzAE6AZcW+OYAn3rurCIbIxFVrUK93lEVS8VkW2BYUAHLLvAD1X1axFpBdwH7AN8AXxPVWeEa10MnAmsA85V1Wca8pBZ4y0vJxmapEsly9q18NVXbrycRMhovFT1EeAREfm9ql7eiGuvBvqq6goRaQm8IiKjgAuA61R1mIjcihmlW8J6saruICIDgauB74lIb2Ag0AfYGnhORHZU1XWNkKluopaXGy8nh+RAl0qTFSts7cbLSYBsEvNeLiLfEpFrwnJCNhdWI7y9tAxL9JX5SCi/l9R8RgPCPuH4kSIioXyYqq5W1U+A6cB+2cjQYBYvtrRQrmxOAjRWl0qWaGJX1ycnAeo1XiJyFXAeMCUs54WyehGR5iLyNjAfGA18BCyJ5XObBXQN212BmQDh+FIsF9yG8jTnxO81WETGici4BVELqqEsW2bZ5BuaJd5xsqApulSSuPFyEiSbDBvHA3uq6noAEbkXmAhcXN+JwbW3Z5j6/DFgl3TVwjqdxdA6ymveaygwFKCqqqpxKduXL888FYrjNJ1G61JJErkNk5yhwalYsh3nFR9l2OBBG6q6BHgROABoJyKR0ewGRGlyZgHdAcLxzbHM2xvK05yTW5Yt869EJ2mapEslhbe8nATJxnhdBUwUkXvCl+J44Mr6ThKRzqHFhYi0Bo4CpgIvYPMYAQwCollkR4Z9wvHnVVVD+UARaRUiFXsBb2bzcA0mchs6TjI0SpdKFjdeToJkk9vwQRF5EdgXc+H9RlU/z+LaWwH3ikhzzEg+pKpPiMgUYJiI/BlzmdwZ6t8J/EtEpmMtroHh/pNF5CGsj2AtcE4ikYZgytaxYyKXdpwm6FJp4sbLSZBs+rxQ1blYCyhrVPVdYK805R+TJlpQVVdh06Knu9YVwBUNuX+jWLYMevZM/DZO5dIYXSpZ3Hg5CeK5DeN4wIbj5A43Xk6CuPGK431ejpM7VqyAZs2gdetCS+KUIXUaLxFpJiKT8iVMQVm/3pTNvxKdBKgoXYpYvtzC5H3cpJMAdRqvMB7lHRHpkSd5Ckc0JsVbXk4CNFWXwoD/iSLyRNjfVkTGisg0ERkuIhuF8lZhf3o43jN2jYtD+QcickwOHqtuPCmvkyDZBGxsBUwWkTeBlVGhqn4rMakKwbJltnZlc5KjKbp0HjbUJPq6uppizREa4cbLSZBsjNcfE5eiGIg6l73l5SRHo3RJRLph2TmuAC4IOT/7At8PVe4FLsOM14CwDZYj9KaaOUKBT8KQlP2A1xv1JNngxstJkGzGeb0kItsAvVT1ORHZBGievGh5xlteTsI0QZeuB34NRC9nR7LMESoi8Ryhb8SumTZHaE5x4+UkSDaJeX+CfcHdFoq6Ao8nKVRB+PJLW7dpU1g5nLKlMboUMs/PV9Xx8eI0VXOSIzTcs+lJrsEDoJxEySZU/hzgYGAZgKpOA7ZIUqiC4MbLSZ7G6NLB2AzMM7BJXPtiLbHEcoSq6lBVrVLVqs5NmVU8ijZ0nATIxnitVtWvo52gEI3L2l7MRMZrk00KK4dTzjRYl1T1YlXtpqo9sYCL51X1/yjmHKER7jZ0EiQb4/WSiFwCtBaRo4GHgf8mK1YBcOPlJE8udek3WPDGdKxPK54jtGMovwAYApYjFIhyhD5NkjlCI9x4OQmSTbThECz89j3gLOAp4I4khSoIbryc5GmSLqnqi9jUQsWdIxRg7VpYtcqNl5MY2UQbrg/TN4zFXBwfBDdEeeHGy0mYitEl8LyGTuLUa7xE5HjgVuAjLGJpWxE5S1VHJS1cXomMl+dhcxKiYnQJ3Hg5iZON2/DvwDdUdTqAiGwPPAmUl8KtWAEbbwwtspolxnEaQ2XoEqTSrXm0oZMQ2QRszI+ULfAxMD8heQrHihWuaE7SVIYugbe8nMTJ2MwQkRPD5mQReQqLVFKsM/itPMiWX9x4OQlRcboEbrycxKnLR/bN2PY84PCwvQBon5hEhWLlSh+g7CRFZekSwNKltt5888LK4ZQtGY2Xqp6eT0EKjre8nISoOF0CN15O4mQTbbgt8AugZ7x+2U2J4sbLSZiK0SVw4+UkTjahdY9jo/b/C6xPVpwCsmIFdOpUaCmc8qYydAlSxsv7vJyEyMZ4rVLVGxOXpNB4y8tJnsrQJYAlS2xuvOblN3uSUxxkY7xuEJFLgWeB1VGhqk5ITKpC4MbLSZ7K0CWwlpe7DJ0EycZ47Qb8EJuKIXJ1aNgvH1as8GhDJ2kqQ5fAjZeTONkYr+8A28WncsgGEekO3AdsiSnqUFW9QUQ6AMOxTusZwCmqujhMVX4DcBzwJfCj6ItURAYBvwuX/rOq3tsQWepl/XoLlfeWl5MsjdKlksSNl5Mw2WTYeAdo14hrrwUuVNVdgAOAc0SkN5ZZe4yq9gLGhH2AY7E5hnoBg4FbAIKxuxTYH8uifamI5HZszFdf2dqNl5MsjdWl0mPJEjdeTqJk0/LqArwvIm9R3U9fZ3ivqs4F5obt5SIyFZv2fABwRKh2LzbFw29C+X0hy/YbItJORLYKdUer6iIAERkN9AcezO4Rs8DzsDn5oVG6VJIsXQo771xoKZwyJhvjdWlTbyIiPYG9sKkgugTDhqrOFZFoGvSuwMzYabNCWaby3OHGy8kPTdalksHdhk7CZDOf10tNuYGIbAo8Cpyvqsusayt91XS3r6O85n0GY+5GevTo0TAhI+PlARtOgjRVl0oGVTdeTuLU2+clIstFZFlYVonIOhFZls3FRaQlZrgeUNX/hOJ5wR1IWEdZtWcB3WOndwPm1FFeDVUdqqpVqlrVuXPnbMRL4S0vJw80RZdKilWrYM0aN15OotRrvFS1rapuFpaNge8CN9V3XogevBOYqqrXxg6NBAaF7UHAiFj5aWIcACwN7sVngH4i0j4EavQLZbnDjZeTBxqrSyXHkiW2blcZsSlOYWjwzIuq+riIDKm/JgdjY1reE5G3Q9klwF+Ah0TkTOAzbFoIgKewMPnpWKj86eF+i0TkclJTR/wpCt7IGW68nALQAF0qLTyvoZMHsknMe2JstxlQRZo+p5qo6iuk768CODJNfQXOyXCtu4C76rtno1m50tZuvJwEaawulRxuvJw8kE3LKz4X0VpsYPGARKQpFN7ycvJD+esSuPFy8kI20YblPxeRRxs6eaAidAlSfV5uvJwEyWi8ROQPdZynqnp5AvIUhhUrQARaty60JE4ZUlG6BKmWlwdsOAlSV8trZZqyNsCZQEegfBQuSsrbLJtsWY7TYCpHl8Ddhk5eyGi8VPXv0baItAXOwyIAhwF/z3ReSeLToTgJUlG6BGa8mjVznXISpc4+r5AU9wLg/7A8hHur6uJ8CJZXfDoUJ2EqRpfAjNdmm5kr3nESoq4+r78BJwJDgd1UdUXepMo3Ph2KkyAVpUvgGeWdvFBXJ8+FwNbYPFpzYmltlpddSht3GzrJUjm6BNby8mANJ2Hq6vOqnOiFlSv9S9FJjIrSJfCkvE5eqCylyoS3vBwnd7jxcvKAGy+wlpcHbDhObvA+LycPuPECb3k5Ti7xPi8nD7jxAm95OU6uUIVly7zl5SSOG6916+Crr7zl5Ti5YOVK0yk3Xk7CuPH68ktbe8vLcZqOp4Zy8oQbL5/LyylyRKS7iLwgIlNFZLKInBfKO4jIaBGZFtbtQ7mIyI0iMl1E3hWRvWPXGhTqTxORQZnu2Wg8o7yTJ9x4+XQoTvGzFrhQVXcBDgDOEZHewBBgjKr2AsaEfYBjgV5hGQzcAhtSVF0K7A/sB1waGbyc4RnlnTzhxitqebnxcooUVZ2rqhPC9nJgKtAVm8jy3lDtXuDbYXsAcJ8abwDtRGQr4BhgtKouCnkVRwP9cyqsuw2dPOHGy2dRdkoIEekJ7AWMBbqo6lwwAwdsEap1BWbGTpsVyjKV5w43Xk6ecOPlLS+nRBCRTYFHgfNVta6ciOnSuWsd5TXvM1hExonIuAULFjRMSO/zcvKEGy8P2HBKABFpiRmuB1T1P6F4XnAHEtbzQ/ksoHvs9G7AnDrKq6GqQ1W1SlWrOnfu3DBBvc/LyRNuvDxgwylyRESAO4Gpqnpt7NBIIIoYHASMiJWfFqIODwCWBrfiM0A/EWkfAjX6hbLcsXQptGgBrVvn9LKOU5M6J6OsCJYF78tmmxVWDsfJzMHAD4H3ROTtUHYJ8BfgIRE5E/gMODkcewo4DpgOfInN2oyqLhKRy4G3Qr0/qeqinEoaJeX1iSidhHHj5R3MTpGjqq+Qvr8K4Mg09RU4J8O17gLuyp10NfCM8k6ecLfh0qXQqpUtjuM0Dc8o7+SJxIyXiNwlIvNFZFKsrPgyAixd6i5Dx8kV3vJy8kSSLa97qD0AsjgzAriyOU5u+OIL6Nix0FI4FUBixktV/wfU7AwuzowAbrwcJze48XLyRL77vBLLCNDogZU+95Dj5AZVN15O3iiWgI0mZQSAJgys9JaX4+SGpUttLq9OnQotiVMB5Nt4JZIRoEm48XKc3PDFF7b2lpeTB/JtvIozI4AbL8dpOm68nDyS2CBlEXkQOALoJCKzsKjB4soIsH49LF/uxstxcoEbLyePJGa8VPXUDIeKJyPA8uXWyezGy3GazsKFtvY+LycPFEvARmFYFBpxngHbcZqOt7ycPFLZxmvePFt36VJYORynHFi4EJo1849BJy+48QI3Xo6TC+bPh86dzYA5TsJU9lvmxstxcse8ea5LTt5w4wWwxRZ113Mcp37mz3ddcvKGG6927Xw6FMfJBW68nDzixsvdHI6TG+bNc+Pl5A03Xm68HKfprFxpi+uTkyfceLmyOU7TiWZy8JaXkyfceLnxcpym8/nntnbj5eSJyjVeq1ZZUl43Xo7TdGbPtnW3boWVw6kYKtd4zQ+zsbjxcpymExmvrmnninWcnFO5xityc2y5ZWHlcJxyYNYs2GgjT8rr5I3KNV6ffmprd3M4TtOZPdtaXZJu8nPHyT2Va7ymTLH1TjsVVg7HKQdmz/YPQSevVK7xmjnTXIabbFJoSRyn9Jk1y/u7nLxSucbrs8+gR49CS+E4pc/XX5sbftttCy2JU0FUrvGaMcONl+PkgmnTYO1a6NOn0JI4FURlGq+vv4ZPPoFevQotieOUPm++aeu99iqsHE5F0aLQAuSNr76C8ePNYC1YYF+Ku+5aaKkcp/R54QWbhHKXXQotiVNBVI7xmj0bDj0U7rsPWoTH3m23wsrkOKWOKrz4IhxxhIfJO3mlctyGUVThypXwyis2oNLD5B2naXzyiUXuHnFEoSVxKozKMV5t2th65UoYNQr69TMD5jhOw5g3D156ybafesrWffsWTh6nIqk84/X66/a12L9/YeVxnFLlqqvg2GNhzRp44AFzv++8c6GlciqMkjFeItJfRD4QkekiMqTBF2jRAlq1gkcfNUN26qkJSOk4FcDBB1sA1PnnwxtvwKBBhZbIqUBKwniJSHPgn8CxQG/gVBHp3eALRaG8v/41dOiQQwkdp4Lo29c+AG++2fq6zj230BI5FUipRBvuB0xX1Y8BRGQYMACY0qCr/OMfFhn1y1/mXEDHKQVEpD9wA9AcuENV/9Lgi3TsCE8+Cf/7H/z859CyZa7FdJx6KRXj1RWYGdufBewfryAig4HBAD0yZc6oqrLFcSqQmAfjaEyH3hKRkarasI9AgMMPt8VxCkRJuA2BdANItNqO6lBVrVLVqs6dO+dJLMcpKTZ4MFT1ayDyYDhOyVEqxmsW0D223w2YUyBZHKdUSefBqJUKXkQGi8g4ERm3YMGCvAnnOA2hVIzXW0AvEdlWRDYCBgIjCyyT45Qa9XowwL0YTmlQEn1eqrpWRH4OPIN1NN+lqpMLLJbjlBruwXDKhpIwXgCq+hTwVKHlcJwSZoMHA5iNeTC+X1iRHKdxlIzxchynabgHwykn3Hg5TgXhHgynXCiVgA3HcRzH2YCo1go2KnlEZAHwaYbDnYCFeRQnF7jMybGNqnpIXQbq0KVi+/8WkzzFJAvkR56861FZGq+6EJFxqlpSaTZcZqfYKLb/bzHJU0yyQPHJkyvcbeg4juOUHG68HMdxnJKjEo3X0EIL0AhcZqfYKLb/bzHJU0yyQPHJkxMqrs/LcRzHKX0qseXlOI7jlDhuvBzHcZySo2KMl4j0F5EPRGS6iAwptDxxRGSGiLwnIm+LyLhQ1kFERovItLBuH8pFRG4Mz/GuiOydJxnvEpH5IjIpVtZgGUVkUKg/TUQG5UN2J7cUQpdy9f7lSJbuIvKCiEwVkckicl6h5BGRjUXkTRF5J8jyx1C+rYiMDbIMD7NxICKtwv70cLxnrmTJO6pa9guWx+0jYDtgI+AdoHeh5YrJNwPoVKPsr8CQsD0EuDpsHweMwqa3OAAYmycZDwP2BiY1VkagA/BxWLcP2+0L/ff3pUHvQUF0KRfvXw5l2QrYO2y3BT4EehdCnnDNTcN2S2BsuMdDwMBQfivws7B9NnBr2B4IDC/0O9XYpVJaXqU4g+wA4N6wfS/w7Vj5fWq8AbQTka2SFkZV/wcsaqKMxwCjVXWRqi4GRgP9k5bdySkF0aUcvX+5kmWuqk4I28uBqdiknnmXJ1xzRdhtGRYF+gKPZJAlkvER4EgRSTfPW9FTKcYrqxlkC4gCz4rIeBEZHMq6qOpcMGUBtgjlxfQsDZWxmGR3Gkcx/Q8LriPB7bYX1uIpiDwi0lxE3gbmYx+EHwFLVHVtmvttkCUcXwp0zJUs+aRSsspnNYNsATlYVeeIyBbAaBF5v466xf4skFnGUpDdqZtS+B/mRUYR2RR4FDhfVZfV0YBJVB5VXQfsKSLtgMeAXeq4Xyn8/7KiUlpeRT2DrKrOCev52Mu3HzAvci2E9fxQvZiepaEyFpPsTuMopv9hwXRERFpihusBVf1PoeUBUNUlwItYn1c7EYkaJ/H7bZAlHN+c2u7YkqBSjNeGGWRD1M1AYGSBZQJARNqISNtoG+gHTMLki6LxBgEjwvZI4LQQwXQAsDRyVRSAhsr4DNBPRNqHSKx+ocwpHYpJlwqiI6GP6E5gqqpeW0h5RKRzaHEhIq2Bo7A+uBeAkzLIEsl4EvC8huiNkqPQESP5WrCInw8xf/BvCy1PTK7tsIitd4DJkWyYH3oMMC2sO4RyAf4ZnuM9oCpPcj4IzAXWYF9vZzZGRuAMYHpYTi/039+XRr0LedelXL1/OZLlEMzV9i7wdliOK4Q8wO7AxCDLJOAPoXw74M2gZw8DrUL5xmF/eji+XaHfp8Yunh7KcRzHKTkqxW3oOI7jlBFuvBzHcZySw42X4ziOU3K48XIcx3FKDjdejuM4TsnhxqvIEJFuIjIiZIP+WERuEpFW4dgRIqIicmas/l6h7KJY2UUi8r6ITArZpk/L8t5HiMhBuX8qx8kvrkfljxuvIiIMfvwP8Liq9gJ6Aa2xbNUR7wHfi+0PxMaIRdf4KXA0sJ+q7opl48428eYRgCudU9K4HlUGbryKi77AKlW9GzbkLPslNjp/01DnM2BjEekSlLQ/Nt1CxCXA2aq6LFxjqareSw1E5FwRmSI2v9CwkGD0p8AvxeYVOzSM3n9URN4Ky8Hh3MtE5F8i8nz4sv1JIn8Nx2kcrkcVQKUk5i0V+gDj4wVqCT9nADvEih8BTsZG1k8AVgOENFNtVfWjLO41BNhWVVeLSDtVXSIitwIrVPWacL1/A9ep6isi0gNL5RQl/dwdy6HWBpgoIk9qyNHoOAXG9agCcONVXAjpMzzXdFc8BAwHdsbS5hwUq5dtypR3gQdE5HHg8Qx1jgJ6Sypb9mZBsQFGqOpXwFci8gKWTDjTdRwnn7geVQDuNiwuJgNV8QIR2QzoAnwQlanq51iOt6OxHGpR+TJgpYhsl8W9jsfyre0DjJdUBuo4zYADVXXPsHRVm3wPaiu35xlzigXXowrAjVdxMQbYJIpqEpHmwN+Bm8LXWZw/AL8J/vw4VwH/DMqKiGwmqQkuCWXNgO6q+gLwa6AdsCmwHJvWPOJZ4Oex8/aMHRsgIhuLSEesg/qtRjyv4ySB61EF4MariFDLkvwd4CQRmQZ8AaxX1SvS1H1NVdO5F27BpkN4S0QmAS8BX9ao0xy4X0Tew/z916nNBfRf4DtRRzNwLlAVOqOnYB3REW8CTwJvAJe7n94pFlyPKgPPKl/EhLEiDwInqur4+urnCxG5jFiHtOMUM65H5YkHbBQxqvoasE2h5XCcUsb1qDzxlpfjOI5Tcnifl+M4jlNyuPFyHMdxSg43Xo7jOE7J4cbLcRzHKTnceDmO4zglx/8DO491Gcrj+EEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Population dynamics runs\n",
    "D3 = Walker_ensemble(random_walk=rw0101, potential=V_harm, n_walkers=500, n_alpha=2, n_beta=1, S = 2.51, pop_dyn=True)\n",
    "Q3 = Walker_ensemble(random_walk=rw0101, potential=V_harm, n_walkers=500, n_alpha=3, n_beta=0, S = 4.505, pop_dyn=True)\n",
    "N_max_Q = 9000\n",
    "N_max_D = 6000\n",
    "N_D = []\n",
    "N_Q = []\n",
    "\n",
    "print(\"Starting population dynamics run for the doublet state...\")\n",
    "while D3.get_number_of_walkers() < N_max_D:\n",
    "    D3.one_time_step()\n",
    "    N_D.append(D3.get_number_of_walkers())\n",
    "    \n",
    "    \n",
    "print(\"Starting population dynamics run for the quartet state...\")\n",
    "while Q3.get_number_of_walkers() < N_max_Q:\n",
    "    Q3.one_time_step()\n",
    "    N_Q.append(Q3.get_number_of_walkers())\n",
    "\n",
    "data = [('1/2', N_D), ('3/2', N_Q)]\n",
    "fig, axes = plt.subplots(nrows=1, ncols=2)\n",
    "for d, ax in zip(data, axes):\n",
    "    spin = d[0]\n",
    "    y = d[1]\n",
    "    x = np.arange(len(y))\n",
    "    ax.plot(x, y, 'r')\n",
    "    ax.set_xlabel('QMC step')\n",
    "    ax.set_ylabel('Number of walkers')\n",
    "    ax.set_title('Population dynamics (S = {0})'.format(spin))\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the results of the simulations with $S = const$ one can conclude that the critical number of walkers for the states with $S = 1/2$ and $S = 3/2$ are around 4000 and 8000, respectively. If we rerun quantum Monte-Carlo simulations with the walker ensembles of the corresponding sizes, the correct energies can be obtained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of walkers at the end of equilibration run is 4945 (doublet)\n",
      "Number of walkers at the end of equilibration run is 7917 (quartet)\n",
      "Starting production runs...\n",
      "Average energy at the end of the production run (S = 1/2):      2.499397\n",
      "Average energy at the end of the production run (S = 3/2):      4.488633\n"
     ]
    }
   ],
   "source": [
    "D3 = Walker_ensemble(random_walk=rw0101, potential=V_harm, n_walkers=5000, n_alpha=2, n_beta=1, S = 2.5)\n",
    "Q3 = Walker_ensemble(random_walk=rw0101, potential=V_harm, n_walkers=8000, n_alpha=3, n_beta=0, S = 4.5)\n",
    "\n",
    "E_D = np.zeros(N_prod)\n",
    "N_D = np.zeros(N_prod)\n",
    "\n",
    "E_Q = np.zeros(N_prod)\n",
    "N_Q = np.zeros(N_prod)\n",
    "\n",
    "for i in range(N_equil):\n",
    "    D3.one_time_step()\n",
    "    Q3.one_time_step()\n",
    "\n",
    "print(\"Number of walkers at the end of equilibration run is {0} (doublet)\".format(D3.get_number_of_walkers()))\n",
    "print(\"Number of walkers at the end of equilibration run is {0} (quartet)\".format(Q3.get_number_of_walkers()))\n",
    "print(\"Starting production runs...\")\n",
    "\n",
    "for i in range(N_prod):\n",
    "    D3.one_time_step()\n",
    "    Q3.one_time_step()\n",
    "\n",
    "    E_D[i] = D3.get_energy()\n",
    "    N_D[i] = D3.get_number_of_walkers()\n",
    "    E_Q[i] = Q3.get_energy()\n",
    "    N_Q[i] = Q3.get_number_of_walkers()\n",
    "    \n",
    "print(\"Average energy at the end of the production run (S = 1/2): {0:13.6f}\".format(np.mean(E_D)))\n",
    "print(\"Average energy at the end of the production run (S = 3/2): {0:13.6f}\".format(np.mean(E_Q)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The run took 953.185147 s for doublet\n",
      "The run took 1433.397774 s for quartet\n"
     ]
    }
   ],
   "source": [
    "# Timing information\n",
    "print(\"The run took {0} s for doublet\".format(D3.time['total']))\n",
    "print(\"The run took {0} s for quartet\".format(Q3.time['total']))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal block: [6]\n",
      "Doublet (QMC): < E > =   2.499397 +-   0.003835\n",
      "Optimal block: [6]\n",
      "Quartet (QMC): < E > =   4.488633 +-   0.004065\n",
      "Doublet (Ref): 2.49781034531\n",
      "Quartet (Ref): 4.49405367447\n"
     ]
    }
   ],
   "source": [
    "print(\"Doublet (QMC): \" + reblock(E_D))\n",
    "print(\"Quartet (QMC): \" + reblock(E_Q))\n",
    "\n",
    "print(\"Doublet (Ref): {}\".format(np.dot(doublet_occ, e5)))\n",
    "print(\"Quartet (Ref): {}\".format(np.dot(quartet_occ, e5)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The doublet energy is in slightly better agreement with the exact result compared to the quartet energy. As shown below, this can be attributed to the timestep error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of walkers at the end of equilibration run is 8064 (quartet)\n",
      "Starting production run...\n",
      "Average energy at the end of the production run (S = 3/2):      4.491615\n"
     ]
    }
   ],
   "source": [
    "rw01005 = Random_walk(delta=0.1, tau=0.05)\n",
    "Q3_1 = Walker_ensemble(random_walk=rw01005, potential=V_harm, n_walkers=8000, n_alpha=3, n_beta=0, S = 4.5)\n",
    "\n",
    "E_Q_1 = np.zeros(N_prod)\n",
    "N_Q_1 = np.zeros(N_prod)\n",
    "\n",
    "for i in range(N_equil):\n",
    "    Q3_1.one_time_step()\n",
    "\n",
    "print(\"Number of walkers at the end of equilibration run is {0} (quartet)\".format(Q3_1.get_number_of_walkers()))\n",
    "print(\"Starting production run...\")\n",
    "\n",
    "for i in range(N_prod):\n",
    "    Q3_1.one_time_step()\n",
    "    E_Q_1[i] = Q3_1.get_energy()\n",
    "    N_Q_1[i] = Q3_1.get_number_of_walkers()\n",
    "    \n",
    "print(\"Average energy at the end of the production run (S = 3/2): {0:13.6f}\".format(np.mean(E_Q_1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal block: [5]\n",
      "Quartet (QMC), tau = 0.05: < E > =   4.491615 +-   0.003630\n",
      "Quartet (Ref):  4.494054\n"
     ]
    }
   ],
   "source": [
    "print(\"Quartet (QMC), tau = 0.05: \" + reblock(E_Q_1))\n",
    "print(\"Quartet (Ref): {0:9.6f}\".format(np.dot(quartet_occ, e5)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}