examples/chemistry/qchem_modelling_basics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VKpLon6Agujb"
   },
   "source": [
    "# Quantum Chemistry: Modelling Basics\n",
    "\n",
    "The aim of this notebook is to provide the basic quantum chemistry knowledge to help solve a chemistry problem with quantum computers. Many parameters can be tweaked, and some have higher impact on the accuracy and the quantum resources (number of qubits, quantum circuit depth and number of measurements) than others. This set of parameters includes the choice of initial state, choice of basis set, and the freezing of molecular orbitals. Some references are provided when needed.\n",
    "\n",
    "The cell below installs Tangelo and PySCF in the current Python environment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    import tangelo\n",
    "except ImportError:\n",
    "    !pip install git+https://github.com/goodchemistryco/Tangelo.git@develop --quiet\n",
    "    !pip install pyscf --quiet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Et892XDtMfV6"
   },
   "source": [
    "## Table of Contents\n",
    "\n",
    "* [Introduction](#1)\n",
    "* [Computation of the Mean-Field](#2)\n",
    "    * [Molecular Information](#21)\n",
    "    * [Restricted vs Unrestricted  Hartree-Fock (RHF/ROHF vs UHF)](#22)\n",
    "    * [Choice of the Basis Set](#23)\n",
    "* [Consideration of Correlation Energy - Potential of Quantum Computing](#3)\n",
    "* [Resource Reduction](#4)\n",
    "    * [Active Space Selection](#41)\n",
    "    * [Considering Symmetries](#42)\n",
    "* [Advanced Usage of PySCF](#5)\n",
    "    * [Extracting Mean-Field Quantities](#51)\n",
    "    * [Outputting Cube Files](#52)\n",
    "* [Closing Words](#6)\n",
    "* [References](#7)\n",
    "\n",
    "## Introduction<a class=\"anchor\" id=\"1\"></a>\n",
    "\n",
    "At the submicroscopic level, experiments have shown that matter behaves according to the principles of quantum mechanics. Quantum effects, modelled via the Schrödinger equation, are essential for studying molecules *in silico*. Quantum chemistry aims at solving the time-independent version of this equation, shown below, for many-electron molecules.\n",
    "\n",
    "$$\n",
    "\\hat{H} \\Psi = E \\Psi\n",
    "$$\n",
    "\n",
    "Where $\\hat{H}$, $\\Psi$ and $E$ are respectively the Hamiltonian of the system, the wavefunction and the eigenvalue associated with $\\hat{H}$. \n",
    "\n",
    "As stated in the postulates of quantum mechanics, the wavefunction $\\Psi$ encodes all the information of the system. One who knows the electronic wavefunction of a molecular system can therefore predict measurement outcomes of an experiment. This is a highly valuable asset for academia and industrial purposes, to help develop new drugs or functional materials.\n",
    "\n",
    "The analysis of natural products is an efficient way of quickly getting inspiration for the design of new materials, while taking into account the million years of evolution that nature had to optimize a biochemical process. For instance, one difficult thing to do in the laboratory is to change the spin state of a compound during a chemical reaction. Nature circumvents this problem by leveraging transition metals to achieve this challenging task: the oxygen fixation process is achieved by the heme biomolecule, where a Fe(II) atom is involved in the spin state change mechanism (see figure below) [[10.1074/jbc.M314007200](https://doi.org/10.1074/jbc.M314007200)].\n",
    "\n",
    "![FeIIPorImO2 system](https://github.com/goodchemistryco/Tangelo-Examples/blob/main/examples/img/FeIIPorImO2.png?raw=1)\n",
    "\n",
    "Fully understanding the mechanism of oxygen fixation, that is to say knowing the rate constants for all elementary reactions, would lead to valuable insights for the design of chemical catalysts. Achieving this goal would take us one step closer to [*making the world cleaner, healthier, and more sustainable*](https://goodchemistry.com/). This is one of the main applications of quantum chemistry. \n",
    "\n",
    "However, the Schrödinger equation shown above cannot be solved exactly for systems beyond one electron, i.e. all relevant chemical systems for industrial applications. There are two equivalent ways of expressing the time-independant Schrödinger equation, whether the phase-space variables are trajectories (first quantization) or fields (second quantization).\n",
    "\n",
    "In the first quantization form, the electronic Hamiltonian $\\hat{H}$ for $N$ electrons is\n",
    "\n",
    "$$\n",
    "\\hat{H} = -\\frac{\\hbar^2}{2 m_e} \\sum_i \\nabla_i^2 -\\frac{e^2}{4 \\pi \\epsilon_0} \\sum_i \\sum_A \\frac{Z_A}{r_{iA}} + \\frac{e^2}{4 \\pi \\epsilon_0} \\sum_i \\sum_{j>i} \\frac{1}{r_{ij}}\n",
    "$$\n",
    "\n",
    "Where $m_e$, $\\nabla^2$, $e$, $Z_A$, $r_{iA}$ and $r_{ij}$ are respectively the electron mass, the Laplacian, the elementary charge, the relative charge of nucleus $A$, the distance between the nucleus $A$ and the electron $i$ and the distance between electron $i$ and $j$. \n",
    "\n",
    "In the second quantization form, the electronic Hamiltonian is expressed as\n",
    "\n",
    "$$\n",
    "H = \\sum_{pq}h_{pq} a_p^{\\dagger}a_q + \\frac{1}{2}\\sum_{pqrs}h_{pqrs}a_p^{\\dagger}a_r^{\\dagger}a_s a_q\n",
    "$$\n",
    "\n",
    "Where $h_{pq}$ and $h_{pqrs}$ refer respectively to the one-body and two-body integrals. \n",
    "\n",
    "For more information about the state-of-the-art theories and methods in quantum chemistry, the reader is invited to read the references at the end of this notebook. The scope of this notebook is to provide clues on how to choose input parameters for modelling molecules on quantum computers with Tangelo. As quantum resources are still scarce, the definition of the problem is a crucial step in an end-to-end quantum chemistry experiment on quantum hardware. To keep a reasonable runtime for this notebook, we work with a simple example, the dioxygen molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "o2_xyz = [(\"O\", (0., 0., 0.)), (\"O\", (0., 0., 1.21))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dfApY4WZMfV8"
   },
   "source": [
    "## Computation of the Mean-Field<a class=\"anchor\" id=\"2\"></a>\n",
    "\n",
    "The computation of the mean-field is usually the first step of all post Hartree-Fock approaches, which we describe in another section below. The mean-field approaches solve the Schrödinger equation by averaging the electron-electron repulsion term in a self-consistent way. In other words, each electron does not explicitely interact with other electrons. Depending on the molecular information, discussed shortly, different flavours of mean-field solutions can be computed.\n",
    "\n",
    "### Molecular Information<a class=\"anchor\" id=\"21\"></a>\n",
    "\n",
    "This notebook focuses on the second quantization formalism. In this case, chemical information is encapsulated in a `SecondQuantizedMolecule` object.\n",
    "\n",
    "When defining a molecule, one must specify molecular properties. One of them is the total charge of the system, i.e. the difference between the number of electrons and the total nuclear charge (default value is 0). This is an important property when dealing with oxydo-reduction reactions where cations (deficit in electrons) and anions (surplus in electrons) are involved. Another important quantity to specify is the spin, defined by the difference between the number of alpha and beta electrons (default value is 0). This allows radical systems to be modelled: a spin of 0 characterizes a restricted-shell molecule while a spin different from 0 defines an open-shell system.\n",
    "\n",
    "In the next cell, the dioxygen molecule data structure is constructed. The ground state of the dioxygen molecule is neutral and is a triplet state (q=0 and spin=2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8 active molecular orbitals\n",
      "12 active electrons\n"
     ]
    }
   ],
   "source": [
    "from tangelo import SecondQuantizedMolecule\n",
    "o2_sto3g = SecondQuantizedMolecule(o2_xyz, q=0, spin=2)\n",
    "print(f\"{o2_sto3g.n_active_mos} active molecular orbitals\")\n",
    "print(f\"{o2_sto3g.n_active_electrons} active electrons\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "X4IzRNI8MfV-"
   },
   "source": [
    "In the above output, there are missing orbitals and electrons. The default implementation of Tangelo freezes the core orbitals; for the O$_2$ molecule, two orbitals, corresponding to the 1p atomic orbitals, are frozen. The next sections elaborate on the impact of the basis sets and frozen orbitals.\n",
    "\n",
    "### Restricted vs Unrestricted  Hartree-Fock (RHF/ROHF vs UHF)<a class=\"anchor\" id=\"22\"></a>\n",
    "\n",
    "The mean-field solution is automatically computed when the `SecondQuantizedMolecule` object is created. With the provided spin, it computes the relevant mean-field -- that is to say the Restricted or Restricted Open-Shell Hatree-Fock solution (RHF or ROHF). Using those methods, the spin-orbitals are treated as pairs, so each molecular orbital can accept up to two electrons. The unrestricted Hartree-Fock (UHF) mean-field calculation, where the spin-orbitals are considered individually, usually results in more accurate energies for open-shell systems. Tangelo supports all three approaches. Further  details about the RHF, ROHF and UHF methods can be found in [Szabo & Ostlund](https://books.google.ca/books/about/Modern_Quantum_Chemistry.html?id=6mV9gYzEkgIC&redir_esc=y).\n",
    "\n",
    "### Choice of the Basis Set<a class=\"anchor\" id=\"23\"></a>\n",
    "\n",
    "The many-electron Schrödinger equation is typically solved in a basis of molecular orbitals which are expressed as a linear combination of atomic orbitals (LCAO). LCAO coefficients, used to construct molecular orbitals, are obtained through the minimization of the total energy, according to the variational principle. For convenience, the basis is usually represented by atom-centered functions for molecular systems, or plane-wave functions in solid-state physics.\n",
    "\n",
    "In early quantum chemistry, basis sets were constructed with decreasing exponential functions ($e^{-\\alpha r}$). This strategy was referred to as *Slater Type Orbitals* (STOs). With the increase of computational resources, the use of Gaussian functions ($e^{-\\alpha r^2}$) was popularized. Using those Gaussian functions is much more efficient than the STO strategy for computing the numerous overlap integrals needed in quantum chemistry. This is the main point of using a *Gaussian Type Orbitals* (GTOs) wavefunction.\n",
    "\n",
    "State-of-the-art basis sets are constructed from many Gaussian primitives. The Pople notation is one of the most explicit notations for basis sets:\n",
    "- The hyphen symbol (-) separates the reference to core and valence orbitals;\n",
    "- The digits represent the number of Gaussian primitives in a function;\n",
    "- There are as many digits as there are functions for each atomic orbital (zeta $\\zeta$);\n",
    "- The letter G refers to the use of Gaussian functions.\n",
    "\n",
    "As an example, the `3-21G` double zeta basis set uses three Gaussian primitives for each core orbital and a set of two plus one Gaussian primitives for the valence orbitals. More information, examples and discussions about the GTOs can be found in this [link](https://chem.libretexts.org/Courses/Pacific_Union_College/Quantum_Chemistry/11%3A_Computational_Quantum_Chemistry/11.02%3A_Gaussian_Basis_Sets).\n",
    "\n",
    "Below are popular basis set examples used in quantum chemistry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "basis_sets = [\n",
    "    \"STO-3G\",       # Simple zeta, minimal basis.\n",
    "    \"3-21G\",        # Double zeta.\n",
    "    \"6-31G\",        # Double zeta with more Gaussian primitives.\n",
    "    \"6-31G(d,p)\",   # Polarization functions (+ 5 d-orbitals for all atoms except H, +3 p-orbitals for H atoms) added.\n",
    "    \"6-311G(d,p)\",  # Triple zeta with polarization functions.\n",
    "    \"6-311+G(d,p)\", # Triple zeta with polarization functions and diffuse functions.\n",
    "    \"cc-pvqz\",      # Quadruple zeta.\n",
    "    \"cc-pv5z\"       # Quintuple zeta.\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nKtrojstMfWA"
   },
   "source": [
    "The wavefunction for a molecule is an approximation when using a finite basis set,\n",
    "\n",
    "$$\n",
    "\\Phi_{\\text{approx}} = \\sum_{i=0}^N c_i |i \\rangle\n",
    "$$\n",
    "\n",
    "According to the variational principle, considering $\\Phi_{\\text{approx}}$ is normalized, the approximated energy is greater or equals to the true ground-state energy.\n",
    "\n",
    "$$\n",
    "\\langle \\Phi_{\\text{approx}} | \\hat{H} | \\Phi_{\\text{approx}} \\rangle = \\epsilon_{\\text{approx}} \\geq  \\epsilon_{\\text{true}}\n",
    "$$\n",
    "\n",
    "As the basis set grows larger, the resulting approximation draws closer to the true wavefunction.\n",
    "\n",
    "$$\n",
    "\\Phi_{\\text{approx}} = \\sum_{i=0}^{\\infty} c_i |i \\rangle \\approx \\Psi_{\\text{true}}\n",
    "$$\n",
    "\n",
    "However, increasing the basis set size has the effect of increasing the number of parameters to optimize, therefore extending the required computational resources in space (amount of memory) and in time (number of operations). As we shall see shortly, it is often a good idea to balance the basis set size to optimize the time to solution versus the target accuracy ($|\\epsilon_{\\text{approx}} -  \\epsilon_{\\text{true}}|$), as the increase in the basis set complexity yields diminishing returns in accuracy gain.\n",
    "\n",
    "A full list of available basis sets supported in [PySCF](https://pyscf.org/) (therefore also in Tangelo) is available at this [link](https://pyscf.org/_modules/pyscf/gto/basis.html). Also, the [Basis Set Exchange](https://www.basissetexchange.org/) initiative hosts a collection of bases for every atom in the periodic table. With this constantly increasing number of bases available to the community, choosing a basis set can be a challenging task. The theory states that an infinite basis is the way of getting the true wavefunction. In practice, it is not possible to perform calculations using an infinite basis, and compute the true wavefunction. We can however select finite basis sets and tailor our approaches to the use cases.\n",
    "- As a rule of thumb, a double zeta (e.g. `6-31G` or `cc-pvdz`) is the minimum basis size required to get semi-quantitative energies.\n",
    "- For molecules with a high electronegativity difference between atoms, adding polarization functions (e.g. `6-31G**` or `6-31G(d,p)`) is recommended.\n",
    "- Diffuse functions (e.g. `6-31G+(d,p)`) should be included when dealing with anions. Diffuse functions are functions with small coefficients (small $\\alpha$ in $e^{-\\alpha r^2}$) that can fit the electronic density far from the nucleus.\n",
    "- Effective Core Potential (ECP) basis sets are needed when dealing with heavier atoms with relativistic and non-trivial effects (for example, the iodine atom is often modelled in the `LANDL2DZ` basis).\n",
    "\n",
    "In the next cell, we demonstrate the effect of the basis on the total energy calculations and time to solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwMAAAH/CAYAAAAL7QFoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAACIYUlEQVR4nOzdd1RUx98G8GcpIkiz0BSsWLCgaCyAXSM2ojGWmBAl2KOvDRtRsYu9/TQYu0YTu0nsGlssJFiCoih2QSMoKk3pzPvHDVdWiruCLLDP55x75M6dnZ0dYfd+d5pCCCFARERERERaR0fTFSAiIiIiIs1gMEBEREREpKUYDBARERERaSkGA0REREREWorBABERERGRlmIwQERERESkpRgMEBERERFpKT1NV4Akqamp+Oeff2BlZQUdHcZoRERERIVNeno6IiMj4eTkBD294nEbXTxeRTHwzz//oEmTJpquBhERERG9R2BgIBo3bqzpauQLBgOFhJWVFQDpl8vGxkbDtSEiIiKidz19+hRNmjSR79uKAwYDhUTG0CAbGxvY2tpquDZERERElJPiNKS7+LwSIiIiIiJSC4MBIiIiIiItxWCAiIiIiCgf/Pnnn3B3d0f58uWhUCjw66+/Kl0XQsDX1xc2NjYwNDRE+/btcefOHc1U9j8MBoiIiIiI8sHr169Rv359rFq1KtvrCxYswIoVK7B69Wr8/fffKFWqFNzc3JCYmFjANX2LE4iJiIiIiPJBp06d0KlTp2yvCSGwbNkyTJkyBd26dQMAbNmyBVZWVvj111/x5ZdfFmRVZewZICIiIiLKRlJSEmJjY+UjLi7ug8t68OABIiIi0L59eznNzMwMTZs2RUBAQH5U94MwGCAiIiIiyoafnx/MzMzko3bt2h9cVkREBABk2aPAyspKvqYJDAaIiIiIiLLh4+ODmJgY+QgJCdF0lfIdgwEiIiIiomwYGBjA1NRUPkxMTD64LGtrawBAZGSkUnpkZKR8TRM4gVhLJSQAsbGAqSlgaKjp2hARERGp4M4dQJVx+yYmQPXqH78+aqhSpQqsra1x4sQJNGjQAAAQGxuLv//+G8OGDdNYvRgMaJlz54AlS4DffgPS0wEdHaBbN8DbG3B11XTtiIiIiHJw5w5Qo4bq+W/fLvCAID4+Hnfv3pXPHzx4gKCgIJQpUwYVK1bE6NGjMXv2bFSvXh1VqlTB1KlTUb58eXTv3r1A65kZgwEt4u8PDB8O6OpKgQAg/bt/P/Drr8APPwBDh2q0ikRERETZU3clnzys/POhLl26hDZt2sjnY8eOBQD0798fmzZtwoQJE/D69WsMHjwY0dHRaN68OY4cOYKSJUsWeF0zKIQQQmPPTrLHjx/Dzs4O4eHhsLW1zffyz50DWrYEcvvfViiAs2fZQ0BERESF0JUrQKNGque/fBlo2DBfq/Cx79c0gROItcSSJVKPQG50dYGlSwumPkRERESkeQwGtEBCgjRHIDU193ypqcC+fVJ+IiIiIir+GAxogdjYt3ME3ic9XcpPRERERMUfgwEtYGoqrRqkCh0dKT8RERERFX8MBrSAoaG0fKjee9aO0tMDPv+c+w4QERERaQsGA1pi7FggLS33PGlpwJgxBVMfIiIiItI8BgNaonlzaR8BhSJrD4GenpT+ww9cVpSIiIhImzAY0CJDh0r7CHTr9nYOQcYOxGfPcsMxIiIiKsRMTD5ufi3FHYi1jKurdCQkSKsGmZpyjgAREREVAdWrA7dvq7azsImJlJ/ei8GAljI0ZBBARERERQxv8PMdhwkREREREWkpBgNERERERFqKwQARERERkZZiMEBEREREpKUYDBARERERaSkGA0REREREWorBABERERGRlmIwQERERESkpYpMMDBnzhy4uLjAyMgI5ubmueZ98eIFbG1toVAoEB0dLad7enpCoVBkOerUqZNjWdOnT8/2MaVKlZLzbNq0Kcv1kiVL5vUlExERERF9VEUmGEhOTkavXr0wbNiw9+YdMGAAHB0ds6QvX74cT58+lY/w8HCUKVMGvXr1yrGscePGKT3m6dOnqF27dpbHmJqaKuV59OiR+i+SiIiIiKgA6Wm6AqqaMWMGAOlb+Nz4+/sjOjoavr6+OHz4sNI1MzMzmJmZyee//vorXr16hW+//TbH8oyNjWFsbCyfX716FSEhIVi9erVSPoVCAWtra1VfDhERERGRxhWZngFVhISEYObMmdiyZQt0dN7/0tavX4/27dujUqVKKj/HunXrUKNGDbRo0UIpPT4+HpUqVYKdnR26deuGGzdu5FpOUlISYmNj5SMuLk7lOhARERER5YdiEwwkJSWhb9++WLhwISpWrPje/P/++y8OHz6MgQMHqvwciYmJ2LZtGwYMGKCUXrNmTWzYsAG//fYbtm7divT0dLi4uODx48c5luXn5yf3VJiZmaF27doq14OIiIiIKD9oNBiYNGlStpNzMx+3bt1SqSwfHx84ODjAw8NDpfybN2+Gubk5unfvrnJ99+3bh7i4OPTv318p3dnZGf369UODBg3QqlUr7N27FxYWFvjxxx9zrW9MTIx8hISEqFwPIiIiIqL8oNE5A97e3vD09Mw1T9WqVVUq6+TJkwgODsbu3bsBAEIIAEC5cuUwefJkec5BxrUNGzbgm2++QYkSJVSu77p169C1a1dYWVnlmk9fXx9OTk64e/dujnkMDAxgYGAgn8fGxqpcDyIiIiKi/KDRYMDCwgIWFhb5UtaePXuQkJAgn1+8eBFeXl44e/YsqlWrppT3zJkzuHv3bpbhPrl58OABTp06hd9///29edPS0hAcHIzOnTur/gKIiIiIiApYkVlNKCwsDC9fvkRYWBjS0tIQFBQEALC3t4exsXGWG/6oqCgAgIODQ5Z9CdavX4+mTZuibt26WZ5n5cqV2LdvH06cOKGUvmHDBtjY2KBTp05ZHjNz5kw0a9YM9vb2iI6OxsKFC/Ho0SO15iMQERERERW0IhMM+Pr6YvPmzfK5k5MTAODUqVNo3bq1yuXExMRgz549WL58ebbXo6KicO/ePaW09PR0bNq0CZ6entDV1c3ymFevXmHQoEGIiIhA6dKl0ahRI1y4cIGTgomIiIioUFOIjMH1pFGPHz+GnZ0dwsPDYWtrq+nqEBEREdE7iuP9WrFZWpSIiIiIiNTDYICIiIiISEsxGCAiIiIi0lIMBoiIiIiItBSDASIiIiIiLcVggIiIiIhISzEYICIiIiLSUgwGiIiIiIi0FIMBIiIiIiItxWCAiIiIiEhLMRggIiIiItJSDAaIiIiIiLQUgwEiIiIiIi3FYICIiIiISEsxGCAiIiIi0lIMBoiIiIiItBSDASIiIiIiLcVggIiIiIhISzEYICIiIiLSUgwGiIiIiIi0FIMBIiIiIiItxWCAiIiIiEhLMRggIiIiItJSDAaIiIiIiLQUgwEiIiIiIi3FYICIiIiISEsxGCAiIiIi0lIMBoiIiIiItBSDASIiIiIiLcVggIiIiIhISzEYICIiIiLSUgwGiIiIiIi0VJEJBubMmQMXFxcYGRnB3Nw817wvXryAra0tFAoFoqOjla5t27YN9evXh5GREWxsbODl5YUXL17kWl5YWBi6dOkCIyMjWFpaYvz48UhNTVXKc/r0aTRs2BAGBgawt7fHpk2bPuBVEhEREREVnCITDCQnJ6NXr14YNmzYe/MOGDAAjo6OWdLPnz+Pfv36YcCAAbhx4wZ27dqFwMBADBo0KMey0tLS0KVLFyQnJ+PChQvYvHkzNm3aBF9fXznPgwcP0KVLF7Rp0wZBQUEYPXo0Bg4ciKNHj37YiyUiIiKiIiUtLQ1Tp05FlSpVYGhoiGrVqmHWrFkQQmi6arkTRczGjRuFmZlZjtd/+OEH0apVK3HixAkBQLx69Uq+tnDhQlG1alWl/CtWrBAVKlTIsbxDhw4JHR0dERERIaf5+/sLU1NTkZSUJIQQYsKECaJOnTpKj+vTp49wc3NT+XWFh4cLACI8PFzlxxARERFRwcntfm3OnDmibNmy4sCBA+LBgwdi165dwtjYWCxfvlwDNVVdkekZUEVISAhmzpyJLVu2QEcn60tzdnZGeHg4Dh06BCEEIiMjsXv3bnTu3DnHMgMCAlCvXj1YWVnJaW5uboiNjcWNGzfkPO3bt1d6nJubGwICAnIsNykpCbGxsfIRFxen7sslIiIiokLiwoUL6NatG7p06YLKlSujZ8+e6NChAwIDAzVdtVwVm2AgKSkJffv2xcKFC1GxYsVs87i6umLbtm3o06cPSpQoAWtra5iZmWHVqlU5lhsREaEUCACQzyMiInLNExsbi4SEhGzL9fPzg5mZmXzUrl1b5ddKRERERB+fOl/euri44MSJE7h9+zYA4OrVqzh37hw6depUUNX9IBoNBiZNmgSFQpHrcevWLZXK8vHxgYODAzw8PHLMExISglGjRsHX1xeXL1/GkSNH8PDhQwwdOjS/XpLKfHx8EBMTIx8hISEFXgciIiIiypk6X95OmjQJX375JWrVqgV9fX04OTlh9OjR+PrrrwuwxurT0+STe3t7w9PTM9c8VatWVamskydPIjg4GLt37wYAebJGuXLlMHnyZMyYMQN+fn5wdXXF+PHjAQCOjo4oVaoUWrRogdmzZ8PGxiZLudbW1lm6dyIjI+VrGf9mpGXOY2pqCkNDw2zra2BgAAMDA/k8NjZWpddJRERERAXDx8cHY8eOlc+fPHmSY0Cwc+dObNu2DT///DPq1KkjLypTvnx59O/fv6CqrDaNBgMWFhawsLDIl7L27NmjNCTn4sWL8PLywtmzZ1GtWjUAwJs3b6Cnp/ySdXV1ASDHmd7Ozs6YM2cOnj17BktLSwDA8ePHYWpqKv8yODs749ChQ0qPO378OJydnfPltRERERFRwVPny9vx48fLvQMAUK9ePTx69Ah+fn6FOhgoMnMGwsLCEBQUhLCwMKSlpSEoKAhBQUGIj48HAFSrVg1169aVjypVqgAAHBwc5Jt4d3d37N27F/7+/rh//z7Onz+PkSNHokmTJihfvjwAYN++fahVq5b8vB06dEDt2rXxzTff4OrVqzh69CimTJmC4cOHy78cQ4cOxf379zFhwgTcunULP/zwA3bu3IkxY8YUZBMRERERkYa8efMmywI2urq6SE9P11CNVKPRngF1+Pr6YvPmzfK5k5MTAODUqVNo3bq1SmV4enoiLi4OK1euhLe3N8zNzdG2bVvMnz9fzhMTE4PQ0FD5XFdXFwcOHMCwYcPg7OyMUqVKoX///pg5c6acp0qVKjh48CDGjBmD5cuXw9bWFuvWrYObm1seXzURERERFQXu7u6YM2cOKlasiDp16uCff/7BkiVL4OXlpemq5UohchofQwXq8ePHsLOzQ3h4OGxtbTVdHSIiIiJ6R273a3FxcZg6dSr27duHZ8+eoXz58ujbty98fX1RokQJDdX4/RgMFBIMBoiIiIgKt+J4v1Zk5gwQEREREVH+YjBARERERKSlGAwQEREREWkpBgNERERERFqKwQARERERkZZiMEBEREREpKUYDBARERERaSkGA0REREREWorBABERERGRlmIwQERERESkpRgMEBERERFpKQYDRERERERaisEAEREREZGWYjBARERERKSlGAwQEREREWkpBgNERERERFqKwQARERERkZZiMEBEREREpKUYDBARERERaSkGA0REREREWorBABERERGRlmIwQERERESkpRgMEBERERFpKQYDRERERERaisEAEREREZGWYjBARERERKSlGAwQEREREWkpBgNERERERFqKwQARERERkZZiMEBEREREpKUYDBARERERaakiEwzMmTMHLi4uMDIygrm5ea55X7x4AVtbWygUCkRHRytd27ZtG+rXrw8jIyPY2NjAy8sLL168yLGsq1evom/fvrCzs4OhoSEcHBywfPlypTynT5+GQqHIckRERHzoyyUiIiIi+uiKTDCQnJyMXr16YdiwYe/NO2DAADg6OmZJP3/+PPr164cBAwbgxo0b2LVrFwIDAzFo0KAcy7p8+TIsLS2xdetW3LhxA5MnT4aPjw9WrlyZJW9oaCiePn0qH5aWluq9SCIiIiKiAqSn6QqoasaMGQCATZs25ZrP398f0dHR8PX1xeHDh5WuBQQEoHLlyhg5ciQAoEqVKhgyZAjmz5+fY3leXl5K51WrVkVAQAD27t2LESNGKF2ztLR8b68FEREREVFhUWR6BlQREhKCmTNnYsuWLdDRyfrSnJ2dER4ejkOHDkEIgcjISOzevRudO3dW63liYmJQpkyZLOkNGjSAjY0NPv30U5w/fz7XMpKSkhAbGysfcXFxatWBiIiIiCivik0wkJSUhL59+2LhwoWoWLFitnlcXV2xbds29OnTByVKlIC1tTXMzMywatUqlZ/nwoUL2LFjBwYPHiyn2djYYPXq1dizZw/27NkDOzs7tG7dGleuXMmxHD8/P5iZmclH7dq1VX+xRERERET5QKPBwKRJk7KdeJv5uHXrlkpl+fj4wMHBAR4eHjnmCQkJwahRo+Dr64vLly/jyJEjePjwIYYOHarSc1y/fh3dunXDtGnT0KFDBzm9Zs2aGDJkCBo1agQXFxds2LABLi4uWLp0aa71jYmJkY+QkBCV6kBERERElF80OmfA29sbnp6eueapWrWqSmWdPHkSwcHB2L17NwBACAEAKFeuHCZPnowZM2bAz88Prq6uGD9+PADA0dERpUqVQosWLTB79mzY2NjkWH5ISAjatWuHwYMHY8qUKe+tT5MmTXDu3LkcrxsYGMDAwEA+j42NVel1EhERERHlF40GAxYWFrCwsMiXsvbs2YOEhAT5/OLFi/Dy8sLZs2dRrVo1AMCbN2+gp6f8knV1dQG8DR6yc+PGDbRt2xb9+/fHnDlzVKpPUFBQrsEFEREREZGmFZnVhMLCwvDy5UuEhYUhLS0NQUFBAAB7e3sYGxvLN/wZoqKiAAAODg7yCj/u7u4YNGgQ/P394ebmhqdPn2L06NFo0qQJypcvDwDYt28ffHx85OFJ169fR9u2beHm5oaxY8fKewfo6urKgcyyZctQpUoV1KlTB4mJiVi3bh1OnjyJY8eOfexmISIiIiJ669kzwNoaSE9XKXuRCQZ8fX2xefNm+dzJyQkAcOrUKbRu3VqlMjw9PREXF4eVK1fC29sb5ubmaNu2rdLSojExMQgNDZXPd+/ejefPn2Pr1q3YunWrnF6pUiU8fPgQgLQHgre3N548eQIjIyM4Ojrijz/+QJs2bfLwiomIiIiIPoBCoXpWkdv4GCowjx8/hp2dHcLDw2Fra6vp6hARERHRO4rE/Vpx7RkgIiIiItJ6mUbKZCsmRq3iGAwQERERERUVY8bkfl3NQT8MBoiIiIiIioqXL3O//vw5YGWlcnHFZgdiIiIiIiKtp2bPAIMBIiIiIqLiRI3VhBgMEBEREREVF8bGwLRpKmfP05yBxESgZMm8lEBERERERCr79tvcr2/cKE0y/vZb6ef3ULtnID0dmDULqFBBCjzu35fSp04F1q9XtzQiIiIiIlJZTEzuByDNG1BxiVG1ewZmz5aWN12wABg06G163brAsmXAgAHqlkhERERERCrZu/f9eUxNVcuHD+gZ2LIFWLMG+PprQFf3bXr9+sCtW+qWRkREREREmqJ2z8CTJ4C9fdb09HQgJSU/qkRERERERNlq2zb35UNPnVKrOLWDgdq1gbNngUqVlNN37wacnNQtjYiIiIiIVNaggfJ5Sgpw7Rpw9SrQr5/axakdDPj6Av37Sz0E6enScKTQUGn40IEDaj8/ERERERGpasmS7NPnzpV2H1aT2nMGunUD9u8H/vgDKFVKCg5u3pTSPv1U7ecnIiIiIqK86tNHWuVHTR+0z0CLFsDx4x/ySCIiIiIiylfp6cCOHUCJEmo/9IOCgehoaY7A/fvAuHFAmTLAlSuAlZW0/wAREREREX0EVasqTyAWAnjxQtoNePVqtYtTOxi4dg1o3x4wMwMePgQGDpSCgb17gbAwae4AERERERF9BKNHK5+npADBwdJ8gb591S5O7WBg7FjA01PadMzE5G16587AV1+p/fxERERERKSqkSOzT1+/HhgxAtiwQa3i1J5AfPEiMGRI1vQKFYCICHVLIyIiIiKiPGvdWhrHrya1gwEDAyA2Nmv67duAhYXaz09ERERERHkRGSktOWplpfZD1Q4GPvsMmDnz7W7DCoU0V2DiROCLL9R+fiIiIiKiYuHJkyfw8PBA2bJlYWhoiHr16uHSpUv5+yS6uoCOjvJhYwPs2gX8739qF6f2nIHFi4GePQFLSyAhAWjVShoe5OwMzJmj9vMTERERERV5r169gqurK9q0aYPDhw/DwsICd+7cQenSpfP3ifbtUz7X0ZFuzB0dgZIl1S5O7WDAzEzaY+D8eWnX4/h4oGFDaYUhIiIiIiJtNH/+fNjZ2WHjxo1yWpUqVfL/iT77LF+LUysYSEkBDA2BoCDA1VU6iIiIiIiKo6SkJCQlJcnncXFxOeb9/fff4ebmhl69euHMmTOoUKECvvvuOwwaNCj/K/b0KbBqlbTRl7ExUL++tJKQmZnaRak1Z0BfH6hYEUhLU/t5iIiIiIiKFD8/P5iZmclH7dq1c8x7//59+Pv7o3r16jh69CiGDRuGkSNHYvPmzflbqXv3ACcnaZOvUqWA334Dzp0DatQArl9XuziFEJm3MHu/9eul5/7pJ2mzMcofjx8/hp2dHcLDw2Fra6vp6hARERFpvXd7Bp48eYLatWtne79WokQJfPLJJ7hw4YKcNnLkSFy8eBEBAQH5V6mePaV5Atu3A48eSXMF4uKAqVOBy5eBQ4fUKk7tOQMrVwJ37wLlywOVKkkBSWZXrqhbIhERERFR4WNgYAADAwP5PDa79fX/Y2Njk6XnwMHBAXv27MnfSp06BRw9KgUEmb/T79cPWLZM7eLUDga6dZOWEyUiIiIiIomrqytCQ0OV0m7fvo1KlSrl7xMlJWU/POfNm6zf0qtA7WBg+nS1n4OIiIiIqFgbM2YMXFxcMHfuXPTu3RuBgYFYs2YN1qxZk79PVKkScOcOULXq27THj6VNvzp0ULs4tTcdq1oVePEia3p0tHKdiIiIiIi0RePGjbFv3z788ssvqFu3LmbNmoVly5bh66+/zt8n6tgR2LHj7fmbN9IKP6mpwNKlahends/Aw4fZryaUlCQFJURERERE2qhr167o2rXrx32SxYvfzhWwsZEmDFerBtjbf1BxKgcDv//+9uejR5WXMU1LA06cAD7GvgpERERERJRJxgReQ0PAzS1PRak8TKh7d+kAgP7935537w58+aW0K/HixXmqS67mzJkDFxcXGBkZwdzcPNe8L168gK2tLRQKBaKjo5Wubdu2DfXr14eRkRFsbGzg5eWFF9mNe8pEoVBkObZv366U5/Tp02jYsCEMDAxgb2+PTZs2fcCrJCIiIiIqOCoHA+np0lGpEvD8+dvz9HRpiFBoKPAxe0WSk5PRq1cvDBs27L15BwwYAEdHxyzp58+fR79+/TBgwADcuHEDu3btQmBgoEo7w23cuBFPnz6Vj+4ZkRGABw8eoEuXLmjTpg2CgoIwevRoDBw4EEePHlXrNRIRERERFSS15gykpEiThF++BMqW/VhVyt6MGTMA4L3fuPv7+yM6Ohq+vr44fPiw0rWAgABUrlwZI0eOBABUqVIFQ4YMwfz589/7/Obm5rC2ts722urVq1GlShUs/q9rxMHBAefOncPSpUvhlseuGyIiIiKij0Wt1YT09YFr1z5WVfIuJCQEM2fOxJYtW6Cjk/WlOTs7Izw8HIcOHYIQApGRkdi9ezc6d+783rKHDx+OcuXKoUmTJtiwYQMyb9wcEBCA9u3bK+V3c3PLdbe5pKQkxMbGykdcXJwar5SIiIiIKJPUVODMGbUfpvZqQh4ewPr1wLx5aj/XR5WUlIS+ffti4cKFqFixIu7fv58lj6urK7Zt24Y+ffogMTERqampcHd3x6pVq3Ite+bMmWjbti2MjIxw7NgxfPfdd4iPj5d7GCIiImBlZaX0GCsrK8TGxiIhIQGGhoZZyvTz85N7O4iIiIiIVBYQIC3xmZz8Ni0mBhg9GtiwQZpg3L+/SkWpHQykpkrP8ccfQKNGWTc6W7JE9bImTZr03iE6N2/eRK1atd5blo+PDxwcHODh4ZFjnpCQEIwaNQq+vr5wc3PD06dPMX78eAwdOhTr16/P8XFTp06Vf3ZycsLr16+xcOFCORj4ED4+Phg7dqx8/uTJkyxbWBMRERERKRk+HFi9GjA2BnR136YLIQUBY8dKP3+sYOD6daBhQ+nn27eVr2WscqQqb29veHp65pqnqoo7mZ08eRLBwcHYvXs3AMjDeMqVK4fJkydjxowZ8PPzg6urK8aPHw8AcHR0RKlSpdCiRQvMnj0bNjY2Kj1X06ZNMWvWLCQlJcHAwADW1taIjIxUyhMZGQlTU9NsewUAwMDAAAYGBvJ5bGysSs9NRERERFps505pGc+2bZXTnz8HrKykyb1qUDsYOHVK3UfkzMLCAhYWFvlS1p49e5CQkCCfX7x4EV5eXjh79iyqVasGAHjz5g309JRfsu5/EVXmOQDvExQUhNKlS8s3887Ozjh06JBSnuPHj8PZ2fmDXgsRERERUbZevgTq18+antEzoCa1gwFNCQsLw8uXLxEWFoa0tDQEBQUBAOzt7WFsbCzf8GeIiooCIK3sk7Evgbu7OwYNGgR/f395mNDo0aPRpEkTlC9fHgCwb98++Pj44NatWwCA/fv3IzIyEs2aNUPJkiVx/PhxzJ07F+PGjZOfa+jQoVi5ciUmTJgALy8vnDx5Ejt37sTBgwc/cqsQERERkVaZNg0wMsqabmwsXVPTBwUDly5JPRRhYcrzFgBg794PKfH9fH19sXnzZvncyckJAHDq1Cm0bt1apTI8PT0RFxeHlStXwtvbG+bm5mjbtq3SvIWYmBiEhobK5/r6+li1ahXGjBkDIQTs7e2xZMkSpb0JqlSpgoMHD2LMmDFYvnw5bG1tsW7dOi4rSkRERET5y9dX+vfuXeDmTelnBwfA3v7tNTUohDrjYwBs3w706yftfHzsGNChgzR3IDIS+PxzYONGtetAAB4/fgw7OzuEh4fD1tZW09UhIiIioncUivu16Gjg22+B338HMoa/p6QA7u7A5s3AfyNiVKXWPgMAMHcusHQpsH8/UKIEsHw5cOsW0Ls3ULGiuqUREREREZHKRo+WegXOnwcSE6UjIAC4dw/4gJUu1Q4G7t0DunSRfi5RAnj9WpqrMGYMsGaN2s9PRERERESq+v13aWnRZs2km3CFAmjaVLoRP3BA7eLUDgZKlwYyNsutUEFaahSQeizevFH7+YmIiIiISFXJydJk4XeZmABJSWoXp3Yw0LKltLQpAPTqBYwaBQwaBPTtC7Rrp/bzExERERGRqlq1AiZNAl68eJv28iUwYYJ0TU1qrya0cqU0NAkAJk8G9PWBCxeAL74ApkxR+/mJiIiIiEhVK1ZIY/YrVgRq1pTSQkMBW1vgA5a1VzsYKFPm7c86OlJgQkREREREBaBaNeDGDWnuwI0bUpqDA9C9O/DfZrrqUDkYiI1VLZ+pqdp1ICIiIiIiVenqSmv6f/55notSORgwN899h+OMHZDT0vJcJyIiIiIiKgAqBwOnTr39WQigc2dg3TppRSEiIiIiIip6VA4G3p2crKsrLW9atWp+V4mIiIiIiAqC2kuLEhERERFR8cBggIiIiIioqPjtN2Dz5rfnjx8DS5cCu3d/UHFqLy2aWW4TiomIiIiIKJ/Nmyft+AtIOw67ugJGRsDTp8D168D06WoVp3Iw0KOH8nliIjB0KFCqlHL63r1qPT8REREREanq1i2gaVPp5z/+kFb2uX4dOHFCChI+VjBgZqZ87uGh1vMQEREREVFepaVJPQEAcPw40LGjtLJP7dpAZKTaxakcDGzcqHbZRERERESUn+rVAzZsAL75Bti5E/jhByk9IgIoV07t4jiBmIiIiIioqJgzR5owXKuWtMa/u7uUfvUq0LOn2sXlaQIxEREREREVoNatgfBw4NEjqZdAV1dKHzDgg4pjMEBEREREVJSULg28eAEcOiSdOzgA9vYfVBSDASIiIiKioiImBvD0BH7/HdD771Y+JUUaLrR5M2BurlZxKs8Z8PUFLl9Wq2wiIiIiIspPo0YBd+8C589La/0nJgIBAcC9e8DIkWoXp3LPwOPHQKdOQIkSUuDx2WdAu3bSORERERERFYDffwf27weaNXub1rQpsGYN0LWr2sWp3DOwYYO0YtEvvwAmJsDo0dLqRV98AWzZArx8qfZzExERERGROpKTAWPjrOkmJtKOxGpSa2lRHR2gRQtgwQIgNBT4+28pEPnxR6B8eaBlS2DRIuDJE7XrQURERERE79OqFTBpkjSBOMPLl8CECdI1NeVpArGDg3RMmAA8fy71Wvz+u3Rt3Li8lExERERERFmsWAF06QJUrAjUrCmlhYYCFSq8XV1IDfm2mpCFhbS86QcucUpERERERO9TrRpw44b0DfyNG1KagwPQvfvbPQfUwKVFiYiIiIiKEl1d4PPPpSOPGAwQERERERUVmzernrd///dmYTBARERERFRUjBmjWj4hPk4w8Po1UKqUuo8iIiIiIqI8y+f1/NVaWhQArKwALy/g3Ll8rQcRERERERUwtXsGtm4FNm0C2rYFKleWAoN+/aR9Boi0SUICEBsLmJoChoaarg0RERFpjWPHgJMnpbX909OVr23cqFZRavcMdO8O/PqrtLHY0KHAzz8DlSpJux/v3QukpqpbomrmzJkDFxcXGBkZwdzcPNe8L168gK2tLRQKBaKjo5Wubdu2DfXr14eRkRFsbGzg5eWFF5k3bXjHpk2boFAosj2ePXsGADh9+nS21yMiIvL6sqkQOncO6NFD2vzP2lr6t0cP4Px5TdeMiIiIir0ZM4DOnYHTp6UhQzExyoeaFEIIkdc6/e9/wPjx0u7I5cpJQcKkSYCRUV5LfmvatGkwNzfH48ePsX79+iw3+Zl1794dycnJOHz4MF69eiUHD+fPn0fLli2xdOlSuLu748mTJxg6dChq1KiBvXv3ZltWQkICYt5pWE9PTyQmJuL06dMApGCgTZs2CA0NhampqZzP0tISOjqqxVuPHz+GnZ0dwsPDYWtrq9JjqOD5+wPDh0sremUOfPX0gLQ04IcfpN9/IiIiKn4Kxf1a+fLA/PnAN9/kS3EfvJpQZKS0stGmTcCjR0DPntKGY48fS/X76y+pByO/zJgxA4D0TX1u/P39ER0dDV9fXxw+fFjpWkBAACpXroyRI0cCAKpUqYIhQ4Zg/vz5OZZnaGgIw0xjQJ4/f46TJ09i/fr1WfJaWlq+t9eCiq5z56RAQIisPWAZ5999B9SrB7i6Fnz9iIiISAskJgLOzvlWnNrDhPbuBdzdATs7aYjQd99JQ4a2bgXatJGClN9+k3ouClpISAhmzpyJLVu2ZPuNvLOzM8LDw3Ho0CEIIRAZGYndu3ejc+fOKj/Hli1bYGRkhJ49e2a51qBBA9jY2ODTTz/F+feMGUlKSkJsbKx8xMXFqVwH0owlS96/sZ+uLrB0acHUh4iIiLSQl5d0451P1O4Z+PZb4MsvpfHRjRtnn6d8eWDy5LxWTT1JSUno27cvFi5ciIoVK+L+/ftZ8ri6umLbtm3o06cPEhMTkZqaCnd3d6xatUrl51m/fj2++uorpd4CGxsbrF69Gp988gmSkpKwbt06tG7dGn///TcaNmyYbTl+fn5ybwcVfgkJUpD77hydd6WmAvv2Sfk5qZiIiIjyXUoKsGwZ8McfQIMGgL6+8nU1v5VUu2fg6VPgxx9zDgQA6SZo2rT3lzVp0qQcJ+dmHLdu3VKpXj4+PnBwcICHh0eOeUJCQjBq1Cj4+vri8uXLOHLkCB4+fIihKg7yDggIwM2bNzFgwACl9Jo1a2LIkCFo1KgRXFxcsGHDBri4uGBpLv8ZPj4+iImJkY+QkBCV6kCaERv7/kAgQ3q6lJ+IiIgo3wUHAw0bAiVKACEhwNWryoea1O4ZSE3N/kZHoQAMDKR6qcrb2xuenp655qlatapKZZ08eRLBwcHYvXs3ACBjXnS5cuUwefJkzJgxA35+fnB1dcX48eMBAI6OjihVqhRatGiB2bNnw8bGJtfnWLduHRo0aIBGjRq9tz5NmjTBuVw2YzAwMICBgYF8Hsu7x0LN1BTQ0VEtINDRkfITERER5buTJ/O1OLWDAXNz6cY/J7a2gKen1DPwvoV0LCwsYGFhoW4VsrVnzx4kJCTI5xcvXoSXlxfOnj2LatWqAQDevHkDPT3ll6z73yDw9y2qFB8fj507d8LPz0+l+gQFBb03uKCiw9AQ6NYN2L8/9+Vz9fSkfBwiRERERAUqNhYYNUrtfQbUDgY2bZLmA3h6Ak2aSGmBgdLKQlOmSHsfLFok9RJ8/726pecsLCwML1++RFhYGNLS0hAUFAQAsLe3h7GxsXzDnyEqKgoA4ODgIK/w4+7ujkGDBsHf3x9ubm54+vQpRo8ejSZNmqD8f7um7du3Dz4+PlmGJ+3YsQOpqanZDkNatmwZqlSpgjp16iAxMRHr1q3DyZMncSw/l1MijRs7VtpjIzdpacCYMQVSHSIiItJGd+5IN9sPH0rr+mdITgYCAqR0ADh1SqXi1A4GNm8GFi8Gevd+m+buLi2n+OOPwIkTQMWKwJw5+RsM+Pr6YvPmzfK5k5MTAODUqVNo3bq1SmV4enoiLi4OK1euhLe3N8zNzdG2bVulpUVjYmIQGhqa5bHr169Hjx49sl06NDk5Gd7e3njy5AmMjIzg6OiIP/74A23atFHvRVKh1ry5tI/Ad9/lvs8AlxUlIiKij8bTU7rpaNpUeZnDN2+ktf3/u0dWldqbjhkaAteuAdWrK6ffuQPUry/V48EDoE4d6WdSTaHYxIJUcv68NFF/3z5pDoGODvD551KPAAMBIiKi4qtQ3K+VKiVNHK5USTn9+XPAykr1FU/+o3bPgJ0dsH49MG+ecvr69dI1AHjxAihdWt2SiYoGV1fpSEiQhueZmnKOABERERWQxETA2DhruhC5T+zNgdrBwKJFQK9ewOHDb5cXvXQJuHUL+G8hH1y8CPTpo3ZdiIoUQ0MGAURERFTAHjwAypTJmm5hIV1Tk9rBwGefAaGh0vyAjKH1nTpJEysrV5bOhw1Tux5ERERERPQ+Dx7kftNfsaI0sfH8eaBVq/cWp1YwkJICdOwIrF4NqLjCJhERERER5Ze2bXMeEiSENGfg5UspX1rae4tTKxjQ15cmDxMRERERkQa8evX+PJaWquUD8J5twbLy8JAmCxMRERERUQEzNc39yJxPBWrPGUhNBTZsAP74A2jUSFrdKLMlS9QtkYiIiIiIVPbPP9LGXzdvSucODtIa540aqV2U2j0D168DDRsCJibA7dtSXTKO/zYFJiIiIiLSavPmzYNCocDo0aPzt+B9+4AmTYCoKKB7d+mIipI2IduzR+3i1O4ZUHFnYyIiIiIirXTx4kX8+OOPcHR0zP/CfX2BadOAKVOU0+fOldK/+EKt4tTuGchw9y5w9Ki08RIgTV4mIiIiIioukpKSEBsbKx9xcXHvfUx8fDy+/vprrF27FqU/xi68d+5kv6FX797SDbqa1A4GXrwA2rUDatQAOncGnj6V0gcMALy91X5+IiIiIqJCyc/PD2ZmZvJRu3bt9z5m+PDh6NKlC9q3b/9xKmVjA1y4kDX9wgXpmprUDgbGjJGWGA0LA4yM3qb36QMcOaL28xMRERERFUo+Pj6IiYmRj5CQkFzzb9++HVeuXIHfx9yQa+RI4LvvgAkTgN9/l46JE6Vdf//v/9QuTu05A8eOScODbG2V06tXBx49Uvv5iYiIiIgKJQMDAxgYGMjnsbGxOeYNDw/HqFGjcPz4cZQsWfLjVWrMGKBcOWDhQmDFCimtenXA3x/o10/t4tQOBl6/Vu4RyPDyJZCprYiIiIiItMbly5fx7NkzNGzYUE5LS0vDn3/+iZUrVyIpKQm6urr582TffCMdGZN2s9uNWEVqDxNq0QLYsuXtuUIh7Xq8YAHQps0H14OIiIiIqMhq164dgoODERQUJB+ffPIJvv76awQFBeVfIJCZQpGnQAD4gJ6BBQukCcSXLgHJydJwpRs3pJ6B8+fzVBciIiIioiLJxMQEdevWVUorVaoUypYtmyU9T6pWzX0ZzwcP1CpO7WCgbl1ps7GVK6WNx+LjgR49gOHDP2gCMxERERERqerdTcxSUoDgYODgQWDsWLWLUzsYAAAzM2Dy5A95JBERERGRdjh9+nT+FzpyZPbpq1cDFy+qXdwHBQPR0UBgIPDsmTRfILMPmMRMRERERER50aGDNH5//Xq1HqZ2MLB/P/D119LwIFNT5TkLCgWDASIiIiKiArdrF/ABOx6rHQx4ewNeXsDcudkvMUpERERERB9Jw4bKE4iFACIigKgoaa8BNakdDDx5Ig1VYiBARERERFTAundXPtfRASwtpTX+q1dXuzi1gwE3N2lZ0apV1X4uIiIiIiLKC1/ffC1O7WCgSxdg/HggJASoVw/Q11e+/tln+VU1IiIiIiL6mNQOBgYNkv6dOTPrNYUCSEvLa5WIiIiIiKggqB0MvLuUKBERERERFU06mq4AERERERFphsrBQOfOQEzM2/N586TNxzK8eAHUrp2PNSMiIiIiIsmiRUC3btKmYs+f51uxKgcDR48CSUlvz+fOBV6+fHuemgqEhuZbvYiIiIiIKEPfvkCLFsDmzYCdHeDqCixYANy8madiVQ4GMu9tkN05ERERERF9JBUqAOPGAX/+Cfz7LzBkCPD330DTpkDNmsCECcDZs2pP8OWcASIiIiKioqRMGaBfP2DPHmnn4WXLgLg4qffA2hrw8lK5KJWDAYVCOt5NIyIiIiIiDSlRAujUCfD3Bx4/Bg4fBmxtVX64WsOEPD2BHj2kIzERGDr07bkaAcgHmTNnDlxcXGBkZARzc/Nc87548QK2trZQKBSIzjzLGcCqVavg4OAAQ0ND1KxZE1u2bHnvc4eFhaFLly4wMjKCpaUlxo8fj9TUVKU8p0+fRsOGDWFgYAB7e3ts2rRJzVdIRERERJRHjRplvyFYDlTeZ6B/f+VzD4+sefr1U/l51ZacnIxevXrB2dkZ69evzzXvgAED4OjoiCdPniil+/v7w8fHB2vXrkXjxo0RGBiIQYMGoXTp0nB3d8+2rLS0NHTp0gXW1ta4cOECnj59in79+kFfXx9z584FADx48ABdunTB0KFDsW3bNpw4cQIDBw6EjY0N3Nzc8qcBiIiIiIjymUKIojUVeNOmTRg9enSWb/wz+Pv7Y8eOHfD19UW7du3w6tUruSfBxcUFrq6uWLhwoZzf29sbf//9N86dO5dteYcPH0bXrl3x77//wsrKCgCwevVqTJw4Ec+fP0eJEiUwceJEHDx4ENevX5cf9+WXXyI6OhpHjhxR6XU9fvwYdnZ2CA8Ph60aXTtEREREVDCK4/1asZpAHBISgpkzZ2LLli3Q0cn60pKSklCyZEmlNENDQwQGBiIlJSXbMgMCAlCvXj05EAAANzc3xMbG4saNG3Ke9u3bKz3Ozc0NAQEBOdY1KSkJsbGx8hEXF6fy6yQiIiIiyg/FJhhISkpC3759sXDhQlSsWDHbPG5ubli3bh0uX74MIQQuXbqEdevWISUlBVFRUdk+JiIiQikQACCfR0RE5JonNjYWCQkJ2Zbr5+cHMzMz+ajNHduIiIiI6H1u3QICA9+ep6cD//wD/Hdfqi6NBgOTJk2CQqHI9bh165ZKZfn4+MDBwQEe2U1m+M/UqVPRqVMnNGvWDPr6+ujWrRv6/zcZIruehI/Jx8cHMTEx8hESElKgz09ERERERdCYMcCJE2/PO3YEPvkEqFgR2L1b7eI0Ggx4e3vj5s2buR5Vq1ZVqayTJ09i165d0NPTg56eHtq1awcAKFeuHKZNmwZAGhK0YcMGvHnzBg8fPkRYWBgqV64MExMTWFhYZFuutbU1IiMjldIyzq2trXPNY2pqCkNDw2zLNTAwgKmpqXyYmJio9DqJiIiISIv98w/QoYP0c0CAtPHYw4fAokVqrSKUQeXVhD4GCwuLHG/C1bVnzx6lITkXL16El5cXzp49i2rVqinl1dfXlyd9bN++HV27ds2xZ8DZ2Rlz5szBs2fPYGlpCQA4fvw4TE1N5aE9zs7OOHTokNLjjh8/Dmdn53x5bUREREREAIDYWOC/e1IcOwa4uQF2dkD37sCkSWoXp9FgQB1hYWF4+fIlwsLCkJaWhqCgIACAvb09jI2Ns9zwZ8wBcHBwkFcTun37NgIDA9G0aVO8evUKS5YswfXr17F582b5cfv27YOPj488PKlDhw6oXbs2vvnmGyxYsAARERGYMmUKhg8fDgMDAwDA0KFDsXLlSkyYMAFeXl44efIkdu7ciYMHD37kViEiIiIirVKpEvDXX0CFCsDevcD//Z+UHhcHGBmpXVyRmUDs6+sLJycnTJs2DfHx8XBycoKTkxMuXbqkchlpaWlYvHgx6tevj08//RSJiYm4cOECKleuLOeJiYlBaGiofK6rq4sDBw5AV1cXzs7O8PDwQL9+/TAzUzdMlSpVcPDgQRw/fhz169fH4sWLsW7dOu4xQERERET5a+RIaXOvSpWAyEigZ08p/exZoGFDtYsrcvsMFFfFcd1aIiIiouKk0NyvHTsGhIQAn38uBQUA8OYNoFAAOcxXzUmRGSZERERERESQJhBnTCLO8AFDhIAiNEyIiIiIiIjyF4MBIiIiIiItxWCAiIiIiEhLMRggIiIiItJSDAaIiIiIiLQUgwEiIiIiIi3FYICIiIiISEsxGCAiIiIi0lIMBoiIiIiItBSDASIiIiIiLcVggIiIiIhISzEYICIiIiLSUgwGiIiIiIi0FIMBIiIiIiItxWCAiIiIiEhLMRggIiIiItJSDAaIiIiIiLQUgwEiIiIiIi3FYICIiIiISEsxGCAiIiIi0lIMBoiIiIiItBSDASIiIiIiLaWn6QoQERERaY07d4C4uPfnMzEBqlf/+PUhrcdggIiIiKgg3LkD1Kihev7btxkQ0EfHYUJEREREBUGVHoG85Cf6AAwGiIiIiIi0FIMBIiIiIiItxWCAiIiIiEhLMRggIiIiItJSDAaIiIiIiLQUgwEiIiIiIi1VZIKBOXPmwMXFBUZGRjA3N88174sXL2BrawuFQoHo6Gila6tWrYKDgwMMDQ1Rs2ZNbNmyJdeyrl69ir59+8LOzg6GhoZwcHDA8uXLlfKcPn0aCoUiyxEREfEhL5WIiIiIqEAUmU3HkpOT0atXLzg7O2P9+vW55h0wYAAcHR3x5MkTpXR/f3/4+Phg7dq1aNy4MQIDAzFo0CCULl0a7u7u2ZZ1+fJlWFpaYuvWrbCzs8OFCxcwePBg6OrqYsSIEUp5Q0NDYWpqKp9bWlp+4KslIiIiIvr4ikwwMGPGDADApk2bcs3n7++P6Oho+Pr64vDhw0rXfvrpJwwZMgR9+vQBAFStWhUXL17E/PnzcwwGvLy8lM6rVq2KgIAA7N27N0swYGlp+d5eCyIiItJSJiYfNz/RBygywYAqQkJCMHPmTPz999+4f/9+lutJSUkoWbKkUpqhoSECAwORkpICfX19lZ4nJiYGZcqUyZLeoEEDJCUloW7dupg+fTpcXV1zLCMpKQlJSUnyeRx3GSQiIireqlcHbt9WbWdhExMpP9FHVmyCgaSkJPTt2xcLFy5ExYoVsw0G3NzcsG7dOnTv3h0NGzbE5cuXsW7dOqSkpCAqKgo2NjbvfZ4LFy5gx44dOHjwoJxmY2OD1atX45NPPkFSUhLWrVuH1q1b4++//0bDhg2zLcfPz0/u7SAiIiItwRt8KmQ0OoF40qRJ2U68zXzcunVLpbJ8fHzg4OAADw+PHPNMnToVnTp1QrNmzaCvr49u3bqhf//+AAAdnfc3xfXr19GtWzdMmzYNHTp0kNNr1qyJIUOGoFGjRnBxccGGDRvg4uKCpUuX5lrfmJgY+QgJCVHpdRIRERFR4ePn54fGjRvDxMQElpaW6N69O0JDQzVdrffSaDDg7e2Nmzdv5npUrVpVpbJOnjyJXbt2QU9PD3p6emjXrh0AoFy5cpg2bRoAaUjQhg0b8ObNGzx8+BBhYWGoXLkyTExMYGFhkWv5ISEhaNeuHQYPHowpU6a8tz5NmjTB3bt3c7xuYGAAU1NT+TDhuEAiIiKiIuvMmTMYPnw4/vrrLxw/fhwpKSno0KEDXr9+remq5Uqjw4QsLCzeexOuqj179iAhIUE+v3jxIry8vHD27FlUq1ZNKa++vj5sbW0BANu3b0fXrl1z7Rm4ceMG2rZti/79+2POnDkq1ScoKEilYUdEREREVDipM8fzyJEjSuebNm2CpaUlLl++jJYtW360OuZVkZkzEBYWhpcvXyIsLAxpaWkICgoCANjb28PY2DjLDX9UVBQAwMHBQV7h5/bt2wgMDETTpk3x6tUrLFmyBNevX8fmzZvlx+3btw8+Pj7y8KTr16+jbdu2cHNzw9ixY+W9A3R1deVAZtmyZahSpQrq1KmDxMRErFu3DidPnsSxY8c+ZpMQERER0UeUlzmeMTExAJDtojOFSZEJBnx9fZVu2p2cnAAAp06dQuvWrVUqIy0tDYsXL0ZoaCj09fXRpk0bXLhwAZUrV5bzxMTEKI3v2r17N54/f46tW7di69atcnqlSpXw8OFDANIeCN7e3njy5AmMjIzg6OiIP/74A23atPnwF0xEREREGuXj44OxY8fK50+ePEHt2rXf+7j09HSMHj0arq6uqFu37sesYp4phBBC05Ug4PHjx7Czs0N4eLg8hImIiIiICg9V79eGDRuGw4cP49y5c4X+vq7I9AwQERERERV2I0aMwIEDB/Dnn38W+kAAYDBARERERJRnQgj83//9H/bt24fTp0+jSpUqmq6SShgMEBERERHl0fDhw/Hzzz/jt99+g4mJibzojJmZGQwNDTVcu5xpdJ8BIiIiIqLiwN/fHzExMWjdujVsbGzkY8eOHZquWq7YM0BERERElEdFdU0e9gwQEREREWkpBgNERERERFqKwQARERERkZZiMEBEREREpKUYDBARERERaSkGA0REREREWorBABERERGRlmIwQERERESkpRgMEBERERFpKQYDRERERERaisEAEREREZGWYjBARERERKSlGAwQEREREWkpBgNERERERFqKwQARERERkZZiMEBEREREpKUYDBARERERaSkGA0REREREWorBABERERGRlmIwQERERESkpRgMEBERERFpKQYDRERERERaisEAEREREZGWYjBARERERKSlGAwQEREREWkpBgNERERERFqKwQARERERkZYqMsHAnDlz4OLiAiMjI5ibm+ea98WLF7C1tYVCoUB0dLTStVWrVsHBwQGGhoaoWbMmtmzZ8t7nVigUWY7t27cr5Tl9+jQaNmwIAwMD2NvbY9OmTWq+QiIiIiKigqWn6QqoKjk5Gb169YKzszPWr1+fa94BAwbA0dERT548UUr39/eHj48P1q5di8aNGyMwMBCDBg1C6dKl4e7unmuZGzduRMeOHeXzzAHJgwcP0KVLFwwdOhTbtm3DiRMnMHDgQNjY2MDNzU39F0tEREREVACKTDAwY8YMAHjvN+7+/v6Ijo6Gr68vDh8+rHTtp59+wpAhQ9CnTx8AQNWqVXHx4kXMnz//vcGAubk5rK2ts722evVqVKlSBYsXLwYAODg44Ny5c1i6dCmDASIiIiIqtIrMMCFVhISEYObMmdiyZQt0dLK+tKSkJJQsWVIpzdDQEIGBgUhJScm17OHDh6NcuXJo0qQJNmzYACGEfC0gIADt27dXyu/m5oaAgIAcy0tKSkJsbKx8xMXFqfISiYqNhAQgMlL6l4iIiDSj2AQDSUlJ6Nu3LxYuXIiKFStmm8fNzQ3r1q3D5cuXIYTApUuXsG7dOqSkpCAqKirHsmfOnImdO3fi+PHj+OKLL/Ddd9/hf//7n3w9IiICVlZWSo+xsrJCbGwsEnK40/Hz84OZmZl81K5d+wNeNVHRc+4c0KMHYGwMWFtL//boAZw/r+maERERaR+NBgOTJk3KdnJu5uPWrVsqleXj4wMHBwd4eHjkmGfq1Kno1KkTmjVrBn19fXTr1g39+/cHgGx7EjI/ztXVFU5OTpg4cSImTJiAhQsXqvdis6lvTEyMfISEhOSpPKKiwN8faNkS2L8fSE+X0tLTpfMWLYDVqzVbPyIiIm2j0WDA29sbN2/ezPWoWrWqSmWdPHkSu3btgp6eHvT09NCuXTsAQLly5TBt2jQA0pCgDRs24M2bN3j48CHCwsJQuXJlmJiYwMLCQuV6N23aFI8fP0ZSUhIAwNraGpGRkUp5IiMjYWpqCkNDw2zLMDAwgKmpqXyYmJio/PxERdG5c8Dw4YAQQGqq8rXUVCn9u+/YQ0BERFSQNDqB2MLCQq2b8Nzs2bNHaUjOxYsX4eXlhbNnz6JatWpKefX19WFrawsA2L59O7p27Zprz8C7goKCULp0aRgYGAAAnJ2dcejQIaU8x48fh7Oz84e+HKJiZ8kSQFc3ayCQma4usHQp4OpacPUiIiLSZkVmNaGwsDC8fPkSYWFhSEtLQ1BQEADA3t4exsbGWW74M+YAODg4yMuA3r59G4GBgWjatClevXqFJUuW4Pr169i8ebP8uH379sHHx0cenrR//35ERkaiWbNmKFmyJI4fP465c+di3Lhx8mOGDh2KlStXYsKECfDy8sLJkyexc+dOHDx48CO2CFHRkZAA/Pbb26FBOUlNBfbtk/Ln0KlGRJp05w6gyoIXJiZA9eofvz5ElGdFJhjw9fVVuml3cnICAJw6dQqtW7dWqYy0tDQsXrwYoaGh0NfXR5s2bXDhwgVUrlxZzhMTE4PQ0FD5XF9fH6tWrcKYMWMghIC9vT2WLFmCQYMGyXmqVKmCgwcPYsyYMVi+fDlsbW2xbt06LitK9J/Y2PcHAhnS06X8DAaICpk7d4AaNVTPf/s2AwKiIkAhMq+RSRrz+PFj2NnZITw8XB7CRFRcJCRIqwapEhDo6ADx8QwGiAqdK1eARo1Uz3/5MtCw4cerD5EGFMf7tWKztCgRFV6GhkC3boDee/oi9fSAzz9nIEBERFRQGAwQUYEYOxZIS8s9T1oaMGZMwdSHiIiIGAwQUQFp3hz44QdAocjaQ6CnJ6X/8ANXEiIiIipIDAaIqMAMHQqcPSsNGcpYzVdHRzo/e1a6TqpLSAAiI6V/iYiIPgSDASIqUK6uwO7d0iThiAjp39272SOgjnPngB49pEnZ1tbSvz16cMM2IiJSH4MBItIIQ0PAyoqThdXl7w+0bAns3/92dab0dOm8RQtg9WrN1o+IiIoWBgNEREXEuXPA8OGAEFl3ck5NldK/+449BEREpLois+kYkcZx503SsCVLAF3drIFAZrq6wNKlHHZFRESqYTBApAruvEkalpAA3Pj1DuqJ9wSkqUDwXhMkJFTnECzKXyYmHzc/EWkEgwEiVajSI5CX/ETv8TroDkKFigGpAKKCbsPQmQEp5aPq1aUvOthDSlSsMBggIioCjN/XI5DH/NoqIQGIjQVMTTmZXSW8wScqdjiBmIioCChZ8uPm1zZcnjV/cK8LoqKPPQNEVDA4AZsKCX9/aVUmXd2sy7P++qu0EzY3wMvduXPShPbffpPaLmPzQG9vTl4nKmoYDBDRx8cJ2FRIvG95VkBanrVePd7U5oTBFFHxwmCAiD4+TsCmQuKXmXfwiU4cUtNyzqOnA/wy0wSuRxmQvovBFFHxw2CAiIi0QmLwHaw6rkIPVRqAY0Bi8G2UrMeAIDPudZH/OImdNI3BgLbgeG0i0nKvI+Kgzrzq1xFxKFnvo1WnyElIeDtHIDepqcC+fVJ+3tzmjPMu8geDqbxjMKANOF6biAilSn3c/MVdbOz7A4EM6elSft6cZY/zLvKOwVT+4dKi2oDjtfOOO28SFXlcnjVvTE2lGy5V6OhI+Smr9827EEKad8FlbnPm7w+0bCkFT+8GUy1aAKtXa7Z+RQ17BohUwZ03SdMYkJKGGRpK37zu35/7nAE9PSkfewWyx3kXecNJ7PmPwQCRqniDT5rEgJQKgbFjpWEsuUlLA8aMKZDqFDmcd5F3DKbyH4MBIqKigjf4pGHNm0vj2b/7LusNmZ6eFAj88ANvwnLCeRd5w2Dq4+CcASIiIlLZ0KHA2bPSUKCMOQQZkzfPnuXE19xw3kXefEgwRe/HYICIPj6OdycqVlxdgd27gfh4ICJC+nf3bvYIvE/GvAu994zL0NMDPv+c32q/q6gEU6tWrULlypVRsmRJNG3aFIGBgZqpiIo4TIiIPj6OdycqlgwNecOqLs67+HBFYRL7jh07MHbsWKxevRpNmzbFsmXL4ObmhtDQUFhaWhZ8hVTAYICICgZv8EnT2ENFhQDnXeRNYQ+mlixZgkGDBuHbb78FAKxevRoHDx7Ehg0bMGnSJM1U6j0YDBARkXZgDxUVEkOHSktfLl0qTXTNvGnWmDEMBHJT0MFUUlISkpKS5PO4XN4/kpOTcfnyZfj4+MhpOjo6aN++PQICAvKnQh8BgwFtwG/DiIgkvMGnQsLVVToSEqSJrqamHHKlqoIMpvz8/DBjxgyV8kZFRSEtLQ1WVlZK6VZWVrh161b+VSqfMRjQBvw2jIiIqFDivIsPU1DBlI+PD8aOHSufP3nyBLVr187/J9IgBgPagjf4REREVMx87GDKwMAABgYG8nlsLuuVlitXDrq6uoiMjFRKj4yMhLW19UerY15xaVEiIiIiojwqUaIEGjVqhBMnTshp6enpOHHiBJydnTVYs9yxZ4CIiIiIKB+MHTsW/fv3xyeffIImTZpg2bJleP36tby6UGFUZHoG5syZAxcXFxgZGcHc3DzXvC9evICtrS0UCgWio6OVrq1atQoODg4wNDREzZo1sWXLllzL2rRpExQKRbbHs2fPAACnT5/O9npEREReXjIRERERFSF9+vTBokWL4OvriwYNGiAoKAhHjhzJMqm4MCkyPQPJycno1asXnJ2dsX79+lzzDhgwAI6Ojnjy5IlSur+/P3x8fLB27Vo0btwYgYGBGDRoEEqXLg13d/dsy+rTpw86duyolObp6YnExMQsm0eEhobCNNN2d4V1cwkiIiIi+jhGjBiBESNGaLoaKisywUDGsk6bNm3KNZ+/vz+io6Ph6+uLw4cPK1376aefMGTIEPTp0wcAULVqVVy8eBHz58/PMRgwNDSEYaaZKc+fP8fJkyezDUgsLS3f22tBRERERFRYFJlhQqoICQnBzJkzsWXLFujoZH1pSUlJKFmypFKaoaEhAgMDkZKSotJzbNmyBUZGRujZs2eWaw0aNICNjQ0+/fRTnD9/PtdykpKSEBsbKx+5bWJBRERERPQxFJtgICkpCX379sXChQtRsWLFbPO4ublh3bp1uHz5MoQQuHTpEtatW4eUlBRERUWp9Dzr16/HV199pdRbYGNjg9WrV2PPnj3Ys2cP7Ozs0Lp1a1y5ciXHcvz8/GBmZiYfxW3NWiIiIiIq/DQaDEyaNCnHybkZh6o7tvn4+MDBwQEeHh455pk6dSo6deqEZs2aQV9fH926dUP//v0BINuehHcFBATg5s2bGDBggFJ6zZo1MWTIEDRq1AguLi7YsGEDXFxcsHTp0lzrGxMTIx8hISEqvU4iIiIiovyi0TkD3t7e8PT0zDVP1apVVSrr5MmTCA4Oxu7duwEAQggA0gYQkydPxowZM2BoaIgNGzbgxx9/RGRkJGxsbLBmzRqYmJjAwsLivc+xbt06NGjQAI0aNXpv3iZNmuDcuXM5Xn93E4uMVY+ePn363rKJiIiIqOBl3Kelp6druCb5R6PBgIWFhUo34arYs2cPEhIS5POLFy/Cy8sLZ8+eRbVq1ZTy6uvrw9bWFgCwfft2dO3a9b09A/Hx8di5cyf8/PxUqk9QUBBsbGxUrn/GbnVNmjRR+TFEREREVPAiIyNzHJZe1BSZ1YTCwsLw8uVLhIWFIS0tDUFBQQAAe3t7GBsbZ7nhz5gD4ODgIK/wc/v2bQQGBqJp06Z49eoVlixZguvXr2Pz5s3y4/bt2wcfH58sw5N27NiB1NTUbIchLVu2DFWqVEGdOnWQmJiIdevW4eTJkzh27JjKr8/JyQmBgYGwsrJSachSXsXFxaF27doICQmBiYnJR3++4ohtmDdsv7xjG+Yd2zBv2H55xzbMm4Juv/T0dERGRsLJyemjP1dBKTLBgK+vr9JNe8Z/wqlTp9C6dWuVykhLS8PixYsRGhoKfX19tGnTBhcuXEDlypXlPDExMQgNDc3y2PXr16NHjx7ZLh2anJwMb29vPHnyBEZGRnB0dMQff/yBNm3aqPz69PT00LhxY5Xz51VsbCwAoEKFCkp7I5Dq2IZ5w/bLO7Zh3rEN84btl3dsw7zRRPsVlx6BDAqRMbietEpsbCzMzMwQExPDN58PxDbMG7Zf3rEN845tmDdsv7xjG+YN2y/vis3SokREREREpB4GA1rKwMAA06ZNU1rRiNTDNswbtl/esQ3zjm2YN2y/vGMb5g3bL+84TIiIiIiISEuxZ4CIiIiISEsxGCAiIiIi0lIMBoiIiIiItBSDASIiIiIiLcVggIiIiIhISzEYICIi+ki4YF/+Y5sqY3t8GLbbWwwGiKjQevPmDV69eoWXL19quipEH0ShUAAA0tPTNVyTou/UqVNK52xTIC4uDgqFgm2hpvT0dPlv8/Xr1wC0OzhgMEAqCQsLw48//ogFCxbg2rVrmq5OkfDy5UvcvHkT586d03RViqTbt2/ju+++w/Dhw7FmzRoA2v1m/bENGzYMy5cvx8GDBzVdlWJh7NixaN68OYYNG4arV6/KNx70Yfbt24dRo0ahXbt2mDlzJsLDw6Gjo923MDNnzkSzZs0QFRWl9W2hjvT0dLm9Fi5ciL59+yI2Nlar/0b520Pvdf36dXTq1AmBgYFITk6Go6OjpqtU6N28eRMeHh6YNm0aTp48qXSNN7TvFxwcjNatW8PKygr9+/fHpEmTAADx8fEarlnxValSJURFRaFPnz7w9PTE1q1bNV2lIm3GjBkYOXIk0tPT0bx5c4wdOzbLewGprmPHjggKCoKbmxuuXr2KOnXqYMuWLVrbazhq1ChMnz4d//77L6KiojRdnSIjLS1NDgRGjx6NiRMn4q+//pJ7B7SWIMrFjRs3ROnSpcWUKVPEmzdv5PTt27eL9evXa7Bmhde1a9eEhYWFmDJlirh8+bKcfuXKFfnn9PR0TVStSLh3756oUKGCGD9+vFL6okWLRMWKFcXNmzc1VLPi6d3fxaCgIPHVV1+Jpk2bisGDB2uoVkXbu226Y8cO0bFjR9G8eXOxZcsWDdWqaEtNTZV/fv36tfD19RWlSpUSU6ZMEWFhYRqsWcHz9vYWZcqUEf/8849o0KCB+P333zVdpSIh89/l2LFjRdmyZcW5c+dE3bp1xV9//aXVn8sMBihH0dHRws3NTQwdOlSkpaXJ6X5+fkJXV1c0aNBArFmzRoM1LHzCw8NF9erVxciRI5XSlyxZIoyNjYWPj4+cps1vPNlJT08X6enp4vvvvxeff/65iIqKkq/5+fkJY2Nj0aBBA1GpUiUREhKiwZoWXxk3XM+fPxf+/v6iZs2aomfPnhquVdGV+W/80qVLYvDgwaJevXril19+0WCtiq7Mn0NCCPHDDz8IS0tLMXPmTJGcnKyhWhWsQYMGCSMjI/mLpurVq4vdu3dnyfduW9FbgwcPFqVKlRJXrlwR6enpwtLSUpw+fTpLPm36jOYwIcpRdHQ0bt++jc6dO8vdaj///DMWLVqEH3/8EY0bN8bmzZvl8dzaTPw39OfQoUOoUKECxo0bJ19bsmQJfH190b17d+zfvx+TJ08GAK0en5gdhUIBhUKBM2fOoHTp0ihbtiwA4N69e/jrr7+wd+9ebNu2DQ0bNoSbmxtCQkI0XOOi791Jh7q6ugCAcuXKoV+/fpg5cyZCQkLkYVqUu3fbM/PEzkaNGmH06NH45JNPsHXrVly/fl0TVSxS3m3PjM+htLQ0ANI8l6lTp2LWrFk4f/58to8pbmrVqoWzZ8+iYcOGAKThfc+fP1fKwzkEym7evCn/LISAjY0Nzp49CycnJ6SlpcHGxgbR0dFKj4mOjtaqz2j+tlCOgoOD8ejRI7Rq1UpOa9OmDQ4fPowBAwZg0qRJqFmzJubNm6f1k2Qz3jTOnz+P9PR02NnZQQiBhIQEPHjwAMePH8f//vc/eHh44Ndff8X48eM1XOPCJz09HfHx8Xj58iUqVKggp1WrVg1r167Fp59+itq1a2P69OlIT0/H8uXLNVzjoi3zJLq7d+8iNDQUcXFx8nUjIyN07twZHh4euHjxIv766y9NVbVIyNyex44dw5YtW/DXX3/h33//lfM4ODjg22+/RXh4uDx/oLjfvH4oIYTcnmvXrsXEiRMxYsQIXLlyRQ4GAGDEiBHw8vLCsGHD8PLly2J7E3zo0CEA0sT0hg0byl9A6evrIygoCIDUZh06dMCwYcM0Vc1CZ9WqVahTpw4ePnwIQPqsnj59OpycnCCEgJ6eHgwMDBAcHAxAasPPP/8c06ZN02CtC17x/KuhPMl4k7G3t4eBgQG2bNkiX7OxsUHjxo0hhEDVqlXRsWNHWFpawtbWVlPVLVTS09PliUgKhQKGhoZYunQpmjVrBnNzc3h5eaFatWq4e/euhmta+Ojo6MDY2BjNmjXDxo0bcffuXfmDvVy5cvJNk42NDZycnNC8eXNNVrfIy2jbCRMmoEuXLmjQoAE+//xzTJ06Vc5jbGyMb775BtHR0di/f7+mqlokZG7PXr16YdasWejcuTO8vLzw+++/y/latGiBgQMHYtasWYiMjCy2N695kXnZx/Hjx2PChAm4du0aLl68iObNm2PFihV49uyZnP///u//UKlSJfzzzz8Ait8iDZcvX8bgwYOVvkRKTU0FAJQvXx4pKSkAgM6dOyMsLAw///yzRupZ2Pz444/w9vbG9u3bUbly5SzXM9rQ3Nwcb968AQB06dIF169fx6JFiwqyqhrHdyGSZfxhZPxrZmaGqlWrYtu2bQgMDFTKm/FGfenSJZQrVw7m5uYFWtfCJuPDx8XFBQ8fPsSGDRuyXBNCwMLCAubm5qhXr16x+8DKi8zf9HXs2BGpqamYOnUq7t27B0D6fcu4aVqxYgXu37+v1GNFqsv8TfS2bduwfft2LF68GLt27ULjxo2xfft2DB06VM5ja2sLPz8/bN++HXfu3NFElQu1zH/H58+fx9GjR3Hw4EGEhIRg165dsLS0xOzZs3HkyBE537Bhw+Ds7IyLFy9qosqFXsbfekREBCIjI/HHH3/g8OHD+PvvvzF58mTMmzcPe/fulfPXqVMHpqam8vtucRre8cMPP2DVqlVISUnBqlWr5CGo+vr6AIB69eohIiICn376Ke7evYvg4GDo6+vLn+PaasOGDRgxYgR27dqF3r17IyIiAleuXMGBAwfkz5WMNnRycsLz58/RuXNn3L17FyEhIdrXhgU/TYEKo1u3bokhQ4aI7t27i++++07cv39fCCHEn3/+KUqWLCnat28vTpw4Ied//vy58Pb2FmZmZiI4OFhT1daopKQk+eeMiUbh4eGiRo0aolatWlkmdaWkpIjvv/9e2NnZiTt37hRoXQuj0NBQsWDBAvk8JSVF/nn8+PGiTJkyon379uLs2bPi5cuX4ty5c2Lo0KHC2NhY/PPPPxqocfFy6tQpMXLkSLFo0SI57dWrV2LdunWievXqYu3atUII6Xf7yZMnomvXrlr7t66KxYsXi//7v/8TAwYMUEoPDAwUXbp0EQMGDJAnyQshxJAhQ8TcuXM1UdUiYcOGDaJ06dLC0dFR3L9/X2ky5+TJk0Xp0qVFZGSknHbjxg0xdepUERcXp4nqfhTTp08X5ubm4ueffxZ79uwRX3/9tahRo4bSAhXLly8XCoVCNGjQQJ5Enfm9VBtFRUWJ2rVriypVqgghhLhz545wdHQU9erVEwqFQjRq1EhMmjRJzj916lShUChEvXr1tLYNGQyQuHbtmihTpozw8vISHh4eolWrVmL48OHyze7+/fuFsbGxqFChgnB3dxd9+vQRbm5uws7OTmm5TG1y69Yt8eWXX4rAwEA5LWMlluDgYFGuXDlRrVo14e3tLW7fvi22bNkiBg4cKMzMzLS2zTK7c+eOsLKyEiVKlFB6U05MTJR/nj17tqhbt67Q0dERZcqUEbVq1RLNmjUTV69e1USVi4309HRx7949YWxsLBQKhfD29la6HhMTI9zd3cWgQYOU0rds2SLOnDlTkFUtUgYNGiTfUDx79kzp2qpVq4SxsbGIiIiQ02JiYsT58+cLuppFQnp6ujh8+LBo3bq1MDExkb+cev36tRBCiKdPnworKytx6NAh+TFxcXHFKlh99uyZcHFxEatXr1ZKmz17trCzsxMTJkwQQgjx5s0bMWPGDPnmVdtuYrOTmpoqzpw5I+zs7ETjxo1FnTp1xLhx48S1a9fE/fv3xdSpU0WdOnXEjz/+KIQQ4p9//hETJkzQ6jZkMKDl7t27J6pWrSomT54sp40ZM0YMHz5cCCFEQkKCEEIKGEaPHi1atmwp3NzcxOzZs8W9e/c0UmdNu3fvnqhYsaIwNzcX3bt3F5cuXRJCSB9gGW8it27dEt26dRNly5YVenp6okqVKuKzzz4T169f12TVC4WXL1+KHj16iO7du4vJkyeL2rVri3HjxsnXMwcEd+7cEb/99ptYvXq1CAgIEM+fP9dElYulM2fOiMqVK4tPPvlEBAQEKF2bOHGiaNWqlUhMTNSq5fXyatq0aUKhUIgVK1aImJgYOf306dOidu3a4sGDB5qrXCGW3TKYqamp4uTJk6Ju3bqiTp06SvvcPHjwQFSoUEEcO3asIKtZoJKSkkTt2rXFxIkTldKjo6OFs7Oz0NPTE2PGjFG6po03sZm9+3t07tw5UatWLdG7d2/x+vVr+XpMTIxo2bKl6NevnxBCeQlRbW1DPU0PUyLNEULgzJkzcHFxwdixYyGEkJfCu3jxIlq3bo0SJUpg9uzZaNKkCRYtWiQvPaitEhMT5fZo27Yt9u7di+nTp2P69Olo1KgRdHV1kZqaipo1a2Lr1q14/fo1bt++jZo1a8LQ0BAmJiaafgkaV6JECVhaWqJDhw5o2bIlTExM5EnqCxcuhIGBAZKTk1GiRAnY29vD3t5ewzUufoQQaNmyJTZs2ABPT0+sWLECSUlJaNWqFV69eoWzZ8+ievXqMDAw0HRVC6WM98oMaWlp0NXVxfTp0xETEwNvb2/ExMSgY8eOMDc3x9y5c2FmZoaKFStqsNaFz5kzZ9CqVSvo6OgorcYESMvctmrVCv/73//wf//3f3BycsLMmTOhr6+P9evXw9LSEm3bttVg7T+OX375BXZ2dmjevDlcXV0RGhqKu3fvyu+DZmZmaNq0KczNzXHu3DmMHDkSK1asAACtn4yeeQn0qlWrwtXVFT///DNSUlJgZGQEQPpbNTU1ReXKleU5AZn/lvX0tPS2WLOxCGna8+fPxY0bN+Tz+fPnixIlSoj58+eLxYsXi969ewtTU9Ms32hp87eFO3fulDdb27dvn2jfvr3o2rWr3EMghPI3FNrcVu/KaIuM7n4hhIiMjBR+fn7Z9hBo67c0+cHT0zPbjXQyZPxf/PHHH6JixYrCyspKtG/fXnz++eeiWbNmcg8Nf38l/fv3F23atJHP322XzDvkjh07VigUCmFkZCS++eYb4e7uLo9F5mZQkqVLl4oaNWoo7cicXdukpaWJU6dOiYYNGwqFQiG+/fZbMW/ePPk9JHO7F3WPHj0SNWvWFMuXLxdCCHH27FlRtmxZMXDgQHHt2jUhhDQs6PPPPxcrVqwQixYtEvXr1xfffvutJqtdqLzbhpl7mjPEx8eLli1bitmzZxd09QotBgMki4uLE3369BFHjx6V00JCQoSVlVW2Oxxqm5w+xPfs2ZMlIEhMTBT//vtvQVavUEtMTFSacJ0h8463c+fOVQoIBg8enKUbnFSTmpoqWrZsKaytrcWFCxdyzJdxQ/vnn38KOzs70ahRI7Fp0yb5enb/Z9rqwIEDwtLSUvTq1UtOyy0gmD59ulAoFGLTpk0iPj5eCKG9QxCyc/nyZeHh4SGaN28uNm7cKKfnFBD88ccfolWrVqJhw4ZyIJB56FBx4evrK6ytreXJ0SdOnBB2dnaiWbNmwsXFRTRq1Eg4ODgIIaRhvDNnzhQuLi5K81G0XUYbZt7FXgjp/ez+/fuic+fOolGjRvx7zITBAAkh3r4Bv/vt1c2bN0X9+vWzjCnWJm/evBHJycnixYsXclpaWprSB//u3bvlgOCvv/4Sw4YNE9bW1hxzLaTfoR49eghXV1fh6uoqHj58KITIeiP1/Plz4efnJ+rVqydq1aol9PX1xd9//62JKhcLycnJ4vPPPxcWFhY5BgRpaWny/8OJEydEpUqVRL9+/ZR6C+mtEydOCBsbG9GnTx85LfPvceb2FEKI4cOHCwMDA7F58+ZitcpNXmXchN2+fVsMHDhQtGrVSvzyyy/y9YzPn8zvscnJyeL48eOiYcOGwtnZudi1Z8Zn7+PHj0Xz5s3FwoUL5Xa6ceOGWLt2rRg+fLiYNm2anFcIaQWwzJ9N2uTdwPHdNly8eLFIT08XaWlpIi0tTWzevFl06tRJNGvWTM5bnHqW8oLBgJZ79w/h3Uj5+++/F40aNVJawk2b3LhxQ/Ts2VN88sknolWrVuKnn35Sup75g3/v3r3Czc1NlClTRpibm/NGVrxdqWrQoEFi2bJlwtnZWTRt2lS+nvFmntGOT548EfXr1xdlypQpViuDFKTMf9OPHz8WTZs2FQ4ODlkC+swfpBkfjMeOHRP29vaie/fubP//ZG7P48ePi1mzZgmFQqG0hOi7gW1YWJj8c8aQoZ9//vnjV7YIyPx7d+DAATFkyBBRtmxZUatWLbF9+/Zs82WsIJacnCxOnTolqlatKtq2bVtwlf6Itm7dKh4/fiyfp6eni0GDBim9T2aH32q/pWobhoWFiR07dsh/02zDtxgMaLGMP4jw8HCxZMkSpWt37twR3t7ewtzcXAQFBWmiehp3/fp1Ubp0aTFq1Cjh5+cnRowYISpVqpRlOcDMNwLt2rUT5ubmvJESQjx8+DDLahh79+4VX331ldKcgQzJycli0qRJwsDAgMuH5oMJEyaI9u3bixYtWghDQ0NRvnz5bHsI5s2bJzZt2iT/Hh86dEg4OjqKJ0+eFHSVC7Xx48eL6tWrizFjxohWrVqJEiVKKA0Zyrh5HTdunOjWrZvSYydNmiRu3rxZkNUtdN4NmMaPHy9sbGzE/PnzxYwZM4SDg4No2rSp0hwCIYT4+eefhUKhkH93U1NTxZ9//ikvN1qUZfS8m5iYiO+//15eHSkqKkpUqlRJaR8Kbe9hzok6bZgZewSUMRjQUhkfXA8fPhTly5dXmrj5zz//iPHjx4sGDRpobSDw/Plz0aJFC6U12B88eCA++eQT4e/vL4TIuhzZhAkThL6+vta22bu2b98u+vfvr7Tmure3t7C2thaOjo6ibt26Yt26dfJ46jdv3ogBAwZwQ7F8sG7dOmFiYiL+/vtv8eTJExEcHCw6deokrKyslAKCpKQkUadOHaWJsUKIbIM1bXbu3Dlhbm4uTp06JYSQ2mfPnj2idOnSSkOGhJD2ZdHV1ZXzUtYb2dDQUGFvby/2798vp12+fFl88cUXwsnJSezcuVNODwkJEX369MnSK1scZHwOL126VPTo0UOYmJiIwYMHi8OHD4vJkyeLfv36iZcvX2q4loXLu0OD2Ib5g8FAMRceHi727t0rJk2aJFasWKH0ARUVFSUsLS3F0KFDs7xZX758WTx9+rSAa1t4BAUFiY4dO2ZZjeXrr78Ww4YNE0Jk/WZh586dDATekbEChhBCrFixQigUCrFs2TJx9uxZMXToUFGuXDmlvRe40or6fvvtNyGEcttNmzZNdO3aVSlfdHS0aNOmjahSpYrSELZnz54JR0dHERAQIL8PaPO3kEOGDJG/dc5oh/3794vy5cuL6OhoOV9SUpLYtGmTUCgUYujQoUpljB8/XumGVpv5+fmJmTNnCiGUhwNaW1srzRMQQvoiqmzZssLR0VH+0kUI6ffZz8+v4Cr9kZ08eVIcPXpU/Pnnn3Las2fPxNGjR0XTpk1Fu3bt5E0Bz549q8GaFi6Z3+PYhvmLwUAxdvXqVVGtWjXh7OwsqlatKiwtLYVCoRCjRo0SDx48EImJieKHH37gMpjZePbsmdLulhk3/v3791caK0zZe/f3KDY2VixbtixLcFWmTBkxf/78gqxasXL58mVRtmxZcffuXaX0iRMnikqVKsnnGWNjf/rpJ6FQKETJkiVFcHCwSE1NFcnJyeLcuXMFWe1CrUuXLvKyhBm/x/fv3xdlypQR27ZtU8p7+/ZtYW1tLRQKhdJO2qGhoQVX4UIsLS1NTJkyRQwePFgpLSIiQjRq1EhMnTpVJCcnK71fdOzYUTg4OIiRI0cWy6Eco0ePFpaWlsLCwkJUqFBBtGvXTrx69Uq+/uzZM3HixAnx+eefi7p163Jc+38y/46wDfMfg4Fi6vbt26Js2bLi+++/lyf/3rx5U8ybN0/o6emJr776qlguy/YxZA6Whg8fLjw9PeXzCRMmKC2LRznLvExlamqqePTokWjWrJk4fPiwBmtV9Hl4eIi1a9cKId5OBA4ODhYODg7C29tb6ff3+PHjYvjw4WLy5MnZfkBq85cBGa99+vTp4osvvpDTU1JSRHR0tPj6669F586dlb4k+Pfff4WHh4c4e/Zstjeu2tyeGUJCQkSFChXEgQMHlNJXr14tdHR0xOrVq+VhafHx8aJ3795iy5Ytcttlbtei3p7nzp0TdevWFX///be4efOmOHnypKhZs6Zo2rRprn+P2n4zm/n/nW34cTAYKGbS09NFamqq+L//+z95q+13P6T8/f3l4RqUu3eHrYwdO1Z88803QghppSWFQsFVg3Lw7pvvu205depUUa9ePaVVIEh9kyZNEg0aNFBKe/36tZg1a5ZwdnYWgwcPFpGRkeLWrVuic+fOSkNa+AGZ1b///itKly4t9w5kuHDhgnBzcxMuLi5i6tSpYt++faJdu3aidevW2S6FSdLffHx8vBg4cKAYMWKEeP78udJ1Pz8/oaenJ/r06SOGDh0qmjdvLho0aCC3Z3EaNrhz507Rs2dP4eXlpZR+7949UalSJaUvmdhbnz224cej3XtXF0MKhQK6urq4cuUKLCwsAEjbumf2zTffoHfv3lixYgWio6M1UMvCL2ObciGE0nl8fDxMTU2xdOlSLF68GJcuXUKTJk00Vs/C4M2bN0hOTlZKS0tLg56eHh48eIDZs2dDCCFvFX/u3DlMmDABy5cvx5YtW1ChQgVNVLvYmDZtGhITE+Hl5SWnGRkZYcSIEejduzcuXLgAOzs7dOrUCf/++y9WrFgh59PT09NElQs1GxsbTJgwAb/88guOHj0qpzs7O2POnDlo06YN1q9fjxkzZiA9PR3Hjh2Djo4O0tPTs7zXajsdHR2UKlUKHTp0wM8//4x9+/bJ76kAMGnSJGzfvh2mpqZ4+PAhqlevjsDAQLk9M94ziroXL15g+/btOHPmDMLDw+X01NRUVK1aFSNHjsT169flz+PMr1uhUBR0dQulzG0YFhYmp2e04ahRoxAcHMw2/FAaDkYoH0VFRYno6GiRlJQkqlWrJsaPH59j3jVr1ghTU1Ot37XwxYsX4ubNm+L27dvyMJaMb/fu378vHB0dlVbDGTlypFAoFMLU1FRcvHhRI3UuTIKDg0WXLl3EmTNn5G3fM76FefDggbCxsZF7qISQfkenTp0qmjRpojS5mFT377//ijt37ijt/bFp0ybRsGFDeaJmhpSUFJGSkiKOHDkizp07x/W1s3Ht2jVx6NAhpcmsQUFBomvXrqJHjx7ZTj58/fq1eP78OYcgZOPs2bNi2bJlYsWKFUobg/n6+goDAwOxc+fO9/YaFsf2vH37tvj222+Fubl5ll6nLVu2iJo1ayp91lBWt2/fFl5eXtm24U8//SRq1KjBNvxADAaKicuXLwtzc3N5DfwBAwaIatWqicuXL8t5MnbhE0J686lXr16x28VRHcHBwcLJyUnUq1dPGBgYiFmzZsljrh8+fChsbW2Fh4eHUhfjihUrhLW1NfcRENI+DObm5mLIkCFKmywJIURMTIxwcnISAwcOzNJF++rVqyzDBUg1u3fvFt26dRM9evRQWmoxMjJSjB8/XjRp0kRMmTJFTs9umAWHsry1detW0aRJE/HNN9+IDRs2KF07cOCAaN68ufjss8/Enj17ciyjOA1lyastW7YIBwcHMWrUqCx71wghxKhRo4Surq5YtGiR0iTr9PT0YruSVebP3YcPH4r+/fuLRo0aiXnz5omEhATx8OFD0b59e9GhQ4di99rzC9vw42MwUAwEBQUJY2NjMWbMGDltz549olSpUsLDw0OEhIRkecywYcNE9+7dRUJCQkFWtdC4ceOGKFu2rBg3bpy4ceOGWLRokVAoFOLRo0ciJSVFuLu7i6FDh2b5hur+/ftaveRqhvj4eNGhQwd5mVUhpAnqQUFBIiwsTKSnp4tDhw4Vy2/4NGX9+vWibNmyYv369eLSpUtyekYPQVhYmJg6daqoW7eu6NSpk3j69GmW3hp6a+PGjcLExETs2rVLad5K5p2Cjxw5Ir755htRrVo1MXHiRKXeAFK2efNmUbJkSbFr1y6lxSn8/f3ljaCEEGLBggWiXr16okePHsLf31/+AqY4yhzkBAUFiaioKHHv3j3Rv39/oaenJypVqiR69+4tunTpIn8WM7hUxjYsGAwGiribN28KMzMzMWHCBCGEcvfqnDlzhEKhEF26dBG///67SE5OFrdu3RI+Pj7CxMREa7/dfv78uWjZsqUYNWqUnJaeni46duwoLly4IG7evCmv3U7ZS0xMFM2bNxdXrlwRqampws3NTTRu3FiYmJiIJk2aKO0iypunvNu/f78wNzdXulEVQogvv/xSlC1bVt6r4cWLF+LPP/8UTZs2FY0bNxZfffWVOHHihNJKTiStSGJnZyevwpShb9++QkdHR3z77bdy2oMHD8TevXuFvb29cHNzE127dhVhYWG84cjk+vXronbt2mLNmjVK6T179hQKhUJ06tRJnDx5Uk4/ffq0WLFihShfvrzo1auX0p4CxUHGe17G78ju3buFtbW1POQs42b23aF9GcG7tspu0i/bsGAwGCjCgoKChJmZmTAyMhJTp06VP/Azf/AvX75c1KpVSygUCmFhYSHq1KkjHB0dtXqX16ioKDF37lxx+/ZtOW3mzJlCoVCIevXqiUqVKgk3NzeuvZ6LiIgIYWFhIY4dOybGjBkj3NzcxNWrV8Xhw4fFuHHjhLW1tdi1a5emq1nkpaWlicTEROHp6SlGjhyp9EH3xRdfiEqVKom2bdsKOzs7pc3bhBBi3759Yvr06WLo0KFc9/4/GTcYfn5+okOHDiIqKkq+5unpKWrUqCFmzJghGjZsqBQQCCH1hgUGBoqVK1eKgICAAq13YXfo0CFhb28v7ty5I9+8jR49WlSpUkWcPHlS1K9fX3Tp0kX88ccfSo+LiooSf//9t7h69aomqv1R/Pzzz+LYsWPyBnW7d+8Wenp6WXauDw0NFZ6ensLFxSXLEDVtlDkQYBsWPAYDRdSVK1dEqVKlhLe3t5gzZ45o2rSpGDdunHyzkDkgCAkJEceOHROrVq0Sf/75J4e5CGkTrAy//PKLUCgUYseOHeLFixfizJkzonHjxmL69OkarGHhlp6eLr788ksxYsQI0bVrV3HkyBH5Wnh4uPDw8BBDhw4Vqamp7BnIozdv3oiKFSuKOXPmyGmXLl0SAwYMEI8ePRLBwcHiiy++EFZWVuLmzZtZHh8fH1+Q1S3UMn4X3dzcRM+ePYUQ0k1IamqqWLx4sYiKihKJiYli7dq1olatWvJShZxnkbtZs2aJcuXKKaWFh4fLc4kePXokatSoIVq3bi0PayuOPSvjx48XNjY2Ys2aNfK8qB9//FH88MMPSvkyfg9v3bolBgwYIBwcHMTmzZsLvL6FEdtQMxgMFEH//vuvMDIykocGxcfHiylTpsgBQUYgUJzHYuanhw8fKk20FkLahdTd3V1DNSoaLl68KEqVKiUUCoX4/fffla55e3uLli1bMhDIBxEREcLa2lr873//E0K8/RDM2KhJCCECAgKEhYWFWLhwoUbqWNR8/fXXom7dujkOKXj16pX47LPPlPZkoJytXbtWGBkZiYCAgCx/8xlDVxcsWCA6depUbIPTH3/8UVhZWYlLly6p9dl748YN8d1334kHDx58vMoVEWxDzSkei/hqkUePHmHz5s2YO3cu5s+fj/T0dJQqVQqTJk3Cp59+irNnz2Ly5MlITk6Gvr4+0tLSNF3lQq9SpUpo2LAhACA9PR2JiYkwNjZGs2bNNFyzwu2TTz7B4cOHAQBr1qzBjRs35GspKSmoUaOGvD8DfTgrKys4Oztj5cqVuHfvHhQKBYQQMDQ0RHp6OgCgTJkyqFu3LurXr6/h2hZeQvryCwDg4uKCyMhIrF69GikpKQAg/wtIe7MkJyejRo0aGqlrUdOpUycYGhpi+fLliI+PV7qmp6eH2NhY/Pnnn6hTpw5KlSqloVp+HBm/U4GBgfjss8/QqFEjef+OjGsi094K76pduzaWLVuGypUrf/S6Fjbvtg/bUHMYDBQhwcHB6NChAy5evAgHBwcA0sYaaWlpKFWqFHx8fLIEBLq6uvINA72fjo4O5s6di4CAAPTq1UvT1Sn0WrRogdOnT+PKlSvw8vLCwIED0a9fP2zevBmjRo2Cvr6+pqtYZGUO5N3d3fHy5Uv4+vri/v37UCgUUCgU0NHRQVxcHLy9vVGiRAm0a9dOgzUunJ4/fw4AchAFSBsvVqpUCfPnz8eGDRuQmpoKfX19CCEQGRmJ3r1749WrVxg5cqQmq14kpKeno0KFCpg3bx727t2LESNG4OHDh/L18PBw9OrVC48fP4afnx+A3G/sipoXL14AAK5du6a06ZwQAgqFAqmpqbhy5UqWICkzbXyfTE9PlzcDYxsWAhrojaAPcPPmTVG6dGkxceJE8eTJkxzzvX79WkydOlW4urqKYcOGcRURNezcuVMMHz5clC1bVly5ckXT1SlSbt26JaZMmSLat28vhg0bprUrVeXV7t27xffffy+fZ+4qHzdunChTpoxo3ry5OHz4sLh27ZrYvXu3aNWqlahTp46ctziOxf5Qv/zyi/jkk0/EoUOH5LSMdoqKihJ16tQRVlZWwt3dXWzfvl1MnDhRtGrVSjRo0EDOx/kCknd/r949j46OFkuXLhWlSpUSlStXFu7u7qJTp06icePGokmTJsWyPRcvXiyGDRsmwsLCxNy5c4WtrW2WzShv374tPD09ucliJpl/d9iGhQODgSIgISFB9OrVSwwfPlwpPTk5WYSHh4tbt24ppb9+/VqMHTtWtG/fXmmXUsrd9evXRe/evbPdl4FUk3lzGFLPvn37hK6urihXrpyYOHGinJ55XPuCBQtE48aNha6urjAwMBANGzYUvXv3lsdlc1+Htw4fPixKly4tatSoITp06KA0yT3jS5JXr16JMWPGiPr16wtTU1PRsmVLMX78eLZnDhITE5X2DMhon7CwMJGYmCjS09PFP//8Izw8PMSnn34q+vfvL1asWFEsd76eMGGCsLCwED///LN49OiRuHTpkvj000+Fm5ub+Pvvv4UQQjx9+lS4u7sLV1dXvi9mg21YeCiEKEb9dcVUamoq2rZti969e2PEiBEAgKNHj+LIkSPYsGEDypYti6pVq+L48eNyt9ubN2/w+vVrWFhYaLLqRU5KSgq7G6nA3b59G4MHD0b9+vVRrlw5bN++HZ07d8bChQsBAElJSTAwMAAARERE4N69e3j9+jXs7e1RpUoVuSs9Y5yttouNjcXAgQNhaWmJ9u3bY+PGjYiPj8eECRPg5uYGAEhOTkaJEiWQlpYGXV1dPHnyBBUqVJDLyEint1auXIk9e/Zg2LBh6N27NwBg+/bt8PPzw+7du1G9evUcH1uc2vPEiRMYNGgQfvrpJ7i6usrpv/32GzZv3oyjR4/C1tYWurq6MDQ0xF9//QV9fX2kp6dDR4ejswG2YWHDT44i4M2bN3j+/DmuXbuG0NBQ7N27F5s3b0bdunUxa9YsGBsbw8/PD+PGjcPixYuRlpYGIyMjGBkZabrqRQ4DAdIEGxsbODo6wsPDAzVq1IC+vj5++uknAMDChQthYGAgB6rW1tawtrZWenx6ejoDgUxMTU0xfPhwJCQkoGPHjihdujSWLFmCBQsWAADc3NxQokQJpKenyzeoVlZW8uOFEMXmxjU/dejQAZGRkdi+fTvKlSuH5ORkDBkyBLNmzco1EABQrNozLCwMRkZGqFOnDgDIN6jdunWDk5MTHjx4gGvXrsHKygpffPEFdHV1Gay/g21YuLBViwBTU1OsWrUKbm5uOHbsGF6+fImFCxeiXbt2sLe3R0pKCnbs2CFPwilOb7pE2sDExAQLFixAyZIlAQBeXl4AoBQQ6OvrIzo6GqVKlcoStPKbsqxatWqV5ed3A4KoqCg8ePAATZs2VbrJyOhhJWU1atRAv379sG/fPkyePBlXr17F6tWr0a9fP634xlb8N6E1ISFBaYK/QqGQez7+/vtvNGjQQOn3Ly0tjTex/2EbFk7F+y+3GGnbti3u37+PPXv24P79+xgyZAjs7e0BSDf/ZmZmsLOzU1o+j4gKr5cvX+LZs2fyecYwoLS0NFhaWmLgwIHw8PDAoUOHMHHiRMTGxqJdu3byzSxl793llDNWU2vVqhXGjBkDY2NjLFy4EDt27MBnn32Gfv36ASheK9x8DOnp6UhPT0f16tVhYmKCK1euoFatWjA2NgYgBaTFfeW6jCCxTZs2uHPnDpYtWyan6+rqIj4+Hlu3bsWxY8eUHscv6N5iGxZOnDNQxCUnJ2PWrFnYsGEDTp8+/d6uWiLSvO3bt2PVqlV4+vQpTExM8Ntvv6FixYry9Yxvz168eIH169dj8+bN+Pfff2FhYYEbN25wONs79uzZgytXrmDOnDkAso5Pz2hPAPjzzz8xb948HDlyBA0bNkRAQADbMxuZ2yyzvXv3YsCAAfD29kZycjJu3LiBL7/8UuuWYl6zZg1GjBiBYcOGoWvXrihRogTmzp2LiIgIXL58md9iq4BtWHgwGCjCtm7diosXL2LHjh04fPgwnJycNF0lInqPjRs3YuTIkZgxYwYqV66MefPmoVy5cjh06BCAt2NnM27GQkJC0KpVK1SvXh1//vkn9PT0OHY2k19//RU9e/ZE6dKlMWDAAMybNw9AzgFBQkICGjZsCFNTU5w/f57t+Y7nz5/LC0+8GxBcvHgRTZs2xQ8//IChQ4fi3r172LRpE86ePQtfX1+0bdtWU9UucEII/P777xg5ciTS0tJgbm6OChUq4MCBA/KGn/w2O3dsw8KDwUARFRoaiqFDh6J06dKYM2eOvAkZERVef/zxBwYNGgQ/Pz98+eWXAIBVq1bh0aNHGD16NMqWLQt9fX157HVMTAz69u2LBw8eIDg4mDeu73jfKkzv3kwkJyfjyy+/xLVr13Dz5k3o6+uzPTPZvn07Fi9ejJkzZ6JTp04AlAOCe/fu4eHDh0qb24WGhuLMmTMYOHBgsZ8zkJ2oqCjExMQgPT0d1apVg46ODn+n1MQ21DwGA0XYs2fPYGBgADMzM01XhYhUsGHDBty5cwdTp06VV/tq3749QkJCUKJECejq6mL8+PHw9PREyZIlER8fj9mzZ2PWrFm8cc1GXFwcJk+ejG+++QY1atSAv78/fvrpJ6WA4N2JrRcvXoSTkxMDq3ccOXIEX331FSwsLFC5cmWMHTtWXoY1pyFD79KGScTvwzbIO7ZhwWMwQERUgB49eoRKlSoBAAYOHIijR49izZo1qFu3LhYsWIBdu3bh4sWLsLOzU3ocb1yzl5iYKK/C9OzZM2zYsCFLQBAdHQ0jIyOUKFFCfhyHILylyr4MvEEjKr4YDBARFYB3v119+fIlfv31V3z66afyjX9CQgLKlSuHTZs2ad2ETFW9fPkSqampsLS0BPC2XTNu7qOiorB27Vps3boVXbt2xeTJk9GmTRv06NEDkydP1nDtC68zZ87I+zKcOXMGS5YsyRIQZFC1p4CIigYGA0REGvLut9NXrlzBkCFDsGbNGi4IkA2uwlRwsgsInj17Ju/LQETFB4MBIqICknGz+u6/gDS5tWfPnkhLS8P+/fs5JOMdXIXp43g3IM08HOj06dNYunQpXr9+jUGDBmHp0qV49eoVQkND2TtAVIzwXZGIKJ9FRESgdOnS8kZiwNtAYOfOnbh69SrmzJkDhUKBxMREHDhwAJs2bcKjR49w5coVeQMnBgSSP/74AzNnzsTatWvlVZiePn2KR48e4d9//5VXYQKkzYtiYmIwbtw4lCtXjoFANjLvy6Crq6sUEGQOqFq3bg0dHR3MmzcPffv2RcOGDXH9+nUA3KWZqDjhJw0RUT7atGkT+vbti2vXriElJUVOVygU2LNnDzw9PWFtbS2nR0dH4+LFizA0NMQ///wjrxrEQOCtsLAwfPnll/jss8/ktH379mHr1q1wcXFB7dq1sWbNGiQmJgKQdit1dHTEtWvXGAi849dff0WfPn2wZs0aTJo0CQDkgCBDRq8VADRu3BgPHjxA48aN8ddff8m/n0RUfHCYEBFRPvnpp58wZMgQLFq0CB4eHjA1NZWvhYaGokuXLvD29sawYcOUHhcfH49SpUopTYQlZVyFKe+4LwMRZYfBABFRPvj333/RrVs3eHl5YdiwYXj+/DlCQkKQkpKCSpUqoXr16vjnn39ynRjMcdhZcRWm/MN9GYgoO/yrJiLKB2lpaUhJScHXX3+NkJAQ9OzZE0ZGRnjw4AFsbGzg7e2Nb7/9NtcyGAhk9W6blClTBv3791f6BvvmzZuoXbs27O3tC7p6RYqJiQkWLFgg78vg5eUFQOrRAoCFCxdCR0dHaV+Gxo0bA5B+vxkIEBVPHJRKRJQPkpKS8PTpU1y6dAkTJ05Ex44dsX//fhw5cgTu7u6YMmUKDh48qOlqFkkZHdgZ/2b+5jo5ORnTp0+HpaUl6tevr5H6FWYvX77Es2fP5POMSe1paWmwtLTEwIED4eHhgUOHDmHixImIjY1Fu3bt5J6CDBy6RlR8MRggIsoHdnZ2aNy4MX777TckJSXB09MTNjY2aNy4MYYPHw4nJydcuHABwNubWlIWERGBpKQkpbTMqzBNmTIFAORVmHbv3o0ePXrgwYMH+PXXX+VVmEiyfft2dOvWDS4uLnByckJYWJjc06KrqwshBMqVK4fBgwejf//+OHDgACpVqoS4uDhMmDBBw7UnooLCYICIKB8YGBhgypQp+N///oc//vgDcXFx8jVbW1uUKVMG4eHhADgcKDtchSl/bdy4EYMGDcLnn3+OBQsWQF9fH0OHDpWvp6eny6sGlS1bFl27dsWzZ8/g4OCAkJAQrhpEpEU4AJCIKJ80a9YMZ86cQatWrbBkyRKYm5ujTp06ePPmDSIjI+Hq6qrpKhZKP/30E7777jssWrQINWvWVNopODQ0FBMnTsTixYuVVmGytrbG1KlTlVZh4ph2CfdlICJ1cDUhIqJ8dvbsWfTr1w8lS5ZEjRo18Pz5c8THx+PKlSu8wXoHV2HKfxs2bMCdO3cwdepUGBkZAQDat2+PkJAQlChRArq6uhg/fjw8PT1RsmRJxMfHY/bs2Zg1axaXDyXSQgwGiIg+goxx7KGhoahcuTLGjRsHPT097iPwjvDwcLi7u+PPP//E48ePP2gVJsqK+zIQkaoYDBARFRAGAlndvXsXrq6u+OWXX7B06VJUr14d48ePx+PHj7F3715s2bIFa9asQZcuXTRd1SKB+zIQkboY/hMRFRAGAlnltApTxhEcHIwLFy6gS5cuHA6kAu7LQETqYjBAREQak7EKk4uLCwBwFaZ8kBE0ZfzLfRmIKDdcg42IiDQqYxUmAFiyZAlu3LgBAPIqTPwGO3vcl4GI8gODASIi0rgWLVrgzJkzuHLlCnr27Ilu3bqhffv2ePr0Kb7//ntNV6/Q4b4MRJRfOIGYiIgKDa7C9H4//fQThgwZgkWLFsHDwwOmpqbytdDQUHTp0gXe3t5K+zIAQHx8vNK+DGxPIgIYDBARUSHHG9e3uC8DEeU3TiAmIqJCjYHAW2lpaUhJScHXX3+NkJCQD9qXgYEAEWXGwYJERERFRFJSEp4+fYpLly5h4sSJ6NixI/bv348jR47A3d0dU6ZMwcGDBzVdTSIqQtgzQEREVERwXwYiym8MBoiIiIoI7stARPmNw4SIiIiKEO7LQET5iasJERERFUFnz55Fv379ULJkSdSoUQPPnz9HfHw8rly5Aj09dvwTkWoYDBARERVR3JeBiPKKwQAREVExwkCAiNTBYICIiIiISEtxAjERERERkZZiMEBEREREpKUYDBARERERaSkGA0REREREWorBABERERGRlmIwQERERESkpRgMEBERERFpKQYDRESULzZt2gRzc3NNV4OIiNTAYICISAt4enpCoVDIR9myZdGxY0dcu3Yt356jT58+uH37dr6Vpw4GIkREH4bBABGRlujYsSOePn2Kp0+f4sSJE9DT00PXrl3zrXxDQ0NYWlrmW3lERPTxMRggItISBgYGsLa2hrW1NRo0aIBJkyYhPDwcz58/l/NMnDgRNWrUgJGREapWrYqpU6ciJSVFvn716lW0adMGJiYmMDU1RaNGjXDp0iUAWb+dzy3vu4QQmD59OipWrAgDAwOUL18eI0eOlK8nJSVh3LhxqFChAkqVKoWmTZvi9OnTAIDTp0/j22+/RUxMjNzzMX369PxrOCKiYkxP0xUgIqKCFx8fj61bt8Le3h5ly5aV001MTLBp0yaUL18ewcHBGDRoEExMTDBhwgQAwNdffw0nJyf4+/tDV1cXQUFB0NfXz/Y51Mm7Z88eLF26FNu3b0edOnUQERGBq1evytdHjBiBkJAQbN++HeXLl8e+ffvQsWNHBAcHw8XFBcuWLYOvry9CQ0MBAMbGxvnVVERExRqDASIiLXHgwAH5Jvn169ewsbHBgQMHoKPztpN4ypQp8s+VK1fGuHHjsH37djkYCAsLw/jx41GrVi0AQPXq1XN8PnXzWltbo3379tDX10fFihXRpEkT+drGjRsRFhaG8uXLAwDGjRuHI0eOYOPGjZg7dy7MzMygUChgbW39IU1DRKS1OEyIiEhLtGnTBkFBQQgKCkJgYCDc3NzQqVMnPHr0SM6zY8cOuLq6wtraGsbGxpgyZQrCwsLk62PHjsXAgQPRvn17zJs3D/fu3cvx+dTJ26tXLyQkJKBq1aoYNGgQ9u3bh9TUVABAcHAw0tLSUKNGDRgbG8vHmTNnci2TiIjej8EAEZGWKFWqFOzt7WFvb4/GjRtj3bp1eP36NdauXQsACAgIwNdff43OnTvjwIED+OeffzB58mQkJyfLZUyfPh03btxAly5dcPLkSdSuXRv79u3L9vnUyWtnZ4fQ0FD88MMPMDQ0xHfffYeWLVsiJSUF8fHx0NXVxeXLl+VgJigoCDdv3sTy5cvzv6GIiLQIhwkREWkphUIBHR0dJCQkAAAuXLiASpUqYfLkyXKezL0GGWrUqIEaNWpgzJgx6Nu3LzZu3IjPP/882+dQJ6+hoSHc3d3h7u6O4cOHo1atWggODoaTkxPS0tLw7NkztGjRItvHlihRAmlpaeo2ARGR1mMwQESkJZKSkhAREQEAePXqFVauXIn4+Hi4u7sDkMb0h4WFYfv27WjcuDEOHjyo9E1+QkICxo8fj549e6JKlSp4/PgxLl68iC+++CLLc6mTF5BWIkpLS0PTpk1hZGSErVu3wtDQEJUqVULZsmXx9ddfo1+/fli8eDGcnJzw/PlznDhxAo6OjujSpQsqV66M+Ph4nDhxAvXr14eRkdH/t3OHOorFABhG/wkOhUFcgwRHcCRgsBiCAnc1AkPQCCyvgOABkBgEGnOfhWdg3JhRm+zuJNNzdJs2dV/aNO12+x+cIsDv4pkQQCHu93uqqkpVVRmPx2maJtfrNbPZLEmyWCyy2+2y3W4zGo3yfD5zOBy+5rdarbxer9R1nX6/n9Vqlfl8nuPx+G2tPxmbJJ1OJ+fzOdPpNMPhMI/HI7fb7euno8vlkrqus9/vMxgMslwu0zRNer1ekmQymWSz2WS9Xqfb7eZ0Ov3l0wP4nT7e7/f7pzcBAAD8f24GAACgUGIAAAAKJQYAAKBQYgAAAAolBgAAoFBiAAAACiUGAACgUGIAAAAKJQYAAKBQYgAAAAolBgAAoFCfmoGSiLPPfgIAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "mf_energies = list()\n",
    "mf_times = list()\n",
    "\n",
    "# Perform a Mean-Field calculation for each basis set.\n",
    "for bs in basis_sets:\n",
    "\n",
    "    # Measure execution time.\n",
    "    start = time.time()\n",
    "    scan_basis_mol = SecondQuantizedMolecule(o2_xyz, q=0, spin=2, basis=bs)\n",
    "    end = time.time()\n",
    "\n",
    "    mf_energies.append(scan_basis_mol.mf_energy)\n",
    "    mf_times.append(end-start)\n",
    "\n",
    "# Create the matplotlib figure.\n",
    "fig, ax = plt.subplots(figsize=(8,5))\n",
    "\n",
    "# Plot the energies.\n",
    "ax.set_xticks(range(len(basis_sets)), basis_sets, rotation=45)\n",
    "ax.set_xlabel(\"Basis set\")\n",
    "ax.set_ylabel(\"Energy / Hartree\", color=\"b\")\n",
    "ax.scatter(range(len(basis_sets)), mf_energies, marker=\"o\", s=50, color=\"b\")\n",
    "\n",
    "# Plot the time to solution.\n",
    "ax_time = ax.twinx()\n",
    "ax_time.scatter(range(len(basis_sets)), mf_times, marker=\"s\", s=50, color=\"r\")\n",
    "ax_time.set_ylabel(\"Time to solution / s\", color=\"r\", rotation=270, va=\"bottom\")\n",
    "\n",
    "# Show the graph.\n",
    "plt.tick_params(axis=\"both\", direction=\"in\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "sd8GEceHMfWB"
   },
   "source": [
    "We see in the previous graph that increasing the basis complexity brings the energy eigenvalue closer to the true Hartree-Fock energy. This trend can be generalized to other methods as well, like [CCSD](https://en.wikipedia.org/wiki/Coupled_cluster), [MP2](https://en.wikipedia.org/wiki/M%C3%B8ller%E2%80%93Plesset_perturbation_theory), [FCI](https://en.wikipedia.org/wiki/Full_configuration_interaction), (...) although the time to solution would be greater for those post Hartree-Fock approaches. The takeaway is that even if an infinite basis is required to converge to the true energy, there are diminishing returns at some point, especially when considering the time required for the calculation. In our case, `6-311G(d,p)` appears as a good choice in term of accuracy and compute time.\n",
    "\n",
    "## Consideration of Correlation Energy - Potential of Quantum Computing\n",
    "<a class=\"anchor\" id=\"3\"></a>\n",
    "\n",
    "Post Hatree-Fock methods are meant to consider the electronic correlation energy. The correlation energy is defined as\n",
    "\n",
    "$$\n",
    "E_{\\text{corr}} = E - E_{\\text{HF}}\n",
    "$$\n",
    "\n",
    "In which $E$ and $E_{\\text{HF}}$ are the energies computed from exact diagonalization of the molecular Hamiltonian within the Born-Oppenheimer approximation (i.e. Full Configuration Interaction (FCI)) and the Hartree-Fock method, respectively. Since it is computationally prohibitive to solve the FCI problem, with the exception of very small systems, it is typical to introduce approximations to include correlation effects in electronic structure calculations. The FCI problem considers the diagonalization of the Hamiltonian in a basis of all possible excited determinants, that is all possible 1-, 2-, 3-, …, N-body determinants for an N-electron system. Hence, the FCI wavefunction is a linear combination of all possible determinants for an N-electron system.  Approximate methods, only include a subset of the excited determinants in the approximate wavefunction in order to capture the most important correlation effects (e.g. many-body perturbation theory, coupled-cluster theory). \n",
    "\n",
    "This correlation energy can be split into two categories: static and dynamical correlation. To consider the former, one should use multi-configuration approaches. The latter can be tackled with post Hartree-Fock methods with a single reference state, like Configuration Interaction (CI), Coupled-Cluster (CC) or perturbative methods (MP2, MP3, ...). Going beyond this explanation is out of scope for this notebook, therefore we encourage the user to read on the subject in [Szabo & Ostlund](https://books.google.ca/books/about/Modern_Quantum_Chemistry.html?id=6mV9gYzEkgIC&redir_esc=y).\n",
    "\n",
    "We saw in the previous section that the `6-311G(d,p)` is a good compromise to compute the O$_2$ ground state. In this basis, the dioxygen molecule is described by 72 bases (spin orbitals) and 16 electrons. To get the highest fidelity versus the true wavefunction with the current best method (Full Configuration Interaction), one would need ${72 \\choose 16} \\approx 3 \\times 10^{13}$ configurations. Even if we somehow managed to represent each configuration as a single byte of data, 4.11 petabytes would still be needed to store the full wavefunction of the dioxygen molecule. This effect is amplified when dealing with relevant industrial molecules like organometallic complexes or proteins. Going through the classical formalism of quantum chemistry is therefore not a viable solution.\n",
    "\n",
    "With the occupation transformation (Jordan-Wigner [[Zeitschrift für Physik, 47, p. 631–651 (1928)](https://doi.org/10.1007/BF01331938)]), this problem vanishes. By leveraging the quantum superposition principle in a straightforward way, it is possible to map this O$_2$ triplet ground state to 72 qubits (empty and filled spin orbitals are mapped to the $|0\\rangle$ and $|1\\rangle$ states, respectively). This fact is behind the emerging enthusiasm of the quantum computing community for chemistry applications. We should note that even if the space limitation is mitigated, no consensus is established for the time complexity, correlated to the number of gates and measurement required to approximate the molecular wavefunction [[arXiv:2009.12472](https://doi.org/10.48550/arXiv.2009.12472)].\n",
    "\n",
    "## Resource Reduction<a class=\"anchor\" id=\"4\"></a>\n",
    "\n",
    "Resource reduction refers to lowering the computational requirements to run an algorithm. This applies to both classical and quantum computing. In conventional computing, this translates to less bits required to store information and less operations to do on these bits. For the gate-model quantum computing, this translates to reducing the number of qubits and the number of gates in quantum circuits.\n",
    "\n",
    "Our starting point is the O$_2$ molecule in the `6-311G(d,p)` basis. We focus on reducing resource requirements before choosing an algorithm to compute the ground state. For approaches relying on choosing a different ansatz or fragmenting a molecule into smaller subproblems (problem decomposition), the user can refer to the relevant notebooks in the [Tangelo-Examples repository](https://github.com/goodchemistryco/Tangelo-Examples)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "36 active molecular orbitals\n",
      "16 active electrons\n"
     ]
    }
   ],
   "source": [
    "o2_6311gdp = SecondQuantizedMolecule(o2_xyz, q=0, spin=2,\n",
    "                                     basis=\"6-311G(d,p)\", \n",
    "                                     frozen_orbitals=None)\n",
    "print(f\"{o2_6311gdp.n_active_mos} active molecular orbitals\")\n",
    "print(f\"{o2_6311gdp.n_active_electrons} active electrons\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WENmx_HFMfWC"
   },
   "source": [
    "The O$_2$ molecule in the 6-311G(d,p) basis would require 72 qubits and inherently a quantum circuit that is considered too deep for current architectures (NISQ devices) at the time of writing. Decreasing the number of qubits required to encode the problem provides the first indirect way to reduce those requirements. The two methods discussed below are the reduction of the active space and the consideration of molecular symmetries.\n",
    "\n",
    "### Active Space Selection<a class=\"anchor\" id=\"41\"></a>\n",
    "\n",
    "Active space selection is a way to reduce the number of qubits by ignoring some molecular orbitals in the post-Hartree-Fock calculation. Applying this method is trivial, but choosing which molecular orbitals to freeze is not. While there is no general algorithm to perform that operation, one can evaluate which orbitals to freeze by first printing the molecular orbital occupancies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "NPArrayWithTag([2., 2., 2., 2., 2., 2., 2., 1., 1., 0., 0., 0., 0., 0.,\n",
       "                0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,\n",
       "                0., 0., 0., 0., 0., 0., 0., 0.])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "o2_6311gdp.mo_occ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "JWSVc370MfWD"
   },
   "source": [
    "One way of picturing this is to only consider frontier orbitals and their neighbors. The correlation lost on the total energy when only considering a subset of the full active space might be small. For example, freezing low-lying occupied molecular orbitals is known as the *frozen-core approximation*, and can be applied because core orbitals do not mix with valence orbitals. Also, there are algorithms to freeze virtual orbitals, like the Frozen Natural Orbitals (FNO) truncation method [[Collect. Czech. Chem. Commun. 2005, 70, 837-850](https://doi.org/10.1135/cccc20050837) and [J. Chem. Phys. 2008, 128, 164101](https://doi.org/10.1063/1.2902285)]. Despite the fact that choosing a suitable active space can be tailored a problem, it can help reducing quantum resources *a priori*.\n",
    "\n",
    "The code cell below results in the consideration of molecular orbitals near the HOMO-LUMO level. The `get_orbitals_excluding_homo_lumo` function helps users to get the labels of these orbitals. In our case, the active space is defined from the HOMO-3 to LUMO+3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8 active molecular orbitals\n",
      "6 active electrons\n"
     ]
    }
   ],
   "source": [
    "from tangelo.toolboxes.molecular_computation.frozen_orbitals import get_orbitals_excluding_homo_lumo\n",
    "\n",
    "# Selecting HOMO-3 to LUMO+3 orbitals.\n",
    "frozen = get_orbitals_excluding_homo_lumo(o2_6311gdp, homo_minus_n=3, lumo_plus_n=3)\n",
    "\n",
    "o2_6311gdp_frozen = SecondQuantizedMolecule(o2_xyz, q=0, spin=2,\n",
    "                                            basis=\"6-311G(d,p)\",\n",
    "                                            frozen_orbitals=frozen)\n",
    "print(f\"{o2_6311gdp_frozen.n_active_mos} active molecular orbitals\")\n",
    "print(f\"{o2_6311gdp_frozen.n_active_electrons} active electrons\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "cnbZ7LW-MfWE"
   },
   "source": [
    "Following this approach, we decrease the number of considered orbitals from 36 to 8. With conventional fermion-to-qubit mappings, the number of qubits can be reduced from 72 to 16. While this is an impressive reduction, configurations involving orbitals below HOMO-4 and above LUMO+4 are not considered here. This results in an energy further away from the true energy as there is nothing telling us that those configurations are irrelevant *a priori*. In the next code cell, we consider the CCSD calculation of O$_2$ with and without those orbitals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy difference: 226.51012414561814 kcal/mol.\n"
     ]
    }
   ],
   "source": [
    "from tangelo.algorithms.classical import CCSDSolver\n",
    "\n",
    "# CCSD energy with all the orbitals.\n",
    "e_o2_6311gdp = CCSDSolver(o2_6311gdp).simulate()\n",
    "\n",
    "# CCSD energy with HOMO-3 to LUMO+3 orbitals.\n",
    "e_o2_6311gdp_frozen = CCSDSolver(o2_6311gdp_frozen).simulate()\n",
    "\n",
    "print(f\"Energy difference: {abs(e_o2_6311gdp-e_o2_6311gdp_frozen)*627.5} kcal/mol.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "FruSKYNbMfWE"
   },
   "source": [
    "This choice results in an energy difference of over 200 kcal/mol, well above [chemical accuracy](https://en.wikipedia.org/wiki/Computational_chemistry) (1 kcal/mol). In this situation, the lower bound with the reduced active space is higher than the CCSD ground state. The takeaway is to be careful when choosing an active space for a molecule.\n",
    "\n",
    "Several approaches have been developed by the community to truncate the active space of a molecule for a given basis set. They follow specific criteria ensuring a controlled impact on accuracy. A built-in method in Tangelo is the [Frozen Natural Orbitals](https://doi.org/10.1063/5.0054647) (FNO) algorithm, that aims at truncating the virtual orbital space. Combined with the frozen core approximation, it constitutes a robust strategy to reduce the quantum resources needed. Here is an example code snippet that performs this operation:\n",
    "\n",
    "```python\n",
    "from tangelo.toolboxes.molecular_computation.fno import FNO\n",
    "\n",
    "# FNO is not compatible with ROHF, so we use UHF.\n",
    "o2_6311gdp_uhf = SecondQuantizedMolecule(o2_xyz, q=0, spin=2, uhf=True, \n",
    "                                         basis=\"6-311G(d,p)\", \n",
    "                                         frozen_orbitals=\"frozen_core\")\n",
    "\n",
    "# The threshold parameter corresponds to %FNO where the truncation happens in the virtual space.\n",
    "fno = FNO(o2_6311gdp_uhf, threshold=0.9)\n",
    "frozen_orbitals_uhf = fno.get_frozen_indices()\n",
    "\n",
    "# The molecular object attributes have to be updated.\n",
    "o2_6311gdp_uhf_frozen = o2_6311gdp_uhf.freeze_mos(frozen_orbitals_uhf, inplace=False)\n",
    "o2_6311gdp_uhf_frozen.mo_coeff = fno.mo_coeff\n",
    "```\n",
    "\n",
    "Beyond the built-in options available in Tangelo, such as frozen core or FNO, users can use external tools to compute frozen orbitals or integrals such as [SCINE autoCAS](https://github.com/qcscine/autocas), as described in this [example script](https://github.com/goodchemistryco/Tangelo-Examples/blob/main/examples/chemistry/autocas.py).\n",
    "\n",
    "### Considering Symmetries<a class=\"anchor\" id=\"42\"></a>\n",
    "\n",
    "Molecular symmetries are reflected in the Hamiltonian structure. The number of electrons in the molecule should be consistent before and after the ansatz has been applied onto the initial state. Electronic configurations not respecting this symmetry can therefore be neglected without loss of accuracy. For our purpose in quantum chemistry, it leads to a smaller number of qubits.\n",
    "\n",
    "The code cell below shows an example of Z2-symmetries qubit tapering with Tangelo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of qubits: 16\n",
      "Number of qubits after tapering: 12\n",
      "4 Z2 symmetries detected\n"
     ]
    }
   ],
   "source": [
    "from tangelo.toolboxes.qubit_mappings.mapping_transform import fermion_to_qubit_mapping\n",
    "from tangelo.toolboxes.operators import count_qubits\n",
    "from tangelo.toolboxes.operators.taper_qubits import QubitTapering\n",
    "\n",
    "hamiltonian = fermion_to_qubit_mapping(o2_6311gdp_frozen.fermionic_hamiltonian, \"JW\")\n",
    "n_qubits = count_qubits(hamiltonian)\n",
    "taper = QubitTapering(hamiltonian, n_qubits,\n",
    "                      n_electrons=o2_6311gdp_frozen.n_active_electrons)\n",
    "hamiltonian_taper = taper.z2_tapered_op\n",
    "n_qubits_tapered = count_qubits(hamiltonian_taper)\n",
    "\n",
    "print(f\"Number of qubits: {n_qubits}\")\n",
    "print(f\"Number of qubits after tapering: {n_qubits_tapered}\")\n",
    "print(f\"{n_qubits - n_qubits_tapered} Z2 symmetries detected\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vYkFnG6CMfWE"
   },
   "source": [
    "## Advanced Usage of PySCF<a class=\"anchor\" id=\"5\"></a>\n",
    "\n",
    "[PySCF](https://pyscf.org/) is one of the classical chemistry backends supported by Tangelo ([Psi4](https://psicode.org/) is also supported). Advanced usage of this open-source quantum chemistry package enables the use of many analysis tools. We show a couple of examples where advanced users can extract more information from a PySCF calculation.\n",
    "\n",
    "A `SecondQuantizedMolecule` object can be exported as a `pyscf.gto.Mole` object. From that point, anything implemented in `PySCF` can be applied to this data structure. The next line shows how to do this operation.\n",
    "\n",
    "```python\n",
    "from tangelo.toolboxes.molecular_computation.integral_solver_pyscf import mol_to_pyscf\n",
    "pyscf_mol = mol_to_pyscf(sec_mol, sec_mol.basis, sec_mol.symmetry, sec_mol.ecp)\n",
    "```\n",
    "\n",
    "### Extracting Mean-Field Quantities<a class=\"anchor\" id=\"51\"></a>\n",
    "\n",
    "The first trick is the extraction of the molecular orbital eigenvalues and coefficients. This can be done with the `analyze` method, and it can be called from the mean-field object to print relevant molecular properties. The output features many sections, described in the table below.\n",
    "\n",
    "| Section | Decription |\n",
    "| :- | :- |\n",
    "| `SCF Summaries` | Self-consistent field energy contributions |\n",
    "| `MO energy` | Molecular orbital eigenvalues |\n",
    "| `MO coefficients` | Linear Combination of Atomic Orbitals (LCAO) coefficients |\n",
    "| `Mulliken pop` | Mulliken electron population |\n",
    "| `Mulliken atomic charges` | Partial charges on each atom and dipole moment |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "**** SCF Summaries ****\n",
      "Total Energy =                        -147.632274661317865\n",
      "Nuclear Repulsion Energy =              27.989538428826449\n",
      "One-electron Energy =                 -260.148728081300931\n",
      "Two-electron Energy =                   84.526914991156602\n",
      "**** MO energy ****\n",
      "                Roothaan           | alpha              | beta\n",
      "MO #1   energy= -20.4249798601204  | -20.4400185053651  | -20.4099412148758  occ= 2\n",
      "MO #2   energy= -20.4243648974171  | -20.4389300607359  | -20.4097997340983  occ= 2\n",
      "MO #3   energy= -1.55251062936318  | -1.61345571763428  | -1.49156554109208  occ= 2\n",
      "MO #4   energy= -0.993242469553747 | -1.09010934626379  | -0.896375592843713 occ= 2\n",
      "MO #5   energy= -0.58154450565406  | -0.616775237024776 | -0.546313774283346 occ= 2\n",
      "MO #6   energy= -0.579365900059497 | -0.71407161146602  | -0.444660188652973 occ= 2\n",
      "MO #7   energy= -0.579365900059497 | -0.71407161146602  | -0.444660188652973 occ= 2\n",
      "MO #8   energy= -0.0715888741631396 | -0.412074886893248 | 0.26889713856697   occ= 1\n",
      "MO #9   energy= -0.0715888741631381 | -0.412074886893249 | 0.268897138566972  occ= 1\n",
      "MO #10  energy= 0.730335506725122  | 0.689507257564998  | 0.771163755885239  occ= 0\n",
      " ** MO coefficients (expansion on meta-Lowdin AOs) **\n",
      "               #1        #2        #3        #4        #5       \n",
      "0 O 1s          -0.70708   0.70711  -0.00060  -0.00323   0.00005\n",
      "0 O 2s          -0.00483   0.00057   0.65237   0.65105  -0.27278\n",
      "0 O 2px         -0.00000  -0.00000  -0.00000   0.00000  -0.00000\n",
      "0 O 2py          0.00000   0.00000   0.00000   0.00000   0.00000\n",
      "0 O 2pz         -0.00311   0.00019   0.27278  -0.27591   0.65237\n",
      "1 O 1s           0.70708   0.70711  -0.00060   0.00323   0.00005\n",
      "1 O 2s           0.00483   0.00057   0.65237  -0.65105  -0.27278\n",
      "1 O 2px         -0.00000  -0.00000  -0.00000  -0.00000  -0.00000\n",
      "1 O 2py          0.00000   0.00000   0.00000   0.00000   0.00000\n",
      "1 O 2pz         -0.00311  -0.00019  -0.27278  -0.27591  -0.65237\n",
      "               #6        #7        #8        #9        #10      \n",
      "0 O 1s          -0.00000  -0.00000  -0.00000   0.00000  -0.00475\n",
      "0 O 2s           0.00000   0.00000   0.00000   0.00000   0.27589\n",
      "0 O 2px         -0.24960   0.66159  -0.10780  -0.69884   0.00000\n",
      "0 O 2py          0.66159   0.24960   0.69884  -0.10780  -0.00000\n",
      "0 O 2pz         -0.00000  -0.00000  -0.00000  -0.00000   0.65105\n",
      "1 O 1s          -0.00000  -0.00000  -0.00000  -0.00000   0.00475\n",
      "1 O 2s           0.00000   0.00000   0.00000   0.00000  -0.27589\n",
      "1 O 2px         -0.24960   0.66159   0.10780   0.69884  -0.00000\n",
      "1 O 2py          0.66159   0.24960  -0.69884   0.10780   0.00000\n",
      "1 O 2pz          0.00000   0.00000   0.00000   0.00000   0.65105\n",
      " ** Mulliken pop on meta-lowdin orthogonal AOs  **\n",
      " ** Mulliken pop  **\n",
      "pop of  0 O 1s            1.99995\n",
      "pop of  0 O 2s            1.84777\n",
      "pop of  0 O 2px           1.50000\n",
      "pop of  0 O 2py           1.50000\n",
      "pop of  0 O 2pz           1.15227\n",
      "pop of  1 O 1s            1.99995\n",
      "pop of  1 O 2s            1.84777\n",
      "pop of  1 O 2px           1.50000\n",
      "pop of  1 O 2py           1.50000\n",
      "pop of  1 O 2pz           1.15227\n",
      " ** Mulliken atomic charges  **\n",
      "charge of    0O =      0.00000\n",
      "charge of    1O =     -0.00000\n",
      "Dipole moment(X, Y, Z, Debye):  0.00000, -0.00000, -0.00000\n"
     ]
    }
   ],
   "source": [
    "o2_sto3g.mean_field.analyze(verbose=5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RLP_CotHMfWF"
   },
   "source": [
    "### Outputting Cube Files<a class=\"anchor\" id=\"52\"></a>\n",
    "\n",
    "Visualization of molecular orbitals is useful for analysis purposes or making figures. Those few lines create a `cube` file, compatible with many chemical visualization softwares ([GaussView](https://gaussian.com/gaussview6/), [avogadro](https://avogadro.cc/), [VMD](https://www.ks.uiuc.edu/Research/vmd/), ...), for the `MO #10`.\n",
    "\n",
    "```python\n",
    "from pyscf import tools\n",
    "tools.cubegen.orbital(pyscf_mol, f\"MO10.cube\", sec_mol.mean_field.mo_coeff[:,9]) # Python indices start with 0.\n",
    "```\n",
    "\n",
    "## Closing Words<a class=\"anchor\" id=\"6\"></a>\n",
    "\n",
    "This notebook glanced over the first steps of modelling a chemical system for quantum computing exploration. One of the challenges of quantum chemistry is the ever-present trade-off between computational resources requirements and accuracy of calculations. We emphasized how the choice of a basis set and frozen orbitals can greatly affect this trade-off.\n",
    "\n",
    "Quantum computing is a nascent technology, and the capabilities of current devices are limited. Efficient modelling of the system can be therefore crucial to sucessfully design experiments on quantum computers, and even simulators.\n",
    "\n",
    "In our [VQE notebook tutorials](https://github.com/goodchemistryco/Tangelo-Examples/blob/main/contents.rst), we show how the modelling choices can greatly impact the quantum computational resources required by this quantum algorithm, and how the algorithm options themselves can further affect both resource requirements and accuracy.\n",
    "\n",
    "\n",
    "## References<a class=\"anchor\" id=\"7\"></a>\n",
    "\n",
    "- McQuarrie, D. A. [Quantum Chemistry](https://books.google.ca/books/about/Quantum_Chemistry.html?id=zzxLTIljQB4C&redir_esc=y), University Science Books 2nd ed. (2007)\n",
    "- Engel, T. [Quantum chemistry & spectroscopy](https://books.google.ca/books/about/Quantum_Chemistry_and_Spectroscopy.html?id=JhUqngEACAAJ&redir_esc=y). Pearson 3rd ed. (2012)\n",
    "- Szabo, A. & Ostlund, N. S. [Modern quantum chemistry: Introduction to advanced electronic structure theory](https://books.google.ca/books/about/Modern_Quantum_Chemistry.html?id=6mV9gYzEkgIC&redir_esc=y). Mineola, Dover Publications 1st ed. (1996).\n",
    "- [Computational_Quantum_Chemistry](https://chem.libretexts.org/Courses/Pacific_Union_College/Quantum_Chemistry/11%3A_Computational_Quantum_Chemistry), Chemistry LibreTexts, visited on Aug 17, 2022."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "Review_quantum_chemistry.ipynb",
   "provenance": [
    {
     "file_id": "https://github.com/AlexandreF-1qbit/Tangelo/blob/notebook_chem/examples/quantum_chemistry.ipynb",
     "timestamp": 1657994939201
    }
   ]
  },
  "gpuClass": "standard",
  "kernelspec": {
   "display_name": "qsdk",
   "language": "python",
   "name": "qsdk"
  },
  "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.10.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "95050af2697fca56ed7491a4fb0b04c1282c0de0a7e0a7cacd318a8297b0b1d8"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}