HF_quantum_dots.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import lib\n",
    "import numpy as np\n",
    "from scipy import special\n",
    "import os\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hartree Fock computation of the ground state of Quantum Dots\n",
    "\n",
    "In this Jupyter Notebook we are going to study the groundstate of an N electron Quantum Dot. Our intention is to study how the chemical potential has jumps when shells are added (1s,2s2p). In both cases, the Self-Consistent Field method is used.\n",
    "\n",
    "### Basis functions\n",
    "First we write the basis we are going to use, i.e., the 3D eigenfunctions of the harmonic oscillator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "lz = 0.5\n",
    "lx = 0.75\n",
    "\n",
    "OMEGA_X = 1/lx**2\n",
    "OMEGA_Z = 1/lz**2\n",
    "\n",
    "def Hermite_pol(x, n):\n",
    "    \"\"\"\n",
    "    Evaluates the Hermite polynomial of order n at x.\n",
    "\n",
    "    Parameters\n",
    "    ==========\n",
    "    x : float or np.ndarray\n",
    "        Values in which to evaluate the Hermite polynomial\n",
    "    n : int\n",
    "        Oder of the Hermite polynomial\n",
    "\n",
    "    Returns\n",
    "    =======\n",
    "    float or np.ndarray\n",
    "    \"\"\"\n",
    "    return special.eval_hermite(n, x)\n",
    "\n",
    "def HO_wf(x, n, omega=OMEGA_X):\n",
    "    \"\"\"\n",
    "    Evaluates the n^th eigenfunction of the Harmonic Oscillator at x.\n",
    "    The units of x are sqrt(hbar/m*omega). \n",
    "\n",
    "    Parameters\n",
    "    ==========\n",
    "    x : float or np.ndarray(N)\n",
    "        Values in which to evaluate eigenfunction\n",
    "    n : int\n",
    "        Number of the eigenfunction\n",
    "    omega : float\n",
    "        Angular frequency of the Harmonic Oscillator\n",
    "\n",
    "    Returns\n",
    "    =======\n",
    "    float or np.ndarray(N)\n",
    "    \"\"\"\n",
    "    #return omega**0.25 * Hermite_pol(np.sqrt(omega)*x, n) * np.exp(-omega*x**2 / 2) / np.sqrt(2**n * np.math.factorial(n) * np.sqrt(np.pi))\n",
    "    return omega**0.25 * Hermite_pol(np.sqrt(omega)*x, n) / np.sqrt(2**n * np.math.factorial(n) * np.sqrt(np.pi))\n",
    "\n",
    "def HO_wf_3D(x, y, z, nx, ny, nz, omega_x=OMEGA_X, omega_y=OMEGA_X, omega_z=OMEGA_Z):\n",
    "    \"\"\"\n",
    "    Evaluates the n^th eigenfunction of the (anisotropic) Harmonic Oscillator at x,y,z.\n",
    "    The units of x,y,z are sqrt(hbar/m*omega). \n",
    "\n",
    "    Parameters\n",
    "    ==========\n",
    "    x, y, z : float or np.ndarray(N)\n",
    "        Position in which to evaluate the wave function\n",
    "    nx, ny, nz : int\n",
    "        Number of the eigenfunction for each cartesian coordinate\n",
    "    omega_x, omega_y, omega_z : float\n",
    "        Harmonic osciallation constant for each cartesian coordinate\n",
    "\n",
    "    Returns\n",
    "    =======\n",
    "    float or np.ndarray(N)\n",
    "    \"\"\"\n",
    "    return HO_wf(x, nx, omega=omega_x)*HO_wf(y, ny, omega=omega_y)*HO_wf(z, nz, omega=omega_z)\n",
    "\n",
    "def index_to_q_numbers(k):\n",
    "    \"\"\"\n",
    "    Returns the quantum numbers nx, ny, nz associated with the basis index k\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    k: int\n",
    "        Index of the basis from 0 to 13\n",
    "\n",
    "    Returns\n",
    "    ----------\n",
    "    nz ,ny, nz : int \n",
    "        Quantum numbers\n",
    "    \"\"\"\n",
    "\n",
    "    #q_numbers = np.array([(0,0,0),(0,0,1),(0,1,0),(1,0,0),(0,1,1),(1,0,1),(1,1,0),(1,1,1),(1,1,2),(1,2,1),(2,1,1),(1,2,2),(2,1,2),(2,2,1)])\n",
    "    q_numbers = np.array([(0,0,0),(0,1,0),(1,0,0),(1,1,0),(2,0,0),(0,2,0),(2,1,0),(1,2,0),(3,0,0),(0,3,0),(2,2,0),(2,2,0),(3,1,0),(1,3,0)])\n",
    "\n",
    "    return q_numbers[k]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now write the functions to calculate the integrands of the two-body integrals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrand_2(R, p, r, q, s):\n",
    "    r1 = R[:,0:3]\n",
    "    r2 = R[:,3:6]\n",
    "    r12 = np.sqrt(np.sum((r1 - r2)**2, axis=-1))\n",
    "    \n",
    "    x1 = R[:,0]\n",
    "    y1 = R[:,1]\n",
    "    z1 = R[:,2]\n",
    "    x2 = R[:,3]\n",
    "    y2 = R[:,4]\n",
    "    z2 = R[:,5]\n",
    "    \n",
    "    n_p = index_to_q_numbers(p-1)\n",
    "    n_r = index_to_q_numbers(r-1)\n",
    "    n_q = index_to_q_numbers(q-1)\n",
    "    n_s = index_to_q_numbers(s-1)\n",
    "    \n",
    "    result = 1/r12\n",
    "    \n",
    "    result *= HO_wf_3D(x1, y1, z1, n_p[0], n_p[1], n_p[2])\n",
    "    result *= HO_wf_3D(x2, y2, z2, n_r[0], n_r[1], n_r[2])\n",
    "    result *= HO_wf_3D(x1, y1, z1, n_q[0], n_q[1], n_q[2])\n",
    "    result *= HO_wf_3D(x2, y2, z2, n_s[0], n_s[1], n_s[2])\n",
    "  \n",
    "    \n",
    "    return result\n",
    "\n",
    "def norm_product(p, r, q, s):\n",
    "    return (np.pi/OMEGA_X)**2 * np.pi/OMEGA_Z\n",
    "\n",
    "def two_body_integrand(p, r, q, s):\n",
    "    f = lambda R: integrand_2(R, p, r, q, s)*norm_product(p, r, q, s)\n",
    "    \n",
    "    return f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### One- and two- body integral computations\n",
    "We first initialize some necessary parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_electrons = 2\n",
    "N_basis = 14\n",
    "N_points = 1000000\n",
    "integrals_file = \"integrals/integrals_QD_14nz0.npy\"\n",
    "\n",
    "normalized_wf = True\n",
    "\n",
    "max_iter_SCF = 500\n",
    "eps_SCF = 1e-4\n",
    "Delta_SCF = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The exact one-body integral function is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analytical_1(p,q):\n",
    "    n_p = index_to_q_numbers(p-1)\n",
    "    I = (p==q)*(OMEGA_X*(n_p[0] + n_p[1] + 1) + OMEGA_Z*(n_p[2] + 0.5))\n",
    "    \n",
    "    return I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The covariance function for the exponents of the gaussians is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_cov(p, r, q, s):\n",
    "    \n",
    "    cov = 0.5*np.diag([1/OMEGA_X, 1/OMEGA_X, 1/OMEGA_Z, 1/OMEGA_X, 1/OMEGA_X, 1/OMEGA_Z])\n",
    "    \n",
    "    return cov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then calculate and store the one- and two-body integrals using lib.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating two electron integrals: 1/14\n",
      "Calculating two electron integrals: 2/14\n",
      "Calculating two electron integrals: 3/14\n",
      "Calculating two electron integrals: 4/14\n",
      "Calculating two electron integrals: 5/14\n",
      "Calculating two electron integrals: 6/14\n",
      "Calculating two electron integrals: 7/14\n",
      "Calculating two electron integrals: 8/14\n",
      "Calculating two electron integrals: 9/14\n",
      "Calculating two electron integrals: 10/14\n",
      "Calculating two electron integrals: 11/14\n",
      "Calculating two electron integrals: 12/14\n",
      "Calculating two electron integrals: 13/14\n",
      "Calculating two electron integrals: 14/14\n"
     ]
    }
   ],
   "source": [
    "# One- and Two-body integrals\n",
    "MC_args = {\"f_cov\":f_cov, \"f_integrand\":two_body_integrand, \"N_points\":N_points}\n",
    "integrals = lib.integral_master(N_basis)\n",
    "integrals.calculate(integrals_file, analytical_1, MC_args = MC_args)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We calculate the overlap matrix S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not normalized_wf:\n",
    "    S = np.zeros((4,4))\n",
    "    for p in range(4):\n",
    "        for q in range(4):\n",
    "            S[p][q] = (np.pi/(ALPHA[p]+ALPHA[q]))**(1.5)\n",
    "    SVAL, SVEC = np.linalg.eigh(S) \n",
    "    SVAL_minhalf = (np.diag(SVAL**(-0.5))) \n",
    "    X = np.dot(SVEC, np.dot(SVAL_minhalf, np.transpose(SVEC)))\n",
    "else:\n",
    "    S = np.eye(N_basis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we run the Self-Consistent Field for Hartree Fock"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "E = 4.7735114 | N(SCF) = 1\n",
      "E = 7.6991238 | N(SCF) = 2\n",
      "E = 8.6360124 | N(SCF) = 3\n",
      "E = 8.7140220 | N(SCF) = 4\n",
      "E = 8.7335401 | N(SCF) = 5\n",
      "E = 8.8241742 | N(SCF) = 6\n",
      "E = 9.0630186 | N(SCF) = 7\n",
      "E = 8.7972728 | N(SCF) = 8\n",
      "E = 8.9477332 | N(SCF) = 9\n",
      "SCF CONVERGED! E = 8.9477332429\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "8.94773324287041"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lib.SCF(N_electrons, integrals, S, max_iter_SCF, eps_SCF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we compute the energies from electron pairs up to 14 electrons (limited by our basis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing energy for N_electrons =  2 \n",
      "E(2)= 8.71269978916053\n",
      "Computing energy for N_electrons =  4 \n",
      "E(4)= 16.202846539646167\n",
      "Computing energy for N_electrons =  6 \n",
      "E(6)= 22.935956737937353\n",
      "Computing energy for N_electrons =  8 \n",
      "E(8)= 32.99520206245607\n",
      "Computing energy for N_electrons =  10 \n",
      "E(10)= 41.63587893171656\n",
      "Computing energy for N_electrons =  12 \n",
      "E(12)= 52.01394107684339\n",
      "Computing energy for N_electrons =  14 \n",
      "E(14)= 64.34117491082426\n",
      "[ 0.          0.          8.71269979  0.         16.20284654  0.\n",
      " 22.93595674  0.         32.99520206  0.         41.63587893  0.\n",
      " 52.01394108  0.         64.34117491]\n"
     ]
    }
   ],
   "source": [
    "N_max = 14\n",
    "N_av = 20 #times we do the SCF to do statistics\n",
    "N_pairs = int(N_max/2)\n",
    "\n",
    "E_i = np.zeros(N_av)\n",
    "E_N = np.zeros(N_max + 1)\n",
    "E_Nvar = np.zeros(N_max + 1)\n",
    "\n",
    "for n in range(1, N_pairs + 1):\n",
    "    N_electrons = n*2\n",
    "    print ('Computing energy for N_electrons = ', N_electrons,'\\r')\n",
    "    for av in range(N_av):\n",
    "        E_i[av] = lib.SCF(N_electrons, integrals, S, max_iter_SCF, eps_SCF, print_E = False)\n",
    "        \n",
    "    E_N[N_electrons] = np.average(E_i)\n",
    "    E_Nvar[N_electrons] = np.std(E_i) / np.sqrt(N_av)\n",
    "    print('E({})='.format(N_electrons), E_N[N_electrons])\n",
    "    \n",
    "print(E_N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we proceed to compute the chemical potential. We don't have a computation for N odd, but we perform a linear interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "energies [ 0.          4.35634989  8.71269979 12.45777316 16.20284654 19.56940164\n",
      " 22.93595674 27.9655794  32.99520206 37.3155405  41.63587893 46.82491\n",
      " 52.01394108 58.17755799 64.34117491]\n",
      "chemical pot [4.35634989 4.35634989 3.74507338 3.74507338 3.3665551  3.3665551\n",
      " 5.02962266 5.02962266 4.32033843 4.32033843 5.18903107 5.18903107\n",
      " 6.16361692 6.16361692]\n"
     ]
    }
   ],
   "source": [
    "for k in range(int(Ne_max/2)+1):\n",
    "    if k != Ne_max/2:\n",
    "        E_N[2*k+1] = (E_N[2*k+2] + E_N[2*k]) / 2\n",
    "\n",
    "mu_N = np.zeros(Ne_max)\n",
    "for k in range(1, Ne_max + 1):\n",
    "    mu_N[k-1] = E_N[k] - E_N[k-1]\n",
    "    \n",
    "print('energies',E_N)\n",
    "print('chemical pot', mu_N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAF7CAYAAADc0IJwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd8VFX6x/HPQxMRARVYUDaJ6GJdFzWsZS1YwIJl13V116CyFhSsa4egohKx/VBcRYy9BPuq2BBZxYIVRNeCnSKIFAsioYbn98eZQBJSJmRmbmbm+3695jUz596ZeYYyz73nnvMcc3dEREQkezSJOgARERFJLSV/ERGRLKPkLyIikmWU/EVERLKMkr+IiEiWUfIXERHJMkr+IiIiWUbJX0REJMs0izqAZGvfvr3n5eVFHYaIiEhKTJkyZaG7d6htn4xP/nl5eUyePDnqMERERFLCzGbWtY+6/UVERLKMkr+IiEiWUfIXERHJMkr+IiIiWUbJX0REJMso+YuIiGQZJX8REZEso+QvIiKSZZT8RUREsoySv4iISJZR8hcREckySv4iIiJZRslfRGR9lZRAXh40aRLuS0qijkgkLhm/qt/6WrZsGQsWLGDZsmWsWrUq6nAyTvPmzenYsSNt2rSJOhSR9VNSAv37Q2lpeD5zZngOUFAQXVwicTB3jzqGpMrPz/f6Lum7aNEi5s2bR4cOHWjdujXNmjXDzJIUYfZxd5YuXcqcOXP4zW9+owMASU95eSHhV5WbCzNmpDoakTXMbIq759e2j7r9q7Fw4UK6dOnCJptsQvPmzZX4E8zMaNWqFVtssQXz58+POhyR9TNrVv3aRRoRJf9qrFixgg033DDqMDLehhtuyMqVK6MOQ2T95OTUr12kEVHyr4HO9pNPf8aS1q66Kgz0q6hVKygqiiYekXpQ8hcRWR+rV4fbZpuBWbjWX1yswX6SFpT8s0S/fv3o0qVLtdsmTpyImTFhwgQAevbsyV577VXtvnfeeSdmxowKA5r69euHmVV700wJyUiLFsFFF8Huu8P8+eEgYMYMJX5JG5rqJwnRoUMHxo4du057s2b6JyYZ6IorYMECeP75dbv+RdKAfpklIVq0aMHuu+8edRgiyffJJ3DzzXDqqbDrrlFHI7JedMgqIhIvdzj7bGjTRgP7JK3pzF9qVN31+tWrV8e9f5MmTWiiLlHJJI8/Di+/DKNGQfv2UUcjst70y5wqaVYDfNKkSTRv3nyd22mnnVbt/nPmzFln38suuyzFUYsk0ZIlcN550L372jK+ImlKZ/6pkIY1wP/whz9w5513rtP+9NNPM2zYsHXaO3bsyHPPPVepbfPNN09afCIpN3w4zJ4NDz8MTZtGHY1Igyj5p0Jh4drEX660NLSnKPk3a9aMsrKyareVt1ccmd+6dWvy89ctDf3BBx9U+x7Nmzevdn+RjPDVV3D99dC3L/zpT1FHI9Jg6vZPhUZQA7xjx44sXLiQFStWrLPtu+++A+A3v/lNyuIRSSv/+he0aAHXXRd1JCIJoeSfCo2gBvh+++3HqlWrqp2L/8QTT9C5c2e22WablMUjkjaefTbchg6Fzp2jjkYkIdTtnwpFRZWv+UPKa4AfeOCB9OrVi379+vHZZ5+x2267sXjxYh5++GGefvpp7rnnHo3MF6lq2TI491zYdls466yooxFJGCX/VCi/rl9YGLr6c3JC4k/hYD8zY+zYsQwbNoz777+fq666ihYtWtC9e3eeeuopjjzyyJTFIpI2RoyAr7+G8eNDt79IhjB3jzqGpMrPz/fJkyfX6zXTpk1ju+22S1JEUpH+rKXR+vbbcMZ/8MHwxBNRRyMSNzOb4u61jsBWP6+ISHUuuCAs2DNiRNSRiCSckr+ISFUvvwyPPgqDBoWlekUyjJK/iEhFK1eGwX1bbgkXXhh1NCJJoQF/IiIV3XorfPopPP00bLhh1NGIJIXO/EVEys2bB5dfHgb5HX541NGIJI2Sv4hIuUsugaVLYeRIMIs6GpGkUfIXEQF46y24996wcl+3blFHI5JUSv4iImVlYZDf5pvDkCFRRyOSdBrwl0I9e4b7iROjjEJE1nHXXTBlCjz0ELRuHXU0IkmnM38RyW4//giDB8O++8Kxx0YdjUhKKPlniXvvvRczW3Nr0aIFW221FYMHD2bZsmVr9ps4cSJmxoQJE6p9ny5dutCvX781z2fMmFHpfSvehqj7VNLBpZfCzz/DzTdrkJ9kjZR2+5tZU2Ao0BfoDMwFSoCh7r4qto8BlwP9gU2Ad4Az3P2TCu+zCXAzcESsaSxwlrv/nJpvkr4ee+wxunTpwuLFi3nyyScZPnw4ixcv5t///neD3nfQoEEcccQRldq6dOnSoPcUSboPPoDRo+GMM2CnnaKORiRlUn3N/2LgDOBE4CNgJ+A+YDlwVWyfi4DzgX7A58BlwEtmto27L47tMwbIAQ4BHLgTeADQxNw6dO/ena233hqAXr168eWXX3LXXXcxcuTIBi3p27VrV3bfffdEhSmSfO5w5pmw2WZw5ZVRRyOSUqnu9t8TeMbdn3H3Ge4+lnDWvhusOes/F7jG3Z9w948JBwobA8fF9tkOOBjo7+5vuvtbwGnAYWa2TYq/T9xKSuDtt+HVVyEvLzxvDHbZZReWLl3KwoULow5FJLVKSmDSJLjmGmjXLupoRFIq1cn/DWA/M9sWwMy2B/YHno9t3xLoBIwvf4G7LwVeIxw4AOwB/Aq8WeF9JwFLKuzTqJSUQP/+sHx5eD5zZnjeGA4AZsyYQdu2bdlss80qta9evZpVq1atc6tJdfuLNFq//BLq9v/xj1BhDItItkh18r+W0D3/qZmtBD4B7nP3UbHtnWL386q8bl6FbZ2ABe7u5Rtjj+dX2KdRKSyE0tLKbaWloT3VysrKWLVqFT/99BN33303TzzxBMOGDaNp06aV9jvooINo3rz5Orc5c+ZU+76nnXbaOvvqAEAarauuCqV8b7kFGnC5SyRdpfqa/7HACYQu/E+A7sBIM5vu7ndV2M+rvM6qtFXdXt0+ACxYsID8/Pw1z/v370///v3XL/r1NGtW/dqTadttt630fODAgZx55pnr7Hfrrbfyxz/+cZ32Pn36VPu+Q4YM4cgjj6zU1qyZykhIIzRtGtx0E5x0EvToEXU0IpFI9a/z9cAN7v5w7PlHZpYLDALuAr6PtXcCvq3wuo6s7Q34HuhoZlZ+9h8bK9CBdXsM6NChA5MnT074F6mPnJzQ1V9de6o9+eSTdOnShQULFjBixAhGjRrFbrvtxgknnFBpv27dulU6aCrXvHnzat83Nze32v1FGhV3OPvsUMhn+PCooxGJTKr7u1oBZVXayirEMZ2Q3HuVbzSzlsDerL3G/xbQmnDtv9wewEZUHgfQaBQVQatWldtatQrtqbbjjjuSn5/PIYccwrPPPku3bt248MILWbJkSeqDEUm1J5+ECRNCt3+HDlFHIxKZVCf/Z4BLzKyPmeWZ2V+A84AnYc21+5ti+xxlZjsC9xIG+I2J7TMNGAfcbma7m9kewO3As+7+eYq/T1wKCqC4GDbYIDzPzQ3PCwqijWuDDTbg+uuvZ/78+YwaNaruF4iks9JS+Ne/4Pe/h9NPjzoakUilutv/LMJ8/lGErvy5wB1AxUm21wEbAreytshP7wpz/AEKCEV+ymcFjAXWvXDdiBQUwB13hMeNqbb/EUccQY8ePbjhhhuqvfYvkjGuvTYMtHn1VdB4FMlyKT3zd/fF7n6uu+e6+4bu3tXdB7v7sgr7uLsPdffO7t7S3feNzfev+D4/untfd28Tu/VVdb/1N2zYMObPn8/o0aOjDkUkOb75JiT/f/wD9tkn6mhEIqfD3yzRr1+/SjX5K+rduzcVZk5WelzV7NmzKz3Py8urdX+RRuG888LZ/vXXRx2JSKOg5C8imW3cOHj66VDJb4stoo5GpFFQ8k+hxnStXyQrLF8epvZ16wbnnht1NCKNhpK/iGSum26CL78MZ//l021EJOVT/UREUmPOnDCf/8gj4aCDoo5GpFFR8heRzHThhbBqFdx4Y9SRiDQ6Sv410Aj25NOfsSTNq6/CQw/BxRfDlltGHY1Io6PkX40WLVqwdOnSqMPIeEuXLq1xrQCR9VJSEkpo9uwJTZuGxyKyDiX/arRv357Zs2fz448/snLlSp2hJpi7U1paypw5c+jYsWPU4UimKCmB/v3XLpdZVgZnnRXaRRqhkhJo2RLMIC8vtf9ULdMTW35+vq/Pqn7Lli1jwYIFLFu2TOvSJ0Hz5s3p2LEjbdq0iToUyRR5edUvn5mbCzNmpDoakVqVH6uWlq5ta9UqMeu+mNkUd691mVUlfxHJDGY1t69endpYROqQzGPVeJJ/jfP8zWx8TdvqMNDdv1rP14qIrJ/NNoMffli3PScn9bGI1KH86lS87YlWW5GfA4FPCMvpxuuPgPpxRSS1Zs0K/adNmlQ+y2/VCoqKootLpAY5OdWf+afqWLWuCn8nu/u78byRmTUDVjQ8JBGReli9Gv75zzC6/4YbYOTIcDCQkxMSf0MvoIokQVFR9df8U3WsWlvyLwEW1OO9ymKvqabfLbMdfvjhzJ07N+owGoXOnTvzzDPPRB2GZJNRo+Dll8NIqVNPhX/9K+qIROpUfkxaWBjNsaoG/CXmM9CgwkB/FpJSX3wB3buHef3PPVfzoD+RLBLPgL8a5/mb2UVm1jnxYYmIJMCqVXDiiWGi9J13KvGL1ENtRX6uAWaZ2Qtm9g8za5mqoERE6nT99fD223DrrbD55lFHI5JWakv+OwIjYvclwDwzu8PM9k5JZCIiNfnwQ7j8cvjb3+Dvf486GpG0U2Pyd/dP3f1iIAc4CHga+Dsw0cy+MbPLzaxriuIUEQmWL4cTToBNNw2D/dTdL1Jvddb29+Aldz8B+A3wT+Ab4FLgSzN73cxOSXKcIiLBFVfA//4XrvO3bx91NCJpqV4L+7h7qbvf7+4HEnoEbgb2BG5PRnASzJ49GzNjhuqTS7Z76y249lo46SQ47LCooxFJW/Ve1c/MWprZccBdwBmAAW8mOrBsYWa88cYbDX6f0aNH061bN1q3bs3OO+/MxIkT1/mcVq1a0bp16zW3RYsWrdk+ZMgQNttsM7p168akSZPWtK9cuZJddtmFKVOmNDhGkQZZsiSM7v/tb+HGG6OORiStxZ38zWw/M7sHmAc8CGwDDAe2dncNAozQY489xqWXXsqjjz7KokWLOO200+jTpw+zqhSJHj9+PL/++uuaW9u2bQF4//33efjhh/n666+55pprGDBgwJrXXH311fTu3Ztdd901pd9JZB2XXAJffgn33ANaDVKkQWpN/ma2nZkNN7NZwATgr8ATwH7u3tXdL3f3b1IRaDq7+eab2XLLLdl4443ZYostGDx4MAB/+MMfAOjduzetW7fmlFPC0Invv/+eI444grZt29KtWzfGjRtX6/s/9thj9O3bl+7du9O0aVNOP/10OnbsyL333htXfF999RU9evSgXbt2HHrooXz1VViX6eOPP+bxxx9n6NCh6/fFRRLlv/+FW26Bc86B/faLOhqRtFfbqn6TgZ0BB14BBgNPuPvSFMWWEb744gsuueQS3nvvPXbYYQd+/vlnPvvsMwA+/PBDzIzx48ez1157rXlNQUEBbdq0YdasWSxdupSjjz661s9YvXo1VSs1ujsffPBBpba//e1vrFy5kq222oqLL76Yo446CoAdd9yRwYMH8+OPP/LKK6+w0047UVZWximnnMJtt91Gy5Yq8SARWrQo1O7fZhsYPjzqaEQyQm1n/hsBQ4A8d+/l7g8q8ddfs2bNcHc++eQTfv31V9q1a8fuu+9e4/5z5szh5Zdf5oYbbqBt27Z06tSJyy+/vNbPOPzww3nwwQeZPHkyK1eu5JZbbmHWrFn88ssva/aZMGEC06dPZ/bs2Zx33nkUFBSs6VHYfvvtueCCCzjggAMYPXo0d999NzfeeCM9evSga9euHH300eyzzz4UaXU0icI558CcOXDffbDhhlFHI5IRapvnv527D3f32akMKNN07dqVkpIS7rjjDjbffHP22msvxo8fX+P+s2eHP+7c3Nw1bVtuuWWtn3HCCSdw4YUXUlBQQKdOnXj//fc54IADaF9hGtQBBxxAy5YtadmyJcceeyx9+/alpKRkzfbTTz+dqVOn8tJLL9GiRQvuvvtuhg8fzvnnn0+fPn145ZVXePHFF3nppZfW949CpP6efjok/UGDYLfdoo5GJGPUZ8DfDmY2xsxmm1mpme0ca7/CzA5MXojp76ijjuKll15i4cKFHHPMMRx55JGUxtZxtCoFSrbYYgsAZlZY6Hn69Om1vr+ZcfHFF/P555/zww8/MHr0aKZNm0bPnj1rfE2TJk3WuVQA4XLBKaecwk033UTr1q358MMP2XPPPWnatCl77rknU6dOjfdrizTMggVhzdPu3eGyy6KORiSjxJX8zWx34D3gj8DzwAaEKX4ALYABNbw0633++eeMGzeO0tJSmjdvTtu2bTEzmjQJf/SdOnXiyy+/XLN/ly5d6NmzJxdddBG//PIL8+bN46qrrqr1MxYtWsS0adNwdxYsWMCAAQNo06YNJ554IhAG7r377rusWLGClStX8tRTT/HAAw9wzDHHrPNeo0ePpmvXrvTu3RuArbbaihdeeIGlS5fyyiuvsPXWWyfqj0akZu5w+unw889w//3QokXUEYlkFnev8wa8Skj6TQmDBFcDu8S2/RWYGc/7RHHbddddPdlq+4z//e9/vvvuu3ubNm28TZs2vssuu/jzzz+/Zvvdd9/tXbp08Xbt2nn//v3d3f27777zPn36eJs2bfx3v/ud33HHHQ749OnTq/2MWbNm+Q477OAbbbSRb7LJJt63b1///vvv12x/+eWXffvtt/dWrVp5u3btfNddd/WHHnqo2vfp1q2b//TTT2vaPvnkE+/evbu3bdvWTzrpJC8rK1vvPwuRuD34oDu4X3NN1JGIpB1gsteRG82r6fqtysxKgb+6+wtm1hRYCeS7+/tmtg/wors3ypE4+fn5nuz15bWG/Vr6s5AGmzMHdtwRtt8eXnsNmjaNOiLJAiUlUFgIs2ZBTg4UFUFBQdRRrR8zm+Lu+bXtU+NUvyqWAzXN9+oE/FLDNhGR+LnDySfDihVhoJ8Sv6RASUkYXhIbisXMmeE5pO8BQF3iHfA3CTjTKo9OK+8y6AdMTGBMIpKtbr8dXnwRrr8eNL5EUqSwcG3iL1daGtozVbxn/pcDrwNTgEcJif8fZjYc2AvQHBwRaZivv4YLLoBevWCAxhBL6lSphF5neyaI68zf3acA+wNLgSLCSP8LgLbAAe7+SdIiFJHMV1YWFu1p1gzuuguqTIEVSaacnPq1Z4K45/m7+9vu/iegHbA1sKm77+Hu7yQtOhHJDiNGwKRJcPPNYdU+kRQqKoJWrSq3tWoV2jNVvZf0dffF7v6Nu/+cjIBEJMt8/DEMGQJ//jMcf3zU0UgWKiiA4mLIzQ2dTrm54XmmDvaD+K/5Y2ZbAUcBv2Xdkf/u7qclMjARyQIrVsAJJ0DbtmGwn7r7JSIFBZmd7KuKK/mb2bHAg4QiPz8BK6rs4oCSv4jUz7BhMHUq/Oc/0LFj1NGIZI14z/yvAsYBJ7n7giTGIyLZ4r334OqrQ1f/X/4SdTQiWSXe5N8ZGKjEX73OnTuTn19rMaWs0blz56hDkHSwdGno7u/cOQzyE5GUijf5v00Y4T8hibGkrWeeeSbqEETSy+DB8NlnMH48tGsXdTQiWSfe5H8W8LiZzQPGu/uSJMYkIpls4kS46SYYODAU9BGRlIs3+X8DvAM8DmBmZVW2u7tvkMjARCQD/fIL9OsXSvded13U0YhkrXiT/2jgeMKgv89Yd7S/iEjdzj8fvv0WXn8dNtoo6mhEsla8yf+vwCXu/n/JDEZEMthzz8Gdd8LFF8Oee0YdjUhWi7fC3zLgw2QGIiIZ7Icf4JRTYMcd4Yoroo5GJOvFm/zvA45JZiAiksHOOAMWLoT774cNNDxIJGrxdvtPAy4zs2cI1/1/qrqDu49JZGAikuZKSsKC6DNnhudHHw077xxtTCICxJ/874rd5wJ9qtnugJK/iAQlJdC/P5SWrm177rnQnk0F1EUaqXiT/3ZJjUJEMkthYeXED6GqX2Ghkr9IIxBX8nf3z5MdiIhkkFmz6tcuIikV74A/EZH4tW9ffXtOTmrjEJFq1Zj8zexTM/t9vG9kZk1ir6n1EoGZdTaz+8xsgZkti71m3wrbzcyGmtl3ZrbUzCaa2Q5V3mMTM3vAzBbFbg+YmQqEizQG06bB4sXQpMrPS6tWUFQUTUySdkpKIC8v/DPKywvPJXFqO/PfFmhZj/ey2Gs2rHGHkKAnxfbtQxhLcBYwv8JuFwHnx9p7xLa9ZGYbV9hnDLALcAhwcOzxA/WIVUSS4eef4cgjoU0buPFGyM0Fs3BfXKzr/RKX8vGiM2eCe7jv318HAIlk7l79BrPVwHLCSP54bQD0cPf3a3jPq4F93f1PNWw34DvgFncvirVtSDgAuMDdb4/1LHwK7OXuk2L77AW8DmxbdXxCfn6+T548uR5fQUTWS1kZHH44vPQSvPwy7L131BFJmsrLWztDtKLcXJgxI9XRpB8zm+Luta4zX9uAv2vX83Pn1rLtz8A4M3sE2I+Q6O8EbvVwFLIl0AkYX/4Cd19qZq8BewK3A3sAvwJvVnjfScCS2D4anCgShSFD4IUX4LbblPilQTReNPlqTP7uPigJn9cVGAjcCFwDdAf+Hdt2CyHxA8yr8rp5wBaxx52ABV6hy8Ld3czmV3i9iKTSI4/ANdfAqafCaadFHY2kuZyc6s/8NV40cVI92r8J8L67D3L3qe5+D3AzcEaV/apearAqbdVdiqi6DwALFiwgPz9/za24uLgB4YvIOj78EP75z7BYzy23hGv8Ig1QVBTGh1ak8aKJFW+Rn0SZS7heX9E04JzY4+9j952Abyvs05G1vQHfAx3NzMrP/mNjBTqwbo8BHTp0QNf8RZJk4cIwwG/TTeGJJ6BFi6gjkgxQPi60sDB09efkhMSv8aKJk+rkPwnYpkpbN6C8g2c6Ibn3At4DMLOWwN7AhbF93gJaE679l1/33wPYiMrjAEQkmVauhGOOge+/h9dfh0666iaJU1CgZJ9MqU7+NwJvmlkh8AiwM3A2MBjWXLu/CSg0s8+AL4AhhAF+Y2L7TDOzccDtZnYqobv/duBZVSIUSaELL4RXXoF774UePaKORkTqIaXJ393fM7M/A1cDlwKzYvejKux2HaFWwK3AJsA7QG93X1xhnwLCWIHyWQFjgTOTG72IrHHffTByJJxzDpx4YtTRiEg91TjPP1Nonr9Igr37LuyzD/zpT/Dii9As1R2IIlKbhs7zr/hGTQnV9nYHNiecmS8kzKl/zd2/a2CsIpIO5s6Fv/wFOncO0/uU+EXSUq3/c80sh3BN/gRgM8L19aWxWzvC1D03s0mErvtHPNO7EkSy1fLl8Ne/hhK+b75Z8+I9ItLo1bawzw2EAXe9CAP19gFau/tG7t7e3ZsRKvIdB8wAioEPzGzXpEctIqnlDmedBW+9BffcA3/4Q9QRiUgD1HbmvyPQ093frmkHd59JmKb3iJltRCjWsxswJaFRiki0Ro+GO+6AQYPC9D4RSWu1lfc9uD5v5O5LCCP1RSSTvP46nH02HHIIXHVV1NGISAKkuryviKSTb7+Fo4+Grl1hzBho2jTqiEQkAWo88zezjvV5I3ef3/BwRKTRWLoU/vzncD9xIrRrF3VEIpIgtV3z/57qF9CpjtfxXiKSTtzDCn3vvw9jx8J220UdkYgkUG0J+2zqTv6HAQclLhwRaRRuvBFKSuDKK+Hww6OORkQSrLYBf7fUtM3M9gWKgD0JI/sLEx+aiETipZdC3f6jjgrLqknWKynRCnuZpl5d9WaWT6jLfwDwGXC0u/8nGYGJSAS+/hqOPRa23z7U72+iMcHZrqQE+veH0tLwfObM8Bx0AJDO4vqfbWY7mNl/CIvsbA38E9hRiV8kg/z6axjgB/DUU9C6dbTxSKNQWLg28ZcrLVWnULqrq7xvV+AK4B/AfOAs4A53X5mC2EQkVdyhXz/49FMYNw622irqiKSRmDWrfu2SHmor7zsamAYcDAwCurr7KCV+kQxUVARPPAHXXQe9ekUdjTQiOTn1a5f0UNuZf3/CaP9vCPX9e5lZTfu6u2vUv0g6euYZuPRS6NsXzjsv6mikkSkqqnzNH6BVq9Au6au25P8ua6f6bZyCWEQk1T77LIza2mUXKC6Gmg/wJUuVD+rTaP/MYpm+Am9+fr5Pnjw56jBEGp+ff4bddoOffoLJk9WPK5IhzGyKu+fXto+q8olko7KycOr2zTfw8stK/CJZprYBf+3X5w3X93UikkKXXQbPPw833wx77x11NCKSYrXN859uZtfGpvvVysyam9lRZvYOMCBx4YlIwj32GFx9NZxyCpx+etTRiEgEauv2Pwi4HrjAzN4FXgc+BBYAy4FNgK7AH4HeQHPgRmBEMgMWkfVQXp915swwqO93v4NbbtEAP5EsVVtt/zeBP5nZn4CTgQLggvLNgAGrCLX9rwTuc/efkhuuiNRb1fqs7jB7Njz+uIZsi2Speo32N7McYHOgJfAD8LW7l9b+qmhptL9kvby8cMZfVW4uzJiR6mhEJMkSPtrf3WcBKuookk5UnzXj9OwZ7idOjDIKSWe1jfbvZmYt6noDM2trZocmNiwRSZiOHatv1/Q+kaxV22j/aUD38idm1sTM5pvZ76vsty3wTDKCE5EG+uEHWLly3YF9qs8qktVqS/5VhwEb0J4wql9EGrvylfp+/RWuvDJc4zcL98XFGuwnksVU4U8kU40cCc8+G+7PPhuGDIk6IhFpJGo78xeRdDV5Mlx0ERx5JJx1VtTRSAKVlMDbb8Orr4aJHCUlUUck6UjJXyTTLFoExx4LnTvD3XerkE8GKS/ZsHx5eD5zZniuAwCpr7q6/Xub2daxx00IxX0ONrNtK+yzVVIiE5H6cw/ZYOZMeO012HTTqCOSBCoKqnfDAAAgAElEQVQsXFurqVxpaWjXEA6pj7qS/5XVtA2rpi2z1wUWSRd33AGPPgrDh8Oee0YdjSSYSjZIotSW/LdLWRQi0nAffQTnnAO9e4fr/ZJxcnKqL9aokg1SX7XV9v88lYGISAMsWQLHHAPt2sEDD0ATDefJREVFlZdpAJVskPWjqX4imeDMM+Hzz2HChJor+knaK7+uf/LJYdBfbm5I/LreL/Wl5C+S7h54AO69Fy69FPbfP+poJMkKCsLQDlBtf1l/6hsUSWeffw4DBsA++8Bll0UdjYikCSV/kXS1bFmYz9+yJYwZA83UkSci8dGvhUi6Ov98+PDDUMJ3iy2ijkZE0khcZ/5mdryZbZDsYEQkTk88AaNGhQOAPn2ijkZE0oy5112fx8xWAz8C9wHF6TQNMD8/3ydPnhx1GCKJM3067LwzbLMNvP46tGgRdUQi0oiY2RR3z69tn3iv+W8PPACcCHxqZq+Y2bFmpuV9RVJpxQr4+9/D44cfVuIXkfUSV/J398/c/V/AFkA/oDnwEDDbzIabWdfkhSgiaxQWwrvvwp13wpZbRh2NiKSpeo32d/fl7v6Au+8F/B6YBlwEfGFm48zsgGQEKSLA88/DDTeEqX1HHx11NCKSxuo91c/MNjSzk4C7gX0IBwDDgM2A8WY2JLEhighz5sAJJ8BOO8GIEVFHIyJpLu7kb2Y7mtm/ge+A24AZwH7uvqO7D3X3HsDVwDlJiVQkW61aBccdF+b1P/pomNcvjUJJCeTlhaUU8vLCc5F0ENc8fzN7E9gNmAOMIIz4n1fNrs8BhYkLT0S46ip47TW4//4wwl8ahZKSyovszJwZnoNq7UvjF++Z/xLgaCDP3a+qIfEDTEVLAYskziuvhOR/4olw/PFRRyMVFBZWXl0PwvNCnf5IGohrnn860zx/SVvz50P37tC2Lbz3HrRuHXVEUkGTJlDdz6cZrF6d+nhEysUzzz/ebv/a1ghdDfzi7ivqE5yI1GL16jDA78cfYdw4Jf5GKCcndPVX1y7S2MXb7f89MLeG2zyg1Mw+MjP1S4okwg03wIsvwsiRYYS/NDpFRdCqVeW2Vq1Cu0hjF+/CPmcDFwBLgf8QEn4n4CigJXAn0BO418xWu7vGvIqsr7fegsGD4W9/WzuCTBqd8kF9hYUwa1Y44y8q0mA/SQ/x1vYfTijx+2ev8AIzM+BpYJq7X2xmDwO/c/ddkxVwfemav6SVn34K1/mbNoWpU8P1fhGRekhkbf8Tgdu9ypFC7Pno2HYIJX812l9kfbjDSSfB3LnwyCNK/CKSNPEm/3axW3U2ATaOPV5EGAAYFzMbbGZuZrdUaDMzG2pm35nZUjObaGY7VHndJmb2gJktit0eMLOa4hNJD7feCk89BddeCz16RB2NiGSweJP/60CRme1YsdHMfg8UxbYDbAXMiucNzWx34FTgf1U2XQScD5wF9ADmAy+Z2cYV9hkD7AIcAhwce/xAnN9FpPF5/304/3w47DA499yooxGRDBdv8j8TaAp8aGbTzOxVM5sGfABYbDvApkBxXW9mZm2BEuBk4KcK7QacC1zj7k+4+8eESwobA8fF9tmOkPD7u/ub7v4WcBpwmJmp/Jmkn8WL4dhjoWNHuPfeMFFcRCSJ4l3S90ugG/Av4H2gLHZ/LrCNu38V2+96d78pjrcsBh5395ertG9JmEUwvsJnLwVeA/aMNe0B/Aq8WeF1kwhVCPdEJJ24w+mnwzffwJgxsNlmUUeU9nr2DDcRqVmdU/3MrDmwP2FE/80N/UAzOxXYGqiuJkCn2H3V8sHzgC0q7LOg4uBDd3czm1/h9WssWLCA/Py1gx779+9Pf02fksbinntC0h82DPbeO+poRCRL1Jn83X2lmT1D6GqP63p+TWLd8lcDe9dREbDq/EOr0lbd/MSq+wDQoUMHNNVPGpWSkjA5fObM0MW/ww5wySVRRyUiWSTea/7TgUT0R+4BtAc+NrNVZrYK2BcYGHv8Q2y/qmfwHVnbG/A90DE2PgBYM1agA+v2GIg0LuVLwZXXhXUPXf4PPxxtXCKSVeJN/iOAwWa2SQM/7yng90D3CrfJwMOxx18Qknuv8heYWUtgb9Ze438LaE04kCi3B7ARlccBiDQ+1S0Ft3SploJLkJISePttePVVyMsLz0VkXfGW980nnFnPMLPXCTX9K3XDu/tpdb2Ju/8M/FyxzcyWAD/GRvZjZjcBhWb2GeFgYAhhgN+Y2HtMM7NxwO2x8QMG3A486+6fx/l9RKIxq4YrZzW1S9zKO1WWLw/PZ85cWx1ZJXdFKos3+R8Wuy8Fqivd64TpdolwHbAhcCuhgNA7QG93X1xhnwLgZtbOChjL2umGIo1Xu3ahhG9VWgquwarrVCktDe1K/iKVxZX83b1zsgJw955VnjswNHar6TU/An2TFZNIUjzxBPz8c6jbX1a2tl1LwSWEOlVE4hfvNX8RaYjXXw+nn3vsAXfcAbm5YaR/bi4UF+vUNAFq6jxRp4rIuuLt9i8feHcCsA9h5P9Z7v6VmR0FfBQrBCQiVX36KRxxRBiBNnZsKOTzz39GHVXGKSoK1/grdv2rU0WkenElfzPbHHiZULv/G0KRnjaxzYcSK7ebjABF0tqcOXDwwbDBBvDCC6rgl0TlnScnnxwG/eXmhsSvThWRdcV75v9/hEsE2wEzgIoFel4BLktsWCIZYNEiOPTQMMDv1Vdhyy2jjijjFRSEqyoAEydGGopIoxZv8j8IGBDr5m9aZdsc1pbeFRGAFSvgr38NXf7PPQe77BJ1RCIia8Sb/Degyvz8CjYmLPQjIgCrV8NJJ8F//wv33Qe9e0cdkYhIJfEm/4+BI4EXq9l2EGGFPxEBGDQoVJwpKoITTog6mqyj7n6RusWb/EcAY8ysjFilPWBrMzsIOBU4OhnBiaSdW26B666DAQPCQYCISCMUb5GfR8ysMzAMGBhrfhhYClzg7s8kKT6R9PGf/8DZZ8ORR8K//x3m8YuINEJxz/N395vM7B7CIjsdCSvwvebu1dQqFckyb7wBxx0Hu+0GY8aEKn4iIo1U3MkfwN0XAc8mKRaR9DRtWijik5sLzzwTKsuIiDRi9Ur+ZtYdyAFaVt3m7o8mKiiRtPHdd6GIT4sWMG4ctG8fdUQiInWKt8JfN+A/hCI/1V3IdEDJX7LLL7+EIj4//qgiPiKSVuI98x9FKOd7AvARsDxpEYmkg/IiPp98As8+qyI+IpJW4k3+fwROdvfHkhmMSFpwDwXkJ0yAe+6Bgw6KOiIRkXqJd0nfH4HSOvcSyQaDB8ODD8KwYdCvX9TRNFolJWEhwyZNwn1JSdQRiUi5eJP/zcDpZpq4LFnu1lvhmmvgtNPCQYBUq6QkLK87c2boKJk5MzzXAYBI42DuXvdOZoXAycASQonfH6vs4u4+PPHhNVx+fr5Pnjw56jAkEzz5ZLjOf/jh8MQT0Kxek2WySl5eSPhV5ebCjBmpjkYku5jZFHfPr3WfOJP/6jp2cXdvlFVNlPwlISZNggMPhO7dw4I9mstfqyZNwhl/VWZh3SMRSZ54kn+8py4bJiAekfT02WfhbP+3v1URnzjl5FR/5p+Tk/pYRGRdcV3zd/fldd2SHahIJObODUV8mjdXEZ96KCpa9xipVavQLiLRqzH5m1lvM9u4Sts6pzxm1tHMBlZtF0l75UV8Fi6E55+Hrl2jjihtFBRAcXG4xm8W7ouLQ7uIRK+2M/8XCBX9ADCzpsBiM6tazWRL4N9JiE0kOitWwNFHw0cfweOPw667Rh1R2ikoCIP7Vq8O90r8Io1Hbdf8q5vWp6l+kvnc4ZRT4KWXQhGfgw+OOiIRkYSKd56/SPYoLIQHHoCrrlIRHxHJSEr+IhWNGgXDh4eKNIWFUUcjIpIUSv6S3SrWoO3YEc48M0zru/XWMFJNRCQD1TXP/0QzOzD2uAlh6d6TzKziRdAuSYlMJNnKa9CWxpatWLAgHAQceaSq94lIRquxwl8cVf0qUoU/ST+qQSsiGaihFf5U1U8y26xZ9WsXEckQNSZ/Ve2TjNelC3z77brtqkErIhlOA/4kO61cCe3arduuGrQikgWU/CX7lBfx+egjOPVU1aAVkayjIc2SfS69FO6/H664Ai67LOpoUqZnz3A/cWKUUYhIY6Azf8kuo0eHbv1TTw0HAVmipATefhtefTVMcigpiToiEYmSkr9kj6efhjPOgMMOC5X8sqSIT3k5g+WxIbwzZ4bnOgAQyV5K/pId3noL/v53yM+Hhx/OqiI+hYVr6xiVKy1V9WKRbFbjL6CZPV+P93F375OAeEQS74svQsneLl3gmWdgo42ijiilVM5ARKqq7fRnU0I5X5H09f33YUneJk1g3LhQvz/L5ORUX8hQ5QxEsldtRX52T2UgIgm3eDH06QPz5oUh7lttFXVEkSgqqryEAaicgUi20zV/yUwrV8Lf/gYffgiPPgo9ekQdUWQKCkL5ApUzEJFy9Rr1ZGYbAVsBLatuc/d3ExWUSIO4h1PdF1+EO+8MZ/9ZrqBAyV5E1oor+ZtZC2A00BeoafW+Rrmqn2Shyy6De++FoUPh5JOjjkZEpNGJt9t/MNAHGAAYcD5wJvAe8DVwVFKiE6mv22+HYcNC+d4sqt4nIlIf8Sb/Y4ErgXtjz19z99tigwI/BfZJQmwi9TN2LAwcCIceCrfdljVFfERE6ive5J8LfOTuZcBKoFWFbcXAcYkOTKRe3n47FPHZddcwwC+LiviIiNRXvMn/B6B17PFsYKcK29oB2VU1RRqX8iI+m28Ozz7b6Iv4lJSE+vpNmqjOvohEI97To/cICf954CngSjPbAFgFXAK8mZzwROowb14o4gNpUcSnvM5++Zz78jr7oNH4IpI68Z75XwfMiD2+CngHuAEYCcwDzkh4ZCJ1+fXXtUV8nnsOtt466ojqpDr7ItIYxHXm7+5vA2/HHv8M9DGz1kArd5+fxPhEqldexGfq1LBa3x//GHVEcVGdfRFpDOI68zezi8ysU8U2d//V3eeb2W/M7KLkhCdSDXc47bTQzT96dFiiN03UVE9fdfZFJJXi7fYfDtT089Qltl0kNYYOhXvuCfP4Tz016mjqpago1NWvSHX2RSTV4k3+tU2YbgusSEAsInUrLoYrr4STTgoHAWlGdfZFpDGo8Zq/me1F5eI9/czswCq7bQgcCUxLQmwilT37LAwYAIccErr707SIj+rsi0jUahvwdwBweeyxA6dXs48DnxNK/dbJzAYRSgFvAywnDCIc5O4fV9jHYp/bH9iEMLPgDHf/pMI+mwA3A0fEmsYCZ8UGI0omeucdOOYY2HnnUMSnefOoIxIRSVu1dfsPI5zZtyJ0++8Te17x1szdt3f31+L8vJ7AKGBPYH9CnYAJZrZphX0uIqwdcBbQA5gPvGRmG1fYZwywC3AIcHDs8QNxxiDp5ssvw6C+zp3DlL7Wret+jYiI1KjGM/9YKd8yADPb0N2XN/TD3P2gis/N7HhgEfAn4JnYWf+5wDXu/kRsnxMJBwDHAbeb2XaEhL+Xu78Z2+c04HUz28bdP29onNKIVC3i85vfRBuPiEgGiHee//JYRb/jgX2BTQklfycCJQ04MNiY0PvwU+z5lkAnYHyFz15qZq8RegtuB/YAfqVyVcFJwJLYPkr+meLXX8MZ/9y58Mor8LvfRR2RiEhGiHeefwdgMmERnwOBzYFewJ3Ae2bWfj0/fyTwAfBW7Hl5LYF5VfabV2FbJ2CBu3v5xtjj+RX2WWPBggXk5+evuRUXF69nqJIS5YXvzaBDB5gyBR55BHbbLerIREQyRry1/a8FOgO93P2/5Y1mdgDwUGz7yfX5YDMbAexF6L4vq7LZq+5epa3q9ur2AaBDhw5Mnjy5PqFJVKoWvl+2LAzs++WXpH5kYWGosJeTE+bbayS+iGS6eOf5H0YYlf/fio2x50Ni2+NmZjcC/wD2d/dvKmz6PnZf9Qy+I2t7A74HOsbGB5S/nwEdWLfHQNJJdYXvV65MWuH78mONmTND0cDyRXa0yp6IZLp4k38boKbq4zNj2+NiZiMJg/f2d/fPqmyeTkjuvSrs3xLYm7XX+N8iLC+8R4XX7UFYVlirC6azFBe+1yI7IpKt4k3+XxDO1KtzbGx7nczsVuCfsff6ycw6xW6tYc21+5uAS8zsKDPbEbiXMMBvTGyfacA4wsj/3c1sD8JAwGc10j+NucPGG1e/LUmF77XIjohkq3iv+d8I3BUb+FcCzCV0zf+d0OV/UpzvMzB2/98q7VcAQ2OPryPUELiVtUV+erv74gr7FxCK/JTPChhLnIWGpBFyh3PPDdf2mzWDVavWbkti4fucnNDVX127iEgmi3eq3z2xIjuXEQrrOGGA3Y/Aue5+X5zvU2c91tjZ/1DWHgxUt8+PQN94PlMaubIyOP10uPPOcACw664wZEhKRuAVFVUeXwhaZEdEskO8Z/64+81mdhuwI2Ge/4/Ax+6+MlnBSYZbtQpOPBHGjAkX2q+6Kkzx65ua47ryYwqN9heRbFPbwj7fAH9x9w/L22KJfmoqApMMt3w5/OMf8OSTcPXVMGhQJGFokR0RyUa1nfnnARukKA7JJkuXwlFHhXK9I0fC2WdHHZGISFaJu9tfJCEWL4YjjoBXX4U77oBTTok6IhGRrFNX8q+ukp7I+vn5ZzjkEHjvPXjwQTjuuKgjEhHJSnUl/yvMbGEc7+PufmIiApIMtXAh9O4NH38Mjz0Gf/lL1BGJiGStupJ/dyCeFfvUQyA1mzsXDjwQvvkGxo5du0SviIhEoq7k/2d3fzclkUhmmjULDjggHAC88AL07Bl1RCIiWS/e8r4i9ffVV7D33rBgAUyYEFfiL1/Rt0mTcK9FdkREEk+j/SU5Pv00dPWvXAmvvAI771znS6qu6Fu+yh5oLr6ISCLpzF8Sb+pU2HffULP/1VfjSvygVfZERFKlxjN/d9eBgdTf22+HAX1t28J//wtbbx33S7XKnohIaijBS+JMnAi9ekGHDvDaa/VK/FDzanpaZU9EJLGU/CUxxo0LBXxyckLiz82t91sUFYVV9SrSKnsiIomn5C8N9+SToWTvdtuFs//OndfrbQoKoLg4HDeYhfviYg32ExFJNI32l4YZMwZOOAF69Ajz+Nu1a9DbaZU9EZHk05m/rL+77oK+fcNc/vHjG5z4RUQkNZT8Zf3cfHNYke+gg+D552HjjaOOSERE4qTkL/U3fDicc05YnOepp2DDDaOOSERE6kHJX+LnDkOGwODBYTneRx+FDTaIOioREaknDfiT+LjDeefBTTeF7v7Ro6Fp06ijEhGR9aAzf6lZ+So7ZtCmTUj8Z58d5t8p8YuIpC0lf6leSQkl/5xA3syJNKGMvF8/osT6hil9ZlFHJyIiDaDkL9UqOecd+q+8hZnk4TRhJnn099GUnPtu1KGJiEgDKfnLusrKKPzhPErZqFJzKRtR+MN5EQUlIiKJouQvlb3/Puy+O7OofjWdmtpFRCR9KPlL8MsvYe5+jx7w7bfktP6p2t1yNitNcWAiIpJoSv7Zzh0efzwsyvPvf8Ppp8Nnn1E0ejNatVhVaddWLVZRNLJ1RIGKiEiiKPlns+nT4bDD4G9/g44d4a234NZboV27sMLe3c0qr7B3dzMtuiMikgFU5CcbrVgBI0bAlVdCkybh8VlnQbPK/xy0wp6ISGZS8s82b7wRuvY/+STU5h85En7726ijEhGRFFK3f7b44YdQlnfvvWHxYhg7Fv7zHyV+EZEspOSf6dzh3nthm23C/UUXwaefwuGHRx2ZiIhERN3+mWzaNBgwAF59FfbcMyzG8/vfRx2ViIhETGf+mWjp0rD07h/+AP/7X1iI5/XXlfhFRATQmX/mefFFGDgQvvkGjj8ebrghTOMTERGJ0Zl/ppg7F/7+dzj44DBl77//hfvvV+IXEZF1KPmnu7KyUJhn223hqafgiitCV//++0cdmYiINFLq9k9n778f5uy/9x4ceCCMGgW/+13UUYmISCOnM/80UTLwDfKazaaJrSav2WxKdromLMIzaxaMGQPjxyvxi4hIXHTmnwZKBr5B/9t2ppSNAJhZ1oX+H50F26+mYNJAaNcu4ghFRCSd6Mw/DRQW561J/OVK2YjCz09Q4hcRkXpT8m/s5s9nVtnm1W6qqV1ERKQ2Sv6N1bJlcO21sPXW5DCr2l1ymn6X4qBERCQTKPk3Nu7wyCNh6t4ll0DPnhT9/SNasaTSbq1YQlH/GdHEKCIiaU3JvzF5+234059CsZ527WDCBBg7loKHDqd4wFRym87GWE1u09kUD5hKwai9oo5YRETSkLl71DEkVX5+vk+ePDnqMGo3cyYMGgQPPQSdOkFREZx4IjRtGnVkIiKSZsxsirvn17aPpvpF6Zdf4JprYMQIMINLLw1L7rZuHXVkIiKSwZT8o7BqFdx9d0j28+dD375w9dXw299GHZmIiGQBJf9UGz8ezj8fPv4Y9toLnn02VOoTERFJEQ34S5VPP4VDD4WDDoLSUnj8cXjtNSV+ERFJOSX/ZFuwAAYOhJ12gjffhBtuCAcCf/1ruM4vIiKSYur2T5bly+Hmm2HYMFiyBAYMgMsvh/bto45MRESynJJ/ormHLv2LL4bp06FPH7j+ethuu6gjExERAdK829/MBprZdDNbZmZTzGzvSAN6990wiO+YY8J0vfHjw4A+JX4REWlE0jb5m9mxwEjgamBn4E3gBTPLSfZnlwx8g7xms2liq8lrNpuS41+AggLYbTf4+mu44w6YOhV69Up2KCIiIvWWthX+zOwd4H/ufmqFti+Bx919UHlboiv8lQx8g/637Vxpid1WLKG4yekUDMoN3f0bb5ywzxMREamPeCr8peWZv5m1AHYFxlfZNB7YM5mfXVicVynxA5SyEYU2PAzuU+IXEZFGLi2TP9AeaArMq9I+D+hUsWHBggXk5+evuRUXFzfog2eVbV6vdhERkcYm3Uf7V71mYVXbOnToQCK7/XOafsfMsi7VtsO67SIiIo1Nup75LwTKqHKWD3Rk3d6AhCrqP4NWLKnU1oolFPWfkcyPFRERSZi0TP7uvgKYAlQdTt+LMOo/aQpG7UXxgKnkNp2NsZrcprMpHjCVglF7JfNjRUREEiadR/sfCzwADAQmAacDJwM7uPvM8v0SPdpfRESkMYtntH/aXvN390fMbDNgCNAZ+Bg4tGLiFxERkXWlbfIHcPdRwKio4xAREUknaXnNX0RERNafkr+IiEiWUfIXERHJMkr+IiIiWUbJX0REJMso+YuIiGQZJX8REZEso+TfAA1dITBdZMP3zIbvCPqemSQbviNkx/eM4jsq+TdANvyjhOz4ntnwHUHfM5Nkw3eE7PieSv4iIiKSdGm7sE+8zGwBkKx6/+0Jywtnumz4ntnwHUHfM5Nkw3eE7Pieif6Oue7eobYdMj75i4iISGXq9hcREckySv4iIiJZRslfREQkyyj5ryczG2hm081smZlNMbO9o44pUcxskJm9Z2a/mNkCM3vGzHaMOq5kM7PBZuZmdkvUsSSamXU2s/tif5/LzOxTM9s36rgSxcyamtlVFf5PTjezYWbWLOrYGsLM9jGzsWY2J/Zvs1+V7WZmQ83sOzNbamYTzWyHiMJdL7V9RzNrbmbXmtn/zGyJmc01szFmlhNhyOulrr/LKvsWx/a5IFnxKPmvBzM7FhgJXA3sDLwJvJCO/yBr0BMYBewJ7A+sAiaY2aZRBpVMZrY7cCrwv6hjSTQzawdMAgzoA2wHnAXMjzKuBLsYOAM4G9gWOCf2fFCUQSVAa+BjwvdZWs32i4DzCX+fPQh/py+Z2cYpi7DhavuOrYBdgKLY/ZHAb4FxaXhgV9ffJQBmdjTh7/K7pEbj7rrV8wa8A9xRpe1LYHjUsSXp+7YGyoDDo44lSd+vLfA14UBnInBL1DEl+PtdDUyKOo4kf8dngfuqtN0HPBt1bAn8jr8C/So8N2AuUFihbUNgMXBa1PEm4jvWsM/2gAO/jzreRH9PIBeYQzhAnwFckKwYdOZfT2bWAtgVGF9l03jCmXIm2pjQS/RT1IEkSTHwuLu/HHUgSfJn4B0ze8TM5pvZB2Z2pplZ1IEl0BvAfma2LYCZbU84mHs+0qiSa0ugExV+i9x9KfAamftbBNAmdp9Rv0exnoyHgGHuPi3Zn5du3SaNQXugKTCvSvs84MDUh5MSI4EPgLeiDiTRzOxUYGvg+KhjSaKuwEDgRuAaoDvw79i2TBnfcC3hIPVTMysj/LYVufuoaMNKqk6x++p+i7ZIcSwpETv5+j/gGXefHXU8CXYF8IO735aKD1PyX39VqyNZNW1pz8xGAHsBe7l7WdTxJJKZbUPoEt/b3VdEHU8SNQEmu3v59e+pZvY7wjXxTEn+xwInAMcBnxAOcEaa2XR3vyvSyJIvW36LmgEPAu2AIyIOJ6Fig2/7Ef7dpoS6/etvIeH6d6cq7R1Z9wg8rZnZjcA/gP3d/Zuo40mCPQg9OR+b2SozWwXsCwyMPd8g2vASZi7waZW2aUCmDFAFuB64wd0fdveP3P0BYATpP+CvNt/H7rPht6i8S3wn4AB3/yHikBJtP6AzMLfCb1EucK2ZJaWHQ8m/nmJniFOAXlU29SKM+s8IZjaScBa1v7t/FnU8SfIU8HvC0Xb5bTLwcOxxpvQGTAK2qdLWjeSteRGFVoSD8orKyOzfuOmEA4A1v0Vm1hLYm8z6LWoOPEJI/Pu5+/d1vCQdjSJ8v4q/Rd8RLtUdkIwPVLf/+hkBPGBm7xJ+WE8HNgdGRxpVgpjZrYRr4H8GfjKz8jOLX9391+giSyx3/xn4uWKbmS0BfnT3j6OJKiluBN40s0LCj+jOhClxgyONKrGeAS4xs+mEbv+dgfOA+yONqoHMrDVhTAqEA5kcM+tO+Dc6y8xuAgrN7DPgC3oIKGQAAAimSURBVGAIYST5mEgCXg+1fUdCAnyMMPXtcMAr/B4tig1wTAt1/V1SZeqtma0Evnf3z5MSUNRTHtL1RhhANQNYTugJ2CfqmBL43byG29CoY0vBd59Ihk31i32vPsCHwDJCkjib2MJemXAjDPa7idCbsRT4hjCeo2XUsTXwe/Ws4f/ivbHtBgwlXNpZBrwK7Bh13In6jkBeLb9H/aKOPZF/l9XsP4MkTvXTqn4iIiJZJpOvh4mIiEg1lPxFRESyjJK/iIhIllHyFxERyTJK/iIiIllGyV9ERCTLKPlL1jGzfmbmZvazmW1SZVuz2LahEcQ1NPbZjbr4lpk1MbObzGyuma02s6ca+H55/9/e2cZoUV1x/PdPCU2l6oJWm1hT2qLR6gfSfmhKC/qhL6RLNbTQmEKFRkraKCFtE4FaJFLMJjWAdomCtr4FpKjgS1RYKooRXSu0NbSN0tqwNVVpu+xCF6yL1NMP544Ms/M8+zzLrqDP+SWTyX05d87c2X3O3HPvnZPue+YgqZhvuyn162cGu+2hJvf38Lf0lbt82Zih6rOgMQjjHzQypwLzjrcS70GmAHPx7+l/Abj6+KpTlSZgEfCeM/45PglccbyVCN5fhPEPGpnNwJzc50Lf9wxSsKLz0/lGM2s3s78MQpsnBCdoMKfNwE/Td/uDYFAI4x80MkvS+ZpqlTL3a0n+nZI6cunMff19SS2S9kjqkbRa0knJVdsm6YCklyXNqHDJ8yU9KemN5FpfLOmo/1VJp0u6RdKrknolvSRpdqFONr0xQdJ9kvYBv+3nXidKapf0X0n7JT2YQh9n5R3452QB/tef6zlNoyxI+vVKek3S0loMmaSLJG1JfXgw9d2FJfUmS3om9et/JD0v6RJJo/HgNwC3JV3f0VfSVknbJH1d0h8k9eKf7UbSKZJWJH17Je2S9ENJyl334tTeJalup6R/p+fdVNBxrqQXU792S9ohaXJ/fZC4Bo8dcmWN9YOgX8L4B43M63g8+9mSPj6I7S7Af6xnANfiseZXAg8AjwKTgZ3AHZIuKJF/EHgcD6x0D7AwtQO4YcIDSjXjhrgZD2xzi6Q5Je2twY3gFGB+JaUlTUz6HUg6/wC4ENgm6axUbTL+zXXwkMifTzKVWI0Hm7kn6dmCu7DXVJFBUjOwJekyHY8weTLwtKSzc/XmABvwoCgzgKl4P4/Gn+83UtWWCvqeC/wCaAW+CmxJL1qPAt8FluIBZTbhAb2uL1H3Jvwb7d8GFgPfTHmZjtNSO2uBrwHTgPuBUdX6IMcLeHCb+ZJOrlEmCKpzvIMdxBHHu30AM/Ef6zH4D/A+4PZUNoxCECPcwFpJO3cCHbn06CT7RKHehpQ/PZc3EjgMLCpeB5hfkL8N6AGaUnohHsTlnJJ6ncCwwn0ur7FfdgB/zeRT3ieAt4BlubwlZf1R0t74dP3LC/nTUv7YQr/NzNV5GdhSkDsl3d+NuXQPsKGKDlnbs0rKtgJvZ3rk8icV9Un5v8QDeZ2e0henencV6q1Iz0e59O8H8Hea/T0Mw19SDgPXprIxZTrGEUetR4z8g4bGzLrwUdnleff2MbKxkH4pndty1+3GR6tn05d7C+lfAx/GR+EAE3H3/e7kVh8m3yHQBpwGfLog/0B/CksagS+KW2dmh3N67sa9DBf110YJE4FDwPqCnptT+YQKupwDfApYU5B7A2jPyY3D++XWAeiW0WFmLxTyJuAvBWsL+auB4bj3IE/R8/FH4IPAmSm9HRgrqVXSlySdVK+S5usq7gJ+LKlWj0EQVCSMfxB4vPsu3GU7GHQX0oeq5JfNff+zQjpzvZ+BG6i3Csd9qfy0gvzr/avMSDw8bFndPdTuos5zBm4sDxT0zOKWF/XMywH8ir73OCknl53/MQDdMsrudxQeY723kL8nV56nq5DO5LJnezc+hfI5/AWtS9KGtCahHq7DXypih0pwzJzQ+4mD4N3AzA5IasE9ADeUVHkTQNJwMzuUy69kvI6VM/F49Pk0wKvpvBc3oHMryO8qpGuJ292d6pXtfPhouma97MX7bnyF8teqyIGvnXi8pDx7Bp3pfBbwpwHoB+V90wWMKnneWd/U1RdmZsAqYJX8uxJfwf/W1uEvBLW284qkVcBVwCP16BAERWLkHwTOzbhxXVJS9vd0fmeleVrNPW6IdPlWIX0ZPnrODNwm4DzgFTPbUXL01HtBMzsI/A6YKukDWX5aCDkOeGoA97EJH/2eWkHPSsZ/F9ABXFBBbmeq9yzeL7MrtANHRuEfqkPvp/DfxqmF/Gn4i8dzdbR1FGbWbWbr8KmdPjsXauB6/IWl6g6VIOiPGPkHAWBmvZIWUz5/vBHYj28XW4S7Xq/GDc9Q8L204nw7vgJ9Fr4AcV8qX46vxn9a0nLcWI7AXwjGm9mlA7zuQnz++hFJN+Pz6dfh97603sbMbKuktcD9kpYBz+Nz6aPxVe/zrOQbAWZmkq4EHpI0HDeUnbgHZBz+0rPMzHokLQBaJa3HdxD0AGOBN82sFZ8y2QtcJmkncBDYbWbVRu8bgW3ASkkfAf6c9J0FtJhZZxXZPki6NenVjntszgW+w5G1DzVjZv+SdBPwk3plgyBPjPyD4Ah34KvdjyIZ3Um44boX3zbWCjw5RHpcCnwZeBjf5rYE+FlOn/24EXwMn/9tA25PcgPWycw24dvxmvD7XAm8CHyxyii9P6bjq9anAA/hW9yuwvu5uLYhr8tj+LqGEfgq+zbg57jrvT1XbwU+Qv8YbvzXp2vtTuVv40Z7JD6FsB3fuleRJNOML7Cbh78QNQM/YmAj7meAz+Lepd+kNlbjWxMHwg30XT8SBHWRbUUJgiAIgqBBiJF/EARBEDQYYfyDIAiCoMEI4x8EQRAEDUYY/yAIgiBoMML4B0EQBEGDEcY/CIIgCBqMMP5BEARB0GCE8Q+CIAiCBiOMfxAEQRA0GP8Hk5jEAAAynRgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1)\n",
    "fig.set_size_inches(8, 6)\n",
    "ax.tick_params(axis=\"x\",direction=\"in\",which='major',labelsize=14)\n",
    "ax.tick_params(axis=\"y\",direction=\"in\",which='major',labelsize=14)\n",
    "ax.set_xlabel('Number of electrons N',fontsize=16)\n",
    "ax.set_ylabel('Total Energy E(N) [meV]',fontsize=16)\n",
    "\n",
    "electrons = np.array([0,2,4,6,8,10,12,14])\n",
    "electrons2 = np.array([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14])\n",
    "E_N_paper = np.array([0,50,105,190,290,390,510,620,750,900])\n",
    "\n",
    "plt.errorbar(electrons2 , 11.61 * E_N, yerr = 2*np.sqrt(N_av) * 11.61 * E_Nvar, fmt='bo', label = 'RHF')\n",
    "ax.scatter(electrons2[:10],E_N_paper,c='r',label = 'UHF')\n",
    "ax.plot(electrons2[:10],E_N_paper,'r')\n",
    "ax.legend(fontsize=16,loc = (0.03,0.8))\n",
    "ax.text(0.3, 685, 'std 95%', style='normal',\n",
    "        bbox={'facecolor': 'white', 'alpha': 0.8, 'pad': 7}, fontsize = 13)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAF7CAYAAADPB8KEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XecVNX9//HXBxABO4JiY1csWDAxBiNWLLGCscWSYI26UdAYv8nPhlFjQE3ssUTRGElEWbuLBWvU2MVodomKUSkqIKioIAgCn98f564Ow87u3d2ZuTN33s/H4z5m594793xmFvYzp9xzzN0RERGR9OiQdAAiIiKSX0ruIiIiKaPkLiIikjJK7iIiIimj5C4iIpIySu4iIiIpo+QuIiKSMkruIiIiKdMp6QDao0ePHl5dXZ10GCIiIkXz2muvfeLuPZs7p6yTe3V1NRMmTEg6DBERkaIxs6ktnaNmeRERkZRRchcREUkZJXcREZGUUXIXERFJGSV3ERGRlFFyFxERSRkldxERkZRRchcREUkZJXcREZGUUXIXERFJGSV3ERGRlFFyFxERSRkldxERkZRRchcRkVQZMwaqq6FDh/A4ZkxllJ2prJd8FRERyTRmDNTUwPz54fnUqeE5wJAh6S07m7l7cUvMo/79+7vWcxcRkUbV1SGpZuvUCTbaqLBlv/ceLF68/P6qKpgyJX/lmNlr7t6/uXNUcxcRkdSYNq3p/YsXw9ZbF7bsSZOa3p8rpkJSchcRkdTo3bvpmntVFYwdW9iyX3qp6bJ79y5suU3RgDoREUmNkSPDYLZM3bqF/cUou1u3ZMrOpuQuIiKpMXgwmMGqq4bHqioYNao4A9qGDAllVVUVv+xsapYXEZHUGD8eliyBhx+GHXcsfvlDhiSTzLOp5i4iIqlRVwc9esCAAUlHkiwldxERSYVvvgk19sGDoWPHpKNJlpK7iIikwvPPw+efw09+knQkyVNyFxGRVKirg86dYc89k44keUVN7mY2xcy8ie2hjHOGmtlkM/vazF4zs52LGaOIiJQf95Dc99gDVl456WiSV+ya+7bAOhnbNoADdwKY2eHA1cBFwA+AF4BHzCyBKQBERKRcvP12mP5VTfJBUZO7u89295mNG7Af8CVwV3TK/wG3uvtN7v6Wu58KzABOLmacIiJSXurqwuPgwcnGUSoS63M3MwOOB25z9/lm1hn4IfBY1qmPATsUOz4RESkfdXWwzTaw/vpJR1IakhxQtyewIXBz9LwH0BH4OOu8j4FeRYxLRETKyKxZ8OKLapLPlOQMdScCr7r7G1n7s9egtSb2ATB79mz69/9u1buamhpqGhfPFRGRivDQQ2FAnZL7dxJJ7ma2FnAAMCxj9yfAEpavpa/F8rV5AHr27InWcxcRqWzjxoXm+EIv6VpOkmqWPw5YCHy7AJ+7LwJeIzTXZ9qTMGpeRERkGV9/DY8+GmrtZklHUzqKXnOPBtKdAIx197lZh68A/mFmrwDPAycB6wI3FDdKEREpB089BfPnw/77Jx1JaUmiWX5XYGNguXVz3L3WzNYEziXcBz8R2M/dpxY1QhERKQt1dWHSmt12SzqS0lL05O7u/yQMkst1/Hrg+uJFJCIi5cg99LfvvTesuGLS0ZQWzS0vIiJl6d//hunTNUq+KUruIiJSlurqoEMH2G+/pCMpPUruIiJSlurqYIcdoEePpCMpPUruIiJSdj74AN54Q03yuSi5i4hI2Rk3LjwquTdNyV1ERMpOXR1suin07Zt0JKVJyV1ERMrKl1+GyWs0cU1uSu4iIlJWHnsMvvlGTfLNUXIXEZGyMm4cdO8eRspL05TcRUSkbCxeHJZ4HTQIOiW5aHmJU3IXEZGy8eKL8Omn6m9viZK7iIiUjbo6WGGFMJ+85KbkLiIiZWPcuLAC3KqrJh1JaVNyFxGRsjBpUtg0Sr5lSu4iIlIWGmelGzw42TjKgZK7iIiUhbo6+P73oaoq6UhKn5K7iIiUvE8/heefV5N8XEruIiJS8h5+GJYuVXKPS8ldRERKXl0drLsubLNN0pGUByV3EREpaQsXwvjxYSBdB2WtWPQxiYhISXv6aZg3T03yraHkLiIiJW3cOOjWDXbfPelIyoeSu4iIlCz30N++117QtWvS0ZQPJXcRESlZ//kPfPCBFoppLSV3EREpWXV1YBaWeJX4cq6Ga2bT2nA9Bwa5+8S2hyQiIhKMGwcDBsDaaycdSXlpbqn79YGHgdkxr9UBOBLo3N6gREREPvoIJkyAiy9OOpLy01xyB7jQ3V+JcyEz6wQc1f6QRERE4MEHw6NugWu95vrcfwd8EPdC7r44es1H7Q1KRESkrg769IHNN086kvKTM7m7+0h3n9Gai0Wv+bj9YYmISCX76it48slQazdLOprykzO5m9ktZrZLMYMREREBePzxMO2smuTbprlm+cOBf5rZZDP7vZltXKygRESkstXVweqrw047JR1JeWouua8NnABMAc4FJpnZc2Z2opmtVozgRESk8ixZEgbT7bsvrLBC0tGUp+b63Oe5+9/cfTegmjBYrjtwIzDDzO4ws33NTBPhiIhI3rz8MsyerSb59oiVmN39A3e/yN23AAYAtwB7AA8CH5nZZXELNLN1zGy0mc02s6/N7E0zG5hx/FYz86ztpVa+LxERKVPjxkGnTrDPPklHUr5aXet291fc/RRgPeBKYC3g9DivNbPVgecBAwYBmwOnArOyTn0CWCdj26+1cYqISHmqq4OBA0Ofu7RNS5PYLCcaWHc0YTa6KmAucFfMl58BzHD3ozP2TW7ivIXuPrO1sYmISHl79114802oqUk6kvIWq+ZuZt3N7GQzexGYBJwDvENI8L3c/cSY5R0IvGxmtWY2y8zeMLNTzJa7i3Gn6Pg7ZnaTma0V9w2JiEj5GjcuPGoVuPZpbuGYFYDBhFr6voQ5498EzgZuc/fpbSivDzCU0Jx/CbA1cE107NrocTxwL6FGXw2MAJ4ysx+6+8I2lCkiImVi3Djo1y/MTCdtZ+7e9AGzT4HVgc+AscBod5/QrsLMFgET3H2HjH0XAQe5e5MTDJrZusBU4HB3vzfzWFVVlffs2fPb5zU1NdSoLUdEpCzNmQM9e8KZZ8LIkUlHU7rM7DV379/cOc31uf8LGA086O7f5CmmGYTaf6a3gNNyvcDdp5vZh8Am2cd69uzJhAnt+r4hIiIl4pFHwj3uugWu/XImd3c/sKn9ZrYysCYwvQ1J/3mgb9a+TQk18yaZWQ/CyPxWzXMvIiLlpa4urNu+7bZJR1L+Yt8KZ2aDzezfwBfAe8BW0f6bzeznMS9zJTDAzIab2cZmdijwK+C66Form9llZra9mVWb2a7AOMKtcvfFflciIlJWFi2C8eNh8GDooKnR2i3uaPkDgQeAT4Azs143GTgmznXc/VXCiPnDgInASMLMd9dHpywhfGl4gDAafzRhdP727j43ThkiIlJ+/vUv+OILNcnnS9z73M8H/ubuJ5hZJ+BPGccmEkbAx+LuDwEP5Ti2ANg77rVERCQd6uqgSxf48Y+TjiQd4jZ+bA7URj9nD6+fQ+iDFxERaTX3kNx//GPo1i3paNIhbnL/EuiR41g1MDsv0YiISMWZOBGmTFGTfD7FTe6PA2dHc8M3cjNbETgFeCTvkYmISEVonJVu8OBk40iTuH3uw4FXCIPbHiY0zZ8FfA9YjTBITkREpNXq6uBHP4J11kk6kvSIu+TrFGAbwhKvexJGte8CvARs18apaEVEpMLNnBnWb9dc8vkVe1U4d/8QOL6AsYiISIV58MHwqP72/NJUASIikphx46CqCrbaKulI0iV2zd3M9gEOBTYAumQddncfmM/AREQk3ebPh8cfhxNOgOUW/pZ2iTtD3RmEgXSDgZUIfe6Z29JCBSgiIun05JOwYIGa5Ashbs39FOBG4BR3X1LAeEREpELU1cGqq8IuuyQdSfrE7XNfFbhLiV1ERPJh6dIwmG6ffaBz56SjSZ+4yf1RYEAhAxERkcoxYUK4DU5N8oXRmmb5+8zMgccI88kvw93fz2dgIiKSXnV10LEj7Ltv0pGkU9zk7sBcwhKtI3Kc0zEvEYmISOrV1cFOO0H37klHkk5xk/utwA7AlcDbwKJCBSQiIuk2ZQo0NMDllycdSXrFTe67EkbK31q4UEREpBI0LhSjKWcLJ+6Auk+AjwsZiIiIVIa6Oth8c9hkk6QjSa+4yf3PwFAz03S1IiLSZl98AU8/rVHyhRa3WX4NoB/wppk9zvKj5d3dz89rZCIikjrjx8PixWqSL7TWrOfeaNMmjjug5C4iIs0aNw569IABmjmloGIld3dXc7yIiLTLN9/AQw/BgQeGe9ylcJS0RUSkKJ5/Hj7/XP3txaDkLiIiRVFXF+aR33PPpCNJv5zJ3cyWmNmP4l7IzDpGr9kmP6GJiEhauIfkvscesPLKSUeTfs31uRvQy8x6t+Ja1v6QREQkbd5+G957D37726QjqQwtDai7r5XX87YGIiIi6VVXFx4HD042jkrRXHI/ro3XnNzG14mISErV1cE228D66ycdSWXImdzdfXQxAxERkXSaNQtefBHO12woRaPR8iIiUlAPPxwG1OkWuOJRchcRkYKqqwvN8VtvnXQklUPJXURECubrr+HRR0Ot3XQ/VdEouYuISME89RTMn6+FYopNyV1ERApm3Lgwac1uuyUdSWVRchcRSakxY6C6Gjp0CI9jxhS37KoquOEGWLoU7r67eGVLM7fCmdktrbiOu/vxcU40s3WAS4D9gFWA94GT3f2Z6LgRlo+tIawj/zIwzN3/24p4REQq2pgxUFMTmsQBpk4NzwGGDClu2fPnF69sCcy96UnlzGwK8Wecc3fv02JhZqsD/waeA64FZgN9gOnu/lZ0zpnAucCxwCTgPGAnoK+7z828Xv/+/X3ChAkxQxQRqRzV1SGhZ1tlFTj66MKW/fe/w9y5y++vqoIpUwpbdiUws9fcvX+z5+RK7oVgZhcBA919xxzHDZgOXOvuI6N9XYFZwG/d/cbM85XcRUSa1qFDuLe8KWuuWdiyP/206f1moYle2idOcm9pbvl8OxAYb2a1wG6ERH4zcJ2HbxkbAr2Axxpf4O4LzOxZYAfgxuUvKSIi2Xr3brrmXozac65Wg95xlyGTdmv1gDozW8vMemdvMV/eBxhK6GffG7ia0P8+LDreK3r8OOt1H2cc+9bs2bPp37//t9uoUaNa+3ZERFJp5EhYccVl93XrFvYXo+xu3ZIpW4JYNXcz6wCMAH4JrJ7jtI4xLtUBmODuZ0fPXzezTQjJ/dqM87Ibk6yJffTs2RM1y4uILG/IELjpJnj22dA8X1UVkmsxBrQ1ljF8OEybFmrsxSpbgrjN8r8mJOA/EpL8SGApMCR6vCTmdWYAb2btews4Lfp5ZvTYC/gg45y1WL42LyIiOSxaBP/5T0io//hH8csfMkTJPElxm+WPAy4kJHeA+9z9fGBz4CMgbrP880DfrH2bAo29M5MJCX7PxoNm1gXYGXghZhkiIhXvscfg88/hiCOSjkSSEDe59yE0py8BFgNdAdz9G+Aq4Bcxr3MlMMDMhpvZxmZ2KPAr4Lroeh5d7ywzO9jM+gG3AvOA22OWISJS8WprYY01YM89Wz5X0idus/wXQJfo5+mE2vfzGdfoHuci7v6qmR0IXAT8DpgWPV6fcdqfCF8eruO7SWz2yr7HXUREmrZgATzwABx6KHTunHQ0koS4yf11YAvg0Wj7vZktINTiRxImponF3R8CHmrmuAMXRJuIiLTSI4+ESWQOPzzpSCQpcZP7VYSmeQhTw24DNM5SPBU4Jc9xiYhIG9XWQo8esPvuSUciSYmV3N398YyfZ5rZj4CNgG7AW1Hfu4iIJOyrr+DBB8MUs52KPU2ZlIw2/eqjpvN38xyLiIi004MPhoVa1CRf2ZpbFW4X4N/uPi/6uVnu/mxeIxMRkVYbOxZ69YKdd046EklSczX3p4EBwCvRz7lWmGmcPS7ODHUiIlIgX34ZBtP98pfQUX+RK1pz97nvxnezye3ezLZb9CgiZWDMmLCwR4cO4XHMmJZekY6yK8EDD8DChWqSl2Zq7u7+TMbPTxclGhEpqDFjoKYm9MlCWLnrxBNhyRL4+c8LW/btt8NJJ4V7sBvLrqkJP2ua0vyorYUNNoABA5KORJIWaz13M3sfOMjd/9PEsX5Anbv3Wf6VhaX13EVaJ9dSnEkqxhKkleCzz2DtteG00+Cyy5KORgopn+u5VwMr5jjWBahqRVwikpBp03If+8MfClv2737X9P7mYpL47rsPFi/WXPIStOZWuFxV/P7A53mIRUQKrHfvpmvuVVVw7rmFLfvmm5suu3fcZaekWbW10KcP/PCHSUcipSDngDozO93MppnZNEJiH9f4PGObTZgDfnyxAhaRths5Erp2XXZft25hfzHK7tZt2X1duxan7LSbPRueeioMpDNLOhopBc3V3N8Hnox+PgaYAMzOOmchYUT9zfkPTUTybcgQ+OgjOPPM8LyqKiTXYgxoayxj+PDvavCHHabBdPlwzz1hUKRGyUuj5kbLPwA8AGDhq+CF7j65SHGJSIE0NoO/8QZ8//vFLXvIkLC5w7bbwquvhp9V22yfsWNhs83ge99LOhIpFbHWc3f345TYRdKhvj7MOb7ZZsnFYAbDhsGbb8Izz7R8vuQ2fTo8+6ya5GVZsQfUmVkf4DCgN9+t7d7I3f34fAYmIoXR0AB9+8KKue5/KZIjjoDf/hauuw523TXZWMrZ3XeH1g81yUumWMndzA4A7iLU9GcR+toztXyzvIiUhPp62HHHpKMIg+l+8Qu48sowDmC99ZKOqDzV1sJWW8HmmycdiZSSWM3ywAjC/PLruPu67r5h1lb0CWxEpPW++CLcV14qfbMnnwxLl8KoUUlHUp6mTYMXXtC97bK8uMm9D3CZu2ePlheRMtLQEB632irZOBr16QP77huS+zffJB1N+bnzzvCoJnnJFje5vw2sWchARKTwGpN7qdTcAYYOhZkzwwxr0jq1tWHSmo02SjoSKTVxk/sZwDnRoDoRKVP19bDaarD++klH8p199oENNwwD6yS+996DCRNUa5emxR0tfwGh5v6Wmf0P+CzruLv7wHwGJiL5V18fau2ldMtUx46h7/2MM0LLQql0GZS62trweNhhycYhpSluzX0JMAl4gTBL3ZKsbWlBohORvHEv3eT5i19Aly5w/fVJR1I+amth++3DLIMi2WLV3N191wLHISIFNnUqzJ1bWv3tjdZcM4z4/sc/4JJLQteB5Pb226EV5qqrko5ESlXcmruIlLlSHEyXadgw+Oor+Pvfk46k9NXWhq6VQw9NOhIpVbGTu5mtZ2ZXmNkEM5tsZv2i/b82s+0KF6KI5EN9fXjs1y/ZOHLp3x9+9KPQNO+aFisn9zCX/C67wLrrJh2NlKpYyd3MtgQagKOA6YQpaDtHh6uA0woSnYjkTUNDGJW+yipJR5Lb0KGhyfmf/0w6ktLV0BA+I42Sl+bErblfDrwFbAgcDGSOtX0BGJDnuEQkz+rrS3MwXabDDw/977otLrfaWujQAQ45JOlIpJTFTe47AZe4+zyWn0f+Y6BXXqMSkbz6+mt4553S7W9v1KULHH88PPAAfPhh0tGUHveQ3HffHdZaK+lopJTFTe7N3erWA1iQh1hEpEDeeguWLCn9mjvASSeF+eZvvDHpSErPa6+FyWs0l7y0JG5yfwU4Lsexw4Dn8xOOiBRC42C6Uq+5QxgXMGgQ3HQTLFqUdDSlpbYWOnWCgw5KOhIpdXGT+x+A/c3sMcKgOgd+bGajgYOAkQWKT0TyoKEhNHlvvHHSkcQzbBh8/DHcc0/SkZQO97BQzF57QffuSUcjpS5Wcnf3Z4ADCQPqbiEMqLsE2Bk40N1fLliEItJu9fWwxRah1lcO9torLIaiGeu+89JLYYlXjZKXOGLf5+7uD7n7JsCmhAF2m7t7H3d/pGDRiUheNDSUR5N8ow4dwnzzzz33XZdCpRs7FlZcEQ44IOlIpBzEvc/9PDNbF8Dd33X3F9x9UnRsHTM7r5BBikjbzZoVllQth8F0mY47LnQl6La4MBjyrrtg3301Na/EE7fmfj6Qa5HIdaPjLTKzC8zMs7aZGcdvbeL4SzFjFJEmlPq0s7l07w4//zncdht8/nnS0STruedgxgw1yUt8cZN7cwtErgEsbEWZk4B1Mrbs+sQTWcf3a8W1RSRLuSZ3CAPr5s+H0aOTjiRZtbXQtSsMHpx0JFIucg6vMbNdgd0zdv3SzLL/aXUFBgH/bUWZi919ZjPHF7ZwXERaob4+THhSjpOebLMNDBgQBtademroi680ixfD3XfD/vvDyisnHY2Ui+bGzg4Ezo1+dpq+z30R8Cbwq1aU2cfMPope+zJwjru/n3F8JzObBXwOPAMMd/dZrbi+iGQot8F02YYNg6OOgiefhD33TDqa4vvnP2H2bDXJS+uYx1h+ycyWAgPc/ZV2FWa2L7AK8DawFuHLw2bAlu7+qZkdAcwHJgPVwAigI/BDd1+u6b+qqsp79uz57fOamhpqamraE6JIqixZEmp7J58MV1yRdDRts3AhrL8+7Lgj3H9/0tEU3wknhGb5WbNC07yImb3m7v2bOyfWXa/unpfGsOzb5qLBcu8DxwBXuPvYjMMNZvYaMJXQ9H9v9vV69uzJhAkT8hGaSCq9916YV76ca+4rrhgS3J/+FO7z7t076YiKZ9EiuPfecPubEru0RmvWc+9mZqeY2V1m9qSZ3WlmQ82sW1sLjxai+S+wSY7j04EPcx0XkeY13iNebrfBZTvppPBYafPNP/44zJmjueSl9eLe594L+DfwZ6A/0A3YFrgWeM3M1m5L4WbWhdAsPyPH8R7AermOi0jz6uvDILQttkg6kvapqgojxW+6KTTTV4raWlh99TBjn0hrxK25/4lwy9vO7r6hu2/v7hsSZqpbHfhjnIuY2WVmNtDMNjSz7YC7gZWA0Wa2cnR8ezOrjkbrjwNmAfe18n2JCGEw3aabpqNJd9iwMLDs7ruTjqQ4vv46jDE46CDo3DnpaKTcxE3u+wJnu/syq7+5+wuEQXGDYl5nfeAOwr3u9xLujx/g7lOBJYR73h8A3gFGR+dt7+5zY15fRDLU15d/k3yjH/8YNtmkcmasGz8e5s7VKHlpm7jLSKwMTM9x7MPoeIvcPWfPkbsvAPaOGY+ItGDePHj//TCNaxp06ABDh8Lpp8Prr8MPfpB0RIVVWws9esDuu7d8rki2uDX3SYSlXptyJOHWNhEpIRMnhse01NwBjj0WunVLf+39q6+grg4OOQRWWCHpaKQcxU3ulwE/M7MnzOwXZravmR1nZo8CPwcuLVyIItIW5TztbC6rrw5DhsDtt4dR5Gn10ENh2l01yUtbxV3P/TbgJKAfcDPwEPBX4HvASe5+e8EiFJE2qa8PE9hUVSUdSX4NHQoLFsCttyYdSeHU1kKvXrDLLklHIuWqNeu5jyKsALclsHP0uJ6731Sg2ESkHRoH06VtPvatt4YddgjzzS9dmnQ0+ffll/Dww3DoodCxY9LRSLlq1X97d1/q7m+5+/PRYwr/a4mUP/fyn1O+OcOGwbvvhkle0qauLtwGpyZ5aY/WzFC3iZmNNrN3zOyr6PFWM9u4kAGKSOt99FHok07TYLpMhxwSVrlL48C62towl/722ycdiZSzuDPU7Qr8BxgMvARcHz3uT5gDfmChAhSR1kvjYLpMK64IJ54IDz4IU6YkHU3+zJkDjz4aau1p606R4or7z+dy4HWgyt2Pdvf/5+5HE1ZueyM6LiIlonFO+X79ko2jkH75SzCDG25IOpL8uf9++OYbNclL+8VN7lsAf4wWevlWNHPcHwmD60SkRDQ0wAYbwBprJB1J4WywQVgt7a9/DX3UaTB2LPTpA/2bXcxTpGVxk/uHQK7ZjTsDH+UnHBHJhzRNO9ucoUPhk0/grruSjqT9Zs+GJ5+Eww4LLRIi7RE3uf8R+L2ZrZe5M3p+PnBRvgMTkbZZtAjeeiu9/e2Z9tgD+vZNx8C6e++FJUu0vKvkR9zkPhBYBXjPzJ42s1ozexp4jzCv/K5m9vdoG12gWEUkhkmTYPHiykjuZqH2/vLL8NprSUfTPrW14YtKJfzepPDiJvedCKu2zQCqgB9FjzOApYRJbTI3EUlI42C6SmiWBzjmGFhppfKuvc+YAU8/HQbSqUle8iHWqnDR2u0iUgYaGsJiI337Jh1Jcay2Ghx5JIweDZdeCmuumXRErXf33WHiIY2Sl3zRnZQiKVNfD5tvXlmriQ0bFkbM/+1vSUfSNrW1oaVliy2SjkTSQsldJGXSPO1sLlttBTvvDH/5S/nNN//BB/D886q1S34puYukyGefwYcfVk5/e6ahQ+H998MMb+XkzjvDo5K75JOSu0iKpH3a2eYcfDCsvXb5DayrrYVttoGNtUqH5JGSu0iKVHJy79wZamrCcqmTJycdTTzvvw+vvqp72yX/lNxFUqS+Hrp3h3XWSTqSZPzyl2HBlb/8JelI4mlskj/ssGTjkPRRchdJkcbBdJV6r/R668GBB4b55hcsSDqalo0dCwMGQFVV0pFI2uRM7ma21MyWxNwWFzNoEVne0qUhuVfiYLpMw4aFgYW1tUlH0rxJk+A//9FAOimM5iaxuRDwYgUiIu0zZQp89VVl9rdn2nXXcJ//ddfBsccmHU1utbWhheXQQ5OORNIoZ3J39wuKGIeItFOlTTubS+N886eeGgarbbtt0hEtzz00ye+8c+hKEMk39bmLpER9fUhsW26ZdCTJO/poWHnl0r0tbuLEsHKfmuSlUGLNLQ9gZp2BfYG+QJesw+7uf8hnYCLSOg0NsNFGIalVulVXhaOOgltugcsugx49ko5oWbW1YVT/IYckHYmkVazkbmbrAs8B1YR++MaxuJl98kruIgmqr1eTfKZhw8ItcbfcAmeckXQ033EPyX333cOkOyKFELdZ/lJgNtCbkNi3A/oAI4F3o59FJCHz58O772owXaYtt4SBA0OCX7Ik6Wi+8/rr4XelJnkppLjJfWfgcmCGQn1VAAAgAElEQVR69Hypu09x9/OAu4E/FyI4EYnnzTfDrXCquS9r2LBwF8EjjyQdyXfGjoVOncJ0uSKFEje5rwlMd/elwFfAGhnHngJ2zXNcItIKjSPlVXNf1oEHwrrrls7AOvcwK92ee4aZBEUKJW5y/xBoHJLyHrBXxrEfAV/nM6ikjRkD1dVhwEt1dXie5nKl/DU0QNeu0EcdZMtYYYUw3/z48fDee0lHAy+/DFOnai55Kby4o+X/CQwE7gduBK4zs62Bb4C9o32pMGYMnHjid1NXTp0ans+fX9j5n++8E047bdlya2rCz0OGFK5cSYf6eujXDzp2TDqS0nPiiTBiROh7v+yyZGOprQ0L3BxwQLJxSPqZe8uT0JlZD6C7u78TPT8VOBzoBowHLnT3otfe+/fv7xMmTMjrNaurQ2ItFVVVoc9QJBd3WGutkDBuvjnpaErTYYfBE0+Ete67dUsmhqVLYYMNwqQ699+fTAySDmb2mrv3b+6cWDV3d/8E+CTj+TXANW0I6ALg/KzdH7t7r+i4RcdrCP36LwPD3P2/rS2rraZNy33siisKV+7//V/T+5uLRwTg44/hk080mK45w4bBXXdB795h3vnevWHkyOK1io0ZA7/5TfhdPf98eK4WOSmkuPe59wTWaKy5Zx3bFPgs+gIQxySWHYCXeZPKGcBvgGOj884DHjezvu4+N+b126V376Zr7lVVcPrphSv36qubLrd378KVKelQyWu4x/XBB2H2vk8/Dc+L2e01Zkwoa/788PyTT9TlJoUXt8/9euAz4JdNHDudMJo+bo/0Ynefmb0zqrX/GrjE3e+J9h0DzAJ+TpH69UeOXPY/IoRmvJEji19ux46FL1fKn+aUb9m554bui0zz54cv7F27Frbs009f9v91Y9nDhyu5S+HETe47AcNyHHsMuLYVZfYxs4+ARYRm93Pc/X1gQ6BXdD0A3H2BmT0L7ECRknvjf7bhw0OTeLGa77LLXWUVmDu3NBe9kNJSXw/rrFN6U6yWklzdW7NnJzcFrLrcpJDiJvc1gC9yHPuSUHOP42VCk/vbwFrAucALZrYlIbEDfJz1mo+Boq6bNGRIMt+oM8udNSsM7rvkkjB9pkguDQ1qkm9Jru62ddYJt8kV0j77wIwZTcckUihxk/uHhClnn2zi2HZAE/90l+fuy8wTZWYvAe8DxwAvNZ6W9TJrYh8As2fPpn//7wYM1tTUUNPYmVXm1lorNNNfdx2cd15I9CLZFi8Os9PtsUfSkZS2XN1tl15a+C9Gl16aTFefVLa4k9jcDZxjZoMyd0bPzwLubEvh7j4P+C+wCdDYD98r67S1WL42D0DPnj2ZMGHCt1taEnuj3/42TGjzpz8lHYmUqv/9DxYuVM29JUOGwKhRYWCsWXgcNao4LXRJli2VK25yvxCoB+rM7CMzeyXqN68DGoDft6VwM+sCbEao+U8mJPg9s47vDLzQluuXu/XXh+OOg7/+FaZPb/l8qTwaTBffkCFhzoilS8NjMZNrkmVLZYqV3N19PmGGuhOBZ4HPgWeA44GB0fEWmdllZjbQzDY0s+0ILQIrAaM9zKZzFXCWmR1sZv2AW4F5wO2te1vpceaZYUWrpGfWktLU0BDuqth886QjEZFSErfPHXf/Brgl2tpqfeAOwjz1swn97APcvXGoy5+ArsB1fDeJzV7Fuse9FG24IRx5JNxwA5x9NvTsmXREUkrq66FvX1hxxaQjEZFSErdZPi/c/Qh3X9fdO7v7eu5+iLu/mXHc3f0Cd1/H3bu4+0B3n1jMGEvR2WfD11/DlVcmHYmUmvp69beLyPJyJncze9/Mvh/9PDl6nmsrgfWW0qtv3zA39rXXwpw5SUcjpeKLL8LtXUruIpKtuWb5Zwj3sDf+3PIKM1Iw55wTVpS65ppwa5zIxKhNS4PpRCRbzuTu7sdl/HxsUaKRnL73vbDq11VXheksV1kl6YgkaZpTXkRyKWqfu7TP8OGhWf4vf0k6EikF9fWw2mphGVERkUyxR8ubWXdgELAB0CXrsLt79lKukmfbbgt77w2XXw6nnJLcutRSGhoaQpO8WdKRiEipibvk617APYR70pviLL9OuxTAuefCzjvDzTfDr36VdDSSFPdQc9dkKCLSlLjN8lcArwPfB1Z09w5ZW8fChSiZdtoJBg4MU9IuXJh0NJKUadPgyy/V3y4iTYub3KuBP7h7QzSZjSTo3HPho49g9OikI5GkaDCdiDQnbnJ/HVi3kIFIfHvsAdttBxdfDN/oq1ZFapxTvl+/ZOMQkdIUN7n/H3CGmW1fyGAkHrNQe58yBe64I+loJAkNDWEZ4FVXTToSESlFcZP7a4S13J8zsy/NbFrWNrWlC0h+DRoE3/8+XHRRWFhGKkt9vSavEZHc4t4KdxlwCqF5/m1gUcEiklgaa++HHgr33BOmp5XKsHAhTJoEBx2UdCQiUqriJvdjCQPqdLtbCTn44LDU54gR8NOfQgdNSVQR3nortNao5i4iucRNB05Yx11KSIcOYc75hgZ48MGko5FiaRxMp5HyIpJL3OR+F7BvIQORtjniCOjTJ9TeXUv7VISGhrB++yabJB2JiJSquM3yjwBXmtlqwHhguYVH3f2pfAYm8XTqFNZ7P/FEePxx2GuvpCOSQquvhy22CL97EZGmmMeo7pnZ0hyHHDDC3PJFn6Wuf//+PmHChGIXW3IWLYKNNoINN4Rn1XmSeuuuG77E3Xpr0pGISBLM7DV379/cOXG/+++Wh3ikQDp3hjPPhFNPDcl9l12SjkgK5ZNPYMYMDaYTkebFSu7u/kyhA5H2Of740O8+YgQ89ljS0UihaNpZEYmjVTdPmVkPMxtsZsdES8BiZl3MTDdhJaxrV/jtb0O/+8svJx2NFErjSHnV3EWkObGSsgWXAh8CdcAthMVkAB4AhhckOmmVk06C7t1h5MikI5FCqa+Hnj1h7bWTjkRESlncGvfZhBnqLgS2IwyiazQOGJznuKQNVl4ZTj8dxo2DN95IOhophIaG0CRv1vK5IlK54ib3E4AL3f0i4N9Zx94FNsprVNJmp5wSFhO56KKkI5F8W7IEJk5Uk7yItCxucl8PeCnHsUXASvkJR9pr9dXDqPm77w7TlEp6vP8+LFigwXQi0rK4yf0jINfK0d8HJucnHMmHX/86DLC7+OKkI5F80mA6EYmrNdPPnmdmO2bsczPbFPgNMDbvkUmb9egBJ58Mt98O772XdDSSLw0NYT2BLbZIOhIRKXVxk/sFhKVenwX+F+27C2iInl+S98ikXX7zmzA96R//mHQkki/19bDxxtCtW9KRiEipi5Xc3X0BsCth6dcXgCeAV4EaYE931/ruJWaddeCEE8IUpR98kHQ0kg/19epvF5F4Yk8+4+5L3P0f7n6ku+/l7j9z99HuvriQAUrbnXFGWCnu0kuTjkTaa968MKBOyV1E4tDMcinWuzcccwzcdBPMnJl0NNIe//1v+KKmwXQiEkfcGeo6m9n5Zva2mc03syVZm2rvJeqss8KqcVdckXQk0h6aU15EWiPuqnCXAsMI67rfCywsWESSVxtvDD/7GVx/fVg5bs01k45I2qK+HlZaCaqrk45ERMpB3OT+U+B8d9es5WXonHNgzBi4+mq48MKko5G2aGgITfId1JEmIjHE/VOxMvBiIQORwtliCzjkEPjzn+GLL5KORlrLPdTc1d8uInHFTe7jgF3yXbiZnWNmbmbXZuy7NdqXueWa+lZiGj48JPbrrks6Emmt6dPhs8/U3y4i8eVsljezPhlPrwH+bmZLgYeBz7LPd/f3W1OwmQ0ATgTqmzj8BHBUxnPdR99OP/gBDBoUBtaddlrov5XyoMF0ItJazdXc3yXMPvc/wsQ1mxBmqns5Y3/mFpuZrQaMAY4H5jRxykJ3n5mxLfdlQlpv+HD49FO48cakI5HW0JzyItJazQ2oO66A5Y4C7nb3p8zsvCaO72Rms4DPgWeA4e4+q4DxVITtt4c99giT2gwdCl26JB2RxNHQAOuvD2uskXQkIlIuciZ3dx9diALN7ERgY5Ztds80nnC73WSgGhgBPGVmP3R33YLXTueeC7vtBrfcEhK8lD4NphOR1srZLG/B/maWa6lXzGwrM9s/bmFm1he4CBiSaz56dx/r7nXu3uDu44B9gb7AoOxzZ8+eTf/+/b/dRo0aFTeUijVwIOy4Y1hQZpFGMpS8b76Bt95Sf7uItE5zzfJHAdcDzdUZ5gJ3mNmJ7n5HjPK2B3oAE82scV9HYBczOwlYKbt27u7TzexDQp//Mnr27MmECRNiFCuNzELtfd994bbb4Be/SDoiac6kSSHBK7mLSGs0N6DuSOBv7j451wnuPgX4K3BMzPLuJ3xZ2Dpjm0BYD35rmhgVb2Y9gPWAGTHLkBbsvTf88Idw8cWwWBMHlzQNphORtmguuW8DPBbjGk8A/eMU5u6fu/vEzA34Cvgs+nklM7vMzLY3s2oz25Vwj/0s4L44ZUjLGmvv774Ld96ZdDTSnIYGWGEF6Ns36UhEpJw0l9xXoenb1LLNic7NhyWEmv0DwDvAaGASsL27z81TGQL85CfQrx+MHAlLlyYdjeRSXw+bbQadOycdiYiUk+aS+ydAVYxr9I7ObRN339XdT4l+XuDue7v7Wu7e2d2r3P1Yd/+grdeXpnXoEO57f/NNuP/+pKORXBoa1N8uIq3XXHJ/jnh96cdG50qZOfRQ2GQTGDEizF8upWXOHPjgA/W3i0jrNZfcrwL2MLMrzWy5RkEzW8HMrgZ2B64sVIBSOB07hhXjXn8dHnkk6Wgk28SJ4VE1dxFprZzJ3d1fBH4D/Ar40MxuM7OR0XYb8CFhjfffuLsWdilTQ4ZAVRX84Q+qvZeaxpHySu4i0lrNrufu7leZ2b+Bs4CDgK7RoQXA08Al7v6vgkYoBbXCCnDWWXDyyfDPf8LuuycdkTSqrw9Tzq67btKRiEi5aXHJV3d/1t33I4yI7xVtq7r7ICX2dDj22JBE9tsvDLSrroYxY4pX/pgxocwkyi5ljYPpvpvvSUQknrjruePuS919VrQtKWRQUlz33APz5sHChaFpfupUqKkpTpIdMyaUNXVq8csuZUuXhuSuwXQi0hbNNstLZRg+PExxmmn+/LCwzOuvF7bsm24KZWWXPXx4GA9QqaZODV+41N8uIm2h5C5Mm9b0/i+/hBtuKGzZX33V9P5cMVUKTTsrIu2h5C707h1qitmqqmDKlMKWXV3ddNm9exe23FLXmNz75VyTUUQkt9h97pJeI0dCt27L7uvWLexPouwOHcKteZWsoQE22ghWXjnpSESkHCm5C0OGwKhRoaZuFh5HjSpOn3d22d27h8Fk77xT+LJLWX29muRFpO2U3AUISXbKlJBYp0wp7mC2zLI//RSOPz5MiTtuXPFiKCULFsD//qfBdCLSdkruUnKuvTasN3/UUWFZ2krz5pvhi45q7iLSVkruUnK6dIG77w5z3x98cO4R9WnV0BAeVXMXkbZScpeSVF0Nd9wRFk+pqamsee/r66Fr1zCgTkSkLZTcpWTttVcYNX/77XDddUlHUzz19bDllqHlQkSkLZTcpaSdfTbsvz+cfjq88ELS0RRH45zyIiJtpeQuJa1DB/j738OtcoceCjNnJh1RYX38McyapcF0ItI+Su5S8lZfHe69F+bMgcMPX34e/DTRYDoRyQcldykL3/teWGTm2WfD+vNppTnlRSQflNylbAwZAqeeCldcAXfemXQ0hdHQAL16Qc+eSUciIuVMyV3KymWXwQ47wC9+ESZ7SRtNOysi+aDkLmWlc2e4666woMpBB4VladNi8WL473/V3y4i7afkLmVn3XVDs/x778Gxx6Zngpt334WFC5XcRaT9lNylLO2yC1x6Kdx3X3hMAw2mE5F8UXKXsvXrX8Nhh4WJbp56Kulo2q+hIcxKt/nmSUciIuVOyV3Klhn89a/Qty8ccQR88EHSEbVPfT1sumlYOEdEpD2U3KWsrbxyaJr/+mv46U9Dn3W50rSzIpIvSu5S9vr2hdGj4ZVXQlN9OfryS5g8Wf3tIpIfSu6SCgcdBGeeCTfcALfemnQ0rTdxYnhUzV1E8kHJXVJjxAjYfXc4+WR4/fWko2kdzSkvIvmk5C6p0akT3HEH9OgBhxwCn32WdETx1dfDqqtC795JRyIiaaDkLqmy1lpwzz3w0UdhLvqlS5OOKJ6GhtDfbpZ0JCKSBkrukjo/+hH8+c8wfjxceGHS0bTMXXPKi0h+JZrczewcM3MzuzZjn5nZBWY23cwWmNnTZrZlknFK+ampCVPT/v738NBDSUfTvA8/hC++UH+7iORPYsndzAYAJwL1WYfOAH4DnApsC8wCHjezVYoboZQzM7j+evjBD+DII+H995OOKDdNOysi+ZZIcjez1YAxwPHAnIz9BvwauMTd73H3icAxwCrAz5OIVcpX166h/90MDj4Y5s9POqKmKbmLSL4lVXMfBdzt7tkzgm8I9AIea9zh7guAZ4EdiheepMWGG8KYMSGBnnRSaa4g19AAVVWw2mpJRyIiaVH05G5mJwIbA79r4nCv6PHjrP0fZxz71uzZs+nfv/+326hRo/IbrKTCvvvCBRfAP/4RJrkpNRpMJyL51qmYhZlZX+AiYGd3X9TMqdn1K2tiHz179mTChAl5jFDS6txzw/S0p50W+uEHDEg6omDhQpg0CQ44IOlIRCRNil1z3x7oAUw0s8VmthgYCAyNfv40Oi+7lr4Wy9fmRWLr0CHU3DfYICwwM2tW0hEFb78Nixer5i4i+VXs5H4/sBWwdcY2ARgb/fwOMBPYs/EFZtYF2Bl4ocixSsqssQbcey98+ikcfnhIqknTtLMiUghFTe7u/rm7T8zcgK+Az6LnDlwFnGVmB5tZP+BWYB5wezFjlXT6/vfhxhvh6afhnHOSjib0t3fuHNZxFxHJl6L2ucf0J6ArcB2wBvAysJe7z000KkmNo4+Gl1+GSy+F7bYL89Anpb4ettgizIsvIpIviU8/6+67uvspGc/d3S9w93XcvYu7D4xq+CJ5c+WVYVDdkUfCeuuFPvnq6nDbXLGMGQNPPAFvvFH8skUk3VRfkIrUuTP87Gfw0kswfXrYN3UqHH88vPsu7LNPYcsfPx4uvhiWLPmu7Jqa8POQIYUtW0TSz7wUZ/WIqX///q5b4aStqqtDUi0lVVUwZUrSUYhIKTOz19y9f3PnqOYuFWvatNzHHn64sGXvt1/T+5uLSUQkLiV3qVi9ezddc6+qCrPaFVJVVdNl9+5d2HJFpDIkPqBOJCkjR0K3bsvu69Yt7E9z2SKSfkruUrGGDIFRo0It2iw8jhpVnAFtSZYtIumnAXUiIiJlJM6AOtXcRUREUkbJXUREJGWU3EVERFJGyV1ERCRllNxFRERSRsldREQkZZTcRUREUkbJXUREJGWU3EVERFJGyb0EjRo1KukQEqH3XTkq8T1DZb7vSnzPkPz7VnIvQUn/o0iK3nflqMT3DJX5vivxPUPy71vJXUREJGXKeuEYM5sNNLEqdtnrAXySdBAJ0PuuHJX4nqEy33clvmco7PuucveezZ1Q1sldRERElqdmeRERkZRRchcREUkZJXcREZGUUXIvEWZ2tpm9amZfmtlsMxtnZv2SjquYzOwcM3MzuzbpWArNzNYxs9HR7/prM3vTzAYmHVchmVlHM/uDmU2O3vNkMxthZp2Sji1fzGwXM6szs4+if8vHZh03M7vAzKab2QIze9rMtkwo3Lxp7n2b2Qpm9kczqzezr8xshpndbma9Ewy53Vr6XWedOyo657fFik/JvXTsClwP7ADsDiwGnjCz7kkGVSxmNgA4EahPOpZCM7PVgecBAwYBmwOnArOSjKsIzgSGAb8CNgNOi56fnWRQebYyMJHw3hY0cfwM4DeE3/e2hN/542a2StEiLIzm3nc3YBtgZPR4ALABML7Mv9i19LsGwMx+SvhdTy9SXKFcjZYvTWa2MvAFcKC7j0s6nkIys9WAfxOS+3nARHc/JdmoCsfMLgIGuvuOScdSTGb2IPCpux+TsW80sKa7D04ussIws3nAKe5+a/TcCH/gr3X3kdG+roQE/1t3vzGpWPMp+33nOGcL4L/A99y9oVixFUqu92xmVcALwI+BRwi/+8uKEZNq7qVrFcLvZ07SgRTBKOBud38q6UCK5EDgZTOrNbNZZvaGmZ0S/fFPs+eA3cxsM/j2D/zuwMOJRlU8GwK9gMcad7j7AuBZQotdJVk1ekzt37eoVeIOYIS7v1Xs8su5SSTtrgbeAF5MOpBCMrMTgY2Bo5KOpYj6AEOBK4FLgK2Ba6JjaR5v8EfCl9Y3zWwJ4e/PSHe/PtmwiqZX9Phx1v6PgfWKHEtizKwzcDkwzt0/TDqeAvo9oaXqL0kUruRegszsCmAnYCd3X5J0PIViZn2Bi4Cd3X1R0vEUUQdggrs39jW/bmabEPqf05zcDweOBn5OaJLdGrjazCa7+18Tjay4svtCrYl9qRTVZm8DVgd+knA4BRMNjj2W8G88EWqWLzFmdiXwM2B3d38/6XgKbHvCFI0TzWyxmS0GBgJDo+crJhtewcwA3sza9xZQ1qOHY7gUuMzdx7p7g7v/A7iCdA2oa87M6LFX1v61WL42nzoZzdTfA/Zw908TDqmQdgPWAWZk/G2rAv5oZkVprVDNvYSY2dXAEcCu7v520vEUwf3AhKx9fwP+R6jRp7U2/zzQN2vfpqRznYRM3YDslqglVE4lYzIhwe8JvApgZl2AnYH/l2BcBWdmKwBjgX6Ev28zW3hJubseuDtr36OELzc3FSMAJfcSYWbXEfqdDwTmmFnjt/t57j4vucgKx90/Bz7P3GdmXwGfufvEZKIqiiuBF8xsOFAL/IBwe9g5iUZVeOOAs8xsMqFZ/gfA/wF/TzSqPIructk4etoB6G1mWxP+TU8zs6uA4Wb2NvAOcC4wD7g9kYDzpLn3TbhD4C7C7WD7A57x9+2LaFBh2Wnpd03Wra1m9g0w090nFSVAd9dWAhuhz62p7YKkYyvy5/A04XaRxGMp8PscBPwH+JrwR/5XRLempnUjDKa7itBCsQB4n9BC0yXp2PL4HnfN8f/41ui4ARcQuma+Bp4B+iUddyHfN1DdzN+3Y5OOvVC/6ybOn0K45bEo8ek+dxERkZSplL4uERGRiqHkLiIikjJK7iIiIimj5C4iIpIySu4iIiIpo+QuIiKSMkrukjpmdqyZuZl9bmZrZB3rFB27IIG4LojKLunJo8ysg5ldZWYzzGypmd3fzutVR+/72DyFmHnt1aPPdZt8X7vQMv49vBfN4JZ5bONCfWZSGZTcJc1WA85MOogy9FPgNMJc8DsCZyQbTrNWB84Hyi65Z+gDHJ90EJIuSu6SZo8Bp2ZMdZl6eVpsZ/Po8Sp3f9Hd38nDNUtCiS5G9BhwbjTPvEheKLlLmo2IHoc3d1Jj82gT+281sykZzxubl08ys4vNbKaZzTWz28ysW9SU+qiZzTOzd83smBxFbm5m/zSz+VHT94Vmtsz/RTPrYWZ/MbOPzGyhmb1tZjVZ5zR2P+xiZneZ2efAyy28133M7EUzW2BmX5jZ/dHSu43HpxCmRwVY0lLTcNTNcXYU30Izm25ml8dJVGY20MyejD7Dr6LPrl8T5x1kZs9Hn+uXZvaKmf3EzKoJi7EA3BTF+m28Zva0mT1nZvub2etmthAYGh1b1cyujeJdaGaTzOx0M7OMcneNrveT6NxPzGx29PtePSvG08zsrehznWNmE8zsoJY+g8hwYF3Ckr8ieaHkLmk2g7A+eo2ZVeXxumcT/hgfA5xHWKf8BuA+4CHgIKAe+JuZbdnE6+8HniAsEnQ78LvoOkBIPISV4wYREu0gwqIrfzGzU5u43hhCkvspcFauoM1snyi+eVHMJxNW6XrOzNaLTjuIMB84hCV5t49ek8tthMVPbo/ivJjQxDymmddgZoOAJ6NYjiSs8b4K8C8z2yDjvFOBewmLcBwDHEr4nKsJv9+Do1MvzhHvpsCfgWuAvYEnoy9SDwHHAZcTFjMZT1h+dmQT4V5NmDP858CFwCHRvsYYh0TXuQPYDxhCWBGse3OfQYY3CAurnGVmq8R8jUjzkp58X5u2fG/AsYQ/xhsT/sB+DtwSHetE1oI8hATqTVznVmBKxvPq6LVPZZ13b7T/yIx9awCLgfOzywHOynr9TcBcYPXo+e8Ii4ps0sR5nwCdst7nlTE/lwmE5XQ7ZezbEPgGuCJj34imPo8mrrdzVP7RWfuHRPu3zvrcjs04513gyazXrRq9v6syns8F7m0mhsZrn9DEsaeBpY1xZOwfnB1PtP9mYCHQI3q+a3Te6Kzzro1+P5bx/N9t+Hfa+O+hE+FLyGLgvOjYxk3FqE1b3E01d0k1d/+MUKs6OrP5uZ0eyXr+dvT4aEa5cwi1zQ1Y3p1Zz8cCKxNq0QD7EJrXJ0fN3p0sjLB/FFgT2CLr9fe1FLCZrUQYdFbr7osz4pxMaCUY2NI1mrAPsAi4JyvOx6Lju+SIZRNgI2BM1uvmAy9mvG4Hwucyqg2xNZri7m9k7duFkPTvyNp/G9CZUPvPlN1y0QCsCKwdPX8V2NrMrjGzH5tZt9YG6WFcw2jgN2YWt8YvkpOSu1SCKwnrSl+Yp+vNyXq+qJn9TfU9f5zjeWPT+FqEBPRN1nZXdHzNrNfPaDlk1iAsN9rUuTOJ34ScaS1CMpyXFWfjOtbZcWa+DuCvLP8eB2e8rvHxwzbE1qip99udsOb2wqz9MzOOZ/os63nj6xp/t38ndHFsR/gC9pmZ3RuNCWiN3xO+NOgOD2m3kr7fViQf3H2emV1MqMFf2sQpXwOYWWd3X5SxP1dyaq+1CWuZZz4H+Ch6/JSQIE/L8fpJWc/jrNs8JzqvqTsHekVlttanhM9u5xzHp36TJMEAAAKzSURBVDfzOghjF55o4njj7+CT6HE9YGIb4oOmP5vPgO5N/L4bP5tWfRbu7sCNwI0W5lXYi/BvrZaQ8ONeZ5qZ3QicAjzYmhhEsqnmLpXiekLyHNHEsanR47cjtaPR0DsUKJbDsp4fQaj9Niaw8cBmwDR3n9DENre1Bbr7V8BrwKFm1rFxfzTQcAfgmTa8j/GE2utqOeLMldwnAVOALXO8rj467wXC51KT4zrwXS26ayvifobwt+/QrP1DCF8sXmrFtZbh7nPcvZbQ9bLcyP8YRhK+kDR7h4dIS1Rzl4rg7gvN7EKa7r99BPiCcDvV+YSm0TMIiaUQToxGbL9KGMF9AmGA3+fR8SsJo9n/ZWZXEpLhSoSEv7O7H9DGcn9H6D9+0MyuJ/Rn/57w3i9v7cXc/WkzuwO428yuAF4h9GVXE0aNn+lN3CPv7m5mw4AHzKwzIRF+QmjB2IHwpeYKd59rZmcD15jZPYQR+HOBrYGv3f0aQpfGp8ARZlYPfAVMdvfmat+PAM8BN5hZT+C/UbwnABe7+yfNvHY5ZjYqiutFQovLpsBRfDf2IDZ3n2VmVwPntPa1IplUc5dK8jfCaPFlREl1MCEx3Um4reoa4J8FiuMAYE+gjnAb2AjgDxnxfEFIcg8T+l8fBW6JXtfmmNx9POF2tdUJ7/MG4C1gp2Zq2S05kjDq+6fAA4RbwE4hfM7ZYwsyY3mYMK5gJcIo9UeBPxGaxl/MOO9aQg17fUJyvycqa3J0fCkhKa9BaOJ/lXBrW07RawYRBrCdSfjCMwj4P9pWY34e+CGhdejx6Bq3EW7da4tLWX78hkirNN7KISIiIimhmruIiEjKKLmLiIikjJK7iIhIyii5i4iIpIySu4j8//bqQAYAAABgkL/1Pb6SCJiROwDMyB0AZuQOADNyB4CZAGu9vkr6UhnQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1)\n",
    "fig.set_size_inches(8, 6)\n",
    "ax.tick_params(axis=\"x\",direction=\"in\",which='major',labelsize=14)\n",
    "ax.tick_params(axis=\"y\",direction=\"in\",which='major',labelsize=14)\n",
    "ax.set_xlabel('Number of electrons N',fontsize=16)\n",
    "ax.set_ylabel('Chemical potential [meV]',fontsize=16)\n",
    "\n",
    "ax.plot(electrons2[1:],11.61*mu_N,color='b')\n",
    "ax.scatter(electrons2[1:],11.61*mu_N,color='b')\n",
    "\n",
    "plt.savefig('chemical.png')\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Different quantum dot (increased size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [],
   "source": [
    "lz = 0.5\n",
    "lx = 5\n",
    "\n",
    "OMEGA_X = 1/lx**2\n",
    "OMEGA_Z = 1/lz**2\n",
    "\n",
    "N_electrons = 2\n",
    "N_basis = 14\n",
    "N_points = 100000\n",
    "integrals_file = \"integrals/integrals_QD_14nz0_a.npy\"\n",
    "\n",
    "normalized_wf = True\n",
    "\n",
    "max_iter_SCF = 500\n",
    "eps_SCF = 1e-4\n",
    "Delta_SCF = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [],
   "source": [
    "def analytical_1(p,q):\n",
    "    n_p = index_to_q_numbers(p-1)\n",
    "    I = (p==q)*(OMEGA_X*(n_p[0] + n_p[1] + 1) + OMEGA_Z*(n_p[2] + 0.5))\n",
    "    \n",
    "    return I\n",
    "def f_cov(p, r, q, s):\n",
    "    \n",
    "    cov = 0.5*np.diag([1/OMEGA_X, 1/OMEGA_X, 1/OMEGA_Z, 1/OMEGA_X, 1/OMEGA_X, 1/OMEGA_Z])\n",
    "    \n",
    "    return cov\n",
    "def Hermite_pol(x, n):\n",
    "    \"\"\"\n",
    "    Evaluates the Hermite polynomial of order n at x.\n",
    "\n",
    "    Parameters\n",
    "    ==========\n",
    "    x : float or np.ndarray\n",
    "        Values in which to evaluate the Hermite polynomial\n",
    "    n : int\n",
    "        Oder of the Hermite polynomial\n",
    "\n",
    "    Returns\n",
    "    =======\n",
    "    float or np.ndarray\n",
    "    \"\"\"\n",
    "    return special.eval_hermite(n, x)\n",
    "\n",
    "def HO_wf(x, n, omega=OMEGA_X):\n",
    "    \"\"\"\n",
    "    Evaluates the n^th eigenfunction of the Harmonic Oscillator at x.\n",
    "    The units of x are sqrt(hbar/m*omega). \n",
    "\n",
    "    Parameters\n",
    "    ==========\n",
    "    x : float or np.ndarray(N)\n",
    "        Values in which to evaluate eigenfunction\n",
    "    n : int\n",
    "        Number of the eigenfunction\n",
    "    omega : float\n",
    "        Angular frequency of the Harmonic Oscillator\n",
    "\n",
    "    Returns\n",
    "    =======\n",
    "    float or np.ndarray(N)\n",
    "    \"\"\"\n",
    "    #return omega**0.25 * Hermite_pol(np.sqrt(omega)*x, n) * np.exp(-omega*x**2 / 2) / np.sqrt(2**n * np.math.factorial(n) * np.sqrt(np.pi))\n",
    "    return omega**0.25 * Hermite_pol(np.sqrt(omega)*x, n) / np.sqrt(2**n * np.math.factorial(n) * np.sqrt(np.pi))\n",
    "\n",
    "def HO_wf_3D(x, y, z, nx, ny, nz, omega_x=OMEGA_X, omega_y=OMEGA_X, omega_z=OMEGA_Z):\n",
    "    \"\"\"\n",
    "    Evaluates the n^th eigenfunction of the (anisotropic) Harmonic Oscillator at x,y,z.\n",
    "    The units of x,y,z are sqrt(hbar/m*omega). \n",
    "\n",
    "    Parameters\n",
    "    ==========\n",
    "    x, y, z : float or np.ndarray(N)\n",
    "        Position in which to evaluate the wave function\n",
    "    nx, ny, nz : int\n",
    "        Number of the eigenfunction for each cartesian coordinate\n",
    "    omega_x, omega_y, omega_z : float\n",
    "        Harmonic osciallation constant for each cartesian coordinate\n",
    "\n",
    "    Returns\n",
    "    =======\n",
    "    float or np.ndarray(N)\n",
    "    \"\"\"\n",
    "    return HO_wf(x, nx, omega=omega_x)*HO_wf(y, ny, omega=omega_y)*HO_wf(z, nz, omega=omega_z)\n",
    "\n",
    "def index_to_q_numbers(k):\n",
    "    \"\"\"\n",
    "    Returns the quantum numbers nx, ny, nz associated with the basis index k\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    k: int\n",
    "        Index of the basis from 0 to 13\n",
    "\n",
    "    Returns\n",
    "    ----------\n",
    "    nz ,ny, nz : int \n",
    "        Quantum numbers\n",
    "    \"\"\"\n",
    "\n",
    "    #q_numbers = np.array([(0,0,0),(0,0,1),(0,1,0),(1,0,0),(0,1,1),(1,0,1),(1,1,0),(1,1,1),(1,1,2),(1,2,1),(2,1,1),(1,2,2),(2,1,2),(2,2,1)])\n",
    "    q_numbers = np.array([(0,0,0),(0,1,0),(1,0,0),(1,1,0),(2,0,0),(0,2,0),(2,1,0),(1,2,0),(3,0,0),(0,3,0),(2,2,0),(2,2,0),(3,1,0),(1,3,0)])\n",
    "\n",
    "    return q_numbers[k]\n",
    "def integrand_2(R, p, r, q, s):\n",
    "    r1 = R[:,0:3]\n",
    "    r2 = R[:,3:6]\n",
    "    r12 = np.sqrt(np.sum((r1 - r2)**2, axis=-1))\n",
    "    \n",
    "    x1 = R[:,0]\n",
    "    y1 = R[:,1]\n",
    "    z1 = R[:,2]\n",
    "    x2 = R[:,3]\n",
    "    y2 = R[:,4]\n",
    "    z2 = R[:,5]\n",
    "    \n",
    "    n_p = index_to_q_numbers(p-1)\n",
    "    n_r = index_to_q_numbers(r-1)\n",
    "    n_q = index_to_q_numbers(q-1)\n",
    "    n_s = index_to_q_numbers(s-1)\n",
    "    \n",
    "    result = 1/r12\n",
    "    \n",
    "    result *= HO_wf_3D(x1, y1, z1, n_p[0], n_p[1], n_p[2])\n",
    "    result *= HO_wf_3D(x2, y2, z2, n_r[0], n_r[1], n_r[2])\n",
    "    result *= HO_wf_3D(x1, y1, z1, n_q[0], n_q[1], n_q[2])\n",
    "    result *= HO_wf_3D(x2, y2, z2, n_s[0], n_s[1], n_s[2])\n",
    "  \n",
    "    \n",
    "    return result\n",
    "\n",
    "def norm_product(p, r, q, s):\n",
    "    return (np.pi/OMEGA_X)**2 * np.pi/OMEGA_Z\n",
    "\n",
    "def two_body_integrand(p, r, q, s):\n",
    "    f = lambda R: integrand_2(R, p, r, q, s)*norm_product(p, r, q, s)\n",
    "    \n",
    "    return f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating two electron integrals: 1/14\n",
      "Calculating two electron integrals: 2/14\n",
      "Calculating two electron integrals: 3/14\n",
      "Calculating two electron integrals: 4/14\n",
      "Calculating two electron integrals: 5/14\n",
      "Calculating two electron integrals: 6/14\n",
      "Calculating two electron integrals: 7/14\n",
      "Calculating two electron integrals: 8/14\n",
      "Calculating two electron integrals: 9/14\n",
      "Calculating two electron integrals: 10/14\n",
      "Calculating two electron integrals: 11/14\n",
      "Calculating two electron integrals: 12/14\n",
      "Calculating two electron integrals: 13/14\n",
      "Calculating two electron integrals: 14/14\n"
     ]
    }
   ],
   "source": [
    "# One- and Two-body integrals\n",
    "MC_args = {\"f_cov\":f_cov, \"f_integrand\":two_body_integrand, \"N_points\":N_points}\n",
    "integrals = lib.integral_master(N_basis)\n",
    "integrals.calculate(integrals_file, analytical_1, MC_args = MC_args)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not normalized_wf:\n",
    "    S = np.zeros((4,4))\n",
    "    for p in range(4):\n",
    "        for q in range(4):\n",
    "            S[p][q] = (np.pi/(ALPHA[p]+ALPHA[q]))**(1.5)\n",
    "    SVAL, SVEC = np.linalg.eigh(S) \n",
    "    SVAL_minhalf = (np.diag(SVAL**(-0.5))) \n",
    "    X = np.dot(SVEC, np.dot(SVAL_minhalf, np.transpose(SVEC)))\n",
    "else:\n",
    "    S = np.eye(N_basis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing energy for N_electrons =  2 \n",
      "E(2)= 4.067343852093808\n",
      "Computing energy for N_electrons =  4 \n",
      "E(4)= 5.996329279611082\n",
      "Computing energy for N_electrons =  6 \n",
      "E(6)= 8.07675198739323\n",
      "Computing energy for N_electrons =  8 \n",
      "E(8)= 10.145576229440753\n",
      "Computing energy for N_electrons =  10 \n",
      "E(10)= 12.252479552241068\n",
      "Computing energy for N_electrons =  12 \n",
      "E(12)= 14.404555831868734\n",
      "Computing energy for N_electrons =  14 \n",
      "E(14)= 16.642144388662526\n",
      "[ 0.          0.          4.06734385  0.          5.99632928  0.\n",
      "  8.07675199  0.         10.14557623  0.         12.25247955  0.\n",
      " 14.40455583  0.         16.64214439]\n"
     ]
    }
   ],
   "source": [
    "N_max = 14\n",
    "N_av = 20 #times we do the SCF to do statistics\n",
    "N_pairs = int(N_max/2)\n",
    "\n",
    "E_i_a = np.zeros(N_av)\n",
    "E_N_a = np.zeros(N_max + 1)\n",
    "E_Nvar_a = np.zeros(N_max + 1)\n",
    "\n",
    "for n in range(1,N_pairs +1):\n",
    "    N_electrons = n*2\n",
    "    print ('Computing energy for N_electrons = ', N_electrons,'\\r')\n",
    "    for av in range(N_av):\n",
    "        E_i_a[av] = lib.SCF(N_electrons, integrals, S, max_iter_SCF, eps_SCF, print_E=False)\n",
    "        \n",
    "    E_N_a[N_electrons] = np.average(E_i_a)\n",
    "    E_Nvar_a[N_electrons] = np.std(E_i_a)/np.sqrt(N_av)\n",
    "    print('E({})='.format(N_electrons),E_N_a[N_electrons])\n",
    "    \n",
    "print(E_N_a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "energies [ 0.          2.03367193  4.06734385  5.03183657  5.99632928  7.03654063\n",
      "  8.07675199  9.11116411 10.14557623 11.19902789 12.25247955 13.32851769\n",
      " 14.40455583 15.52335011 16.64214439]\n",
      "chemical pot [2.03367193 2.03367193 0.96449271 0.96449271 1.04021135 1.04021135\n",
      " 1.03441212 1.03441212 1.05345166 1.05345166 1.07603814 1.07603814\n",
      " 1.11879428 1.11879428]\n"
     ]
    }
   ],
   "source": [
    "for k in range(int(Ne_max/2)+1):\n",
    "    if k != Ne_max/2:\n",
    "        E_N_a[2*k+1] = (E_N_a[2*k+2] + E_N_a[2*k]) / 2\n",
    "\n",
    "mu_N_a = np.zeros(Ne_max)\n",
    "for k in range(1, Ne_max + 1):\n",
    "    mu_N_a[k-1] = E_N_a[k] - E_N_a[k-1]\n",
    "    \n",
    "print('energies',E_N_a)\n",
    "print('chemical pot', mu_N_a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAF7CAYAAADc0IJwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4VGX2wPHvoRMpUaSJhoCIAhaU8BMQlBVBV2V1XRQlSBOD7EpTbGDBgm0VBEWaCgrBgrt2RVBAXBQlkcUFwUoRqYKAQEBIzu+PdyaZhCRMwszcmcz5PM88mbn3zp2TEHLufct5RVUxxhhjTPwo53UAxhhjjIksS/7GGGNMnLHkb4wxxsQZS/7GGGNMnLHkb4wxxsQZS/7GGGNMnLHkb4wxxsQZS/7GGGNMnKngdQDhdvzxx2tycrLXYRhjjDERkZmZ+auq1i7umDKf/JOTk8nIyPA6DGOMMSYiRGTdkY6xZn9jjDEmzljyN8YYY+KMJX9jjDEmzljyN8YYY+KMJX9jjDEmzljyN8YYY+KMJX9jjDEmzljyN8YYY+KMJX9jjDEmzljyN8YYY+KMJX9jjDEmzljyN8aYUlIt/rUx0cqSvzHGlMKoUTBsWF7CV3WvR43yMipjgmPJ3xhjSkgVdu6EcePyLgCGDXOvd+60FgAT/cr8kr7GGBNqIjB2rHs+bpx7AAwZ4raLeBebMcEQLeOXqCkpKZqRkVHi9+3evZutW7dy8ODBMERlAlWsWJE6depQo0YNr0MxpkRUoVxA+2lOjiV+4z0RyVTVlOKOsTv/QuzevZstW7bQoEEDqlatitj/5rBRVbKysvjll18A7ALAxAxVSE3Nv23YMLvzN7HB+vwLsXXrVho0aEBCQoIl/jATERISEmjQoAFbt271OhxjgqIKQ4fCyy9D1arw22+uyT9wDIAx0czu/Atx8OBBqlat6nUYcaVq1arWxWJihgisW+eeT5gAiYl5YwASE+3O30Q/S/5FsDv+yLKft4klu3bBkiXQpg307u22+QcB2q+yiQWW/I0xpoTuvx+2boX33ss/4M8Sv4kV1udfhnXs2JFu3boVui8lJYU+ffoA0KdPH1JSCh8YWvAc06dPR0QOezRp0iTk8RsTjVauhPHjIS0NWrXyOhpjSsfu/E2pzJ8/P9+4iCpVqngYjTGRoQqDBkHNmjB6tNfRGFN6lvxNqbRu3Zpq1ap5HYYxETV7NixYABMnQq1aXkdjTOlZs78Ji6lTp3LGGWdQpUoV6tatS7du3di1a5fXYRlTanv2wK23wtlnw403eh2NMUfHkn+YxdKqX4cOHTrsUVQFyOzs7CKPe+ihhxgwYAAXXHABb775JhMnTqRmzZrs2bMnUt+KMSH38MOwYQM8/TSUL+91NMYcHWv2D6NRo9wiH/7pP/7FPxITo2/lr8zMTCpWrFjovr/97W+HbUtMTMz3eurUqfTv35+dO3fy8MMPM3ToUMaMGZO7/6qrrgptwMZE0Pffw5NPQq9ecN55XkdjzNGz5B8mgat+gbsA8K/6NWSI2x9N04KaNWvGSy+9dNj2AQMGFHr8okWL8g34S05OBuDzzz8nKyuLvn37hiVOYyJN1f2frVwZHnvM62iMCQ1L/mESDat+VahQgezs7EL3ZWdnU6FC3j9/QkJCodP9qlevXuj7zz777EIH/G3fvh2A+vXrlyZkY6LOu+/CBx/AmDFQr57X0RgTGtbnH0aBFwB+kawAVrt2bTZv3lzovk2bNlGnTp2Qf2Yt3xDoTZs2hfzcxkTa/v3ugr15c7j5Zq+jMSZ0LPmHkb+PP1AkF/3o0KEDmZmZuSvm+X3xxRds2bKFDh06hPwz27ZtS9WqVXnxxRdDfm5jIu2f/4Q1a9wgvyKGxBgTkyz5h4k/8fv7+HNyIr/qV69evUhKSuL8889n2rRpzJ8/nwkTJtC1a1fatWvHxRdfHPLPTExM5J577mHMmDEMHjyYOXPm8NZbb5GWlnbYRYgx0WztWjfC/+qr4cILvY7GmNCyPv8wEXGj+gP7+CO96le1atVYtGgRI0aM4M4772THjh3UrVuX7t27M3r0aMqVC8+131133cVxxx3HuHHjmDx5Msceeyznn39+keMHjIlGt97q6vY/8YTXkRgTelLUPO6yIiUlRTMyMkr0nlWrVtGsWbOQfH7BUf3RNso/moTy527M0Zg3D7p0cSV8R4zwOhpjSkZEMlW18AVbfKzZP8wKJnpL/MZEtz/+gMGDoUkTd/dvTFlkzf7GGBNg/HhYvdot11u5stfRGBMedudvjDE+GzfC/fdD165w6aVeR2NM+FjyN8YYn9tvh4MHD6/PYUxZY8nfGGOARYsgPR1uuw1OPtnraIwJL0v+xpi4d+gQDBoESUlw111eR2NM+NmAP2NM3Js8Gb7+Gl5/HRISvI7GmPCzO39jTFzbtg3uvhsuughs5WkTLyz5G2Pi2ogRsGePm+JndThMvIho8heR8iLyoIisEZH9vq8PiUiFgGNEREaJyEYRyRKRhSLSosB5jhWRGSKyy/eYISKJkfxeYkHHjh3p1q1boftSUlLo06cPAH369Cl0Od/CzjF9+nRE5LBHkyZNQh6/MeH25Zfw/PMwdChYcUkTTyLd538H8A+gN/A/4EzgReAA8KDvmNuBW4E+wLfAvcA8ETlVVX/3HTMLSAL+DCjwHDAD6BqR78Iwf/58qlatmvu6SpUqHkZjTMnl5LhleuvWhXvu8ToaYyIr0sm/HfCOqr7je71WRN4GzgV31w8MBR5V1X/5tvUGtgI9gMki0gy4BGivqp/5jhkAfOq7QPg2ot9RnGrdujXVqlXzOgxjSm3aNFi6FGbMgBo1vI7GmMiKdJ//f4A/ichpACLSHLgQeN+3vxFQD5jrf4OqZgGLcBcOAG2BPcBnAeddDOwNOMZ4aNq0aZx33nkcd9xxHHfccXTq1ImvvvrK67CMyfXbb3DnndC+PaSmeh2NMZEX6Tv/x4DqwDciku37/NGq+qxvfz3f1y0F3rcFaBBwzDYNWI5QVVVEtga8P3rE0LJ+hw4dOmxbUas+Zmdn5zu+fPnyiO/7WrduHX369KFx48b88ccfzJw5kw4dOvDNN9/QsGHD8ARvTAncdx/s2AFPPx21/x2NCatIJ//uQC9cE/5KoCUwTkTWqOrzAccVzDhSYFthGangMQBs27Yt32C2tLQ00tLSShd9SY0aBTt3ulqhIi7xDxsGiYluXxTJzMykYsWKhe7729/+dti2xMT84yunTp1K//79ARgV8L3l5OTQuXNnTjvtNNLT0xlh66Maj339NUyYAAMHQsuWXkdjjDcinfz/CTyhqq/4Xv9PRBoCdwHPA5t92+sBPwe8rw55rQGbgToiIv67f99Ygdoc3mJA7dq1ycjICPk3ckSqLvGPG+dejx3rEv+4cTBkSNS1ADRr1oyXXnrpsO0DBgwo9PhFixblG/CXnJyc+3zlypWMHDmSzz//nK1bt+Zu/+6770IXsDGloOoG+R13HDzwgNfRGOOdSCf/BCC7wLZs8sYerMEl987AUgARqQJ0AG7zHfM5UA3X9+/v928LHEP+cQDeEslbHWTcuLyLgCFD8loCwqxChQpkZxf8cTvZ2dlUqJD3z5+QkFDodL/q1asX+v6zzz670AF/u3btokuXLjRo0ICxY8eSlJRElSpV6Nu3L/v37y/ld2JMaLz8Mnz6KUyd6i4AjIlXkU7+7wB3isgaXLP/2cAtwEuQ23f/FDBSRFYD3wF34wb4zfIds0pE5uBG/t+Ia+6fDLwbdSP9/RcA/sQPEUv84Fo91q5dW+i+TZs2UadOnZB/5uLFi9m4cSOffPJJvrn/O3fuDPlnGVMSv/8Ow4dDSgr06+d1NMZ4K9Kj/QcBrwPPAquAJ4GpwMiAYx4HxgATgAygPtAlYI4/QCqwHDcr4EPf8+vDHXyJ+fv4Aw0b5rZHQIcOHcjMzOSXX37Jt/2LL75gy5YtdOjQIeSfmZWVBUDlypVzty1atIgNGzaE/LOMKYkHH4RNm+CZZ6Cc1TY1cS6i/wVU9XdVHaqqDVW1qqo2VtURqro/4BhV1VGqWl9Vq6jqBaq6osB5dqhqT1Wt4Xv0VNXourX0J35/H39Ojvs6blzELgB69epFUlIS559/PtOmTWP+/PlMmDCBrl270q5dOy6++OKQf2a7du1ISEigf//+zJs3j+eff56ePXtywgknhPyzjAnW6tWu0a1fPzj3XK+jMcZ7tqpfuIi4Uf2Bffz+MQCJiRFp+q9WrRqLFi1ixIgR3HnnnezYsYO6devSvXt3Ro8eTbkw3P7Ur1+f2bNnM3z4cLp27cqpp57KlClTeOihh0L+WcYEQxUGD4ZjjoFHHvE6GmOigxQ1j7usSElJ0ZKO9l+1ahXNQlXoO4bm+XstpD93Y3zeeMOt1jd+PAwa5HU0xoSfiGSqauELtvhYz1e4FUz0lviNiZh9+1wv2xlnuHn9xhjHmv2NMWXWY4/BunXwySdQwf7aGZPL7vyNMWXSTz+55H/ddXD++V5HY0x0seRvjCmThg1zd/v//KfXkRgTfawhzBhTZvjH037wAbz9Njz6KDRocOT3GRNv7M7fGFMmjBrl7vb373czbJs2hQ0bom4NLWOigiV/Y0zMC1xHq1Mn+P57aN7cVfPbuTNiRTWNKZH0dKhSxbVWJSe715Fizf7GmJjnr6G1ezdMm+a2vflmRNfRMqZE0tMhLQ0OHHCv161zrwFSU8P/+Xbnb4wpE1Sh4DpWlvhNtBo50tWhCLRvn9seCUXe+YvI3FKe8++q+kMp31smdezovi5c6GUUxpRtTz8NCxbk3zZsmF0AmOi0fn3JtodacXf+F+FW1KtegkcnoEYY4zUlMGrUKEQk91GvXj0uv/xyvv7669xj1q5di4jw7rvvHvb+hQsXIiKsWJG3rlJycnK+c/ofM2fOjMj3ZExhVq+GW291zwcP9mQdLWNKJCmpZNtD7Uh9/jeo6pfBnEhEKgB/HH1IJpRq1qzJnDlzAJfo7733Xjp37syqVas47rjjSnXOHj16MKhAkfQmTZocdazGlMahQ9C7N1Sq5PpKn3rKk3W0jCmR0aNdH39g039CgtseCcUl/3RgWwnOle17z/ajisiEVIUKFWjTpg0Abdq0ITk5mbZt2zJnzhx69OhRqnPWr18/95zGeO2xx+DLL+HVV+Hqq/MSvf8CwBK/iUb+QX0jR7qm/qQkl/gjMdgPimn2V9XrVXVNsCdS53pVXRea0MqG9HRYssTVFo/0VI7CnHXWWQD8/PPPYfuMdevWcd1113H88ceTkJDAmWeeyaxZs8L2eSZ+LVvm5vFfey1cc42to2ViS2qqG6Sak+O+RirxQ/ED/m4HZqjqpsiFU7Z4PZWjMOt9o0kaNWqUb3tOTg6HDh3Kty07O7vQc6hqvmNFhPLlywOwdetW2rZtS0JCAk888QQnnXQSK1asCOvFholPBw5Ar15QuzZMmOB1NMbEluKa/R8FRovIR8BLwBuquj8yYZUNxU3liGTy9yfqdevWcfPNN9OyZUuuuOKKfMcUfF2cMWPGMGbMmNzXDRo0YMOGDQCMHTuWXbt2kZmZSf369QHo1KnT0X4Lxhzm3nthxQp47z0o5fAVY+JWccn/dKA30APXl/+7iLwGvKSqn0YiuFjn9VQOgO3bt1OxYsXc17Vq1WLp0qVUrlw533Fjx46lffv2+bZlZmZy0003HXbOnj17MmTIkNzXlSpVyn0+f/58LrnkktzEb0w4/Oc/bsGetDS49FKvozEm9hSZ/FX1G+AOEbkTN+3veuBaoJ+IrANexHUL/BSRSGNQUpJr6i9se6TUrFmTjz76iOzsbJYvX87w4cPp0aMHixcvply5vCEfTZo0ISUlJd979+zZU+g569ate9ixftu3b6d169ah+waMKWDPHje6PzkZnnjC62iMiU1HrPDnG8g3T1V7AXWBvsBPwD3A9yLyqYj0D3OcMWn0aDd1I1Akp3KAG+2fkpLCueeeS1paGhMmTGDJkiXMnj07LJ9Xq1YtNm2yYSImfG67DdasgRdfhOrVvY7GmNhUovK+qrpPVV9S1YuAJGA80A6YHI7gYl1qKkyZAv4W9oYN3WuvBvuBa7Jv0aIFjz32WFjO36lTJz788EO2bNkSlvOb+PbhhzBpEtxyC3To4HU0xsSuEtf2F5EqItIDeB74ByDAZ6EOrKxITYU2beCCCyI/laMwIsKIESNYtmwZH3/8ccjPP2zYMBITE+nQoQMvvvgi8+fPZ8KECTz++OMh/ywTX377Dfr1c6v1PfSQ19EYE9uCTv4i8icRmQZsAWYCpwKPAE1U1a7BY0j37t055ZRTwpKQa9euzeLFizn77LMZOnQol19+OVOmTCEpkgMdTJk0aBBs3QovveSWQTXGlJ5oMUWvRaQZ0AtIBRoAe4HXgRdV9ZOIRHiUUlJSNCMjo0TvWbVqFc2aNQtZDLawT3BC/XM3Zcfrr7vqffff76b4GWOKJiKZqlr4qGyf4or8ZABnAwosAEYA/1LVrJBGaYwxxdi8GW66CVJS4K67vI7GmLKhuHn+xwB346bzbYhQPGWS3fEbUzqqbi7/nj2uuT+gZIUx5igUN8/f2l+NMZ6aPh3eeQfGjAHrETImdEoy4K+FiMwSkQ0isk9EzvZtv19ELgpfiMaYeLR2LQwZ4mbKBBSUNMaEQFDJX0TaAEuB/wPeByrjpvgBVAIGhiU6Y0xcysmBvn3d8+nToVyJJyUbY4oT7H+px4CFuOl9fycv8QNkAMWOKjTGmJJ4+mk3VmbsWFfG1xgTWsEm/9bA06qajRv9H2gbUCekURlj4tbq1XDnnXDZZa6ojzGRkJ7uLjTLlXNf09O9jii8ihvtH+gAUFRZjXrA7tCEY4yJZ4cOQa9ebg2MqVNB5MjvMeZopae7WSX+JdjXrXOvwfuqrOES7J3/YuBmkXz/Ff0tAH1wXQKmKB075lX6McYU6ZFHYOlSV7/fVoU2kTJyZF7i99u3z20vq4JN/vcBbYFM4DZc4r9OROYA5wMPhCc8czRGjRqFiOQ+6tWrx+WXX87XX3+de8zatWsREd59993D3r9w4UJEhBUrVuRuS05OzndO/2PmzJkR+Z5M2fXVV/DAA3Ddda6anzGRsn59ybaXBUE1+6tqpohcCDwJjMYN+BsOfAF0UtWV4QvRHI2aNWsyZ84cwCX6e++9l86dO7Nq1SqOO+64Up2zR48eDBo0KN+2Jk2aHHWsJn7t3++a+2vXhmee8ToaE2+SklxTf2Hby6pg+/xR1SXAeSJSHagN7FDVnWGLzIREhQoVaNOmDQBt2rQhOTmZtm3bMmfOHHr06FGqc9avXz/3nMaEwr33wsqV8P77UMprUmNKbfTo/H3+4MadjB7tXUzhVuLZs6r6u6r+ZIk/Np111lkA/Pzzz2E5/y+//ELfvn1p1KgRVatWpWnTptx3330cPHgwLJ9nYt+nn8ITT8CAAfDnP3sdjYlHqakwZQo0bOgGmTZs6F6X1cF+UII7fxE5GbgKOInDR/6rqg4IZWBlRno6LFkCBw64+SOjR3v6G7Xe14nVqFGjfNtzcnI4dOhQvm3Z2dmFnkNV8x0rIpQvXx6Abdu2cfzxx/PUU0+RmJjI6tWruf/++/n111+ZMGFCKL8VUwbs2QN9+kCjRu4CwBivpKaW7WRfUFDJX0S6AzOB8sBvwB8FDlHAkn9B/vkjBw641x7NH/En6nXr1nHzzTfTsmVLrrjiinzHFHxdnDFjxjBmzJjc1w0aNGDDBrf2U8uWLWnZsmXuvvPOO4+qVaty0003MW7cOCpUCPp608SB4cNhzRr45BOoVs3raIyJH8H+JX4QmAP0U9VtYYynbClu/kiEkv/27dupGLAUWq1atVi6dCmVK1fOd9zYsWNp3759vm2ZmZncdNNNh52zZ8+eDAkotl6pUqXc5zk5OYwdO5bnnnuOtWvXsn///tx9GzZsINnKtRmfOXNg8mR3AdChg9fRGBNfgk3+9YG/W+IvoSiYP1KzZk0++ugjsrOzWb58OcOHD6dHjx4sXryYcgEF05s0aUJKSv4qzXv27Cn0nHXr1j3sWL8nn3ySu+66ixEjRtChQwcSExNZsmQJgwcPznchYOLbjh1www3QogU8+KDX0RgTf4JN/kuAJsBHYYyl7ImC+SMVKlTITdTnnnsuVatWpVevXsyePZvu3buH/PNmz57NtddeywMP5JV+CKwrYAzAzTfD1q1uud4qRdUONcaETbCj/QfhKvz9VUSOCWdAZcro0W6+SCCP54/07NmTFi1a8Nhjj4Xl/FlZWYd1KaSX9SLZpkRmz4aXX3bT+845x+tojIlPwSb/n3AFfV4HdovIHwUeB8IXYgzzzx/xJ8MomD8iIowYMYJly5bx8ccfh/z8nTt3ZtasWUycOJEPP/yQ1NRU1q5dG/LPMbFp82YYOBBat4a77vI6GmPiV7DN/pOA63GD/lZz+Gh/U5TUVLdCCbg1SqNA9+7dGTVqFI8//jiTJ08O6bnvv/9+tm/fzogRIxARunXrxtixY7nyyitD+jkm9qjCjTfC3r3w0ktgEz+M8Y6oFlyht5CDRHYBD6jqk+EPKbRSUlI0IyOjRO9ZtWoVzZo1C10Q/kV9oiT5R6uQ/9xNVHn+eejfH8aOhaFDvY7GmLJLRDJVtfBR2T7BNvvvB5YffUjGmHi0dq1L+B07wuDBXkdjjAm24e1F4BpstH/p2B2/iUOqrlRqTo6r4icC06dDuRIXFTfGhFqwyX8VcK+IvIPr9/+t4AGqOiuUgRljYteoUbBzp2viHz/eVfC76CKYNs3tM8Z4K9jk/7zva0PgskL2K2DJ3xiDqkv848a5Yj6vveZq93/0kSvq428RMMZ4J9jkb6OwjDFBEXF3/IcOgX8tpzVrYMgQt90SvzHeCyr5q+q34Q4k2qgqYn+lIiaYWScmthRc1sISvzHRw4beFKJixYpkZWV5HUZcycrKyrcAkYltzzzj+vcDDRvmmvyNMd4rMvmLyDcickawJxKRcr73FNtFICL1ReRFEdkmIvt977kgYL+IyCgR2SgiWSKyUERaFDjHsSIyQ0R2+R4zRCQx2FiPpE6dOvzyyy/s27fP7kjDTFXZt28fv/zyC3Xq1PE6HBMCCxa4Jn5w0/pyctzrcePsAsAELz0dkpPd7JDkZPfahE5xzf6nASVZckN876la5AEuQS8G/oMbOLgNaAxsDTjsduBWoA/wLXAvME9ETlXV333HzAKSgD/jBhs+B8wAupYg3iLVqFEDgI0bN3Lw4MFQnNIUo2LFitStWzf3525i17p1cM01UKsWdOsGTz2VNwYAIDHRmv7NkaWnQ1paXtfRunXuNXhaHb1MKbLCn4jkAAdwyTVYlYHWqvpVEed8GLhAVc8rYr8AG4FnVHW0b1tV3MXBcFWd7GtZ+AZor6qLfce0Bz4FTis4PqE0Ff6MMSW3bx+cd54b3Pfll3DKKfkTvY3yN8FKTi58QdSGDV3BKFO8YCr8FXfnX9pl3zYVs+9KYI6IvAr8CZfonwMmqLsKaQTUA+b636CqWSKyCGgHTAbaAnuAzwLOuxjY6zsm7gYnGuM1VejXD5Yvh3ffhaZNDz/GEr8J1vr1JdtuSq7I5K+q4VhzqzHwd2As8CjQEnjat+8ZXOIH2FLgfVuABr7n9YBtGtBkoaoqIlsD3p9r27ZtuevZA6SlpZHmbz8Kka5du7JpU3HXPPGjfv36vPPOO16HYSLs8cfh1VfhkUfg0ku9jsbEuqSkwu/8k5IiH0tZFel1tcoBGQEXFstE5BTgH7jk71ewq0EKbCusK6LgMQDUrl2bcDf7b9q0KeyfESsCL7RMfPjgA7c8b/fucMcdXkdjyoLRo/P3+QMkJLjtJjQiPdVvE66/PtAq3OA9gM2+rwXv4OuQ1xqwGagjAZPwfc9rc3iLgTEmjL77Dq67Ds48063aZ037JhRSU2HKFNfHL+K+Tplig/1CKdLJfzFwaoFtTQF/A88aXHLv7N8pIlWADuT18X8OVMP1/fu1BY4h/zgAY0wY7d4NV14JFSvCm2/CMcd4HZEpS1JT3eC+nBz31RJ/aEW62X8s8JmIjAReBc4GBgMjILfv/ilgpIisBr4D7sYN8JvlO2aViMwBJovIjbjm/snAu/FYidAYL+TkQM+e7s7/o4/c6GxjTOyIaPJX1aUiciXwMHAPsN739dmAwx7H1QqYABwLfAF0CZjjD5AKjCdvVsDbwM3hjd4Y4zdqFLzzDjz9NHTs6HU0xpiSivSdP6r6HvBeMfsVGOV7FHXMDqBnqGMzxhzZv/8NDz7opvb94x9eR2OMKY2gkr+IlAdaA22AE3B35r/i5tQvUtWNYYvQsGHDBk466STWrFlDsrWvGg/973/Qqxe0aQPPPmsD/IyJVcUO+BORJBF5AjdKfzEwBjct7zpcc/0s4GcRWSQi1waOwDfBERH+85//HPV5Jk2aRNOmTalWrRpnn302CxcuPOxzEhISqFatWu5j165dufvvvvtuatWqRdOmTVm8eHHu9oMHD3LOOeeQmZl51DGa2LZjhxvgV6MG/OtfULmy1xEZY0qruIV9nsANuOuMG6h3PlBNVY9R1eNVtQKuIl8PYC0wBfiviLQKe9Qmn9mzZ3PPPffw2muvsWvXLgYMGMBll13G+gLlsObOncuePXtyHzVr1gTgq6++4pVXXuHHH3/k0UcfZeDAgbnvefjhh+nSpQutWtk/azw7dMjN49+wAd54A044weuIjDFHo7g7/9OBjqp6lqo+oqqLVTXfCt2quk5VX1XVXkB9IB04N4zxxqTx48fTqFEjqlevToMGDRgxYgQAZ511FgBdunShWrVq9O/fH4DNmzfzl7/8hZo1a9K0aVPmzJlT7Plnz55Nz549admyJeXLl+emm26iTp06TJ8+Paj4fvjhB1q3bk1iYiKXXnopP/zwAwArVqzg9ddfZ9SoUaX7xk2ZcccdblT/pElwrv0PNybmFVfe95KSnEhV9+JG6psA3333HXfeeSdLly6lRYsW7Ny5k9WrVwOwfPlyRIS5c+fSvn373PekpqZSo0b5DJnaAAAgAElEQVQN1q9fT1ZWFt26dSv2M3Jycg5belhV+e9//5tv29VXX83Bgwc5+eSTueOOO7jqqqsAOP300xkxYgQ7duxgwYIFnHnmmWRnZ9O/f38mTpxIlSolWdzRlDUzZsCYMTBoEPTt63U0xphQiHSRn7hToUIFVJWVK1eyZ88eEhMTadOmTZHH//LLL8yfP58nnniCmjVrUq9ePe67775iP6Nr167MnDmTjIwMDh48yDPPPMP69evZvXt37jEfffQRa9asYcOGDdxyyy2kpqbmtig0b96c4cOH06lTJyZNmsQLL7zA2LFjad26NY0bN6Zbt26cf/75jLbamnFn6VK48Ub405/gySe9jsYYEyrF9fnXKckjkkHHksaNG5Oens7UqVM54YQTaN++PXPnzi3y+A0bNgDQsGHD3G2NGjUq9jN69erFbbfdRmpqKvXq1eOrr76iU6dOHH/88bnHdOrUiSpVqlClShW6d+9Oz549SU9Pz91/0003sWzZMubNm0elSpV44YUXeOSRR7j11lu57LLLWLBgAR9++CHz5s0r7Y/CxJjNm+Gvf4V69eC111wlP2NM2VDcnf9m3Cj/YB421a8YV111FfPmzePXX3/lmmuu4YorrmCfb8WKghMkGjRwixeuC1jSas2aNcWeX0S44447+Pbbb9m+fTuTJk1i1apVdCym+kq5cuUO6yoA113Qv39/nnrqKapVq8by5ctp164d5cuXp127dixbtizYb9vEsD/+gG7d3Aj/N9+EgOtIY0wZUFzyHxzEYy6uvK5N8SvCt99+y5w5c9i3bx8VK1akZs2aiAjlyrkffb169fj+++9zjz/xxBPp2LEjt99+O7t372bLli08+OCDxX7Grl27WLVqFarKtm3bGDhwIDVq1KB3796AG7j35Zdf8scff3Dw4EHefPNNZsyYwTXXXHPYuSZNmkTjxo3p0qULACeffDIffPABWVlZLFiwgCZNmoTqR2Oi2KBBsHgxTJsGLVt6HY0xJuRUtcQP4ALgP0AOsBRXfrdU5wr3o1WrVhpuxX3G119/rW3atNEaNWpojRo19JxzztH3338/d/8LL7ygJ554oiYmJmpaWpqqqm7cuFEvu+wyrVGjhp5yyik6depUBXTNmjWFfsb69eu1RYsWeswxx+ixxx6rPXv21M2bN+funz9/vjZv3lwTEhI0MTFRW7VqpS+//HKh52natKn+9ttvudtWrlypLVu21Jo1a2q/fv00Ozu71D8LExsmTlQF1Tvv9DoSEy1mzlRt2FBVxH2dOdPriExxgAw9Qm4ULaTptygikoKry98JWA3co6r/DvH1SEilpKRoRkZGuD+DcH9GrLCfRWz79FO48ELo3NnV7i9f3uuIjNfS0yEtDfYFTPROSLAldqOZiGSqakpxxwQ12l9EWojIv3GL7DQB+gKnR3viN8YE7+efXT9/48Ywa5YlfuOMHJk/8YN7PXKkN/GY0DhSed/GIjIDWI6r6z8IOFVVX9KSNBkYY6JaVpYr3bt/P7z1FiQmeh2RiRYFCoUecbuJDcVN9ZsErAIuAe4CGqvqs6p6MFLBGWPCT9XN5V+2zDXxnnaa1xGZaJKUVLLtJjYUt6pfGqDAT7j6/p2LWbdHVfXiEMdmjImAMWNc0n/oIbj8cq+jMdFm9OjC+/yt5ldsKy75f4lL/gDVIxCLMSbC5s6F2293ff2+JSeMycc/qG/kSNfUn5TkEr8N9ottxdX2L7oGrTEm5v34I1x7LbRo4ebz24LcpiipqZbsyxqr7W9MHPr9d7jiCpfw33oLqlXzOiJjTCQVN+CvVAU9S/s+Y0xk5ORA796werWr2X+EpSOMMWVQcXf+a0TkMRFpfKSTiEhFEblKRL4ABoYuPGNMKAROzH3oIXjjDXjiCejUybuYjDHeKW7A38XAP4HhIvIl8Cluvv824ABwLNAY+D+gC1ARGAuMCWfA0ah+/fqkpBRbTClu1K9f3+sQTAGjRsHOnTB2LLz9Ntx3HzRrBr/95nVkxhivFDfg7zPgPBE5D7gBSAWG+3fjFvM5BGQCDwAvqmpc/jl55513vA7BmEKpusQ/bhxs3+5W6KtbF1atgl273H4b6GdM/Cnuzh8AVV0MLAYQkSTgBKAKsB34UVX3FfN2Y4yHRNwd/x9/wMSJbtuePTBkiNtuid+Y+HTE5B9IVdcDVtTRmBjz++/5X1vij20dO7qvCxd6GYWJZcWN9m8qIpWOdAIRqSkil4Y2LGNMqEyfDjNn5t82bFj+QYDGmPhS3Gj/VUBL/wsRKSciW0XkjALHnQZYp7cxUeibb1xpVoBBg9w0vyFD3BgAuwAwJn4V1+xfsFFQgONxo/qNMVEuKwu6d4dKlaBPH5fw/WMAwK3cZ03/xsSnEvX5G2Nix9ChsGIFfPghdO6cl+j9FwCW+I2JX1be15gy6JVXYMoUuPNO6NLl8ERviT92pafDkiXwySeQnOxeG1NSlvyNKWN++MH187drBw884HU0JpTS092/7YED7vW6de61XQCYkjpSs38XEWnie14OV9znEhE5LeCYk8MSmTGmxA4ccP38FSrAyy9DRRuhU6aMHAn7ClRW2bfPbbdV90xJHCn5F3bf8FAh22zMsDFR4I474Kuv3Ep9SUleR2NCbX0RVVaK2m5MUYpL/s0iFoUx5qi9+aYb0T90KPzlL15HY8IhKck19Re23ZiSKK62/7eRDMQYU3rr1kHfvtCqFTz6qNfRmHAZPdr18Qc2/SckuO3GlIQN+DMmxh08CNddB9nZ8OqrULmy1xGZcElNdbM4/P/GDRu619bfb0rK5vkbE+PuuQc+/9wl/pNt+G2Zl5oKU6e651bb35SW3fkbE8M+/BAeewwGDIBrrvE6GmNMrLDkb0yM2rgRrr8ezjgjr2SvMcYEw5K/MTEoO9s1/+7dC6+9BlWreh2RMSaWBJX8ReR6EbFhRMZEiYcecv29zz4Lp512xMONMSafYO/8XwR+EZEnReTUcAZkjCneggVw//3Qqxf07u11NMYLCxfaYD9zdIJN/s2BGUBv4BsRWSAi3UXEiocaE0Fbt7rm/qZNYcIEr6MxxsSqoJK/qq5W1WFAA6APUBF4GdggIo+ISOPwhWiMAcjJcXf6O3a4fv5q1byOyBgTq0o04E9VD6jqDFVtD5wBrAJuB74TkTki0ikcQRpj4IknYM4ceOopOPNMr6MxxsSyEo/2F5GqItIPeAE4H3cB8BBQC5grIneHNkRjzGefwYgRcPXVbk6/iQ7p6ZCcDOXKua+2tK6JFUEnfxE5XUSeBjYCE4G1wJ9U9XRVHaWqrYGHgSFhidSYOLVjhyvfm5TkKruJeB2RAZfo09Lcugqq7mtaml0AmNgQ7FS/z4DlwBXAGCBJVbur6icFDn0P1wJgjAkBVejXDzZtcuV7a9b0OiLjN3Jk/gV2wL0eOdKbeIwpiWBr++8FugFvqWpOMcctw5YCNiZknn4a3nrLVfBr3drraEyg9etLtt2YaBJU8lfVzkEedwCwpYCNCYHMTLjtNujaFYZYZ1rUSUpyTf2FbTcm2gXb7F+nmMfxIlIp3IEaE09274bu3aFOHZg2zfr5o9Ho0ZCQkH9bQoLbbky0C7bZfzOgxexXEVkFPK6qM44+LGPil6obOLZ2raviVstG0USl1FT3deRI19SflOQSv3+7MdEs2OQ/GBgOZAH/BrYA9YCrgCrAc0BHYLqI5KiqjXc1ppSee84N7nv4YWjf3utoTHFSUy3Zm9gUbPJvgBvtf6Wq5rYAiMhI4C2guqpeJCKvALcAlvyNKYX//Q8GD4bOneGOO7yOxhhTVgU7z783MDkw8QP4Xk/y7QdX8tdG+xtTCnv3un7+xESYMcMVjjHGmHAI9s9Lou9RmGOB6r7nu4DipgLmIyIjRERF5JmAbSIio0Rko4hkichCEWlR4H3HisgMEdnle8wQkaLiMyYmDBoEq1fDzJlQt67X0RhjyrJgk/+nwGgROT1wo4icAYz27Qc4GQhqlquItAFuBL4usOt24FZgENAa2ArME5HqAcfMAs4B/gxc4ntuAw1NzJoxw43qv/tu6GQrZBhjwizY5H8zUB5YLiKrROQT3+j+/wLi2w9wHDDlSCcTkZq4cQE3AL8FbBdgKPCoqv5LVVfguhSqAz18xzTDJfw0Vf1MVT8HBgCXi8ipQX4/xkSNb7+FgQPh/PPh3nu9jsYYEw+CXdL3e6ApMAz4Csj2fR0KnKqqP/iO+6eqPhXEKacAr6vq/ALbG+FmEcwN+OwsYBHQzrepLbAH+CzgfYtxVQjbYUyUCxw5k5UF11wDVarArFlQIdghuKZIHTu6hzGmaEf8UyMiFYELgVWqOv5oP1BEbgSaANcXsrue7+uWAtu34GYc+I/ZFjj4UFVVRLYGvD/Xtm3bSElJyX2dlpZGWlpa6b8BY47CqFGwc6cr1ysCt9wCX38NPXpAgwZHfLsxxoTEEZO/qh4UkXdwTe1HVbXa1yz/MNBBVf8o7mMLvrXAtsIKDhU8BoDatWuTkZFR0lCNCTlVl/jHjXOvzzsPJk1yz2vXdvutkp8xJhKCbWRcQ2hW62sLHA+skLy/cuWB80XkJsA/qr8e8HPA++qQ1xqwGagjIuK/+/eNFajN4S0GxkQNEXfHD+4CwH8RcPPNeS0BxhgTCcEO+BsDjBCRY4/y894EzgBaBjwygFd8z7/DJffchYREpArQgbw+/s+BargLCb+2wDHkHwdgTNQRgXvuyb9t/HhL/KGSng5LlsAnn0BysnttjDlcsHf+Kbg767Ui8imwiQLN8Ko64EgnUdWdwM7AbSKyF9jhG9mPiDwFjBSR1biLgbtxA/xm+c6xSkTmAJN94wcEmAy8q6q2oqCJallZcOaZ+bcNG2Z3/qGQnu7WRDhwwL1et869BivBa0xBwSb/y31f9wGtCtmvuOl2ofA4UBWYgCsg9AXQRVV/DzgmFRhP3qyAt8mbbmhMVDp0yCX+jRvhz3+G995zid/f/G8XAEdn5EjYty//tn373HZL/sbkF1TyV9X64QpAVTsWeK3AKN+jqPfsAHqGKyZjQk3VJfoffnDz+d97L/8YgMRES/xHa30Rw5GL2m5MPLNZxcZEwD//Cc88A7fe6p77E73/AsAS/9FLSnJN/YVtN8bkF/TSISJSRUTSRGSmiHwgIk18268SkVPCF6IxsS093a3Qd+218Pjjhyd6S/yhMXo0JCTk35aQ4LYbY/IL6s5fRE4A5uNq9/+EK9JTw7f7UnzldsMRoDGx7KOPoG9fV3Fu+nRbqS+c/P36N9zgBv01bOgSv/X3G3O4YJv9n8S1EjQD1gKBBXoWAFaR3JgC/vtfuOoqOO00eOMNqFzZ64jKvtRUmDrVPV+40NNQjIlqwSb/i4GBqvqDiJQvsO8X8krvGmOAtWvdiP7ERPjgA/fVGGOiRbDJvzIF5ucHqI5b6McYA+zY4RL//v2u2d9q9htTCgXrXVv965AKNvmvAK4APixk38W4Ff6MiXtZWfCXv8BPP8G8edCixZHfY0LLmvvLgIIrYPnnyiYmun3mqJWkvO8AEXka+D/ftiYichdwo2+/MXEtO9v1OX/2Gcyc6ebzG2NKKHAFrGHD8hL/uHFuuxa2rpspqWCL/LwqIvWBh4C/+za/AmQBw1X1nTDFZ0xMUIXBg93AvnHj4OqrvY7ImBglAk8+Cdu25V8Ba8gQK4oRQkEX+VHVp0RkGm6RnTrAdmCRqv4WruCMiRWPPQbPPgu33eYuAowxJaAK330H8+fDxx+7vpvt2/MfY4k/pEpU4U9VdwHvhikWY2LSjBlw111w3XXw6KNeR2NMjFi/3iX6+fPdY+NGt/2kk6BrV/d67ty8420FrJAqUfIXkZZAElCl4D5VfS1UQRkTK+bNg3794MILYdo0K+JjTJG2bIEFC/KS/Y8/uu21a7v/QP5H48Zwyy2uKpa/qd9WwAq5YCv8NQX+jSvyU9hPXQFL/iauLFvmivg0bw7//rcV8TEmn5074ZNP8pryV65022vUcCUvBw1yyf700w9P5omJ+fv4bQWskBMNYuSkiHwENAXuAv4HHCh4jKp+G/LoQiAlJUUzMjK8DsOUMWvWQLt2UKkSfP45nHCC1xEZE0bBzLnfuxcWL85ryv/qK8jJgapVoX17l+g7dYKzz4YKQdx32jz/UhORTFVNKe6YYJv9/w+4QVVnH31YxsS2X3+FSy5xRXw+/tgSf1HS02HkSNe1m5RkdfZjVlFz7qtXh86d85rxlyyBgwddYm/TBu65xyX8c88tXbOYrYAVVsEm/x3AvnAGYkws2LfPFfFZt8719zdv7nVE0Sk9HdLS3M8L3M8rzbf0l10AxJDAOfc5OdCzp+uPX7zYJfmHHnJJuVUrd0Fw4YXuLv+YY7yO3BxBsM3+twB/Av6iwbwhilizvwmV7Gz429/g7bdh9mz33BQuOdkl/IIaNnTrHpgYsn69K1zx5Zd522rVctNbOnWCCy6AY4/1Lj5zmFA2+1cFWgBfi8iHuJaAQKqqj5QiRmNigqobn/TWW/D005b4j2T9+pJtN1FG1fVpPfus+6UveM+3bZs1w8e4YJP/gwHPC6tWroAlf1NmPfIITJwIt98ON9/sdTTRLymp8Dv/pKTIx2JKYOdOePFF98v+7bfuDn/4cJfsp03LO87m3Me8YGclVz3CIyEs0RkTBaZPdwPXUlPdRYA5stGjIaHAX4WEBLfdRKHly92gjAYNYOhQ14z/0kvw889w4IBL/EOGuH7/IUPy1903MSnY2v6HTe0zJh58+CHceKPr2nzhBSviEyz/oD4b7R/FDhyAf/0LJkxwq1FVrQo9esDf/w7nnJN3nM25L5OKHPAnIl2Az1X194BtCaq6r8BxdYBuqvpsWCMtJRvwZ0orM9ONZWrSBBYtcrVJjIl569fDpEnw3HOuOb9JE5fw+/QpeuCezbmPKcEM+CvuPuYDXEU//8nKA7+LyDkFjmsEPF3qKI2JQj/9BJdeCscfD++/b4nfxLicHFcn/4oroFEjtxJVu3auaevbb10TfnEj9m3OfZlTXLN/Yf+69i9uyjx/EZ+DB93iYlbEx8Ss335zg1YmToTvv3d19O+8EwYMsNGXca5EC/sYU9bt2weXX+7GOX30ETRrduT3GBN1vvrKTdObNQuysuC881ylvr/9zRahMIAlfxPnArsuDx2Ca691tUz+9S/399KYqFNU//v+/a761LPPulK7CQlw/fUwcCC0bOldvCYqWfI3cSuwZDm4MU/vvAN//jP89a+ehmZM4Qqrs9+vn2vS//Zb12d16qluKl7v3lCzptcRmyh1pOTfW0Qu8j0vhyvm009ELgk45sSwRGZMGAWWLAdXy2TqVPe8aVMbzGyiUOAvrSpcfLGrOLVmjftl/etf3RXshRfaL685ouKm+uWU4DyqquVDE1Jo2VQ/UxT/4mT+CwCAwYPhqafsb6eJUnv3QteusGBB3rb/+z94/XU46STv4jJR5Win+h2pqp9V+DMxTQS6dMm/zRK/iUpbtrglchs2zJ/4wfXvW+I3JVRk8lfVAyV5RDJoY0Jh6VI37TmQVSw1UWXlSrjhhrwSiR06QLdu+Y+xX1pTClas1MSlH35w1fsOHYL+/a1kuYki/hX1Lr0UTj8dXn7ZXQCsXu3u/F9/3ersm6Nmo/1N3Nm2zY3oBzcTasoUK1luosAff8Crr8KTT7qFdurWhQcfhJtucqUmwersm5ApcsBfWWED/kygvXvdYOivv3Y3V23bxkfJ8vR0d/N44IC7ebRFdqLIb7+5K9Dx42HjRmjeHG691S2yU6XK4cdbnX1zBMEM+LM7fxM3Dh2C7t0hI8MV8WnX7vBjyuLf0PR0t1rrAd/InHXr3GuwCwBPrVnjRpg+/7y7Ku3UyT2/+OLifxGtzr4JAevzN3FB1U2Bfu89t4LplVd6HVHkjBzpyhYH2rfPbTce+OILuPpqt5res8/CVVfBsmWunvQll1gyNxFhd/4mLjz4oCviM2KE60KNJ+vXl2y7CYPsbHj7bdefv3ixq7x3220waBA0aOB1dCYOFZn8ReT9EpxHVfWyEMRjTMg9/zzcdx/06gUPPeR1NJGXlOSa+gvbbsJs7163qt7YsfDjj2453XHjXEneatW8js7EseLu/I/DlfM1Jma9/75bvfTii+G55+KzRXX0aNfHH9j0n5DgtpujVNTgu02b4JlnYNIk2LEDzj0XHn3U9TdVsAZX470ifwtVtU0kAzEm1JYudV2rZ53lFjurWNHriLzhH9Rno/1DrLBFdq6/3hXm+eYbOHjQJftbb3WjS+PxytNELbsENWXSDz/AZZe5qdLvvQfVq3sdkbdSUy3Zh1TBRXYuuwxuvNENpKhQwTW1DBvmBvUZE4VKNM9fRI4BTgYOm3yqql+GMK6QsXn+8WfrVnejtXMnfPaZW6XPmJDbssWV2v3Pf/K2tW3rBvb5i/IY44GQzfMXkUrAJKAnUNTqfVG5qp+JL3v3wuWXu1op8+db4jchlp3tpuRNnQpvveWKRwRavNia901MCHae/wjgMmAgIMCtwM3AUuBH4KqwRGdMCRw6BNdcA5mZ8Mor0MZGrZhQ2bABHngATj7ZzcVfuNBN07v++vzHWY19EyOCTf7dgQeA6b7Xi1R1om9Q4DfA+WGIzZigqbr5+++/7+qm/OUvXkdkYt7Bg/Dmm64/v2FDN1/0lFNc/f0NG9zCOjNm2CI7JiYFO+CvIfA/Vc0WkYNAQsC+KcDzuNYAYzxx//1uPv/dd7upfcaU2o8/unmh06fD5s1Qvz7cdZebm9+4cd5xtsiOiWHBJv/tgL8ixQbgTMA/yiUROCbEcRkTtKlTXfLv08e1zEa79HRXWnf9+rxl2m0kvsf274c33nBJf/58KFfO3fH37++W1i1sbv6oUfnn+fsvACzxmxgQbPJfikv47wNvAg+ISGXgEHAn8Fl4wjOmeO++CwMHum5Y/9K80cy/yI6/4I4tsuOxlStdwn/pJVeMJznZ1YLu2ze4sru2yI6JUUFN9RORNkCyqr4iIolAOnAJbvDff4GrVfXHsEZaSjbVr+z64gv405/cCqgLF8ZGtdTk5MJL7TZsCGvXRjqaOLV3L7z2mmsy+vxzV/3pyivdPP1OndxdvzExLGRT/VR1CbDE93wncJmIVAMSVHXrUUdqTAl9/72b0le/viviEwuJH2yRnbArbq37zEx3lz9rFuzeDaeeCk884RZ9qF3bm3iN8Uiw8/xvB15S1c3+baq6B9gjInWB3qr6eJhiNCafLVtcMz/AnDmuil+ssEV2wqiwcrt//7u7stq0yS2bW6WKmw96441w3nnWTG/iVrB9/o8AC4HNhew70bffkr8Juz173B3/pk2wYIGbeRVLbJGdMClYbvfqq90PetUqt79lS5gwAXr0cKPxjYlzwSb/4i6PawJ/hCAWY4p18KC7afvqK1dc7dxzvY6o5PyD+my0f4iJuCb81ath/Hj3ADj9dHjhBUhJsbt8YwIUmfxFpD35i/f0EZGLChxWFbgCWBXMh4nIXbhqgKcCB3DjCO5S1RUBxwhwH5AGHAt8AfxDVVcGHHMsMB7wl3J5GxjkG49gyiBVN3//gw/cqP7LL/c6otKzRXZCLDvbDeB74AGX/AN9/bUlfWMKUdyw1k7AQ76HAjcFvPY/RuAuIIYG+XkdgWeBdsCFuKmCH4nIcQHH3I4rGDQIaA1sBeaJSOC6bLOAc4A/42YdnAPMCDIGE4Puuw+mTXNfb7zR62hMVMjJcUn/jDNcc36FCm5ufiCrtmdM4VS10AduoZ7KuBX8coD2vteBDynq/cE8cIWDsoGuvtcCbAJGBhxTFfgdGOB73Qx3MXJewDHtfdtOLfgZrVq1UhPbJk9WBdUbblDNyfE6GuO57GzV119XPf1094vRvLnqq6+qDh7sXg8Z4n5RhgzJ/9qYOAFk6BHyb5HN/qqa7UvMiEhVVT0QsiuOPNVxrQ+/+V43AuoBcwPiyBKRRbjWgslAW2AP+QsLLQb2+o75NgxxGo+8844r4nPppTBxorXgxjVVN9jjvvtcc/5pp8HLL7vBfeXLwzffWLldY4IU7Dz/A76KftcDFwDH4Ur+LgTSj+LCYByuSNDnvtf1fF+3FDhuC9Ag4Jhtvqsbf3wqIlsD3m9iVOC07CVLoHt3aNXKte5WrOhtbMYjqq6U4333uel6TZvCzJlw7bUu6ftZuV0Tazys9R1UKSsRqQ1k4BbxuQg4AegMPAcsFZHjS/rBIjIG11z/N18rQ6CCnXRSYFthnXgFjwFg27ZtpKSk5D6mTJlS0lBNhIwalddF+913blBfpUpwwQVwjK0eEX9UXQWn//s/t0zj7t3w4ouuJG9qav7E72fldk2s8Nf6XrfO/a77a32np0fk44OtY/kYUB/orKr1VfVsVa2PuwCo59sfNBEZC1wHXKiqPwXs8tcRKHgHX4e81oDNQB3frAD/+QSozeEtBtSuXZuMjIzcR5q/kLqJKoHTtNPSXBGffftg1y43xS9cY7bS013J3XLl3NcI/b8zxVF11ZvatHFXgNu3u+l6q1e7anyFLbJjTKwZOTJ/wQ9wr0eOjMznH2lQgK91fStwYxH70oAtwZzHd/w4XJJuVsg+/4C/EQHbqgC7OXzAX7uAY9phA/5iXk6O6oABboyW/xHOsVozZ6omJOT/vIQEt914ICdHde5c1bZt3T9Gw4aqU6eq/vGH15EZE3oi+f/4+B8iR31qghjwF+ydfw2gqOrj63z7j0hEJgB9cXf9v4lIPd+jmu9CRIGngDtF5CoROR2YjhvgN8t3zCpgDjBZRNqISFvcQMB3VdUG+8WwX391i/UECmeXrdcX3ibAggVw/vnQpQts2ACTJrm+n/79bbCHiYxINwMWVdM7QrW+g03+3+ESdmG6+/YH4++4EVIWz7IAAB66SURBVP4f4+7w/Y/hAcc8DowBJuDGGdQHuqjq7wHHpALLcbMCPvQ9vz7IGEwU2rjR9e2vWJF/ezinadsiO1Hgk0+gY0e48EL46SdXgvf7711Fp0qVvI7OxAsv+t9Hj3a1vQNFstb3kZoG3M04fXFz/d8DegB/wl0MvIWbDtg7mPN48bBm/+i3Zo1q48aqFSvmb+oP9zTthg0Lb3Vr2DD0n2UK+PRT1QsvdD/wevVUx49XzcryOioTr7z6YzBzpvsMEfc1RH2OHM08/wIXCNN8FfbuxVXVU1z//A5gqKq+GOJrEhMnvvsOLroIfv8devd2o/ojNU3bFtkJMy1ked0lS9yUvXnz3HKMY8e6u/yqVb2L0xivmgE9rPUd9LBZVR0vIhOB03Hz/HcAK1T1YLiCM2XbihUu8efkwMKFcNZZkZ2mbYvshFHB5XW/+MIVbVi3DmrXdovwDBx4eLOnMX6RnAMfh2ttF7ewz0/AX1V1uX+bL9Evi0RgpmzLzHRju6pUcYn/tNPc9khP07ZFdsJAA+ZtbtgA+/e7+foA553npvFVq+ZtjCa6+fvg/c1y/j54CM9/2DhsBhQtYjSViOQAbVT1y8iGFFopKSmakZHhdRgmwOLFrlzvccfBxx9D48ZeR2SOyr59sGqVa8pZuTLvUfBOauBAN6DPCu+YI0lOLvxOvGFDWLs2PJ/pYbW9UBORTFVNKfYYS/4mkj7+2BVrO/FE9/zEE72OKE4U1v9e0iS8f78rtFMwya9Zkzclo1Il14zTogU0bw733JP3/pwcS/wmOOXKFT7NR8T9HpliBZP8j9Tnb2thmpB5913o1s2VZveP9zIRULD/XdXNoUxMdPsKOnDAjcQsmOR//DHvD2+FCnDqqZCS4kZqtmjhHk2auH3+zwg0bJjV2o9Vkb4rjsM++Eg7UvK/X0R+DeI8qqq9QxGQKZtmz3ZLrrds6bp8a9XyOqI4Edj/Di75DhvmXt98c/7k7n98/z1k+5bbKF8eTjkFzjwTrrsuL8mfckrR8/D9iX/cuLxV9vyv/THYBUDsiHT/O8RlH3ykHanZfzMQzIp9qqpR2XNrzf7ee/FF6NcP2rVzd/81a3odUZxRhaFDYfz4vG21armFcg76JuuUKwcnn5yX3Fu0gNNPd800lSuX/DNL2tpgopcX/e9QpvrgI836/LHk77Vnn4V//MNN6XvzTVudL+I2bYLp0+H5512zvV/XrnkJvkUL14Qf6rn2oRhnYLxn/e8xJ5jkH2x5X2NK7IknXOLv2hXeeSe4xG+r7IVAdrabWvfXv8JJJ8GIEW6wXqDGjeHhh92dVMuW4SmyY8vrhkec1aA34WHJ34ScqmvZve02V9flX/9y8/mPxOPlrWPfunWuel5yslsK97PP4JZb3DK4v/zi+t9zctzXcePCu3CCCY94rEFvwsKa/U1IqcLtt7u7/j594Lnn3JixYHjVtRjTDh6Et9+GqVNh7ly37eKL3Wp4Xbu6QXnW/152WP+7CcJR9fmXFZb8Iycnxw0gnzjRNfePH+9aJoNlXYsl8N13rh9/+nTYutUVTOjXzz0aNjz8eOt/LxvsP4kJgvX5m4g5dAj69nWJ//bb4emnS5b4wboWjygry92BdezoBug9+WTeFIq1a+H++wtP/GD97+Fi/e8mRlnyN0ftjz/cHP6XXoIHHoBHHy1dbrGuxSL8738weDA0aAA9e8LPP7vBej//DG+8AZddFnzfigkd6383sexIa/7G+qNVq1bBLoFsSiErS/Wyy9zS108+efTnC9Py1rHn999Vp05VPfdc98OtVEn12mtVP/5YNTvb6+iMaplbA96UHUCGHiE3Wp+/KbU9e+CKK2DBAtfcP2CA1xHFiKL631UhI8MN3nv5ZfcDbt4cbrwRrr/eyiJGG+t/N1HK+vxN2Px/e/cdJ1V97nH88whioaqASLFwsdB82VLEoBTRVaPRRIMJtli4lqDXEhsxYo8aECJBglcjCUGxaxQpQVEpsSVeRNEoQlYpIkUEBYTd3/3jOePODjPbmNmzM/N9v17zGubMKc+ZWeY551fXrvVG5TNn+gh+Svw1NGxY5S52Ifhsd8cfDwcfDN/9rhcbn3aad9WbP99H51Pir5n6rINX/bvkserG9hfZysqVnvjfeQcefRR+8pO4I8oTyePshwA//jFceKHPlAdw6KFehPKzn2kM5LrQHPAiNaZif6mV5ct9qN6PPoInn/QbVqlCebn3jU6eOGfKFO+el9Czp3fZO+SQ2MIsCJoDXgRQP39AyT+bPvkE+veHpUt9XJl+/eKOqAEJAT79dOsZ8t57z+vuE9q397H0p0+vWKZ57rNDdfAigOr8ZRsl/44uXAi9e/sN67RpBZT4U5NFdRfDIfjVz/TpMHKkN8br1ctHy9tzTzjuOLjqKpg8GZo188EPxo6FWbNg9Wq/QOjWrfI+C3WYXfWBF2mwVOcvaSWPCLtggRf1r1njw8T36hV3dFlS1bC3N97oVzrJd/Hz5/vzF19U7KN1a7+TP+OMihnyundP30CvmOa51xzwIg2akr9sJbld2v33+29pIjfutFOBjAybfJJlZd66/sYbvftChw4wejSsWlWx/i67eFI//fTKc963bVvzY5r5hUUi8Zv5M/jyvP9QkwwdWjkJg78eOjR3yT+xX9XBi1RLdf6S1kMPwXnnVa4qbdzYlxfEb+nmzfDii3DllX43n9CkCRx2WOUE3707tGuXveRcDOPsq/5dJDY1qfPXnb9s5aWXvCo79Td6y5bc3rjl3JYt8PLLMGmSd1VYtQqaN6+8zoYNtZ+UoLbiGGe/vlul77ln+pb3qn8XaRDU4E++tWKF1+n36+d5Mp3S0vqNaZuVlXlR/kUXeUv7o4+GiRPhmGP8AuCssyqvf8UVhdf4TmPQi0gKJX+hvNzr9g84AB55xG8QmzVLv26zZnmQG8vL4dVXfX7hDh2gb1+fdahvX3jiCfj8c098L78Mf/iD18GXl/vzqFGF1/q+qvr3XBk0CMaN8z72Zv48blweFxuJFBYV+xe5efN8kLm5c+Goo3yAua5dvYX/s89WLgFo3BgGDGig1dPl5X4Sjz4Kjz/u3fF23NFnvBs40Ecjatq08jbF0vguU3FNrotxBg1SshdpqKqb+SffH5rVL71160K46qoQGjUKoXXrEMaPD6G8vPI6DX7ysPLyEObODeHyy0Po2NFnVNthhxBOPjmEhx/2k6zJPqp6nQv1/cHGNfuciMSCGszqpzv/IvTMMzBkiI/Yd/758Nvfpu+WHsuNW3Ut4RMz3z36qD9KS72FfkmJn8iJJ0KLFjU/Xn03vlP/dxFpANTVr4iUlsKll3ry79HDB5474oi4o0qSadCdli197uBEwl+0CLbf3hvt/fSncNJJXlSfD+IYfx40Br1IEVFXPwG8S/uoUT6GDcBdd/kssdtvH29clSQPugMwYoTPYT9xoif/m2/2RgdHHw033AAnn+wD7+Qb1b+LSAOg1v4Fbu5cH7PmV7/ySXnee8//3aASP/id/t13+0h7o0ZBo0ae+M18jvv77/cpBV94wcfLz1bi1/jzIlKEdOdfoNasgWuv9d5VHTvCU095yXmDa8S+ZAlMnerT3E6fXnncfPCEX5shdGtD9e8iUqR0519gQoAJE2D//eGBB3zMmgULvJS8QST+TZt8WN2rr4YDD/Qrk/POg9mz4ZRTfFa8ZLffnrs+9+r/LiJFSsm/gHzwgRftn3kmdO7sjeKHD888YE+9+fhjGDPGG+bttpsHOXIktGnjDRDmzfOuBy1aeLF+fQ26E2f9++LFfo6LFyvxi0i9U7F/AdiwAe64A+6800uQx471sflzPUR9Rl9/7UPqTpnijw8/9OX77ANnn+3d8vr23fqqpFUrOPZYr6P4/e+hUyd/natBdzT+vIgUKSX/PJGp+/u0aXDxxbBwod9ADh8Ou++e44OmW2/Bgopk/8orXry/006e5IcM8YTfpUvVSXzffb3RX6IovrQUVq70ooxcUP27iBQp9fPPA8OGwT/+4fn1k0/8hrhzZ8+L8+d7zrzvPi9Nz+pB0/W5b9XK31u7FmbMqEj4n3zi23Xv7on+2GOhd28fYrem4ugDr/7vIlJgatLPX8m/gQvB28BNnbr1e2bed/+aa2qXY2t00Msv9/r25s1h/Xovol+3Dg4/3Pvbz5njM+a1aOED/icSfqdOdT+u5oAXEdlmGuSnAJj5HX86e+xRMXBP1g/6ne94kl+3zpclnufOhUMP9X6EJSXwve9lb9AA1cGLiNQLJf88kChRT7VsWZYPtGULvPaaF+PfdVflKf0SOnb0bgS5oDp4EZF6oeTfgK1f7zk4k20pYf/Wp59WHmRn7VofXa+sLP36S5Zk4aAZJOraVQcvIpJTSv4NUFkZjB/vOXD5cmjXzp9Tde2auQF+Rps2wauvVjTUe/ddX96xow+tW1IC/fp5I7tEUX+yZs3qcNBa0Bj0IiI5p+TfwMyY4aPyzZvnbeueespzdGpr/65d4fvfr2EO/uijimT/0kterN6kCRx5pI+TX1IC3bpV3tmAAfDss5WL/hs39uUNYqhAERGpKyX/BuL9933Cneee8x5vkyb5jbiZJ/lQHrDtKpJu6utK1q+vPMjOwoW+vEsXOPdcT/Z9+kDTppkDeuIJdYMTESlQ6uoXs5Ur4aabvJ9+06aeay+9NKXrXnV97kPwDv+JZD9rFnzzjTeW69evohtely4xnaWIiNQXdfVrwDZtgtGj4ZZb/EZ98GDP41tNYJc8z/2DD1buc3/88X4nP3UqLF3q6/fs6WPil5TAEUfADjvU96mJiEgDp+Rfz0LwEvVrrvH5bo4/3ke07dYtwwZV9bmfPNnv/hOD7BxzjDfcExERqYKSfz16/XVvzDd7NvTo4TfsxxxTxQaffeaD9194Yfo+97vv7l31GutrFBGRmlPWqAelpXD99d5+rm1bn7793HO9O30lmzf7CHqJfvf//GfVO16xQolfRERqTZkjh9at82l2hw/319df76PiNm+etFJpaUVDvRkz4Msv/aqgVy9vXX/ssT4zXhx97kVEpCAp+edAWZm3zbvhBi+5HzQIbr89GqJ+40aY9kpFwk8M3N+pEwwc6HX3/ftDy5YVO1SfexERyaK8Tv5mdjHwK2AP4F3gf0IIr+b6uFX1uZ8+Ha68Et55xxvbP/tM4Lu7fAhPR8l+5kzYsMFb4R91FFxwgSf8Aw7InMjV515ERLIob/v5m9lAYAJwMTArev4F0C2EUJpYL9v9/Gf2GcaKd5bTa+0LtC/7hKWNOjGn5XF83bknj7W9hMmTocde6xg78CV6fTkFmzoFFi3yjffbzxN9SYkn/p13zlpcIiIiULN+/vmc/F8D5oUQLkha9iHweAjhusSybCb/UB54rM1FnLD6zzRlw7fLv2JnruM2dtlxI+d1mEan0lnY5s0+ak///hWD7HTunJU4REREMinYQX7MrAlwKPC7lLemAb1ydtztjF5rX6iU+AGa8jVX8zs6blwCTQ/00fcSg+w0aZKrcEREROokL5M/0BpoBHyWsvwz4OjkBZ9//jmHHVZxATR48GAGDx5c5wO3L/sk/XKW+nS37dvXed8iIiL1IV+Tf0JqnYWlLmvTpg3ZrPNf2qgTHctK0y9X4hcRkTywXdwB1NFKoAxol7K8LVuXBmRNKA/MaXkcX7FTpeVfsRNzWh5HKM/P9hMiIlJc8jL5hxC+Ad4CBqS8NQCYk6vj2nZG257teH7Xs/i00Z6UY3zaaE+e3/Us2vZsl3mKXRERkQYkn4v9RwB/MbPXgdnAhUB7YGwuD9pn5rBK/fo7Aqel9PsXERFpyPI2+YcQJpnZbsCv8UF+5gPHhxD+k+tjpyZ6JX4REckneZv8AUIIY4AxccchIiKST/Kyzl9ERETqTslfRESkyCj5i4iIFBklfxERkSKj5C8iIlJklPxFRESKjJK/iIhIkVHyFxERKTJK/ttg3LhxcYdQL4rhPIvhHEHnWUiK4RyhOM4zjnNU8t8GxfBHCcVxnsVwjqDzLCTFcI5QHOep5C8iIiI5ZyEU9hz0ZvY5kKvJfloDK3O074akGM6zGM4RdJ6FpBjOEYrjPLN9jnuFENpUtULBJ38RERGpTMX+IiIiRUbJX0REpMgo+YuIiBQZJf86MrOLzWyRmW00s7fMrHfcMWWLmV1nZm+Y2Zdm9rmZ/c3MesQdV66Z2fVmFsxsdNyxZJuZ7WFm46Pvc6OZvWdmR8UdV7aYWSMzuyXp/+QiM7vVzBrHHdu2MLMjzexZM1sS/W2ek/K+mdkwM1tqZhvMbKaZdY8p3Dqp6hzNbHszu9PM5pnZV2a2zMwmmtmeMYZcJ9V9lynrjovWuSpX8Sj514GZDQRGAbcDBwNzgBfy8Q8ygz7AGKAX0A/YAvzdzHaNM6hcMrPvAxcA8+KOJdvMrBUwGzDgBKArMARYEWdcWXYNcAlwKXAAcFn0+ro4g8qCZsB8/Hw2pHn/auBK/Pv8Dv6dTjez5vUW4bar6hx3Bg4BbouefwR0Aqbk4YVddd8lAGZ2Kv5dLs1pNCEEPWr5AF4D7k9Z9iFwR9yx5eh8mwFlwIlxx5Kj82sJLMQvdGYCo+OOKcvndzswO+44cnyOzwHjU5aNB56LO7YsnuN64Jyk1wYsA4YmLdsJWAf8d9zxZuMcM6zTDQhAz7jjzfZ5AnsBS/AL9MXAVbmKQXf+tWRmTYBDgWkpb03D75QLUXO8lGhN3IHkyDjg8RDCi3EHkiMnA6+Z2SQzW2Fmb5vZL83M4g4si2YBfc3sAAAz64ZfzE2ONarc2gdoR9JvUQhhA/AKhftbBNAiei6o36OoJONh4NYQwoJcHy/fik0agtZAI+CzlOWfAUfXfzj1YhTwNjA37kCyzcwuALoAZ8YdSw51Bi4G7gF+CxwE3Bu9VyjtG+7EL1LfM7My/LftthDCmHjDyql20XO636IO9RxLvYhuvoYDfwshfBp3PFl2E7AqhHBffRxMyb/uUkdHsjTL8p6ZjQB+APwghFAWdzzZZGb740XivUMI38QdTw5tB7wZQkjUf//LzPbF68QLJfkPBM4Cfg68i1/gjDKzRSGEB2KNLPeK5beoMTABaAWcFHM4WRU1vj0H/7utFyr2r72VeP13u5Tlbdn6Cjyvmdk9wM+AfiGEj+OOJwcOx0ty5pvZFjPbAhwFXBy93iHe8LJmGfBeyrIFQKE0UAW4G/hdCOGREMI7IYS/ACPI/wZ/VVkePRfDb1GiSPxAoH8IYVXMIWVbX2APYFnSb9FewJ1mlpMSDiX/WoruEN8CBqS8NQBv9V8QzGwUfhfVL4Twftzx5MjTQE/8ajvxeBN4JPp3oZQGzAb2T1m2H7mb8yIOO+MX5cnKKOzfuEX4BcC3v0VmtiPQm8L6LdoemIQn/r4hhOXVbJKPxuDnl/xbtBSvquufiwOq2L9uRgB/MbPX8R/WC4H2wNhYo8oSM/sDXgd+MrDGzBJ3FutDCOvjiyy7QghfAF8kLzOzr4DVIYT58USVE/cAc8xsKP4jejDeJe76WKPKrr8B15rZIrzY/2DgCuDPsUa1jcysGd4mBfxCZk8zOwj/Gy01s5HAUDN7H/g38Gu8JfnEWAKug6rOEU+Aj+Fd304EQtLv0dqogWNeqO67JKXrrZltBpaHED7ISUBxd3nI1wfegGoxsAkvCTgy7piyeG4hw2NY3LHVw7nPpMC6+kXndQLwf8BGPElcSjSxVyE88MZ+I/HSjA3Ax3h7jh3jjm0bz6tPhv+LD0XvGzAMr9rZCLwM9Ig77mydI7B3Fb9H58Qdeza/yzTrLyaHXf00q5+IiEiRKeT6MBEREUlDyV9ERKTIKPmLiIgUGSV/ERGRIqPkLyIiUmSU/EVERIqMkr8UHTM7x8yCmX1hZrukvNc4em9YDHENi47doAffMrPtzGykmS0zs3Ize3ob97d3dN7nZCnE5H23ij7XQ7K971xL+ntYGI1yl/xel1x9ZlIclPylmLUErok7iDx0KnAZPp7+EcDV8YZTpVbAjUDeJf8knYHz4g5CCouSvxSzacCQpOFCC16WJivqGj2PDCHMDSH8Owv7bBAa6GRO04BfR+P2i2SFkr8Us1uj56FVrZQofk2z/CEzW5z0OlF8faGZ3WFmy81snZlNMLOdo6LaqWa23sw+MrOzMxyyq5m9ZGZfR0XrN5tZpf+rZtbazO4zsyVmtsnM3jezwSnrJKo3jjSzx8zsC+C1as61xMzmmtkGM1trZk9HUx8n3l+MDycLUFZd0XNUjXJdFN8mM1tqZsNrksjM7CgzmxF9hl9Fn12PNOudYmazo8/1SzN73cxOMrO98clvAO6PYv02XjObaWazzOxEM/uXmW3Ch+3GzFqY2ego3k1m9oGZXW5mlnTcPtH+TorWXWlmn0ffd6uUGC8zswXR57rGzN40s1Oq+wwiQ/G5Qy6p4foi1VLyl2K2DJ/PfrCZ7ZXF/V6H/1ifDfwGn2t+LPAU8DxwCjAP+JOZdU+z/dPA3/GJlSYCN0T7ATwx4RNKnYAn4hPwiW3uM7Mhafb3VzwJngpcmyloMyuJ4lsfxXwR0AOYZWYdotVOwcdcB58S+fBom0wm4JPNTIzivAMvwv5rFdtgZicAM6JYzsBnmGwOvGpmnZLWGwI8iU+KcjZwGv45741/vz+OVr0jQ7z7Ab8H7gWOBWZEF1rPA78AhuMTykzBJ/S6LU24o/Ax2n8O3Az8JFqWiHFQtJ+HgeOBQcDjwK5VfQZJ3sYnt7nWzJrXcBuRqsU92YEeetT3AzgH/7Hugv8AfwE8GL3XmJRJjPAEG9Ls5yFgcdLrvaNtX0xZ78lo+RlJy3YBtgA3ph4HuDZl+/uBdUCr6PUN+CQu+6ZZbyXQOOU876nh5/Im8GFi+2jZPsBmYETSslvTfR5p9tc7Ov5ZKcsHRcsPSvnczkla5yNgRsp2LaLzG5n0eh3wZBUxJPZ9fpr3ZgLliTiSlv8wNZ5o+f/iE3m1jl73idYbn7Le6Oj7saTX/6zD32ni76ExfpGyBfhN9F6XdDHqoUdNH7rzl6IWQliN35WdlVy8vY1eSHn9fvQ8Nem4a/C71U5s7dGU148AzfC7cIASvPh+UVSs3ti8h8BUYDegW8r2T1UXsJk1xRvFTQohbEmKcxFeynBUdftIowT4BngiJc5p0ftHZohlX+C/gL+mbPc1MDdpu1745zKuDrElLA4hvJ2y7Ej8ouDhlOUTgCZ46UGy1JKPd4AdgN2j128AB5nZvWZ2tJntXNsgg7erGA9caWY1LTEQyUjJX8Tnu1+NF9lmw5qU199UsTxd3fdnGV4nit7b4glqc8rjsej93VK2X1Z9yOyCTw+bbt3l1LyIOllbPFmuT4kzMW95apzJ2wE8wNbn+MOk7RLPn9YhtoR057srPsf6ppTly5PeT7Y65XViu8R3+2e8CuV7+AXaajN7MmqTUBs34RcV6qEi26xB9ycWqQ8hhPVmdgdeAnB3mlU2AphZkxDCN0nLMyWvbbU7Ph998muAJdHzKjyBXpZh+w9SXtdk3u410Xrpej60i45ZW6vwz653hveXVrEdeNuJv6d5P/EdrIyeOwDz6xAfpP9sVgO7pvm+E59NrT6LEEIA/gj80XxciWPwv7VJ+AVBTfdTamZ/BH4JPFebGERS6c5fxI3Bk+utad77T/T8bUvzqDV3rxzF8tOU16fjd8+JBDcFOAAoDSG8meaxrrYHDCF8BbwFnGZmjRLLo4aQvYCX63AeU/C735YZ4syU/D8AFgPdM2w3L1pvDv65DM6wH6i4C9+pFnG/jP82npayfBB+4fGPWuyrkhDCmhDCJLxqZ6ueCzVwG37BUmUPFZHq6M5fBAghbDKzm0lff/wCsBbvLnYjXvR6NZ54cuGCqMX5G3gL9PPxBohfRO/fg7fGf9XM7sGTZVP8gqB3COFHdTzuDXj99XNmNgavT78JP/fhtd1ZCGGmmT0MPG5mI4DX8br0vfFW79eENGMEhBCCmV0CPGNmTfBEuRIvAemFX/SMCCGsM7PrgHvN7Am8B8E64CBgYwjhXrzKZBVwupnNA74CFoUQqrp7fwGYBYw1szbAu1G85wN3hBBWVrHtVsxsXBTXXLzEZj/gTCraPtRYCGGFmY0Crq/ttiLJdOcvUuFPeGv3SqKk+0M8cT2Kdxu7F3gpR3H8CBgAPIt3c7sVuCUpnrV4EpyM1/9OBR6MtqtzTCGEKXh3vFb4eY4FFgA/qOIuvTpn4K3WTwWewbu4/RL/nFPbNiTHMhlv19AUb2U/FbgLL3qfm7TeaPwOvSOe/J+IjrUoer8cT9q74FUIb+Bd9zKKtjkBb2B3DX5BdAJwBXW7454NHIqXLk2P9jEB75pYF3ezdfsRkVpJdEURERGRIqE7fxERkSKj5C8iIlJklPxFRESKjJK/iIhIkVHyFxERKTJK/iIiIkVGyV9ERKTIKPmLiIgUGSV/ERGRIvP/BhUp/QNh+pIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1)\n",
    "fig.set_size_inches(8, 6)\n",
    "ax.tick_params(axis=\"x\",direction=\"in\",which='major',labelsize=14)\n",
    "ax.tick_params(axis=\"y\",direction=\"in\",which='major',labelsize=14)\n",
    "ax.set_xlabel('Number of electrons N',fontsize=16)\n",
    "ax.set_ylabel('Total Energy E(N) [meV]',fontsize=16)\n",
    "\n",
    "electrons = np.array([0,2,4,6,8,10,12,14])\n",
    "electrons2 = np.array([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14])\n",
    "E_N_paper_a = np.array([0,25,50,75,105,120,155,170,215,250,290,327,350])\n",
    "\n",
    "plt.errorbar(electrons2 , 11.61 * E_N, yerr = 2*np.sqrt(N_av) * 11.61 * E_Nvar, fmt='bo', label = 'RHFc')\n",
    "ax.scatter(electrons2[:10],E_N_paper,c='b',label = 'UHFc',marker = 'x')\n",
    "ax.plot(electrons2[:10],E_N_paper,'b')\n",
    "\n",
    "plt.errorbar(electrons2 , 11.61 * E_N_a, yerr = 2*np.sqrt(N_av) * 11.61 * E_Nvar_a, fmt='ro', label = 'RHFa')\n",
    "ax.scatter(electrons2[:13],E_N_paper_a,c='r',label = 'UHFa',marker='x')\n",
    "ax.plot(electrons2[:13], E_N_paper_a, 'r')\n",
    "\n",
    "ax.legend(fontsize=15, loc = (0.03,0.66))\n",
    "ax.text(0.3,550, 'std 95%', style='normal',\n",
    "        bbox={'facecolor': 'white', 'alpha': 0.8, 'pad': 7}, fontsize = 13)\n",
    "plt.savefig('energias.png')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAF7CAYAAADPB8KEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8nHW5///XlXXaJjMFks5AWQpCC1oQsYqyg6JyAEW/oh4VwS2en+J24Ci4sQgqynpQVPQocBDUIgoosnkARRAIi7TsSFm6pStp0jbNdv3+uO+U6TSZ3Elm5p7MvJ+PxzySe5m5r0naXPP53J/P9TF3R0RERCpHTdwBiIiISGEpuYuIiFQYJXcREZEKo+QuIiJSYZTcRUREKoySu4iISIVRchcREakwSu4iIiIVpi7uACaipaXFZ82aFXcYIiIiJfPQQw+tcvfWfOeUNLmb2enA+4A5wCbgH8Dp7r5whPMvBz4N/Je7n597fNasWbS3txcxYhERkfJiZi+Odk6pu+UPAy4DDgCOAPqBO8xs29wTzez9wJuApaUMUEREZLIracvd3d+ZvW1mJwCdwIHATVn7dwEuAd4O/LmUMYqIiEx2cQ+oaw5jWDu0w8zqgGuBc9z9ybgCExERmaziTu6XAI8C92XtOwtY7e4/jickERGRyS220fJmdiFwEHCQuw+E+w4FTgL2jfIaK1euZN68eZu329raaGtrK3ywIiIik0gsyd3MLgI+BBzu7s9nHToc2B5YZmZD+2qB88zsS+6+Y/brtLa2arS8iIhIjpIndzO7hCCxH+buT+Ucvgy4LmffrQT34H9WgvBEREQmvVLPc/8RcAJwHLDWzDLhoW5373b3FcCKnOf0Acvd/elSxioiIjJZlXpA3WcJRsj/BViW9Ti1xHGIiIhUrFLPc7fRz9rqObOKEIqIiEjFinsqnIiIiBSYkruIiEiFmdSrwhXS4KDz3h/fG2sMH5i3Ix/Zf5dYYxARkclPyT3L9Cn1sV17wZJOrn94iZK7iIhMmJJ7qKbGuPITb47t+l/69SO0v7h29BNFRERGoXvuZSKdSrBi3SbcPe5QRERkklNyLxOZZILegUHWrO+NOxQREZnklNzLRCaZAGD5up6YIxERkclOyb1MpFNBcu9QchcRkQlSci8Tm1vunZtijkRERCY7Jfcy0drciJm65UVEZOKU3MtEfW0NLU2NdHQquYuIyMQouZeRTDKhlruIiEyYknsZSScTGlAnIiITpuReRjKpRrXcRURkwpTcy0gmmeCVDX309A3EHYqIiExiSu5lJJ3UXHcREZk4JfcykkkNzXVXchcRkfFTci8jKkErIiKFoOReRmaoW15ERApAyb2MJBN1TKmvVQlaERGZECX3MmJmZFKa6y4iIhOj5F5m0knNdRcRkYlRci8zmWRCo+VFRGRClNzLTDqVYEVXD4ODHncoIiIySSm5l5lMMkHfgLN2Q2/coYiIyCSl5F5mNNddREQmSsm9zKRTmusuIiITo+ReZja33DXXXURExknJvcy0Njdipm55EREZPyX3MlNfW0NLUyMdmg4nIiLjpORehjLJhFruIiIybkruZSidVAlaEREZPyX3MpRJqQStiIiMn5J7GcokE7yyoY+evoG4QxERkUlIyb0MpbWuu4iITICSexnKpIbmuiu5i4jI2JU0uZvZ6Wb2oJmtM7OVZnaTmc3NOl5vZueZ2WNmtt7MlpnZNWa2cynjjJtK0IqIyESUuuV+GHAZcABwBNAP3GFm24bHpwL7AeeGX98D7ATcYmZ1JY41NipBKyIiE1HShOnu78zeNrMTgE7gQOAmd+8Ejsw55zPA48BewIIShRqr5sY6pjbUqgStiIiMS9z33JvDGNbmOScZfs13TkUxMzKa6y4iIuMUd1f3JcCjwH3DHTSzBuACglb94tzjK1euZN68eZu329raaGtrK1KopTUjqbnuIiIyPrEldzO7EDgIOMjdt5rQHd5jvxqYDrx7uNdobW2lvb29qHHGJZNM8OALVdNZISIiBRRLcjezi4APAYe7+/PDHK8DrgX2Bg5z99UlDjF26VSCFV09DA46NTUWdzgiIjKJlPyeu5ldAnwYOMLdnxrmeD3wG2AfguS/vMQhloVMMkHfgLNmQ2/coYiIyCRT0pa7mf0IOAE4DlhrZpnwULe7d4ct9vnAm4BjAc86p9PdN5Yy3jhtnuve2UNLU2PM0YiIyGRS6pb7ZwlGyP8FWJb1ODU8viPB3PYdgIdyzvlgiWONlea6i4jIeJV6nnvem8fu/gKgG8y82nLvWKe57iIiMjZxz3OXEbQ2N2KmErQiIjJ2Su5lqr62hpamRjq0eIyIiIyRknsZyyQTarmLiMiYjXjP3cxeGsfrOXC0uy8cf0gyJJ1MsHjthrjDEBGRSSbfgLodgZuBlRFfqwb4KNAw0aAkkEk10v7imrjDEBGRSWa00fJnu/sDUV4onKN+wsRDkiGZZIJXNvTR0zdAor427nBERGSSyHfP/ZvAy1FfyN37w+csmWhQEkgnNdddRETGbsTk7u7nuvuysbxY+JyOiYclAJnUq1XqREREohoxuZvZL8zskFIGI1vaXIJWLXcRERmDfN3yHwTuNLNFZnaWme1eqqAkoBK0IiIyHvmSexr4FPAC8A3gaTO7x8w+bWapUgRX7Zob65jaUMvyTpWgFRGR6PLdc+9291+6++HALILBctsCPwWWmdm1ZnaUmakQTpGYGZlkQi13EREZk0iJ2d1fdvfvuPtrgbcAvwDeBvwRWGJm5xcxxqqWVpU6EREZozG3ut39AXc/GZgJXATMAL5c6MAkkEklNFpeRETGZMxLvoYD6z5GUI1uF6ALmF/guCSUTiZY0dXD4KBTU6PVcEVEZHSRkruZbUswev5jwJsJasjfAXwd+L27q2lZJOlkI30DzpoNvbQ0NcYdjoiITAL5Fo6pB44hSOhHEdSMfwI4Hbja3ZeWJMIqt3mue2ePkruIiESSr+W+HJgOrAF+Blzp7u0liUo2y57rPnemZiCKiMjo8iX3vwFXAn90974SxSM5VKVORETGasTk7u7HDbffzJqA7YClSvrF19rciBl0aMS8iIhEFHkqnJkdY2YPA53Av4C9w/0/N7MPFym+qldfW0NLU6Na7iIiElmk5G5mxwE3AKuAr+Y8bxFwYuFDkyFBlTqVoBURkWiittzPAH7p7u8ALs45thCYW9CoZAtplaAVEZExiJrc9wJ+E37vOcfWEtyDlyLJpNQtLyIi0UVN7uuAlhGOzQJWFiQaGVYmmeCVDX309A3EHYqIiEwCUZP77cDpZjY9a5+bWSNwMvDngkcmm6WTWtddRESii5rcvw5kgKeBnxN0zZ8GPArsCJxZjOAkkEm9WqVORERkNFGXfH0B2I9gidcjgQHgEOAfwP4qRVtcKmQjIiJjEXlVOHdfDHyyiLHICLJL0IqIiIxmzOu5S+k1N9YxtaGW5Z2a6y4iIqOL3HI3s3cBxwM7AYmcw+7uhxYyMHmVmYWFbNRyFxGR0UWtUPcV4GaCJWCnEdxzz34MFitACaSTCd1zFxGRSKK23E8Gfgqc7O6abB2DTCrBA4vWxB2GiIhMAlHvuSeB+Urs8UknE6zo6mFwMLdAoIiIyJaiJvdbgbcUMxDJL5NspG/AWbOhN+5QRESkzI2lW/73ZubAbQT15Lfg7s8XMjDZUnYhm5amxpijERGRchY1uTvQBZwLnDPCObUFiUiGNSOrBO3cmamYoxERkXIWNblfARwAXAQ8BYyrb9jMTgfeB8wBNhFUuDvd3RdmnWMES8y2AdsA9wOfc/fHx3PNSqEqdSIiElXU5H4YwUj5KyZ4vcOAy4AHAQPOBu4ws9e6+9BQ8K8ApwAnEdSy/xZwu5nNcfeuCV5/0mptbsQMOlRfXkRERhE1ua8COiZ6MXd/Z/a2mZ0AdAIHAjeFrfYvAd9z99+F55wIrAA+TDAdryrV19bQ0qR13UVEZHRRR8v/N/BZMyt0udrmMIahAXq7Eqw+d9vQCe6+EfgrwW2BqpZJJli+TiVoRUQkv6gt922AucATZnY7W4+Wd3c/YxzXv4Rg2dj7wu1M+DW3l6ADmJn75JUrVzJv3rzN221tbbS1tY0jjMkhnUzw8poNcYchIiJlLmpy/3rW97OHOe4Eg+AiM7MLgYOAg4YpjpNbqcWG2Udrayvt7e1jueyklkk18uALqlInIiL5RV3PvWaUx5imwZnZRcC/A0fkzI9fHn7N5DxlBgW45z/ZZZIJOjf20dOnQoEiIjKyki/5amaXEAyOO8Ldn8o5vIggwR+ZdX4COBi4t2RBlql0Uuu6i4jI6Eqa3M3sR8DHCVrta80sEz6aILhxD1wMnGZm7zOzuQRz7LuBa0oZaznKrlInIiIykhGTu5kNmNmbo76QmdWGz9kvz2mfJRgh/xdgWdbj1Kxzvg9cCPwIaAe2B95RzXPch6iQjYiIRJFvQJ0BGTPbeQyvZflOcPe8x8NzHDgzfEiWdErd8iIiMrrRRsv/foyvp/VIi6i5sY6pDbUs79RcdxERGVm+5P7xcb7monE+T0ZhZmSSCbXcRUQkrxGTu7tfWcpAJJp0MqF77iIiklfJp8LJxGRSCY2WFxGRvJTcJ5l0MsGKrh4GBzW8QUREhqfkPslkko30DThrNvTGHYqIiJQpJfdJRoVsRERkNEruk4xK0IqIyGiU3CeZzS13JXcRERnBiFPhzOwXY3gdd/dPFiAeGUVrUyNm0KFueRERGUG+IjZHEL3inIZul0hdbQ0tTY1quYuIyIjyFbGZVcI4ZAwyyQTL16kErYiIDE/33CehdDKhbnkRERnRaAvHbMXMZgCJ3P3u/lJBIpJRZVKNPPjCmrjDEBGRMhUpuZtZDXAO8Blg+gin1RYqKMkvk0zQubGPnr4BEvX6sYuIyJaidst/CfgccAHBmu3fIUj2i4B/AZ8uSnQyrKG57ipkIyIiw4ma3D8OnA2cF27/3t3PAPYClgA7FyE2GYHmuouISD5Rk/tuQLu7DwD9wBQAd+8DLgY+UZzwZDgZVakTEZE8oib3Tl4dRLcUmJN1rA7YtpBBSX7plJK7iIiMLOpo+UeA1wK3ho+zzGwjQSv+XODh4oQnw2lurGNqQy3LOzXXXUREthY1uV9M0DUPcAawH/CrcPtF4OQCxyV5mBmZZEItdxERGVak5O7ut2d9v9zM3gy8BpgKPBnee5cSSicTGlAnIiLDGnMRGwhWiQGeK3AsMgaZVIIHFqmQjYiIbC3fqnCHAA+7e3f4fV7u/teCRiZ5pZMJVnT1MDjo1NRY3OGIiEgZyddyvwt4C/BA+P1IK79ZeEyl0kook2ykb8BZs6GXlqbGuMMREZEyki+5Hw48EX4/luVfpQQ2F7Lp7FFyFxGRLeRb8vXurO/vKkk0Elk6q5DN3JmpmKMREZFyEqmIjZk9b2avH+HYXDN7vrBhyWhUglZEREYStULdLGCkvt8EsEtBopHIWpsaqTG0rruIiGwlanKHke+5zwNeKUAsMgZ1tTW0NDWq5S4iIlvJNxXuy8CXw00HbjKz3pzTphDUlf91ccKTfDKpBMvXqQStiIhsKd9o+eeBv4Tfnwi0AytzztlEMKL+54UPTUaTTiZ4afWGuMMQEZEyk2+0/A3ADRDUMgfOdvdFJYpLIkgnG1WlTkREthK1tvzHix2IjF0mmaBzYx89fQMk6lVDSEREApFry5vZbsAHgJ15dW33Ie7unyxkYDK6obnuyzt7mNUyLeZoRESkXERK7mb2HmA+wej6FQT32rOpel0Msue6K7mLiMiQqC33cwjqy3/E3XMH1UlMMllV6kRERIZEnee+G3B+IRK7mR1iZjea2RIzczM7Ked4k5ldamaLzWyjmT0dTsuTHOms+vIiIiJDorbcnwK2K9A1m4CFwFXhI9eFwNuBE4BFwCHAz8xslbv/b4FiqAjNjXVMbailQ3PdRUQkS9SW+1eAr4WD6ibE3W9296+5+3XA4DCnHAD8r7vf6e4vuPtVwD+A/Sd67UpjZmSSCXXLi4jIFqK23M8kaLk/aWbPArmTq93dDy1QTPcAx5rZz939ZTM7ANgX+EGBXr+ipJMJlaAVEZEtRE3uA8DTxQwkyxeAnwAvmVl/uO/z7v7H3BNXrlzJvHnzNm+3tbXR1tZWmijLRCaVUCEbERHZQtQiNocVOY5snwcOBN4NvEhwz/18M3vB3W/JPrG1tZX29vYShlZ+0skEK7p6GBx0amos7nBERKQMRC5iUwpmNgX4LnC8u98U7n7MzPYFTgVuGfHJVSqTbKRvwFmzoZeWppFW5RURkWoSeclXM5tpZheaWbuZLTKzueH+L5lZoQa71YePgZz9A2OJtZpkNB1ORERyRK1Q9zrgbwRJ9j7gDUBDeHgX4M3AhyO+VhOwe7hZA+wctszXuPtLZnY38D0z6ybolj8U+BjBiH3Jkc4qZDN3ZirmaEREpBxEbQ1fADwJ7Aq8D8i+uXsv8JYxXHMe8Ej4mAKcFX5/dnj8Q8CDwK8IlpM9Dfgm8MMxXKNqZJegFRERgej33A8C/t3du80sd/mxDiAT9YLufhdbfjjIPb4c0Cp0EbU2NVJj0KFueRERCUVtuQ9XbGZIC7CxALHIONTV1tDS1KiWu4iIbBY1uT/AyK3pDwB/L0w4Mh6ZVILlKkErIiKhqN3y3wbuMLPbgGsIlnh9u5l9EXgvwVx0iUk6meCl1RviDkNERMpEpJa7u98NHEcwoO4XBPfMvwccDBzn7vcXLUIZVUYlaEVEJEvkIjbu/ifgT2a2OzADWO3upSpJK3lkUgk6N/bR0zdAoj53vKOIiFSbSC13M/uWme0A4O7Pufu9Q4ndzLY3s28VM0jJb0ZzUJlOhWxERASiD6g7A9hxhGM7hMclJprrLiIi2aIm93wrkmwDaKh2jDJZVepERERGvOduZocBR2Tt+oyZHZNz2hTgaODxwocmUaVVX15ERLLkG1B3KPCN8Htn+HnuvQQlYr9Q4LhkDJob65jaUKtueRERAfJ0y7v7We5e4+41BN3ybxnaznok3H0/d7+vdCFLLjMjk0yoW15ERICIU+HCBC9lLJ1MqFteRESAMcxzN7OpwCcIuuu3BVYDdwFXuLvKo8Usk0rwwKI1cYchIiJlIOo89wzwMPDfBEu2TgXeRLAM60Nmli5ahBJJOplgRVcPg4MedygiIhKzqN3t3yeY8nawu+/q7m91910JloKdDpxXrAAlmkyykb4BZ82G3rhDERGRmEVN7kcBp7v7Fqu/ufu9BCPqjy50YDI2GU2HExGRUNTk3gQsHeHY4vC4xCitQjYiIhKKmtyfBk4Y4dhHgacKE46Ml0rQiojIkKij5c8HrgoHzl0DLAMywIeAtzNy4pcSaW1qpMagQ93yIiJVL+o896vDqXBnAz/POtQB/Ie7X1OM4CS6utoaWpoa1XIXEZExred+uZn9HJhDMM99DfC0uw8WKzgZm0wqwfJ1WsNHRKTaRU7uAGEif7JIscgEpZMJXlqtekIiItUucllZM9vDzK40s2fMbH349Qoz272YAUp0mWRC3fIiIhKt5R4u/3ozsBH4E8G99jRwLPBBM3uXu99drCAlmkwqQefGPnr6BkjU18YdjoiIxCRqt/wFwCPAO929e2inmTUDt4XH5xU+PBmLobnuyzt7mNUyLeZoREQkLlG75V8LnJed2AHcvYug9OzrCh2YjF0mqbnuIiISPbkvBhpGONYALClMODIRmVQjoCp1IiLVLmpyPw84y8xmZu8Mt88AvlPowGTsZiRVX15ERKLfcz8UaAb+ZWb/4NUBdW8Jvz8sHHQH4O5+YqEDldE1N9YxtaFW3fIiIlUuanI/CBggKDu7S/gg3AY4OOtcLSgeEzMjk0yoW15EpMpFLT+7a7EDkcJIJxPqlhcRqXKRi9jI5JBJJehQCVoRkaqm5F5h0mG3/OCg7o6IiFQrJfcKk0k20j/orF7fG3coIiISEyX3CpNJBdPhNKhORKR6KblXmKEStEruIiLVq+TJ3cwOMbMbzWyJmbmZnTTMObPN7Hoze8XMNpjZw2a2V6ljnYyGWu6a6y4iUr3iaLk3AQuBLxKsMrcFM9sV+DuwCDgCmAt8A+jOPVe21trUSI1Bh6bDiYhUrRHnuZvZINEL0ri7R50zfzPB8rGY2RXDnHIucJu7n5K17/mIcVS9utoaWpoa1XIXEali+RLy2ZS42pyZ1RCsEf89M7sFeCPwAnC+u/+mlLFMZplUguWa6y4iUrVGTO7ufmYJ4xgyg6Db/mvAN4HTCLrmf2Vm6939jzHENOmkkwleWr0h7jBERCQmUWvLl8rQGIAb3P3C8PtHzWwe8Dlgi+S+cuVK5s2bt3m7ra2Ntra2kgRazjLJBA8sWhN3GCIiEpPIyd3MGoCjgDlAIuewu/u3CxDPKqAfeCJn/5PAh3JPbm1tpb29vQCXrSyZVILOjX309A2QqK+NOxwRESmxSMndzHYA7gFmEdyHt/BQ9j35CSd3d+81swcJPkBkmw28ONHXrxbprHXdZ7VMizkaEREptahT4X4ArAR2Jkjs+wO7EYxsfy78PhIzazKzfc1s3/D6O4fbO4enfB/4oJm1mdnuZvZpglb7j6Jeo9plkprrLiJSzaIm94OBC4Cl4fagu7/g7t8CrgP+ewzXnAc8Ej6mAGeF358N4O5/ANqAU4EFwOeBj7n7n8ZwjaqWSTUCqlInIlKtot5z3w5Y6u6DZrYe2Cbr2P8BJ0e9oLvfxavd+iOdcwVwRdTXlC1ld8uLiEj1idpyXwy0hN//C3hH1rE3A8oiZaQ5Uc+0hlp1y4uIVKmoLfc7gUOBPwA/BX4U3jPvA94Z7pMyMrSuu4iIVJ+oyf0bwLYA7v5jM6sDPghMJRgAd3ZxwpPxSicT6pYXEalSUevBryKYgz60fSlwabGCkonLpFTIRkSkWkW6525mrWY2e4Rjs82sZbhjEp+hbvnBwZIuDyAiImUg6oC6y4BTRjj25fC4lJFMspH+QWf1+t64QxERkRKLmtwPAm4d4dhtwIGFCUcKJZMKpsNpUJ2ISPWJmty3ATpHOLaOYB68lBHNdRcRqV5jmee+/wjH9geWFSYcKZShlrvmuouIVJ+oyf064GtmdnT2znD7NOC3hQ5MJqa1qZEagxVK7iIiVSfqPPezgUOAG81sObAEmAlkgH8Q1IeXMlJXW0NLU6Na7iIiVSjqPPcNZnYocAJwJME99ucIBtNd7e79xQtRxiuTSrB83aa4wxARkRKL2nLH3fuAX4QPmQTSyQQvrd4QdxgiIlJiUe+5yySUSSbULS8iUoVGbLmb2fPAe939n2a2CMhX6szd/TUFj04mJJNK0Lmxj56+ARL1tXGHIyIiJZKvW/5ugjnsQ9+rjukkkz3XfVbLtJijERGRUhkxubv7x7O+P6kk0UhBZZKvznVXchcRqR66517BMqlGQCVoRUSqTeTR8ma2LXA0sBOQyDns7n5GIQOTiVMJWhGR6hQpuZvZO4DfASP17Tqg5F5mmhP1TGuo1Yh5EZEqE7Vb/kLgEeD1QKO71+Q8NBS7TKVTCXXLi4hUmajd8rOAL7v7giLGIkWQSSbULS8iUmWittwfAXYoZiBSHJlkgg6VoBURqSpRk/t/Al8xs7cWMxgpvBnJoFt+cFBlCkREqkXUbvmHgL8A95jZeuCVnOPu7rsUNDIpiEyykf5BZ/X6XlqbG+MOR0RESiBqcj8fOJmge/4poLdoEUlBZVLBdLiOdT1K7iIiVSJqcj8J+Lbmsk8+2XPd585MxRyNiIiUQtR77g78tZiBSHEMtdw1111EpHpETe7zgaOKGYgUR2tTIzWmErQiItUkarf8n4GLzCwF3AKszT3B3f+vkIFJYdTV1tDS1Ki57iIiVSRqcv99+PWT4WOIAxZ+VZW6MpVJJdQtLyJSRaIm98OLGoUUVTqZ4KXVG+IOQ0RESiRScnf3u4sdiBRPJpnggUVr4g5DRERKJPKSrwBm1gK8BdgOuMnd15hZAuh198FiBCgTl0kl6NzYR0/fAIl63T0REal0kUbLW+AHwGLgRuAXBIvJANwAfL0o0UlBaF13EZHqEnUq3OkEFerOBvYnGEQ35CbgmALHJQWUSWquu4hINYnaLf8p4Gx3/66Z5fbrPge8prBhSSFlUkHZWc11FxGpDlFb7jOBf4xwrBeYFvWCZnaImd1oZkvMzM3spDznXh6ec2rU15etqVteRKS6RE3uS4C5Ixx7PbBoDNdsAhYCXwQ2jnSSmb0feBOwdAyvLcNoTtQzraFW3fIiIlViLOVnv2VmB2btczObDZwC/DrqBd39Znf/mrtfBww7wt7MdgEuAT4M9EV9bRlZOpVQt7yISJWImtzPJFjq9a/As+G++cCCcPt7hQrIzOqAa4Fz3P3JQr1utcskE+qWFxGpEpGSu7tvBA4jWPr1XuAO4EGgDTjS3Qu5vvtZwGp3/3EBX7PqZZIJOtZtijsMEREpgchFbNx9APjf8FEUZnYowQeIfaOcv3LlSubNm7d5u62tjba2tuIEN8kNdcsPDjo1NTb6E0REZNIaU4W6Ejgc2B5YZrY5AdUC55nZl9x9x+yTW1tbaW9vL3GIk1MmmaB/0Fm9vpfW5sa4wxERkSKKWqGuwczOMLOnzGyDmQ3kPPoLFM9lwD4ELfehx1LgIuBtBbpGVRqaDqdBdSIilS9qy/0HwOcI1nW/Hhj3zVszawJ2DzdrgJ3NbF9gjbu/BKzIOb8PWO7uT4/3mgLpZNBaX97Zw9yZqZijERGRYoqa3N8PnOHu5xbgmvOAO7O2zwofVxLcb5ciyKRUglZEpFpETe5NwH2FuKC738WWtelHO39WIa5b7VqbGqkxdcuLiFSDqPPcbwIOKWYgUlx1tTW0NDVqrruISBUYseVuZrtlbV4KXGVmg8DNwJrc8939+cKHJ4WUSSXULS8iUgXydcs/B3jWthFUqjtjhPNzV4uTMpNOJnhx9fq4wxARkSLLl9w/XrIopCQyyQT3P7867jBERKTIRkzu7n5lKQOR4sukEqzr6Wdj7wBTGtTRIiJSqUYcUGeBY81spKXbdrKOAAAZCUlEQVReMbO9zezY4oQmhaZCNiIi1SHfaPkTCFZny3eTtgu41sz+vaBRSVFkkprrLiJSDfIl948Cv3T3RSOd4O4vAP8DnFjguKQIMqmgSp1a7iIilS1fct8PuC3Ca9xBUHVOytxQt7zmuouIVLZ8yb0ZWBvhNdaG50qZa07UM62hVt3yIiIVLl9yXwXsEuE1dg7PlUlgaF13ERGpXPmS+z1Eu5d+UniuTAKZZELd8iIiFS5fcr8YeJuZXWRmDbkHzazezC4BjiBYb10mgUwyQce6ca/YKyIik0C+Ijb3mdkpwAXAR8zsNuDF8PAuwJHAdsAp7v6PokcqBTHULT846NTURF6cT0REJpG8S766+8Vm9jBwGvBeYEp4aCNwF/A9d/9bUSOUgsokE/QPOqvX99La3Bh3OCIiUgSjrufu7n8F/mpmNUBLuHu1uw8UNTIpiuwqdUruIiKVadTkPsTdB4EVRYxFSiCTenWu+9yZqZijEREpvnU9ffhgfNef0lBLQ12+IW6FFzm5S2VQCVoRqQZr1/dyw6NLmP/QYh5fui7WWC4/4Y2843WZkl5Tyb3KtDQ1UGMqQSsilWdg0Pnrsyu5rn0xtz/RQe/AIK/bIcmp75jN1Ib40t2cTOnrvCm5V5m62hpamho1111EKsbzK7uZ/9Birn94MR3rNrHN1Ho+vP/OHD9vR163Q3XeflRyr0KZVELd8iIyqXVv6udPjy1lfvti2l9cS43BYXNmcOaxO3LEXjNorKuNO8RYKblXoXQywYur863kKyJSftyd+xetYX77Ym5esIyNfQPs1jqNr75rT96338zNs4FEyb0qZZIJ7n9+ddxhiIhEsuSVjfzuocVc99BiXlqzgabGOo57ww68/407sd/O0zFTQa5cSu5VKJNKsK6nn429A0xpqO6uKxEpTz19A9z6+HKue2gx9zy3Cnd4627b8aW378G75mZiHSA3GeinU4XSWdPhdm2ZFnM0IiIBd+exxZ38tv1lbvznUrp6+pk5fQqfP2IPjn/jjuy07dS4Q5w0lNyr0Oa57p1K7iISv5Vdm/jDI0uY/9DLPNPRTWNdDUfNzXD8vJ14627baR2McVByr0KZVFB2dkWXRsyLSDz6Bga586kVzH9oMXc+tYL+QWffnaZz7nvncsw+O5CaUh93iJOaknsVSme13EVESunp5V3Mb3+ZPzy6hFXdvbQ0NfKJg3bl+DfuyB7p0hd7qVRK7lWoOVHPtIZazXUXqRKDg85Lazbwysa+WK7v7ixcuo7r2l/mn4s7qasx3rbXDI5/404cOqeV+trS1l2vBkruVWpoXXcRqSyDg86LazawYEknC5d08tjiV3h8yTq6NvXHHRp7Zpr5xtF7cdwbZtLSpFUpi0nJvUplkgl1y4tMcu7Oi6uzE3knC5d20tUTJPKGuhr22j7Je96wA3vPTDGjOb4iL+lkgr22b9ac9BJRcq9SmWSC+xetiTsMEYnIPehaX7CkkwWLOzcn9HVDiby2hr22b+bdrw8S+d47ppidblaXd5VScq9SQ93yg4OuaSYiZcbdeXnNRhYs6eSxJa+wMEzo2Yl8z+2bOeb1O7DPzBRzZwaJvNRrhkv5UnKvUplkgv5BZ/X6Xlqbde9LJC7uzuK1YSJfHLTGFyzppDMc/FZfa+yZSXL0Pjuwz44p9lYilwiU3KvU0HS4jnU9Su4iJZKdyIe61Rcs6eSVDa8m8jmZZv5t7wx7z5weJPJMU9WvcCZjp+RepTKpV+e6z51ZnesdS2Vzd5Z29vBMRxfPdnTx9PJunl3RxQur1jPo8cTUPzhIT98gAHU1QSI/am6GuTODFvmcTLMSuRSEknuVymTVlxeZzNydFV2beKaji6eXd/FsRzfPrAi+dmdN/5rR3MjsdDPv2XdmbF3aBsxqmcY+OyqRS3GVPLmb2SHAqcAbgR2Aj7v7FeGxeuAc4CjgNcA64E7gNHd/qdSxVrKWpgZqDM11l0llVfcmnlnexTMdXTyzojtskXdtHmgGsN20BvZIN/H/9pvJHulmZqebmZ1uYvrUhhgjFymtOFruTcBC4KrwkW0qsB9wLvAokAIuAG4xs33cPf4qDBWirraG1ubGqp7r7u78c3En62Kq2gXQWFfDtMY6mhrrmNZYR3Oijsa6mqqfC7x2fe9WCfzZFd2sWd+7+ZzUlHrmpJs59vU7hAk8SOLbqTiKSOmTu7vfDNwMYGZX5BzrBI7M3mdmnwEeB/YCFpQmyuqQSSaqtlt+Zdcmvvb7Bdz+REfcoWyltsaY1lC7OeE3JcLk3xBuN9bSlKh79UNBw6sfDIaOTwufO62hjtoynuq4rqePZzu6eKajO0zgwfcruzZtPqe5sY490k2847Vp9kg3MydM4q3NjVX/IUhkJJPhnnsy/Lo21igq0IxkghdXr487jJJyd/742DK+dcNC1vcOcNpRe/KmWdvEFAts6h+ke1M/6zf10x0+1m/qZ/2mgWC7p5/1vcH+jnU9rN80QFdPH+t7BxiIOCpsakNtVu9ALdMagg8CDXU1GPEkx3U9fTy3optlWT1HUxtq2WNGE4fObmVOupk90k3MTjezfSqhJC4yRmWd3M2sgaBb/iZ3X5x7fOXKlcybN2/zdltbG21tbSWMcHLLJBPc//zquMMomdXdm/jmDQu5ecFyXr/TdC44/vXsPqMp7rDGxd23+GDQ1RN+KOjtp3vTQPBhoSfrw0K4v7unj/WbBlj6Sg+9A4OxxT+lvpa37rZdeE88SOIzp09RQSWRAinb5G5mdcDVwHTg3cOd09raSnt7e0njqiSZVIJ1Pf1s7B1gSkNlj9r984JlfOMPC+nq6ecr75pD28G7UTeJy3KaGYn6WhL1tVqAQ0S2UpbJPUzs1wJ7A4e5e/U0L0sonTUdbteWaTFHUxxr1/dyxo2Pc+M/l7L3zBTnH/965mS0ZrSIVLayS+7hdLhfA3MJEvvymEOqWJvnundWZnK//YkOTr9+AZ0beznlyNn8x2Gv0SIaIlIV4pjn3gTsHm7WADub2b7AGmApMB94E3As4GaWCc/tdPeNpY63kmVSQXdupc1179zQx1k3Pc71jyxhr+2TXPWJN/PaHZKjP1FEpELE0XKfR1CYZshZ4eNK4EzgPeH+h3Ke93HgiiLHVlXSFVil7s6nVnDa9Y+xqruXL7xtD04+fHctsCEiVSeOee53Qd75NxouWyLNiXqmNdRWRCGbdT19nPPHJ/ht+2LmpJv5+cfexN47qma+iFSnsrvnLqWVTiVY0TW5k/tfn1nJV3/3GB3revjsYa/hi2/fQzW7RaSqKblXuUwyMWlb7t2b+vnOzU9yzf0v8ZrWaVz/2QPZd6fpcYclIhI7Jfcql0kmuH/RmrjDGLN7n1vFf133GEs7N/KZQ3bjy0fOJlGv1rqICCi5V710KkHHuh4GB31SVAdbv6mf8255iqvue5FdW6Zx3X+8lTfusm3cYYmIlBUl9yqXSSboH3RWr++ltbm8K53d//xq/uu6x3h57QY+ceCu/Nc751R8ZT0RkfFQcq9yQ9PhOtb1lG1y39g7wPdvfYor7n2BnbaZyq8//Rb23227uMMSESlbSu5VLpN6tUrd3JnlN3XsoRfXcOr8x1i0aj0nvnUXvnrUnkxt0D9bEZF89FeyymXKtJBNT98AF97+DD/72/PskJrCNZ/anwN2b4k7LBGRSUHJvcq1NDVQY+VVgvaRl9Zy6vx/8q+V6/nw/jvztX/bi6ZG/VMVEYlKfzGrXF1tDa3NjWUx131T/wCX3PEsP7n7X2SSCa76xJs5ZHZr3GGJiEw6Su4SFLKJueW+YHEnp8x/lGc6uvnAvB35xjGvJZmojzUmEZHJSsldSCcTPLl8HXc/szKW6z+4aA0/vvtftDQ18MuT3sThe86IJQ4RkUqh5C7s2jKN257o4MRfPBBbDO/bbyZnHPM6UlPVWhcRmSgld+HLR87mnXMzuMdz/WSijj3SzfFcXESkAim5C4n6WvbbeZu4wxARkQKpiTsAERERKSwldxERkQqj5C4iIlJhlNxFREQqjJK7iIhIhVFyFxERqTBK7iIiIhVGyV1ERKTCKLmLiIhUGCX3MnL55ZfHHUIsqvV9Q/W+92p931C9771a3zfE896V3MtItf7jr9b3DdX73qv1fUP1vvdqfd+g5C4iIiIFYB7XUmAFYGYrgRfjjqOAWoBVcQcRg2p931C9771a3zdU73uv1vcNhX/vu7h7a74TJnVyFxERka2pW15ERKTCKLmLiIhUGCV3ERGRCqPkHjMzO93MHjSzdWa20sxuMrO5ccdVamb2NTNzM/th3LEUm5ltb2ZXhr/vHjN7wswOjTuuYjOzWjP7tpktCt/3IjM7x8zq4o6tkMzsEDO70cyWhP+mT8o5bmZ2ppktNbONZnaXmb0upnALKt97N7N6MzvPzB4zs/VmtszMrjGznWMMuSBG+53nnHt5eM6pxYxJyT1+hwGXAQcARwD9wB1mtm2cQZWSmb0F+DTwWNyxFJuZTQf+DhhwNLAX8HlgRZxxlchXgc8BXwD2BL4Ybp8eZ1BF0AQsJHh/G4c5/hXgFILf+5sIfve3m1lzySIsnnzvfSqwH3Bu+PU9wE7ALRXwAW+03zkAZvZ+gt/50mIHpNHyZcbMmoBO4Dh3vynueIrNzFLAwwTJ/VvAQnc/Od6oisfMvgMc6u4Hxh1LqZnZH4HV7n5i1r4rge3c/Zj4IiseM+sGTnb3K8JtI/jD/kN3PzfcN4UgwZ/q7j+NK9ZCy33vI5zzWuBxYB93X1Cq2IpppPdtZrsA9wJvB/5M8G/g/GLFoZZ7+Wkm+L2sjTuQErkcuM7d/y/uQErkOOB+M/uNma0ws0fN7OTwj36luwc43Mz2hM1/2I8Abo41qtLaFcgAtw3tcPeNwF8Jeu+qTTL8WtF/78KeiWuBc9z9yVJcc7J3hVSiS4BHgfviDqTYzOzTwO7ACXHHUkK7AZ8FLgK+B+wLXBoeq/TxBucRfHh9wswGCP7+nOvul8UbVkllwq8dOfs7gJkljiVWZtYAXADc5O6L446nyM4i6LX6cakuqOReRszsQuAg4CB3H4g7nmIysznAd4CD3b037nhKqAZod/eh+8yPmNkeBPeeKz25fxD4GPBhgq7YfYFLzGyRu/9PrJGVXu79UBtmX8UKW7JXA9OBd8ccTlGFg2VPIvj3XjLqli8TZnYR8O/AEe7+fNzxlMBbCUoyLjSzfjPrBw4FPhtuN8YbXtEsA57I2fckMOlHDEfwA+B8d/+1uy9w9/8FLqTyBtTlszz8msnZP4OtW/MVKauLeh/gbe6+OuaQiu1wYHtgWdbful2A88ysaD0WarmXATO7BPgQcJi7PxV3PCXyB6A9Z98vgWcJWvSV2pr/OzAnZ99sKmuNhJFMBXJ7pAaorkbGIoIEfyTwIICZJYCDgf+KMa6SMLN64NfAXIK/d8tHeUoluAy4LmffrQQfcH5WrIsqucfMzH5EcM/5OGCtmQ19ou929+74Iisud38FeCV7n5mtB9a4+8J4oiqJi4B7zezrwG+ANxBMDftarFGVxk3AaWa2iKBb/g3AfwJXxRpVgYUzXnYPN2uAnc1sX4J/2y+Z2cXA183sKeAZ4BtAN3BNLAEXUL73TjBLYD7BVLBjAc/6e9cZDiyclEb7nZMz1dXM+oDl7v500YJydz1ifBDcZxvucWbcscXws7iLYHpI7LEU+X0eDfwT6CH44/4FwmmplfwgGEx3MUEvxUbgeYJemkTcsRX4fR42wv/pK8LjBpxJcIumB7gbmBt33MV+78CsPH/vToo79mL+zoc5/wWCqY9Fi0nz3EVERCpMNd3rEhERqQpK7iIiIhVGyV1ERKTCKLmLiIhUGCV3ERGRCqPkLiIiUmGU3KXimNlJZuZm9oqZbZNzrC48dmYMcZ0ZXrusi0eZWY2ZXWxmy8xs0Mz+MMHXmxW+75MKFGL2a08Pf677Ffq1iy3r38O/wspt2cd2L9bPTKqDkrtUshTw1biDmITeD3yRoBb8gcBX4g0nr+nAGcCkS+5ZdgM+GXcQUlmU3KWS3QZ8PqvEZcUr0II7e4VfL3b3+9z9mQK8Zlko0wWJbgO+EdaYFykIJXepZOeEX7+e76Sh7tFh9l9hZi9kbQ91L/+HmX3XzJabWZeZXW1mU8Ou1FvNrNvMnjOzE0e45F5mdqeZbQi7vs82sy3+L5pZi5n92MyWmNkmM3vKzNpyzhm6/XCImc03s1eA+0d5r+8ys/vMbKOZdZrZH8Lld4eOv0BQGhVgYLSu4fA2x+lhfJvMbKmZXRAlUZnZoWb2l/BnuD782c0d5rz3mtnfw5/rOjN7wMzebWazCBZiAfhZGOvmeM3sLjO7x8yONbNHzGwT8NnwWNLMfhjGu8nMnjazL5uZZV33sPD13h2eu8rMVoa/7+k5MX7RzJ4Mf65rzazdzN472s8g9HVgB4Jlf0UKQsldKtkygjXS28xslwK+7ukEf4xPBL5FsE75T4DfA38C3gs8BvzSzF43zPP/ANxBsFjQNcA3w9cBgsRDsHrc0QSJ9miCRVd+bGafH+b1fkWQ5N4PnDZS0Gb2rjC+7jDm/49gda57zGxmeNp7CeqAQ7As71vD54zkaoKFT64J4/wuQRfzr/I8BzM7GvhLGMtHCdZ4bwb+ZmY7ZZ33eeB6goU3TgSOJ/g5zyL4/b4vPPW7I8Q7G/hv4FLgncBfwg9SfwI+DlxAsIjJLQTLz547TLiXENQJ/zBwNvD/wn1DMX4kfJ1rgX8DPkKwCti2+X4GWR4lWFDlNDNrjvgckfziLrivhx6FfgAnEfwx3p3gD+wrwC/CY3XkLMxDkEB9mNe5Angha3tW+Nz/yznv+nD/R7P2bQP0A2fkXgc4Lef5PwO6gOnh9jcJFhTZY5jzVgF1Oe/zoog/l3aCJXXrsvbtCvQBF2btO2e4n8cwr3dweP2P5ez/SLh/35yf20lZ5zwH/CXnecnw/V2ctd0FXJ8nhqHX/tQwx+4CBofiyNp/TG484f6fA5uAlnD7sPC8K3PO+2H4+7Gs7YfH8e906N9DHcGHkH7gW+Gx3YeLUQ89oj7UcpeK5u5rCFpVH8vufp6gP+dsPxV+vTXrumsJWps7sbXf5mz/GmgiaEUDvIuge31R2O1dZ8EI+1uB7YDX5jz/96MFbGbTCAad/cbd+7PiXETQS3DoaK8xjHcBvcDvcuK8LTx+yAix7AG8BvhVzvM2APdlPe8Agp/L5eOIbcgL7v5ozr5DCJL+tTn7rwYaCFr/2XJ7LhYAjUA63H4Q2NfMLjWzt5vZ1LEG6cG4hiuBU8wsaotfZERK7lINLiJYT/rsAr3e2pzt3jz7h7v33DHC9lDX+AyCBNSX85gfHt8u5/nLRg+ZbQiWGh3u3OVE70LONoMgGXbnxDm0dnVunNnPA/gftn6Px2Q9b+jr4nHENmS497stwTrbm3L2L886nm1NzvbQ84Z+t1cR3OLYn+AD2Bozuz4cEzAWZxF8aNAMD5mwsp5vK1II7t5tZt8laMH/YJhTegDMrMHde7P2j5ScJipNsJZ59jbAkvDraoIE+cURnv90znaUdZvXhucNN3MgE15zrFYT/OwOHuH40jzPg2Dswh3DHB/6HawKv84EFo4jPhj+Z7MG2HaY3/fQz2ZMPwt3d+CnwE8tqKvwDoJ/a78hSPhRX+clM/spcDLwx7HEIJJLLXepFpcRJM9zhjn2Yvh180jtcDT0AUWK5QM52x8iaP0OJbBbgD2Bl9y9fZhH11gv6O7rgYeA482sdmh/ONDwAODucbyPWwhar6kR4hwpuT8NvAC8boTnPRaedy/Bz6VthNeBV1vRU8YQ990Ef/uOz9n/EYIPFv8Yw2ttwd3XuvtvCG69bDXyP4JzCT6Q5J3hITIatdylKrj7JjM7m+Hv3/4Z6CSYTnUGQdfoVwgSSzF8Ohyx/SDBCO5PEQzweyU8fhHBaPa/mdlFBMlwGkHCP9jd3zPO636T4P7xH83sMoL72WcRvPcLxvpi7n6XmV0LXGdmFwIPENzLnkUwavyrPswceXd3M/sccIOZNRAkwlUEPRgHEHyoudDdu8zsdOBSM/sdwQj8LmBfoMfdLyW4pbEa+JCZPQasBxa5e77W95+Be4CfmFkr8HgY76eA77r7qjzP3YqZXR7GdR9Bj8ts4AReHXsQmbuvMLNLgK+N9bki2dRyl2ryS4LR4lsIk+oxBInptwTTqi4F7ixSHO8BjgRuJJgGdg7w7ax4OgmS3M0E919vBX4RPm/cMbn7LQTT1aYTvM+fAE8CB+VpZY/mowSjvt8P3EAwBexkgp9z7tiC7FhuJhhXMI1glPqtwPcJusbvyzrvhwQt7B0JkvvvwmstCo8PEiTlbQi6+B8kmNo2ovA5RxMMYPsqwQeeo4H/ZHwt5r8DbyToHbo9fI2rCabujccP2Hr8hsiYDE3lEBERkQqhlruIiEiFUXIXERGpMEruIiIiFUbJXUREpMIouYuIiFQYJXcREZEKo+QuIiJSYZTcRUREKoySu4iISIX5/wHO1S62Zhd0YwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1)\n",
    "fig.set_size_inches(8, 6)\n",
    "ax.tick_params(axis=\"x\",direction=\"in\",which='major',labelsize=14)\n",
    "ax.tick_params(axis=\"y\",direction=\"in\",which='major',labelsize=14)\n",
    "ax.set_xlabel('Number of electrons N',fontsize=16)\n",
    "ax.set_ylabel('Chemical potential [meV]',fontsize=16)\n",
    "\n",
    "ax.plot(electrons2[1:],11.61*mu_N_a)\n",
    "\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.03367193, 2.03367193, 0.96449271, 0.96449271, 1.04021135,\n",
       "       1.04021135, 1.03441212, 1.03441212, 1.05345166, 1.05345166,\n",
       "       1.07603814, 1.07603814, 1.11879428, 1.11879428])"
      ]
     },
     "execution_count": 202,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mu_N_a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}