biphenyl_mm_energy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Biphenyl MM Energies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adapted from: https://github.com/choderalab/fragmenter_data/blob/master/wbo-manuscript-figures/biphenyl_example/biphenyl_mm_energy.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:100% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bb718514771d4afab4bc007c1a217a1b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from functools import partial\n",
    "\n",
    "import numpy as np\n",
    "from scipy.optimize import minimize\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from openforcefield.topology import Molecule, Topology\n",
    "from openforcefield.typing.engines.smirnoff import ForceField\n",
    "from openeye import oechem, oedepict\n",
    "\n",
    "from simtk.openmm import app\n",
    "from simtk import openmm as mm\n",
    "from simtk import unit\n",
    "from simtk import openmm, unit\n",
    "\n",
    "from fragmenter import torsions, chemi\n",
    "\n",
    "from autograd import grad\n",
    "from autograd import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Function definitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# From Josh - constrained minimzation \n",
    "\n",
    "def compute_dihedral(xyz, indices):\n",
    "    \"\"\"Copied from mdtraj, except without calls to a non-python distance library,\n",
    "    and assuming we dont have to deal with periodic stuff\"\"\"\n",
    "    \n",
    "    a,b,c,d = indices\n",
    "    b1 = xyz[b] - xyz[a]\n",
    "    b2 = xyz[c] - xyz[b]\n",
    "    b3 = xyz[d] - xyz[c]\n",
    "    \n",
    "    c1 = np.cross(b2, b3) # bc x cd\n",
    "    c2 = np.cross(b1, b2) # ab x bc\n",
    "\n",
    "    p1 = np.sum(b1 * c1)\n",
    "    p1 *= np.sum(b2 * b2) ** 0.5\n",
    "    p2 = np.sum(c1 * c2)\n",
    "\n",
    "    return np.arctan2(p1, p2)\n",
    "\n",
    "def unflatten(flat_xyz):\n",
    "    N = int(len(flat_xyz) / 3)\n",
    "    return np.reshape(flat_xyz, (N, 3))\n",
    "\n",
    "def set_positions(xyz, simulation):\n",
    "    simulation.context.setPositions(xyz)\n",
    "        \n",
    "def get_energy(xyz, simulation):\n",
    "    set_positions(xyz, simulation)\n",
    "    U = simulation.context.getState(getEnergy=True).getPotentialEnergy()\n",
    "    return U/U.unit\n",
    "\n",
    "def get_forces(xyz, simulation):\n",
    "    set_positions(xyz, simulation)\n",
    "    F =  simulation.context.getState(getForces=True).getForces(asNumpy=True)\n",
    "    return F/F.unit\n",
    "\n",
    "def fxn_to_minimize(flat_xyz, simulation):\n",
    "    return get_energy(unflatten(flat_xyz), simulation)\n",
    "\n",
    "def jacobian_of_fxn(flat_xyz, simulation):\n",
    "    return - get_forces(unflatten(flat_xyz), simulation)\n",
    "\n",
    "\n",
    "def dihedral_constraint(flat_xyz, torsion_indices, target_dihedral=0.0):\n",
    "    \"\"\"equality constraint: want the output of this function to be 0\"\"\"\n",
    "    xyz = unflatten(flat_xyz)\n",
    "    return compute_dihedral(xyz, torsion_indices) - target_dihedral\n",
    "\n",
    "def form_eq_constraint(torsion_indices, target_dihedral):\n",
    "    \"\"\"adapted from scipy documentation\n",
    "    https://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html#sequential-least-squares-programming-slsqp-algorithm-method-slsqp\n",
    "    \"\"\"\n",
    "    \n",
    "    fun = partial(dihedral_constraint, torsion_indices=torsion_indices, target_dihedral=target_dihedral)\n",
    "\n",
    "    eq_cons = {'type': 'eq',\n",
    "               'fun' : fun,\n",
    "               'jac' : grad(fun)}\n",
    "    return eq_cons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Molecule generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create openeye molecules and add conformers with torsion angle at 15 degree intervals\n",
    "smiles = ['c1cc(ccc1c2ccncc2)O', 'c1cc(ccc1c2cc[nH+]cc2)O', 'c1cc(ccc1c2ccncc2)[O-]', 'c1cc(ccc1c2cc[nH+]cc2)[O-]']\n",
    "names = ['neutral', 'cation', 'anion', 'zwitterion']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create openeye molecules\n",
    "oemols = [chemi.smiles_to_oemol(sm, add_atom_map=True) for sm in smiles]\n",
    "for i, mol in enumerate(oemols):\n",
    "    mol.SetTitle(names[i])\n",
    "    \n",
    "conf_mols = []\n",
    "dihedrals = []\n",
    "# Find torsion to constrain\n",
    "for mol in oemols:\n",
    "    for bond in mol.GetBonds():\n",
    "        if bond.IsRotor():\n",
    "            rot_bond = bond\n",
    "            break\n",
    "    bond_map = (rot_bond.GetBgn().GetMapIdx(), rot_bond.GetEnd().GetMapIdx())\n",
    "    dih = torsions.find_torsion_around_bond(molecule=mol, bond=bond_map)\n",
    "    conformers = chemi.generate_grid_conformers(mol, dihedrals=[dih], intervals=[15])\n",
    "    conf_mols.append(conformers)\n",
    "    dihedrals.append(dih)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, 8, 9, 4], [0, 8, 9, 4], [0, 8, 9, 4], [0, 8, 9, 4]]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dihedrals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Molecule with name 'neutral' and SMILES '[H]c1c(c(c(c(c1c2c(c(nc(c2[H])[H])[H])[H])[H])[H])O[H])[H]',\n",
       " Molecule with name 'cation' and SMILES '[H]c1c(c(c(c(c1c2c(c([n+](c(c2[H])[H])[H])[H])[H])[H])[H])O[H])[H]',\n",
       " Molecule with name 'anion' and SMILES '[H]c1c(c(c(c(c1c2c(c(nc(c2[H])[H])[H])[H])[H])[H])[O-])[H]',\n",
       " Molecule with name 'zwitterion' and SMILES '[H]c1c(c(c(c(c1c2c(c([n+](c(c2[H])[H])[H])[H])[H])[H])[H])[O-])[H]']"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "off_molecules = [Molecule.from_openeye(mol) for mol in conf_mols]\n",
    "off_molecules"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "topologies = [Topology.from_molecules(mol) for mol in off_molecules]\n",
    "forcefield = ForceField('openff-1.0.0.offxml')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Molecule with name 'neutral' and SMILES '[H]c1c(c(c(c(c1c2c(c(nc(c2[H])[H])[H])[H])[H])[H])O[H])[H]'\n",
      "Forces for molecule 0\n",
      "\n",
      "ProperTorsions:\n",
      "atoms:   0  8  9  4  parameter_id: t43  smirks [*:1]~[#6X3:2]-[#6X3:3]~[*:4]\n",
      "\n",
      "\n",
      "Molecule with name 'cation' and SMILES '[H]c1c(c(c(c(c1c2c(c([n+](c(c2[H])[H])[H])[H])[H])[H])[H])O[H])[H]'\n",
      "Forces for molecule 0\n",
      "\n",
      "ProperTorsions:\n",
      "atoms:   0  8  9  4  parameter_id: t43  smirks [*:1]~[#6X3:2]-[#6X3:3]~[*:4]\n",
      "\n",
      "\n",
      "Molecule with name 'anion' and SMILES '[H]c1c(c(c(c(c1c2c(c(nc(c2[H])[H])[H])[H])[H])[H])[O-])[H]'\n",
      "Forces for molecule 0\n",
      "\n",
      "ProperTorsions:\n",
      "atoms:   0  8  9  4  parameter_id: t43  smirks [*:1]~[#6X3:2]-[#6X3:3]~[*:4]\n",
      "\n",
      "\n",
      "Molecule with name 'zwitterion' and SMILES '[H]c1c(c(c(c(c1c2c(c([n+](c(c2[H])[H])[H])[H])[H])[H])[H])[O-])[H]'\n",
      "Forces for molecule 0\n",
      "\n",
      "ProperTorsions:\n",
      "atoms:   0  8  9  4  parameter_id: t43  smirks [*:1]~[#6X3:2]-[#6X3:3]~[*:4]\n"
     ]
    }
   ],
   "source": [
    "for i, top in enumerate(topologies):\n",
    "    print('\\n')\n",
    "    print(off_molecules[i])\n",
    "    molecule_force_list = forcefield.label_molecules(top)\n",
    "    for mol_idx, mol_forces in enumerate(molecule_force_list):\n",
    "        print(f'Forces for molecule {mol_idx}')\n",
    "        for force_tag, force_dict in mol_forces.items():\n",
    "            if force_tag != 'ProperTorsions':\n",
    "                continue\n",
    "            print(f\"\\n{force_tag}:\")\n",
    "            for (atom_indices, parameter) in force_dict.items():\n",
    "                if atom_indices == tuple(dih):\n",
    "                    atomstr=''\n",
    "                    for idx in atom_indices:\n",
    "                        atomstr += '%3s' % idx\n",
    "                    print(\"atoms: %s  parameter_id: %s  smirks %s\" % (atomstr, parameter.id, parameter.smirks) )\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create OMM systems\n",
    "\n",
    "#integrator = mm.VerletIntegrator(1*unit.femtosecond)\n",
    "platform = mm.Platform.getPlatformByName(\"Reference\")\n",
    "\n",
    "systems = [forcefield.create_openmm_system(top) for top in topologies]\n",
    "simulations = [openmm.app.Simulation(topologies[i], systems[i], mm.VerletIntegrator(1*unit.femtosecond)) for i in range(4)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_energies = []\n",
    "all_angles = []\n",
    "for i, mol in enumerate(off_molecules):\n",
    "    energies = []\n",
    "    angles = []\n",
    "    for j in range(24):\n",
    "        simulation = simulations[i]\n",
    "        angles.append(torsions.measure_dihedral_angle(coords=mol.conformers[j]._value, dihedral=dih))\n",
    "        simulation.context.setPositions(mol.conformers[j])\n",
    "        state = simulation.context.getState(getEnergy=True)\n",
    "        energy = state.getPotentialEnergy()\n",
    "        energies.append(energy._value)\n",
    "    energies = np.asanyarray(energies) - min(energies)\n",
    "    all_energies.append(energies)\n",
    "    all_angles.append(angles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABtDElEQVR4nO2dd3hc5ZW432960WhGvctykSsGA6bY9N57CSEkBBJCAglhU2FTNrv5bTabniyEUAIhlJDQe2+mGbDBuMhNlpvKjOpoiqbf7/fHnZFlW7JlW6NR+d7n0aOZO3PvPTM69+jc850ipJQoFAqFYvJgyLUACoVCoRhdlOFXKBSKSYYy/AqFQjHJUIZfoVAoJhnK8CsUCsUkw5RrAYZDcXGxrKury7UYCoVCMa5Yvnx5p5SyZNft48Lw19XVsWzZslyLoVAoFOMKIcTWwbarUI9CoVBMMpThVygUikmGMvwKhUIxyVCGX6FQKCYZyvArFArFJEMZfoVCoZhkTGjDv2VFK6/d8yFbVrTmWhTFJEPpnmIsMy7y+PeHLStaue3af5FKpjBbTdxw92XULajMtViKScCWFa3c/hVd90wWpXuKsceE9fgbl20nlUiBhGQiReOy7bkWSTFJ2LB0K8l4CqlBSumeYgwyYQ3/jIU1GEz6xzOZjMxYWJNjiRSTB9H/yGhWuqcYe0xYw1+3oJIzrl8EwGX/cZq61VaMGpuW6x6+wSj4xt2XKt1TjDkmrOEHmHl0LQBOty3HkigmC53b/WxYug1PWR5aSlJcU5BrkRSK3ZjQht9T7gLA7w3mWBLFZOHDJ1cjDIKTrz0SULqnGJtMaMOfX+LEYBTq4lOMCqmkxkdPrWbOsVOpO0QP7/i9gRxLpVDszoQ2/AajgfySPPy+UK5FUUwC1r7TRKAjzNGXzO+/2+xRTodiDDJh8/gzeMpd+H3q4lNknw8eX0V+iZO5x09DGARGs1HpnmJMMqE9fgBPWZ4K9Siyjt8bZO07mznywoMwmgwYDAJPWR69XnW3qRh7THzDX+7C7w0ipcy1KIoJzEdPrUZqkqMuOqh/m6fcRY/y+BVjkElg+PNJxJKE/ZFci6KYoGiaZOmTq5l5dC3FNZ7+7RmnQ6EYa0x4w1+gUjoVWWbDB1vpaQ1w9CUH77S9oNxFb3sILaXlSDKFYnAmvOFXufyKbPPB4ytxemzMP3n6Tts95S60pEawqy9HkikUgzN5DL9K6VRkgWBnmNVvbuKI8+dhsuycJKecDsVYZcIb/rxCB0aTQV18iqzw0TNr0JIaR18yf7fXPGUZp0PpnmJsMeENv8EgcKuUTkUWkFLy4ROrmXpYFWXTinZ7XXn8irHKhDf8oHte6uJTjDSbljXTsbWHRYN4+wAOtw2zzaR0TzHmmByGv9ylSucVI84Hj63E5rJyyGkzB31dCKF0TzEmmdAtG7b4o2wPxrDPLKH3lQ1omsRgEHvfUaHYC+HeCCtf28jRF8/HYjfv9npzIIo3nMBzcCX+Ld05kFChGJqse/xCCKMQ4lMhxHPp54VCiFeFEBvTv7PSsNwbivNSUw+rO/roPKQa85RCQt0qrU4xMix7di3JeIqjL909zNMWivFsYw8ft4UQF8wnbJ3Q/pViHDIaoZ5vA2sHPL8FeF1KWQ+8nn4+4rSG4mSaNEgB9jmlKrtCMSJIKfnw8VXUHlRO1azS3V7f2B3d8cQg0MrzSSVVEZdi7JBVwy+EqAbOAe4ZsPkC4P704/uBC7Nx7so8C5mojkEIImvb1SKbYkTYurKNtsZOjrp48EXdYCzZ/1gAfQ0+Ap2qjkQxdsi2x/8H4AfAQHenTErZBpD+vbvLBAghviaEWCaEWNbR0bHPJy7Ps3BMlZ5Od1iBjeimTmX4FSPC0sdXYbGbOezs2bu9FktqtITiGNNOx0GxuNI9xZgja4ZfCHEu0C6lXL4/+0sp75JSLpRSLiwpKdkvGWrTs3bzXBbMVpVWpzhwoqEYn760jkPPmoXNadnt9fXdEVISDinLA6C4yAmoXH7F2CKbHv8xwPlCiC3AI8DJQogHAZ8QogIg/bs9WwI4LUYAgnFNFXEpRoRPXlhHPJJk0S4N2UCP/Td09lHmNFPntgJgSDsfSvcUY4msGX4p5a1SymopZR1wBfCGlPIq4Bng6vTbrgaezpYMRiFwmg2E4imVT60YET54fBUV9cXUzi/f7bW2UJyeaJK5xQ5caacjbjJgdVqU7inGFLko4PolcJoQYiNwWvp5VkiEV+Cgk2DUr3qjKw6Y5rU+mht8HH3JfITYvR5kTWcfFqNgRoEdU2wVBlL0htuU7inGHKNi+KWUb0kpz00/7pJSniKlrE//zkp1SyK8gt6N12JLfEagr4v8wj4CHao3umL/Wfr4KsxWEwvPnbvba32JFE3+KLMK7cjIZwQav4IdHz3dy3GXSGX4FWOKCduyIRH8GGQcB+30UYwz34uWkgQ6w7kWTTEOifUlWP78Wg4+rR5HOm4/kPVdETQJ84odJIIfgYzhxEcfxbg8QdUWXDGmmLCG3+w6AjDhoB0NC/l18wC1yKbYPz57ZT3RUHyPi7qVeRYK7GYQeqWu7nSUUlhdTagrTDKRGm2xFYpBmbiG37kAR/nXcaSThuxV+u25MvyK/eGDx1dRUlfAtMOrdnutORgnEE8xt9gBQCK4FAAHHUQoorC2HimhV3n9ijHChDX8AOb8Y3DgA0CotDrFftLW2MmWFa0cffHQi7o2k4FpHhvJaBOJ4PsYLFVp3RPYVF9+xRhjQht+o6Wy3+OPmwxY7GaVVqfYZ175ywcIg6Bs+u7DVsLxFFv8UWYX2TEaBNGOf4Aw4yi/oV/3TKqIS7EfSCnZvKKF1+75kC0rWkf02BPa8AtTIWaRwizi/bn8qlGbYl9o+rSFFS9vQGqS+7/77G4X4NquPiQwt9iBlgoS7X4Ka8FZmJwH7bjbdOl3mz3ewGiLrxjHfPzMGv70xUd44U/v8ufrHh1R4z+xDb8QmGyVOETvDsOvvC7FPtDwdlP/41QiReOy7f3PtfSibo3LgttqItb1FGgR7CVfSN9t6j2mooAj36Z0T7FPbFy6DQApd9e9A2VCG34gHWvtIBhP4SnLUwtsin2iqNoNgDCA0WxkxsKa/te29cYIJzTmljiRUiPS8TAm5wJMjnkIgx2TOQ+riBBK6C1DlO4p9oW8Qj1ZQBjEbrp3oEz4CRFGSyUOuZ2WxEGUl7sIdIRIJVIYzcZci6YYB9hdes+dY69YwGFnzaFuQWX/a2s6+3CaDUxxW0kElqDFt+Os/Fb/60ZLJc5UN8G4R28ZosKMin3AYDRgMArOvPEY6o+o2Un3DpQJbfgjTcsIf7IRa6lGNE/DXeHS0+o6whRW5udaPMU4IBOeOeOGxTjd9v7tgViSbYEYh5fnYRSCUMfDGMylWDynAbruRT7pwTa1hZBlCp5yF9tWe3PyGRTjE783iKc8n9OuO2rEjz1hQz2RpmW0/O5Cgm+8gmH5FgCs5fptu18tsimGSY83iMVuxpG/c7Xu2q4IAn1RV0/hfA9b8eUIYe7XvfA7n2HauIFgNI6n3EW4J0I8msjNB1GMO3p8QQoqXFk59sQ1/BveB02vlLQFdU/LlI6ZqUU2xXDRvS7XTvn7KSlZ29lHrdtKnsXYn8JpK7oUgPCK50HTp3DZgl6SGHBV6HeYKs6vGC7+tgCecmX49wn7zMVg1CNZtj7d8It0vFbl8iuGS8bwD2SLP0okqTEvncIZ634aa8FZGMzFACR6dqTdZXTPmh7MopwOxXDQUhq97SEKlOHfN+zTFlJ8yc8AcB8UQ6ARFQKby6ouPsWw6fHufru9prMPl8VITb6VWNdTSK0PW8mVAGjREJE1b/Q7HQUH6WFFY4Gj/3gKxd4IdIbRUjI3Hr8QoloI8T0hxNNCiI+FEEuEEH8WQpwjhBjz/zScB50KgNnmxi5Cei6/SqtTDJNkPEmwM7zTxeePJmkJxplT7EAgiXT+A5PzEMyOgwAIfPAIWjSI56SvAeA067n8WnpMY6/K7FEMg4xz6qnIThLKkMZbCHEfcC8QB/4X+DxwA/AacCbwrhDi+KxINUKYCqrAYESGLThEZzqXX03iUgyPTCvlgbfbDZ19GIA5RXYSgXfRYtuwF38BAKlp9L51D7apC8lbeAEA5lAMA0kiGjgL7Er3FMOip03Xk4Ky7Hj8e0rn/K2UcvUg21cDTwghLEBtVqQaIYTRhKmgEi2o4ZAt+OP1eMpdNK/N2phfxQSi3+tKX3xJTbK+q486jw2H2Uhvx0MYTCX9KZx9q18j0bGFwvNuwVykXxqiz4HT4CcUz1OV44phk2ktM+qhniGM/sDX41LKxpEXaWQxF9WSCkRxyGZCCb1tQ6i7j2Q8mWvRFGOcjJEuSN9uN/mjRFOSecUOktHNegpnyecQBrP+/jfvwuSpJO/QczA4PBhsLrSQEQfthBIpCsqU4VcMD783iNVhxp5vzcrxh/T4hRCrADnYS4CUUu4+kWIMYiqqJda6Sh/IIiEvXbjl94UorvHkVjjFmCYTlsl4XWs6writRqpcFsLNO6dwxpobiKx/l6ILf4Qw6v8ITMW1aIEIdq2Z9thBeMpdbFrenJsPoxhXDJZGPJLsKdRzblbOOMqYi2vRQgHsqXYwgqVsR290ZfgVe8LvDeJw27DYzXRFEnjDCRZVuZBamFj3U1g9Z/ancPa+dQ/CbCf/mKv69zcXTSHWuhwHbfQlNSorXESCMWJ9cawOS64+lmIc0JPFHH7Yc6hna+YHvcHg/PRPJL1tXGAungKALdwGgNGjl92rRTbF3hiYw9/Q2YdRwKwix24pnKlgJ8GPHsd19GUYnQX9+5uLa0j5e/r78jsq3P3HVSj2xGD1IyPJXlMyhRCXAx8BlwGXAx8KIS7NmkQjjCm9yGYLdgEgVVqdYpj40zn8iZTGhq4I0zw2bEaIdD6MyXEIZud8AHrffQCZjOE58as77W8qqkUm4thjel9+S4kayKLYO8l4kmBXHwXl2esnNpwmbT8CjpBStgMIIUrQUzofy5pUI0gmu8IUBjNRItKBw616oyv2Tk9bgKmHVrLcGyKuSSryLCQC76HFtuGYciMAMhmn9+37cMw9CUvFzJ32z+ieLdQGVjCkm7ypYUCKPZFJI/aUu/CG4rSG4lTmWSjPG7nw4HAMvyFj9NN0MY4qfo35JQizHRk24xDdBONuvUWuMvyKPRDrixMJxrBOL+JTXxiA95oDWO1vU2Aqxuo5HYDQ8mdIBdpxn/SH3Y5hLs7cbQahCFJ2M0KoMKNiz2Sc0kRFPk9t6EICJgHn1ReNmPEfjgF/SQjxshDiy0KILwPPAy+OyNlHASEEpqJqPZefNjWJSzEsMvqhVXr6t2lS0tYH9mI9hVNKif/NuzGX1eOYc+JuxzAV6YMzDGETNhGkL6WRV+RUuqfYIxnHIJhn7U+rTEloDcVH7Bx7NfxSyu8DdwEHA4cAd0kpfzBiEowCO3L5WwjGUxQow6/YC5mLrzDP3L/NQIpSsRpb8WUARJs+JrbtMzwnfRVh2P1SMlgcGPNL0EIGHHQSUrqnGAYZ/XCnm0oKwCigcpRDPUgpHxdCvJp5vxCiUErZPWJSZBlzcS2Rxvex004sJcmvzKevN0o8ksBiN+/9AIpJh79Nb67m8dihN87cIjPlPT+g0lPXn8Lpf+MuDA4PrqOGznUwFdWiBTuwyxaC8Rl4yl34mrpG5TMoxid+bxCnx4YpPSVwYUUe1S7riMb4h5PVc70QwgesBJYBy9O/xw2molpkLIIjobfItacrMZXnpRiKHm8QIUBL59vPt75BkfwUW4nelyfRtZ3wihfIP/YqDFbnkMcxF9WS6g3jwEsoniK/LA+/N4iUg9VGKhQ7cvhD8RROs4GFFa4RNfowvBj/94B5Uso6KeU0KeVUKeW0EZUiy/RnV4T1tDpzcbo3usquUAyB3xskvySPcEpiFCC77sfkOHhHCufb94IQuI+/Zo/H0Q1/Dw7ZTlJCfpWbWF+C6AjGaxUTC783SEF5PsF4CpclO7PBh2P4NwF9WTn7KLEju0JvkWtw62P0lMevGIpMAU0wnsJpDCPjWzHnHwuAFg0TeO9h8hacg7mwao/HMRdPAU3DFs3cbe6oHFcoBmMn3TME6PPeTSK8YkTPMZwY/63A+0KID4FYZqOU8qYRlSSLmNLVu9ZgBIFGKh3XVxefYih6vEEq60sIRoPYk2tBQMT3Vyz5i+lb9ilapBf3yV/b63FMaafDHmoDO5iKdhRxVdQXZ/UzKMYf0bCeRuwud9ERT1IVf5q+0H3gteCuvwezc8GInGc4Hv+dwBvAUvT4fuZn3GC052NweCBsxC789KUkeYUOlU+tGBQp5U5elwNf+oUk8cCH+N+8B2vdodimHr7XY5nTKZ22oL5YjEu/21S6pxiMjDPqqMxHInDgBTSQCRLBj0fsPMPx+JNSyu+M2BlzhLmoFi3oxyG9hOLlKpdfMSR9vVES0STuShdtmg0n7YAAYSbVIki0b6LsmjuG1TkxMwzIHIxiJE7S4kAYhNI9xaBk9CIzo9mJDzCAMGN2HTFi5xmOx/+mEOJrQogKIURh5mfEJBglTMU1aL19OGglGE/o+dRqcVcxCD3pVM5M9pcDH9biy3DX30Pw/dcxusvJO2x4zWszw4BkyIBDdBFOarhLVBGXYnD8Xl33hFuf0eygHUf5DSMa5oHhGf4rScf52RHmGVfpnKB7/MleP3baCcU13MrjVwxBRi9M6ewvBz6c5Teg9dqIrHsb9wnX9PfcHw7molq0QByH1CvH3WXK6VAMTk9bOo3Ypgdj8szgqPj6iBp92PPM3QqAdPrmrj97TecUQtiEEB8JIT4TQqwRQvxnenuhEOJVIcTG9O+CvR1rJDAX1UIyiSPuQyJwVruJhuJEQ7G976yYVGQMv0hnfzlFL8JURO9bf0WYbbiPvWpPu++GqaiWVG8IBz5C8YQeZlSGXzEIA9OIrSKExZqdBIA9efz3CiGWCiF+KYQ4UQgxrCrfAcSAk6WUhwALgDOFEEcDtwCvSynrgdfTz7NOJrsi05ffnu51nemEp1Bk6PEGMZoMJMxGBBpOqxUt3EPww0dxHXkpxryifTpe/zAgzUdfEjyV+aqISzEoPTslFbRjtOw5XXh/2dMglrOAE4G3gIuApUKIJ9Lx/r0OWZc6GatqTv9I4ALg/vT2+4EL91f4faG/iCvYCeycVqdQDMTfFsRdpldO2kUPJmsFve8+iExEcZ/01b0fYBcyumdP5/I7q90kokn6eqMjKrdi/OP3pQ1/LIlDtmCwjrLhB5BSRqWUL0kpvy2lXAh8Fz0T6DYhxEd7O7gQwiiEWAG0A69KKT8EyqSUbenjtwGlQ+z7NSHEMiHEso6Ojn37VIOQ6ZRoC6QXT9INkDKLKQpFhh5fkIJMAY30YjBW0LvkPuyzT8BaOXufj5epI7GFdMNvLVVFXIrdkVLib8t4/Emc+DBaqrNyrn3qqy+l3Cyl/LOU8nzg2GG8PyWlXABUA0cKIQ7ah3PdJaVcKKVcWFJSsi9iDorBbMPoLsMYTGCmj4TZpHqjKwbF7w3iqXARiidx0ErfRxtJ+dtwzD1xv463426zBwBjoZ6xoXRPMZCwP0IilsRV7SElBQ58JNoCdL/0JyJNI5tPM2TcXggRRA/NiPTvgcSATUKIH0kpX9/bSaSUfiHEW8CZgE8IUSGlbEsvILfvee+Rw1w8hVSwGQfthFMF5JfkKa9LsRNaSqO3PYSnwkVzQqMq3Ebw7dcA6H72f7FNOwL7tIX7dMzMMCBrr353mRn/qXRPMZCMPtjS64+WjhZ8z94HUiJMFiq//dg+695Q7CnG75JS5g/43f8DlAPXA38can8hRIkQwpN+bAdOBdYBzwBXp992NfD0iHySYWAuqk3n8rcRiiVUEZdiNwKdYbSkhrPKjURg7WiB9CKsTCWIbHh/n4+ZGQZkCKaw0UvcZMBgMijdU+xEfxpxkX5HaNnWDFoKpLbfujcUw2nL/JVBNv+3lPIz4P/2sGsFevHXSuBj9Bj/c8AvgdOEEBuB09LPRwVTUQ2pQC8O6UvnU+eprB7FTmQuPkuZ7nU5TOl2DcKAMJqxz1y8X8fVhwHF9PbMCQ13aZ5K6VTsRE9bWh/y9PVHhyFtm4TxgHRvMIaTonmpECIqpXwIQAjxZ8AKIKW8c6idpJQrgUMH2d4FnLJ/4h4Y5uJakBJ73EfMKiitdLPu3S1IKYdVfq+Y+GQMv9HjgLhGHr0AuE+4lryFF+73rba5uJbIpvdw4CUYm6XuNhW74fcGMZqNxE1GTESxCkHUaKLw7O9hn3XsiIV5YHiG/2LgGSGEBpwFdEspbxgxCUYRU39ffi9YIa/aTTySIBKI4UgX6ygmN5kFV+m0QDyKLRgiBhSdfysG29ADV/aGqagWGY1iT/loS0g8ZS62rWobIakVEwG/N4inLI9QIoWTdmRQYi6aQuFZN4/4ufZUuZvpyWMHvgr8AAgA/zUee/XAwOwKPT00czuvPC9FBr83iMVuJmoQ2OhFBDWMruIDMvowMJffR0oK8ms8+H0hNE0VcSl0evpz+BM4aEULRPtniYw0e4rxZ3ryLAfeBDzAOYzTXj0ApoJKMJiwBfwAGD12QP/CFQpgl3bMXrTeRH8e/oFg3qVy3FGVTyqRItw9rmccKUYQf3rkYqYVeLLHj7mkLivnGjLUI6WcmpUz5hBhMGIqrMLgDyBIIZ1qEpdiZ/zeIAUVLkKxOC68pHrD2OsP3PBnwoz2YDeUgrlEbwDX4wviKj6wuwnF+CeTRuyudtOsCexJHzIaGRGnYzD2FOrZY4GWECJ/XwqyxgrmolpkIIaDLuImIwaj6o2u2EFPWwB3uYtQQuKQXlK9Pfr4xAPE6HBjcHiw+f0AGPL1u02lewpIpxGnJPYqNwD2Pv3O0FyUHcO/p8XdS4QQvwJeQg/vdAA2YAZwEjAFvYXDuMJcXEts+yfY8RJKVKgiLkU/yXiSYFcfrhoPfin03jpSjojhh3Qdib8DIzE0h1pfUuzAn07ltKTvBG0hva51pHRvV/YU6vm3dMvkS4HL0PPyI8Ba4E4p5btZkSjLmIqmoPWFcWg+ulWLXMUAMjUdjopdvK4RuvhMxTVEt+mV4zFDPiaLUemeAtiRTWZw2yGu4QgE0MiB4QeQUvYAd6d/JgSZRTZ7zEvYIPTS/DW+HEulGAvsGMCix9zt6d46I+nxh1e9nO7LX4unTOXyK3QyDoB0mDHEw5h7QyRHIJtsKPapSdtEoD+trs+LRJBfW0CvL6R6oyv6vS6Rr1dO2nr9CJMVo7t8RI6fGQZkT3oJxlOqiEvRj98bxOowE0HgoAMZGJlssqGYdIbfVJxuzxzSc/ntlfkkYknCPZFciqUYA2SMcMpmxkwYoz+CqagWYRiZyyQzDMge9RFNGXCnB7IoFD39qZxxHLSR8oezFuaBSWj4jXnFCIu9P5ffnJ6rqlrkKvzeIE6PjXBKw4GXVCA2ogU0/QWE6Vz+vNoCettDaCltxM6hGJ/oacT56QEsPlKB3twa/vQwlBtHazZuthFCYC6qxdqje/wGt8rlV+j0e12xOA58pPyBEb34+ocBBbv13+UutJQk0BkesXMoxid+bxB3ZT6RlAF7rG1Es8kGYzge/xVAJfCxEOIRIcQZYpx3NDMV1WLoCWAhSMpuBpThV6QvvnIXoYSGI+lFxqIjWjmZGQZk600vGqfb7/aqDrGTmkwacV61B9DXHyF7GT0wDMMvpWyUUv4ImAk8DNwLbBNC/Oe47dlTXEuqN4QDHzEhMZoMKq1OoRv+KjcJzdA/H3ekLz5zUS3WrnZAg3z9blOFGSc3mTRia0U+APZ0L7FstWuAYcb4hRAHA78Ffg08jp7bHwDeyJpkWcRcVIuMx7GnfITicdwqrW7SEw3HiQRjOGv1iKY9HYc3jXCTLHPxFOgNYaOHlFXPpla6N7nZdQCLLZDOJssvy9o599qWWQixHPADfwVukVLG0i99KIQ4JmuSZZH+vilxH11mSYFKq5v0ZP7+1tJ05WTQD4y8x68PAwrikD4ishiL3aR0b5LTf8eXZ0VEE1h7umEEs8kGYzj9+C+TUjYN9oKU8uIRlmdUyFzM9j4vcbsRd42bLR9vz7FUilySMb6GgrTX1dONMb8Ug8UxoufpHwaU9BGMz1RFXIr+v3/CYsQW9SJ7I1iK67N6zuEY/gsHWcvtBZZLKVeMuESjQH+L3FA7FIGzpoDe59aiaRKDYVyvWyv2E793xyB0Y18YU3d3Vgpo+u82Y168FoMeZlTrS5OanrYAzgI7famEnkbcE8A8py6r5xzOvcRC4OtAVfrna8CJwN1CiB9kT7TsYbDlYXAWYgvsSKtLJTWCXSqtbrLS4w0iBCRMAgftyEBfVrIq+ivHw15S0oB7SoHy+Cc5/TMgYkkcKR8yHstqRg8Mz/AXAYdJKb8rpfwu+j+CEuB44MtZlC2rmItqsPboHfAyiyrqApy8+NuC5JfkEYrH0gU0wex4/JlhQMEuAJw1HgIdIVKJ1IifSzE+yBj+cNKAPTqyjQGHYjiGvxaID3ieAKZIKSNAbPBdxj7m4lrMHT4MJBEuvTeLyqeevPRkvK6ExJ4Y2XbMA8kMA7L2dgJgLctDSlQR1yTG7w2SP6UAiejP4R/pbLJdGU6M/2FgqRDi6fTz84B/CCGcQEPWJMsypqJatBVB7LKdpKUaUPnUkxm/L0jFnDKiKSP2SHa9Lj2XX08mMKYXk3vSJfuKyUUmjdherbcCtwV1hyCnHn+6QvdvwHXoKZ29wNellP8lpQxLKb+QVemyiLmoFrQUDs1HRCYwW039C3yKyYWUEn9bENeUdA5/KLuVk+biWoxd7ZjoQzotgAozTlb604jTA1jsvdnJJtuVvfXjl0KIp6SUh6NP4ZowZG6lbHEf3QlNtcidxIT9ERKxJPa0x20LdiPS7RWygamoFhnuwyHbSZr1fy5K9yYnmb+78OijOK2d7Vn39mF4Mf6lQogjsi7JKDOwL39f0oSnQhn+ycquA1hs3R2Yi2vJVkuq/jqSpI9wKonVaVG6N0npadOjDJrdjBU/9IzMjOe9MRzDfxK68d8khFgphFglhFiZbcGyjbmwGoTAGmpHInCptLpJS7/XlW9DkMTS4cOUpSHXMMDpiPoIJ4S625zE+L1BhEEQEykc0ocWDGV1AEuG4SzunpV1KXKAMFsxuStwBNJpddUeejvCpJIaRtOkG1MwqelJD7pOWg3YIx1ovQHM87J38fUPA+rzEnWZ8FS5VRHXJMXvDZJf7CSUSOJMZi+bbFeG051zK1ADnJx+3Dec/cYDpuLa/lx+a2keUpMEVVrdpMPvDWI0GYiSwKH5kPF4Vjsj9g8DSk+By69Td5uTlUwacSgpcIxSDj8MbxDLfwA/BG5NbzIDD2ZTqNHCXFSDuX3ntDp1AU4+MgU0obiGPZbJ6MleHnVmGJAtnctvr3QT7OojGU9m7ZyKsYnfG8RTV0BKGrGF0oY/i05HhuF47hcB5wNhACllK+DKplCjhamoFoO/G4vsRXMaAZXLPxnp8QXxVObTlzJhj2R/CAboumft1C90S7ojqF8VEE4qpJT6uM9MK/BQl55Nll+a9XMPx/DHpZQSkADpwq0Jge7VSexaOwmTXjKvYq2TD39bgPy6Qr1yMuQDdjRTyxbm4lrMvm2ApsZ/TlIyacS28nQacU9nVrPJBjIcw/8vIcSdgEcIcR3wGnB3dsUaHTLZFY64l0hKw+owq4tvkqGlNHrbQzhrPADYAl0Y3eUYLPasntdcVIuIRbDLLlK63Ve6N8nYdQCLtaMlq9lkA9lrVo+U8jdCiNPQJ27NAn4qpXw165KNAhmvzhbx0ukwqklck5BAZxgtJbGU6dFLW5dvVBbX+tszJ70kbLrHp+42JxeZbDKZZ8akhTF0+jDPHp0kyuGkc5I29BPC2A/E5CkHoxlbyEei0ISn1qMM/yQj8/c2pisnze3NmGefkvXzZhaP7TEfAftMHPk2pXuTjMzfO2nWcER8yERiVJwOGF5Wz8VCiI1CiF4hREAIERRCTIimNsJgxFxYjT3dItdVV6i8rklGxuvSHAZssgvR2zMqWRWZMKOtz0s4adKLuJTuTSr83iBGs5E+qeFIZD+bbCDDifH/CjhfSumWUuZLKV1Syr22ERRC1Agh3hRCrBVCrBFCfDu9vVAI8Wr6n8mrQoiCA/0QB4KpqBZrt/6lOyrzCXaGSare6JOGjLFNGFPYUz5AYh6FOKvB7tKHAYU60TDinlqoPP5Jht8bxFOWRzgxsA9/3aiceziG3yelXLsfx04C35VSzgGOBm4UQswFbgFel1LWA6+nn+cMc/EULO3b9Mcl6d7o7SqtbrLg9waxOsyENYkjnumFPjq32+aiGmx+vYAwT7UMmXT0eAMUTCkgLk3Y0h1hTUU1o3Lu4Rj+ZUKIfwohPp8O+1wshNjrkHUpZZuU8pP04yCwFn104wXA/em33Q9cuH+ijwzm4hrMPW0YZByDywyoXP7JRE9bID39yIytb/QqJ/Xz1GLragXAXu4i7I8SjyRG5dyK3OP3BnFNLQT0PvyjkU2WYTiGPx+9TcPp6ENYzgPO3ZeTCCHqgEOBD4EyKWUb6P8cgEGrFYQQXxNCLBNCLOvo6NiX0+0TpqJaBBK71kHKpgEqrW4y4fcGcU8rRMOAPdSOsNgx5peMyrlNRbWY27cAYCzSy2NUnH9ykEkjtlfpA1js3aOTTZZhOOmc1xzICYQQecDjwM1SysBwixOklHcBdwEsXLhQHogMe6K/U2LCR9ymLzcowz958HuDlC2aQgCw+jswF08ZlQIa0HXPFOnFJMOIPEO/PKV1haNyfkXuyKQRZwawWNq3Y56xeNTOP5ysnplCiNeFEKvTzw8WQvx4OAcXQpjRjf5DUson0pt9QoiK9OsVQPv+iT4y9KfVRdoIawZsLqsy/JOEZDxJsKtvR+VkV9voel3FtQjAnmwnadV9G9W2YXLgT2eTCbcFg4xj6mwZVd0bTqjnbvQGbQkAKeVK4Iq97ZQe2/hXYK2U8ncDXnoGuDr9+Grg6V33HU0MzkKE1Ykt5COSslBQka9utycJGSObGcBi9TaN2sIu7Fw5Hk9ficrpmBxk1hFTNg2H5kOMUjZZhuEYfoeU8qNdtg2njeAxwBeBk4UQK9I/ZwO/BE4TQmwETks/zxn9nRKDnUgMuKcVqYtvktD/d84zYpYBjJHAqHpdmWFAtj4v4ZSJvEK70r1JQsa5jBv0MDOMXjYZDK9yt1MIMZ0dTdouBdr2tpOU8l1gqGBp9ksj9wFz8Y5c/rxaD80fbc2xRIrRoN/rsmo4YvrFN5qGX5itGN3l2ELtxEoseKrdyvBPEnraAlidFsKakdKIntk1GoWDGYbj8d8I3AnMFkK0ADcDX8+mUKONqagWa4eey5+ymwl1R2j8eHuOpVJkm4yRjRrEgD78o2f4M+ez9upZa5YSF80NPrasaB1VGRSjj98bxFPtJqpZsId8ejaZq3jUzj+cCVxNUspTgRJgtpTy2PQkrgmDuagWq78ZgG2bewC48+uPqwtwgtPTFsDpsRFOmbH3tYEQo1ZAk8FcVNOfy9/eESLY3cefr3tU6d4Ex+8N4pleBIAtMLrZZLAPIxSllOF0IdaEw1w8BWMqhiEawFSoF1CkEikalymvfyKTmX6UlCZswXZM7goMZtuoymAqqsXSvhkAY6HenjcZV7o30fF7g+SlB7BYRzmbDCbI7NwDJTP82pH0YS7RMzyEwcCMhaPr/SlGF714S7+9tvX6MI1Sg6yBmItrsfZ1gEz1T+IyGIXSvQnMjjRivRW4tX0LplHq0ZNBGX7AXKhf8HnJVkrm5OMqdlIzt5S6BZU5lkyRTXq8QZzVHkAfgjHaXhfoYUaDTGFPdTL79BoQcPjZc5TuTWAyacSGAitCprD0to1aV84MwyngcgghfiKEuDv9vF4IsU8tG8Y6BpsTQ14RtrCPmMnKrMV1dLdOiM7TiiGIhuNEgzHMZbqXbe3YMqpZFRn6B7LEfQiXgbJpRYT8kVGXQzF69GduOcGmdWCQqTEZ6rkPiAGL0s+bgf+XNYlyhLm4FnvAR1KaqZhXRrCrj94OVUU5UekfwOI2YpRRzNGenHj8/cOAIm2EEgaq55TSsjanxeyKLJNJI05YRH9H2NF2OoZj+KdLKX/FjsrdCEPn549bzEW1WHr0P0LhTL1Jl7oAJy4Zw685wJ70IRj9VE7YMQzIFuqgL2WlanYpve0hgl19oy6LYnTwt+nRhKgwYuvz6tlkhdWjKsNwDH9cCGFnRwHXdPQ7gAmFqagWa6eey++o1jvmtaxThn+i0pO++BImA/Zopg//6C/ugq57tl4fGkaK5+jNapXuTVx6vEGcRQ4imhVb2IsxB9lkwzH8PwNeAmqEEA+hD0/5QTaFygXm4lpsQT13OiaSFNd4aFYX34TF7w0iBETQc/iF1Ykxb/QKaAZiLq7F2tUCgGuK3pmzea0vJ7Ioso/fG6RgRjESA3Z/+6gv7MLw2jK/IoRYjj5FSwDfllJ2Zl2yUcZcVIsl0oNBxghEElTNLlUX3wTG7w2SX+kmJq3YAu2jXkAzEHNxLdY1HwCQMBsorMpXHv8Exu8NUnB4NUnA2tWCubRu1GUYTlbPM+hDWN6SUj43EY0+DBjIkmgnGI9SNaeUruZeIsEJF9VSoN9uF9Snc/j9o19AMxBTUS32kN7+KhgNUTW7lJZ12Rs+pMgtPd5g/wAWa+fWnOjecEI9vwWOAxqEEI8KIS4VQoxuQGoUMBdWgdDjvaG4RvVsFWudyPi9QfKm6JWTlo5tObndzmAuqsWUCGPSggSjIapnl9KxtYdoOJ4zmRTZIZNGbCm2AmALecem4ZdSvi2lvAGYhj4R63JyPDwlGwiTBZOnAlvYRyhpplotsk1YpJT4vUFslXrlpL23BfMoV04OpH8YUNxHMB6nKq17reuV7k00MtlkwmXAmurGmIqNTcMPkM7quQS9K+cR7BiWPqEwFdciE5KIZidoNZFf4qRZpXROOML+CIlYEnOhGSGTWPs6RrUX+q5khgGZ4n10RK1Yp+rNu5TuTTwyhj9lN+SsIywMY3FXCPFP4Cj0zJ7b0WP9WrYFywWB6qNprToTEDzX2E3l8dNo+WyvowcU44yetswAFgP2ZDtCajmN8QshCE0/BX/eXKRm5I2uCAWHVqm7zQlIJo04bjKR79ezyQx5RaMux3Ard6dLKb8upXxjohp9gO7Sg5FC/0pSUpJ/cCW+zV3Eo4kcS6YYSTJeV9JqxB7V2zGbR7kd824y1SxK654gJSWli+uUxz8B8XuDCKMggg1byJezbLIhDb8Q4uT0QwdwgRDi4oE/oyPe6FJm9GHQEiAlghSVTj9aSuJt7Mq1aIoRJGP4Y0YztrAXk6cSYbLkVKaSVAuGlO5gGEhSJHrxbuoiGR/OlFPFeMHvDVJQV4iGCZvfm7O1pT15/Cekf583yM+EatKWodjSyOEvXI8j5cXFdqaXrgFUrHWi4fcGMdrMRLFjD/hyGubJUGxs4rAXr0fIBBV8QK1nE1pSo005HROKHm8Q9ww9jdjauT1n2WRDGn4p5X+kH/6XlPKagT/Az0dHvNHFVn0invaVVPS+RYgqimYcic1lVbHWCUaPN0BhfREgsHa35nRhN4O14ggKfJ9RlFhNH2XUHLwAUFllEw2/N4iz1gOALZCbVuAwvBj/44Nse2ykBRkL2MpPBpMZd3cDGhaChnlUzSpRF98Ew+8Nkj9NX1Czdm/PSTvmXbFVnwiAp6+BXmZQXH80VqdFNQqcQGTSiK1lehmULdiWM6djyKweIcRsYB7g3iWmnw9MuAIu0KdumQtrKehqhnrwhfXWDR88thItpWEwqrk1EwG/N0jN4jqigD2U26rdDJlhQAW9m9A8RrqjSSpnlah+UROITBqx0WPClApiToRy5nTsyZLNQo/le9g5vn8YcF3WJcsR5uJaLF0d2OjEGwpSPbeMRDRJ+5aeXIumGAG0lEZvewhLiRnIVE7mrmo3g8HmxJhXREFPMwC+cJzq2aW0ru9AS03YRLpJxY40YqOewy8E5lFux5xhSI9fSvk08LQQYpGU8oNRlCm3mG2kOjooiDbQHjqKg9KtG5obfJRPH/18W8XIEugIo6UkBrcJa6ITg5YYEx4/oOdzb9qMfb4XbwiqZpcSjyTo3OandGphrsVTHCD9acQWM46e3GaTDSd28akQ4kYhxJ+FEPdmfrIuWQ6INC2jb/WryGScvFWfEEhaya/xYLaaVJx/grBjAIsJe7QNYcvD4My9UY00LSPhayTV00n+tpX4esP9bUNUl9iJQUb3okYbtmBuevRkGI7hfwAoB84A3gaqgWA2hcoVkQ3vg6bfVuf7VgPQGUtSUV+sYq0ThB1j7yzY0w2yctWOeSCRDe+DlAC421cRxomr1oPRbFS6N0Hwe4NY3DaSWLH1tOU0qWA4hn+GlPInQFhKeT9wDjA/u2LlBvvMxf23Xu7utSBT+EJRvU3u2nZk+sJUjF/8viAIiBns2AK5K6DZFfvMxWDS1x3cnbrT0RXXqJhRpDJ7Jgg93gAF6bGutp7tue0PNYz3ZPoV+IUQBwFuoC5rEuUQ+7SFVH77MYyuEmwlZbjFFrxBP1VzSokEY/S0BnItouIA6WkL4Ch3oWHC2t0yJhZ2Ia17Nz4EwkBJiUCQxBfqo2pOKS3rlNMxEfB7g+TV6a3AbTnOJhuO4b9LCFEA/AR4BmgAfpVVqXKIfdpC8g4/n2R7G0WygfaIgapZ6ViruuUe9/i9QTz1aa8rkNt2zLvimHUctqmHI3pCeNiEN9hL9exSwv6ofqeiGNf4vUEclQ5gHBh+KeU9UsqedF/+aVLKUinlX0ZDuFxhr1+MTETxhLeTkGYcUzwIg1C33BMAfQCLBxg7OfwDsdcvJt7SRKG2lo6oicrMQCCle+OaTBqxudiCQYthiXTndurbUC8IIb6zpx2llL8beXHGBvb6owEo6GiBPOhKapROLVQe/wSgpy1IQbmVKLrXNRbaNQzEXr+Inpf/SEFwO5vcJux1boTQ7zYPOmlGrsVT7CeZNGKRb8Ye82KwOTE4C3Imz5768btGTYoxhjGvCEvlbNjShHlqCF9QUj2nlMaPtudaNMUBkIwnCXX3YSo0Y04GMCWjmItyU0AzFLZpR4DBhNvXAm7oSUlKphQoj3+csyON2Iw9sBlzcV1Os8n2VMD1n6MpyFjDXr+IwPv/oECuwxeeS/XsUpY/t5ZgVx+uIkeuxVPsB35fCIBwyogz5MVUUIEwmnMs1c4YbE5sUxaQ2rIRy8wA3mCKqjmlbFnRmmvRFAfAhqXbAIgKOwUBb877Q+01xi+EmCmEeF0IsTr9/GAhxI+zL1puycT5C/vW0xO3Ua5m8I57Vr/ZCEDUYMfQ2UnSVpljiQbHVr+IRPNWCrW1+MIxqmaX0tMWJOyP5Fo0xX6wZUUrL9/5AcJiJGW2Y+1pwVyU2xDjcLJ67gZuJZ3WKaVcCVyRTaHGArZ0nN/T3YLEgKVWj8epKsrxy2evbgTAXOTEGW4mGPPkVqAhsM9cDFqKgvAG/Ak75XPKADUXYrzSuGw7WlLDVOQEwB5szXlSwXAMv0NK+dEu2/Y6Fijd2qE9c6eQ3lYohHhVCLEx/Tt3qxt7weQqwVIxk/ytDQD4ZZLCqnzl8Y9TpJT0tAYwuqwYbGZckW24p8/OtViDYp92BBiMeDq3AwJrOvdb6d74pCqdmWUu1g2/LZj7bLLhGP5OIcR0QAIIIS4FhjOB/G/AmbtsuwV4XUpZD7yefj5msdcvhqb1OGULvkC3XsG7riPXYin2g44tPfS2hzjqS4cDekZPyfyDcizV4BhseVhrD8HVpPtMfhnHU+5SC7zjlL7eKACzTtONvS2U23YNMDzDfyNwJzBbCNEC3Ax8fW87SSmXAN27bL4AuD/9+H7gwuEKmgts9YuQ8SiFsQZ8EUH17FI6tvYQDcdzLZpiH2lY0gTAzFP0oer2MeB17Ql7/SLkto245DZ8AT/Vs0tVOvE4peHtJvIK7Uw/oQAhk1gjXZgKq3Iq03AKuJqklKcCJcBs4ETg2P08X5mUsi193DagdKg3CiG+JoRYJoRY1tGRGy/bXr8YAE/PZiKaneK5eqy1db26AMcba5Y0UTGjGM2ue1+5rpzcG/aZiyGVxBNpoD1qpHJOKR1buon1Jfa+s2LMkEpqrHtvC3OOnUYonsAW7cDiKcl5NtmQhl8IkS+EuFUIcZsQ4jSgD7gaaAQuz7ZgUsq7pJQLpZQLS0pKsn26QTHll2Auqye/eZ3+vMYDqEW28UYkGKPpkxbmnjCNYLQPYyqC1ZDEmMMCmr1hn3akHufv3kxUs1E8twwpoW2DCjWOJ7Z81kpfIMrcE6YRSghs4TbMxVNzLdYePf4H0KdwrUKfuPUKcBlwoZTygv08n08IUQGQ/j3mLah95iLs6z/CIOMERJC8QodaZBtnrH9/C1pSY+7xUwnGU9giY9vbBzDYXVhr5pO/TU8uMNfkA6pf1HijYUkTBpOBWYumEEpasQXaMJeMbcM/TUr5ZSnlncDngYXAuVLKFQdwvmfQ7xpI/376AI41KtjrFyOiYTzJDbSHI1SnuyVOBiJNy+h+6U9EmpblWpQDYs2SJhxuG1MOriSUNGEPeseE17U37PWLsW34GKOMEjCEcbhtkyKdeKLoHeiGf/rh1VjyLERkHrbe3Kdywp4Nf38wUUqZAjZLKYfdIlAI8Q/gA2CWEKJZCPEV4JfAaUKIjcBp6edjmv44f+96OmNOKueU4m3sIhnfa0bruEVKSWjFC7T+4WK6n/1fWv946bi9CLWUxrp3NzP7mDqMJgOhpAPrGLn49oa9fhGGZAxPfD3t6UKuiZ7ZE2n8kNY/XEL3s78c13oH0N3Si7exi7nHTyMcTyExpPtD5b4V+J569RwihMg0oBeAPf1cAFJKmb+nA0spPz/ES6fsu5i5w+QuxVw2A7d3A6liEwXzykklNdoau6hJL/aOZ7RYmHjremIta4i3NBBrWUu8pQEtsmP2gEzqE6Ls0xbmUNL9Y9tqL6HuCPNOmEYipRHHQSSvEn/xFIpzLdxesM04CoQBT+96NlvnUTWvjHce+IRUIoXRbMy1eAdMKtRNrKWBeMva9O8GYs1rQNOdKpmIEVn/zrjUO9iRSTb3+Kls9es+c8pkzXkqJ+y5V8/416wRwl6/iLyGd+AgMFTaAL2YZrwY/kjTMiLr38NSNh0MRuKtDcSa9Qst0bmlf+SfsOVhrZpL3sKLMFid+N+8C1JJkNqYGViyrzQs2YwwCGYdU0dTuuVBd+VRvGEQ5IXilOflZtj1cDDa87HWHER+63q0UhMFB5WTSqTwNnX1z4gY6+i69w6molqEMOzkYKT8O8qBjHlFWKrnkXfYOYQ+fUHXOyShz17CfcI1GB2enH2G/WXN202UTClAK87j/Q2dgKDxiJuZ43SR62Yhe/L4FWns9YuxvvsAVq2bkMmI1WkZN7fckaZltPz+IkgNSAMUAnPJVCzV83AddRmWqrlYq+ZgKqxBGHZE/5wLzia49J8Elz9Fxz9+iNFZgGPOCTn4FPtPw5JNTF1QidNtZ3OjV99oMKIhaR3jhh903cv76BlYAIYKXdaWte3jwvBHNn2s6542ICxqNGOpmIlj5jFpvZuLpWouJveOzxNpWkZkw/vIeJieV/9M86/OpuIbf8dSNn7aUsf6EjR+vJ3Flx9CayiOlt6uGYx4E1Zl+McD9vpFCMATWkOH+RAqZ5WMm+yKwJK/7TD6QpB/3JcovuinGKzOve5rn7YQ+7SFFJz+Tdr+cjWtt19J8WU/x3PCtdkVeoTw+4K0rOvg3JuPA6A30gsyD6SGwZCi1LQVmJdbIfeCfeZibK//BVuqnbDJgsVuonldO0fmWrBh0P38r3cYfWHAfcK1FF/80/651kOR0TsAx9yTabvrGpp/dQ7lX71r3DgeGz/cSjKeYt4J07AZtwIFut5pyTGhd8Op3J30mDzlmEun4W5fTyCVT+X8clrXd6CltL3vnEMijUsJLn8GhABhRJisuI68bFhGfyDm4ilUf+85HHNPpvOf/077P36ITI39QqKGJZsB9BzqeIruRB7Vbc9Qv/ovnMAtFCbfz7GEe8c2/SgQAk9wLe0xI5UzS8bF3WbvkvuJrFsCBiMYjAiThbyFF+7V6O+KfcZR1PzwJUwFFbTefiX+t+/NksQjS8OSzVidFqYdXk1faBMgqNzyPEd9+s0xoXfK8A8Te/1iXJs/ASB/XhnxSILObf7cCrUHYs0NtN3xJczFtVR84wEKz/shld9+bL8Xygy2PCq+/jc8p95A4J37ab3tSlLhnhGWemRZu6SJgsp8yqcXsaErCBiYtuZvzGh5gCJDI2bXEbkWca8YHW6s1fPxeNcRSrmoOLiClvUdaNrYHb4e+uRZOv55C46DTqXq5scpPPfAdM9cVLuT49HxyC1j2vGQUtKwpIlZi6dgMhtpis3AKnuY89Z/UcSGMaF3yvAPE3v9IvJbP9Nv18r1i26shnsSnVtpvf3zCKuTym/+A+dBp1J45k37fOF5Q3E+8YbwhvTeRMJgpPjin1L6xT8QaVxK86/PIe5rzMZHOGASsSQbPtzKvOOnAbC+s5si1mDv6MRSdjDu+nswOxfkVshhYqtfRN7m5QDkH1RGLByna7s/t0INQd/6d/H+7UZsUxdS/tW7sM84ep91b1e9g50dj94lf9Mdjz5/Fj7BgdO6voPe9hBzj5tGLKmxvS+f6uSbGFJJnDOuGRN6pwz/MLHXL8KU6MMVayJs0dPpxuItdzLYQev/XYFMxKj85j8wF9UMa79YSsMbjrO2s4/3mwM8sa6dJzd08mFrgGc2du50EeYvuoKqbz9Gqs9P86/OoW/t29n6OPtN48fbiUeSzD1+Gl2RJP6EnSksResLQ9RF0jd+6jAcMxeT71uFkCkMZXp4cSy2DYlu+4y2v1yNpWQqFTc8gMGy90l1Ukr6EilagjFWd4R5uam7X++e3kXvdnM8fnX2mHQ81rydTuM8biqbesJoGKlp1zutGgxjY1FeLe4OE1NBJeaSOjxda/HZT6NiRtGY8/i1SJC2275A0u+l8tuPYq3U+817Q3FaQ3Eq8yy4rUZ6osn+n+5InJ5onL7kDh/AQAIzAaAQEKSkxtauLZTnzex/Tyb22nbHl8bkom/D201Y7CZmHFnDh94OBAmq2n1Egb61bxHZ+MEBhR9GE9uMozCmYrgimwi7CzGYDLSsa+fQM2flWrR+4u1NtN12JUanh8pvPdKffpnRvYo8My6Lie5okp5IQte/SIyeaJKYNlD34oAZEGhS8nFLC+fO3Hk+bf6iKzCXTsN759hc9G1Y0kTt/HJcxU7WNzSRl+zA9MozAHQ9+XOsUxbkXO+U4d8H7PWLcW37jO1V51G+sJq1zzYgpczp0OQMWiJK251fJtayhoqv39+vWN5QnKc3drIjJLxDVhNRXGyhhG3ks5V8tuOxJMi3u+lMuHmj7xpSWAADG/0as8uTuK07VCYTe/XedwOd//x34q3rcR1xEZHGj7DPXJwz5ZZSsmZJE/VH1WK0GGnsiVLBcuSaTZk3IFOJcVOUZnR4sFTNw93ZgNdxFuXTisZU25Ck30vr/30OiaTqW//E5KkA0rq3oZMdKRA7dM9KABdbqMronmjFY5WEtQLein+LFCZA0By28dImLyfXlWE17fgHYZ9+JNU/fJG2O3Zkm1lrDtb/pjnUvVB3H9tWtXHGNxYTjCXxRu3M6nwdkUoBILXkmNA7Zfj3AXv9YtzP/B4A1xw34Qei+H1BCsr3WMScdaSWwve3G4lseI/Sq/8P50Gn9r+2xrsdTdrRLzqNcj6kXjyHxyJx2Qsw2adhsk3HaFuE0ToFYbAC4Ayv4IQNP6ZdzkWQYF3qSh5taOOkumKmF9j7j2+w5VFx/X10Pf0L/K/eTuDdBwAQJkvOPGrvpi56WgOc+tWjaA7EiGg2DjM1Em/ZCMIAQiCMZr318TjBPnMx+dtWsL32QsqPrGHDC+vGhNOR6vPTetsVpELdVN38+E659stbmtHI6IpGJe8zy/AcHpsRp70Mo206RtsMTLYzMFiqEMJIIrwCueHHtMt5lLCKbjGHlYFreLShmTNmVFDi2NHOWHc8nu13PDAYQYIwmXOme2vf3YyUeibZ+i4/AFNNHehmf+zonTL8+4CtfhFO/02YUmEo0S+4lrXtOTX8Uko6HrmV8KfPU3zJf5J/1GX9r23rjdIYsKAPT9MwkmBBQQ91dXcixJ7/9GbnAqbN/D41wY8x2mcztf0+3gmeziubTcwPhlhUXYzRoH8HwmCk+KKfkOjcSvjT53S5cuhRN7y9o1T+445WzKQoa2knGA1SfMUv0foCOfUK9wf7jEW4l/8nAK45+YQe7CPQEcZdmpczmbR4hLY7ribu20TljQ9hm7Kg/7XPvL1sC9uBFCAxkuTggjB1dfft8Z/VQL0zu37IDKOLksbf8178ap5YZ+CYmnzmFbv6j5FxPFr+cAnRxqUAyFTuWoysebuJ/BInVbNLeGd1E8U0YW8LEDJaKDjtBhwHnTom9E4Z/n3AXFiFuagad+9a+orrEELP7DnopNxVFHY//2sC7/4dz+nfwnPK9f3bt/ZGeGlTJy62Mp/78DOdUrGG6tLv79XoZzA7F/RnIFTmH8u57Q/wQct6VnVegC+0ndNnVOGy7Ojs4Tnl64RXvpwuGBM582waljRRNbsEZ7GTLS091LCE6IcfYK09BPdxV+fcS94f7PVH4whsw5QM9jsdzevac2b4ZSqJ96/XE236iLJr/4Jj9vH9r61o6+GDtihVvEM9T9PJvAG6t/fvfqDeAcyY+wsKtv2Bd7rn8M72I2kLhDixrhyzUQ/9CIORogt/TMvvLgAtBQZjTnQvlUix/v0tLDhjFl3RFL1JJ4fzGeHlr5G/6AqKzh87k2ZVVs8+Yq9fTH7LKvxaASXTi3Ka2eN/66/0vPA7XIs+T9EF/96/fYs/wkubushnC2cUL2fWzG9yeGUJ02Z+f79TyYQQ5JV9iRNnHcNi0+10RxM8uqaZren+N6BXXFbd/ASWqnn6BZiDXOtwb4TNK1qZe/x0mnqCJDEzJbiJRPtmPKdcPy6NPoDRWYC1chaenjX02Xe0bsgFUkraH/oufateoeTyX+A6fMd4jk/auvigLUo1Szh1ip0Zs75/4LpnsFFcdwtn1Vk4SDzIpl6NRxu20RXZoV/2aQupvOlfGPPLECYr5oLRb4rQ9GkL0VCcucdNZV17KwbiVDQ1IpMxPCdfN+ry7All+PcRe/1i8ts+Q2Kk/OjanC2yBZc9ReejP8Z58JmUXvnrfoO22R/m5aYu3GzijNLVFNbchCXvUBzl141I/rDZOZ/5c3/AWfkPYZPbeaHJz9Lt7WjpRm/26UdQ/d2nMZdOxXvfDaRCXQd8zn1h/XtbkJpk7vFTWd/hxYEP5/KPMHoqyDvsvFGVZaSx1x+Dq2UVAa2A4hm5W+Dteur/EVz6TwrO/i7uE67p376spZMP2+LU8Ban1rlxFJ2H2blgxHTPVngWi+d+jpOsfyaaCPP4Wi/rOnd0infMPIaqf3sCZArvfTcgU6ObstvwdhNGs5EZR9fS6E9RLj8m+dG7OOadgqW8flRl2RvK8O8j9vpFuDtWAeCc5aCnLUh4gNebbSJNy/D9/WZ8f7sR2/SjKLv2DoRRD9009YR4pakHDxs4s7yRwupvZMXDNZjyqZr+E86ramUqL/FpR4pn1m0hnNCXsAy2PMq/ciepUBe+v9+MlKNXZbomPdi6eFYJrZE8amLvEW/8DM+J1+Z8zumBYp+5GLd3JRID5UdVj/pQlkjTMlpvuxL/q7eTf9zVFJ7zvf7XPmrp4GNfglre4JSpZdgLz8yKDEZrLTPn/Ixzi16nkDW8uS3EG00tJNNpa5ay6ZRc8b9EG5fS/eLvsiLDUDQsaWLGEdW0J5NENTs1XZ+RCnXhOeXroyrHcFCGfx8xF9XgtNuxR1qhVDe4o+V5RZqW0fqHSwgufQSkRsGZN2Ow6FkTjd1BXtnsx8N6zqhoxVP5lQM+35YVrbx2z4dsWdG622tCCFxll3Pq7EUcZb6X9gj8a802Pt2+jg83vkGH3UDxRT+lb/Wr9L5x1wHLMhwyg61nHzuVxu5OJAYqN3yEsNjJP/aLoyJDNrHPOAp3h14IlDfLSXdLgL7e6KicO9Plta/hDRAGXEdcjBACKSUfNvtY7ksyhdc4ZfoU7AUnH9C59qR3AMJgoXTKTZwz1cIc8Rjr/QYea9jM+ra1fLjxDQLT5+I6+nJ6Xvw9fevfPSBZhkvHth7at/Qw9/hprPc1YyaA56PXsFTNxT7r2FGRYV9Qi7v7gX3msbh9q+itOxrQDf/Mo7M/0Sm8+lVkMpZ+ZiC2bSXOuSeysSvA61uDFLKWM6p6cJddtc/HTiU1en1Bulp66W7pZdOyZpY/vxYpJSaLiRvuvoy6BbvHTU2O2SyYezMlm2/njcBpLO0oA1x8Fkxw1kFH49xwFp1P/T9sM47aKesjG2xd2Upfb5R5x+updB7Nh2nZG+Qfe9W47Oe+K8a8IvKKynCEt0GpfvfSsr6d+iOzPysh+NFjO3V5jTR+iG36kSzd7mVFJ0wVr3DS9NlY84/a52MnEyl6WgN0tfSyYelW3v77crSUxGQ1cuM9lw+qdwD2gpM5zjmH0sZ7eC96BW+0udF1L8lZp34Z8+ZP8P3tRmr+/TVMrpID+PR7J9MQcOZxU3mmM0BN3+torZvxfPEPY3JdSRn+/cBev4j8Tz7CW3cGhfXFo1I+L7XUjtYIwoAw6fnAGzr9vLEtTBGrObOmD1fJ5bvtu2VFKxs/3k7VrBKcBXa6mnXj3tXc2/+4xxtESw4otRH981lIRJOsWdI05AVoMDqpnv59pjc8xqpYGWAghYlWfzsLr/od2//nVLx/vZ6aW1/FaM9e6mvDks0YTAZKD6/mg20B5nmfBC2J+6SxtbB2INjrjyHfu5qeaXqr6Za12Tf8qXAP4U+fB4Sue0YztvrFvL+tjZVdgmniZU6acTAW12G77ZvRvYoZRTg8ad0bqH8tvfT6ggwWDUzGUrxw23tcd/tFmK2DmyqjpYI5c26hffVjrE2cABhJAa1hPwu+chfNvzqL9vtvouKGh3aaNTHSNLy9ibJphfjtKVJYqFz3Dsb8ElwLL8raOQ8EZfj3A/vMxbhf+isAZUdX0vxuc9bP2f3Cb4lvXUHB6TchbHnYZy5mq2sab20LU8xKzqjVcBXvrmRNn7Zw+7X/2smoZ8grtFNY5aZ2fjmHnjmLomo3hVVuiqo9+L0B7vz6EyTiSZCw5KFPcJc4WXz5IRiMu19AQgjqyubSsC3eX+2b5yzD6Cyg7Jq/0PL7C+l4+PuUXfuXrHlADW83Me2wKraEWkHaKVn6FM75p2MpnZaV8+UC+8zFuN97Fe/0s/HUZ79tiNQ0fPffRKrPT8mVvyEV6sJWv5hPjOWs6hJMFy9yYv0RWPLm77bvhg+3cef1j6Gldrfq7tI8CqvdzDiyhqIqd7/u9fVGeeAHz5NMrxdt/HAb/3vh37jolpOYd8L0QWUUwsyMioNZvy2OhgUQlOSXYS2bS/Gl/0nHI7fgf/0OCk67cUS/mwzRcJxNy5o5/ouHsb7diyMVx/HJi7jP/QHCbM3KOQ8UZfj3A3NRLR6tB6ElcNQ7WPdgN7G+BFZHdhYP+xreoufF3+M6+nIKL7gVIQQN7V28vT1CKSs4c4oFZ9Fpu+0X6unjkZ+8tMPoC1h4zhxOuuYIiqrdWB1D90YvqnZzwz2X0bhsO4WV+Sx9YjWP/+INlj6xikt+dApTF1Tttk918TzOZg2bfJ/SGJvN8nYzdcVJnNOPoPDcH9L9zC+wzzoO97H7HoraG92tAdoaOznveyfQ6JeU9i3D3L0N95d+P+LnyiX2GUfjflJftCw/qpKWpYPHwUcK/2u307f6VYov+3+4j/0CUkre2bKdNd1G6g0vcEL9MZids3fbz9fUxQM/eH6H0RdwxPnzOPUrR1JQmT+kBw/0692MhTXEIgme+J83uOebTzH3hGlc9MOTKK7x7LZPdfE8zparWdu2nqbkEXzWaWBKiST/uKvpW/8OXU//D/YZR2ObevhIfTX9bPhgK6mkxrTjpvFuzMCMtn9iMFlxH/elET/XSKEWd/eTvBmLcHWvR5bYkRLaNnRk5TyJ7ha8992ApWIWJVf8El84wYuNPt5ujlPGp5w5NW9Qo795RQu/vewBupoDGE0GhFFgtpo45nMLqJxZskejn6FuQSWnfvUoDjt7Dt+4+1Ku/s25hHsi/OmLj/Dwj14i2BnebZ/q4nmcMO8iTiteTlQz8cL6DSQ1ScHp38Q++wQ6H/0xsda1I/LdDCQz2LpsUSVhzUV5w2tYa+Zjr899efxIYnQVU2hOYkjFcMx00r65m3g0O/USkQ3v0/XML8k77HzcJ36FtlCMJ9Y2s6bHzEzD85ww84RBjf7y59fyuyseIpVIYjQbMaR1b/GlB1M6tXCPRh926F3dgkpmLZrC9x//Eud953gaP9rO/174N168/b1BP3NNyUGcNv9CFtmfxxvz8Oam9QCUfuF3mAoq8N77dVJ9vSPz5QygYUkTdpeVWGkciYGSjx7HddRlGF3FI36ukUIZ/v3EMXMxbt8qIrZCMIis3HLLVALfX69HJmOUX3cP7XETT2/oZEtAA1IcVlmBs2DnroRSSt66fxm3XfMvjGYjNz98Jd+873Oc/c1jhlygHQ5CCBacMYtbnrmGU75yJJ+8sJZfnH8fSx76hNQgYaSami9wbN7rdCbcvLqhAYSg7Mu3YbC78P31erR4337JMRQNbzdRXOOhjU6MWh/Fa57DffLXxuTC2oHirD8GV9c6ZIkdLSVp29g54udI9rbjvffrmIvrKL3qd3jDcZ7e0EV71IQgyeya0zA7dq5YT8SS/Ou/XuXBW16gak4pP3jyy3zz3ss56wB1z2Q2cvI1R3Drs9cw/5QZvPKXpfzvhX9j1RuNu6UKC2Hk4JlXc5DpRTYG3SxvbsLocFN+7Z0ke9pof+g7I5perGn60JXZx9SxsTeKJ9yAs3vTmCvY2hVl+PcTe/1i3O0rSRmsuGeVZCWls/PJ/0d08zJKr/odlrIZfNa6eUCnQ0lH3LXT+yOBKPfd/AxP/+Zt5p0wje/+8ypq5pbt5EEdKFaHhXNvPo4fPHE1U+aX8+Qv3+R3n3uQTct3XucQQjBnxpc51PoyW/oK+WDLekz5JZR9+Xbi3o10/uvHByxLhlhfgo0fbWPOSdPZEsqjrG0JFqd7p4rSiYRj5qK001EERjHiyQVSS+G77xtokSDl192DwZbH0q2bkenumhJoDcd22qdzu58/fvEffPDoSk6+9ghu/OvleMpcI6p7njIXX/rVudzw18sw28zc++2nufuGJ+nYtvMkOGF0sHj2edQaPuDjDgcbO1qxTT2Mogv0nlaBd+4/YFkyNDf4CHb1UXfSNHqShZSvewnH3JOwVIydltmDoQz/fmIqqqUw0gJAyZGlI37xhT59nt437sR9wrW4Dr+AbX4/m0MO9MsuiZEkJXzW//7ta3z85vIHWbOkiQu+fyLX/P587Pm2EZVpIKVTC7n+L5dwze/Ppy8Q5bYv/5OH/v1F1ry9qT8HWxjMHDHzEqYZ3+GzHjcNbVtwzD6egtNvIvD+wwQ/fnJEZNn40TaS8RRFx5SSwE75yqdwn3DtPs93HS/Y0kWEmsFC/qyRn8Hb/dyviWx4j5LP/xJr1RxWtGzHG3MjSDGY7q18fSO//dyD9LQE+Mr/Xch5/3Y8RlP2TEv9kbV8/9EvcsH3TqDp0xb+98L7eeFP77Lxw239ume0lHBK/SEUsY43tyfwBnrwnPINHPNOpvOx/yDWvGZEZGlY0oQQkKpLIGSSsrXP4Tn5+r3vmGPU4u5+IoSgsGoa5mgPxhl5bH54LamEPpnrQIm3b6b9gZux1h1K8cX/QXsoxMtNfly0cgj30MMMSsUaqoq+j5SS9x9dyZO/fBNXoZ1v3fe5EfGuhoMQgoNPrWfW4jpev+dDXr/3Y5Y92wACzNYduf8nzz6WUMMq3mmbQ76tnapzv0+k8QPaH/4e1ikLsJROPSA5GpY0YXWYCeRHsQY7KOpYifvY+0boU449TK4SilN6eKfkiFJalo+c4Q+vfo2el/5A/uIryT/6c6xrb+MDn4kq8QEz5eN0pBuuVRV9n2QixXO/f4e3H1hO7UHlXP2bcymsco+YLHvCaDZy4tULOfTs2Tz72yW8eveHvHr3hwgDA+pOpnPG1B6e2tzNi5tiXDLHTtmX/sS2X5yipxf/8GUMNucBybHm7SamHFLJtpiNYt+HOAtKsY+hoTBDoTz+A8A58wTc7SuRxXmkEime/s3bQ1YbDhctHsF7z1fBYKT8K3fhT0qe2+jFQg9nV4eYNeum/qZXmpjLg7e8wGM/f42ZR9Xy3Ue/OGpGfyBWh5mzbzqW479wqL5BZmbebgPAbK3kjBlVOGnj5c0B/LEIZdf8GWE062sYidgejr5npJSsXdJE/YnTaYmVUL7xJfKPugxjXuFIfLQxS0HtXCyRTsw1eTSvbd8t1LY/JLqb8d3/LSzV8yi+/L/Z0t3BW80pSljF6fWzmT6g4VooMI3bvvxP3n5gOcddeSjf+vsVo2b0B+IuyeOqX57N0Zfo6aRSg2QsSeOy7frrBQs5o7KDpDTw/PqNJO0FlH35dhLtm+j4160HdO7ejhDNDT7qzpxGn8ynvOFZCk4eH40AleE/AGz1i3C3rybmKMXgMPPOw5/y568+ekDGv/PRHxNvXkPZ1f9HzFXGs+u3gIxzRkUzBaWn9ze96myt5ndXPMSnL63n7JuO5au3X0Rewd5nnGaTQ06bqWdsCEDqRS2RgN5SIM81jzNrE0CS5zdsJekqo+xLfyS2fSXe+79F90t/ItK0bJ/P2bqhA78vRMmJRUhhomLDc3hO/trIfrAxiGPmcbh9q0jme9BSGn/+6qNsXtGy38eTyTjee74GqSTlX70bbzTGK1v6cLOZs6YXY8ub2697jZ+4+c1lD+Dd1MXVvzmXi289GdMI3OkeCEddeFB/tpCUejV9JumgvPx0Tir8jN5UAS+vX4lt5jEUnPVvBJf+i65nf7Xfurf2Hb1aV8wUGJMhyrpWkHfkxSP3obKIMvwHgLmkjoKQnkZom1YE6J5uw7tN+3W8wNJ/EnjvIQrOuAnj7JN5du1GYpqZ00rWUF6hd5bcsqKV+77zDL+94kEiwRg33HMZp113FAZD7r2MugWV3HDPZZxz07GcfO0RNK9t5/dfeJj2Ld0AFBcfwymlTYQ1Fy+uW4V13mk4Dz2P8CfP0P3s/9D6x0v3+QLMDF3pKzWT52+ktKJqpylQExX7jEW4O1ZjKCvE4LSgJTWe+tVbpNKFT/tK5xP/RWzLJ5R+8Q8EnaW82NiFjU7OnmLE6dYrcjctb+a2a//JXTc8gafMxXce+QILzhgbi5gZ3Tvrm8ew4IyZrHh5A3ff+ER/L6PpUy7mCOc7tMTKWdK4gsKzvoOlah49L/6O7md/uX+6t6QJT7Ubn6GIsqbXKTzmKgzm7K2rjSTK8B8AQgjKClwgNZz1RWTu8D5+umG3TIO9EWtZS8c/btHn+p79fZ5bt4ZAKo+TCpZSU6NP1Vr//hb+dPUjrHx1I1pS43M/O40ZR9SM9Mc6IDJZHOf92/F84+7L6OuN8vsrH2bd+1v016vPZnH+R/gSFby58RMsVXP0HaVEJmJE9rGpVsOSJmoX19IjKinf8DwFp9wwwp9obGJyl5If0sM79hlFGIyCbau83P3NJ4kE9y10FvrkGXrfugf3yV9Dm3MKz25owUCEs6p7cRfpdRCfvbqB26/5F5s+bsZgFFz4wxMprRtb4bS6BZWcfv3RXP2b87jiv06n8aPt/OGqh2nf3I0QgsNmXspsy3usDZbzWes6HPNO0neUGjIZJ7Lh/WGfKxlPsv79rUy9cDpJ4aC86SXcx385Ox8sCyjDf4C464/F6d/M9NMcnH3TsVz6k1NIRBP84cqHhx131aIhvPdch8HuouSaO3h5wxo6EkUcm7+E6XVXArBh6Vb+9t1nkdqOSshs5G+PJNMXVvOdf3yBgnIXd33jCd5+YDlSSg6afgnzbe/T2FfFhsqDEf1ekiS04nmSPcMLlYW6+9i6so3S08pAatQEV4/JTojZoqqiGKTG3IvL+dbfruBzPzuNjR9t509f+gfdLcMrVIr7NuF78DvYpi7Ecc4PeXZdI0lp4IzybZSUnoyUkg8eW8kDP3x+p/z3rSvbsvWxRoSjLprPDX+9jEggyu+/8DDr3tuCEGaOm306VYbPWNrhoWvGYoQp3VJBakQ3L0OL7V6UOBibljUTjyQwz7FhDfuYUl2NKT+7jeBGEmX4DxD7zGNwt68iZC3llK8cyTGXL+Dmh79AXoGdO657TM9y2QNSStof/A6J9iZKr7mDt3xtNMfKONLxJnOnf4F4JMnj//06d1z3GHaXFZNFr4Q0mY3MWDi2vP3BKKxy8+0HP89BJ07nqV+9xT//4xVSCY3Fs85miukTVmjziH7lbgrPv5WCM79Non0T235xCqEVL+z12JnB1okqNwVty6g45qpxsbA2UrhnHk9eTyP2WiN1Cyo5+pKDuf6Oi/H7Qvz+yof3apy1eB/eu7+KMFoovOYvPL+xkZCWxynFa6isPJve9hB33/gk//rPVymfUYzJquuecZzo3rTDqvm3f1xFYWU+d93wBG/9fTlGo4vTZx2Cm+0siddh/cYDFJ73Q1xHXU7f6lfZ/j+nE922cq/HXrOkCUuRg15nNeWbXqTglOz0AcoWyvAfIOaSqXgCG4mb3CxZ/xbNnWsorvHw7QevZOqhlTz07y/y4m3vDVkt2Pv2vYQ+eYai82/lE6OTTeEy5lvf4rBZn2PrZ15+c9nfefeRFRx/1WHc8sw13PjXA6+EHG2sDgtf/v35nPa1o/nwydX8+auPEvZLTp29iCLRyDvJmXxSu4DI4vOpufU1zEW1eO+6lvaHf7DHCt+GJU0UHlpBxFpK5fbXyT9ibHZCzBb2+kXYAq20y1q2d+h56TOPnsLND34eq8PM7df+i09fWj/k/h2P3Eq8bR3FV9/GK952upNFnOD5mKm1F/LJi+v41UX30/jxdi665SS+88hV3HjP+NO9wsp8bvr7Fcw/eQZP//otHvnpy5iMJZw1vQAzfbwSKWXNtENJnHM9ld9+DC3eR/Ovz6HntTuQ2u4V6aA7aw1vNzH1ojqkMFEXb8SaCVmOE8RoTkfaXxYuXCiXLdv3VffRYum//otPZ1wHSAQpji7u4eCa+WhJjcd+/hofPrmaw86azRU/P2NA5oGk990H6Xzkh1inHo7vkv/Hxz2l1Jvf5/gZZ/DKnZ/yxn0fU1Dh4vM/P3PMxfL3l09fWsc/fvwyeYV2vvKnC4nkbeTV9jr0VCCNw93bOazmUPzP/xr/q7djLqun/No7sNYctNNxNi3bzl+uf5yZN81BmzOT89v/QcW5/z7YKScszZ1reG6zC2k0AxozHds5pOYgihxOwj0R/vrtp9myopWzbzqWU796ZP/dkNRSdD39P/hfvY28RVfw6eFXsSVaxSLXO0wvOpcnfvE6K17ewJSDK7jyv88cc7H8/UHTJC/f8T6v/GUpdQsqufb357M58B7v+XW9EmgsLmpjbuE0Oh7+PuHPXsA++wTKrv4TJnfZTsf65IW1PPDDFzj4t4sxiTgXVcTJO2j3flljASHEcinlwt22K8N/4Lzz1t9Z7ToFhEHPJRMCcypEaWgd5YHNiM+WEVy7nvx8jaISAzIWQAv70atwoXX2xaw59ifUGj9hvnEej/z0bdo2dnL0JfO54PsnYnNOrArU7Wt8/PWmp4gEoxzz3/NoKZoNGNG/D4FDdDMrP8C0SB/hh75DKtxD8YU/xn3iV0nEU6x+axMP3/oiKSmZfdu5lHg/4IIzzhnTTbGywYcb3+CT3nowmPr1DsCW6KI8uI5S/2Zi779P3NtCQRG43RIt0osWCQD6t73+mFvZPudyDrW/S57vUP71s1fp641y5o2LOenLR2S1AjcXrHh5PQ//+CWcHjuLfz6dRsssBuqex+BlTkGK6m1r6X3spwirnbIv/gHn/NOJBGOsfmsTj/zkJYzFTqb+6nxmrLiLU7/8k6z2+j8QhjL8qnJ3BKiwGmhIxdEMJgxakqmf3kMkv5qOmmNoqVkINZfhWrSezvUr8G5by0H1DizhrUQ3fkB77fE0LL6VkuAnmNY5+NNfnsThsfPV2y9i3vETp4/8QGrmlfFvj3yB+25+hg9ub6H232eCAQRJZhkaaOvL41NZz6daCsvRt1Oz8hHkYz9lxX3388qaM4gk9GpLz2FFaA4X+a0bMLquzvGnGn1Ko30YtCQaYNCSzH/jhyRshXTWHMv2qqPYUngMYsrl5PvW0LH2M9wdjcxbUIzs2kR0w3tsXvBVts+5nGkdL9D0XiUfP/k0lTNLuP7OS6iaVZrrj5cVFpwxi6IaD/fe9BRLft9C5fdm6PNlSFGnrcGXKuEDrRphKcK26E7qV96BdseXWNO+kHc2HEdK6iaz5sRCkBqWbsOYNfp7IicevxDiTOCP6P9q75FS/nJP7x/rHn/3S3+i6cPn6Sk/nALvcuoOOY7Cs7+DRpT2rs/Y4u+hOVpIl1YPwkAqGKUi1YS1dSXNsy7BHthG6sG3WbXMxqFnzuKSH52C02PP9cfKOolYknu//TRbvCHsc0qJrG0nuknPVDKXuXCfNAP3cVMx5tkw9XiZuu5hyja/QVfxl3n+XylmfnsGYsY0Tox5mXHygc15HY/sqndT5h1FwenfwmjPJ5HspLlrNdt6+2iJVRBAHw2a6g4yzbaFeMsW2urPo2j7u3Te+xmNWwo55dojOOMbizBZJr4/GOgMc8d1j+KXYjfds00vxn3yDPKPqkWYTdjbNjB93d8o8G6go+omXnigjZm/Oh5nuJlFlQcxdeHYddDGTKhHCGEENgCnAc3Ax8DnpZRDpr+MdcMfaVpG6x8vRaYSCKOZym8/hn3azt+1lBrh0HrWNG5k6acFGOtrMDqt6fz1JB3/t4Tzv7iQQ8/cvb/5RObVuz/khf97N3OnzfyTZrDosoMpqMjHU+7CYImyvmU5a3ssdMlpGFJRyje9TF6il42zP09x2/tcev4Vuf4YOWE4egcgtQTd/pV8sqqZVVsqsM6qwGAx9ete8L73+dwNx1N3yPhYsB0pXvnLB7z45/f7de/gU+s59ooFeMpdeMpcxFMdNLSsYl2whBDlmKM9VG18BuFws3n6hUzZ9BRnX/aNXH+MPTKWDP8i4GdSyjPSz28FkFL+z1D7jHXDD/pFGNnwPvaZiwe9+AYSDcW48643MR83H2E0IFMaNYG1nHfKqaMk7dhhy4pW/nzdo/0N7obKGJFS4utuYFVbC0199WgmPfffkIpxum01U+efO9qijwn2Re9Ab5987wNLyTt9PsKg615dZC1nH690byjd07QEW32fstoboDk1q39NxZCKcbp9zZjWvbEU468Ctg943gwcteubhBBfA74GUFub3WHSI4F92sJhXXgAtjwrUwqNtKR7icikhrb/bVbGNXULKrnh7h2j9oZKExRCUF40j/Kiebz52n2sKzgdhAFNGNnu8zF195Gvk4J90TuA4hoP00sF3oQGJl33EluA47Mm4phluLpnMJiZWnEkUyvgtdf+zsYCPZFDM5jGre7lwvAPVmGz222HlPIu4C7QPf5sCzXaLFgwj+W/eA3LjBLijR2c8++Tz+PKULegcp/ywuvKStgQ2bGYXlNWtvedFP0cvnA+f1W6B+y77k0vK2TTBNA9FerJIVtWtO7V21AMzuZVz7Hd56OmrGxM32qPVZTu7T/jSffGUozfhL64ewrQgr64e6WUcsiROBPV8CsUCkU2GTMxfillUgjxTeBl9HTOe/dk9BUKhUIxsuQkYVdK+QKw9y5cCoVCoRhxxl/JmUKhUCgOCGX4FQqFYpKhDL9CoVBMMpThVygUiknGuGjLLIToALbmWo5dKAbG9uzDnVHyZp/xJrOSN7uMBXmnSCl3mwk5Lgz/WEQIsWyw/NixipI3+4w3mZW82WUsy6tCPQqFQjHJUIZfoVAoJhnK8O8/d+VagH1EyZt9xpvMSt7sMmblVTF+hUKhmGQoj1+hUCgmGcrwKxQKxSRDGf5hIIS4TAixRgihCSEWDtheJ4SICCFWpH/+MuC1w4UQq4QQjUKIPwkhBhtAM6rypl+7NS3TeiHEGWNB3l3k+5kQomXAd3r23mTPNUKIM9MyNQohbsm1PIMhhNiS/vuuEEIsS28rFEK8KoTYmP5dkEP57hVCtAshVg/YNqR8udaFIeQdP7orpVQ/e/kB5gCzgLeAhQO21wGrh9jnI2AR+sSxF4GzxoC8c4HPACswFdgEGHMt7y6y/wz43iDbh5Q9x7phTMsyDbCkZZyba7kGkXMLULzLtl8Bt6Qf3wL8bw7lOx44bOD1NJR8Y0EXhpB33Oiu8viHgZRyrZRy/XDfL4SoAPKllB9I/S//d+DCbMm3K3uQ9wLgESllTEq5GWgEjsy1vMNkUNlzLBPoMjRKKZuklHHgEXRZxwMXAPenH99PDv/mUsolQPcum4eSL+e6MIS8Q5FzeXdFGf4DZ6oQ4lMhxNtCiOPS26rQh8hnaE5vyzWDDbqvYuzJ+00hxMr07XTm9n4o2XPNWJVrVyTwihBiuRDia+ltZVLKNoD079KcSTc4Q8k3lr/zcaG7ORnEMhYRQrwGlA/y0o+klE8PsVsbUCul7BJCHA48JYSYxzAHyh8I+ynvUHJlXd6dhNiD7MAdwM/T5/858Fvg2tGWcR8Yq3LtyjFSylYhRCnwqhBiXa4FOgDG6nc+bnRXGf40UspT92OfGBBLP14uhNgEzET/j1494K3VQOtIyDng3PssL7pcNQOeZ+TKurwDGa7sQoi7gefST4eSPdeMVbl2QkrZmv7dLoR4Ej3U4BNCVEgp29LhvvacCrk7Q8k3Jr9zKaUv83is664K9RwAQogSIYQx/XgaUA80pW9Lg0KIo9PZMV8ChvLCR5NngCuEEFYhxFR0eT8aS/KmL/AMFwGZrIlBZR9t+QbhY6BeCDFVCGEBrkCXdcwghHAKIVyZx8Dp6N/rM8DV6bddzdjQ0YEMJd+Y1IVxpbu5XFkeLz/of8RmdO/eB7yc3n4JsAZ9xf4T4LwB+yxE/8NvAm4jXSWdS3nTr/0oLdN6BmTu5FLeXWR/AFgFrES/YCr2Jnuuf4CzgQ1p2X6Ua3kGkW9aWkc/S+vrj9Lbi4DXgY3p34U5lPEf6KHTRFp3v7In+XKtC0PIO250V7VsUCgUikmGCvUoFArFJEMZfoVCoZhkKMOvUCgUkwxl+BUKhWKSoQy/QqFQTDKU4VdkBSFE0YAuhd5duhZa9uN454+1zpdC7866eu/vHNaxHkvXguy6/ctCiNtG4hz7gxDitVx27VRkB1W5q8gKUsouYAHo7WqBkJTyN3vbTwhhlFKmBjneM4yxwqiRIt3mwyilbMriOQR6bYa2j7s+ANwA/PfIS6XIFcrjV4waQohT0g3tVqWbWFnT27cIIX4qhHgXuEwIcZMQoiHd7OqR9Hv6PV8hxBQhxOvp118XQtSmt/9N6LME3hdCNAkhLh1CjqfSzcrWDGhYhhAiJIT4byHEZ0KIpUKIsvT26ennHwsh/ksIERrkmEYhxK/T71kphLg+vb1CCLEkfaezekAjv4F8gQFVs0KIa4QQG4QQbwPHDNheIoR4PH2Oj4UQxwzY/qoQ4hMhxJ1CiK1CiOL0HclaIcSf0QsMa4QQ3x8g438OOPZVQoiP0nLemalIR/9n+/m9/W0V4wtl+BWjhQ34G/A5KeV89LvNbwx4PSqlPFZK+Qh67/VDpZQHA18f5Fi3AX9Pv/4Q8KcBr1UAxwLnAr8cQpZrpZSHo1cr3ySEKEpvdwJLpZSHAEuA69Lb/wj8UUp5BEP3WPkK0Jt+zxHAdeny/CvRK6cXAIcAKwbZ9xhgOfSX/f9nettp6L3cM/wR+H36HJcA96S3/wfwhpTyMOBJoHbAPrPQv6tD04/r0fv0LAAOF0IcL4SYA3wOvZHbAiCF/s8IKWUPYB3wHSkmACrUoxgtjMBmKeWG9PP7gRuBP6Sf/3PAe1cCDwkhngKeGuRYi4CL048fQB/YkeGpdDijIeOxD8JNQoiL0o9r0I1hFxBnR2Ot5eiGN3O+C9OPHwYGC1mdDhw84C7DnT7ux8C9QghzWrYVg+xbAXSkHx8FvCWl7AAQQvwTvfEfwKnAXLFjOFp+ugfPsehtOpBSviSE6Blw7K1SyqUDZDwd+DT9PC8t48HA4cDH6WPb2blhWztQif4dKSYAyvArRovwPrx+DvqEo/OBn6Rj4HtiYN+R2IDHu7XDFUKciG5AF0kp+4QQb6HfjQAk5I4eJin27foQwLeklC8Pcs7j0T/TA0KIX0sp/77LWyIDZIChW/Ya0nJHdjn+nsZkDvxeBfA/Uso7d9n/W8D9UspbhziGLS2jYoKgQj2K0cIG1AkhZqSffxF4e9c3CSEMQI2U8k3gB4AH3TMdyPvoXTBBD0m8uw9yuIGetNGfDRw9jH2WoodWGHDeXXkZ+Ebas0cIMVPoXTGnAO1SyruBv6KP69uVtUDme/kQOFHoWVFm4LIB73sF+GbmiRBiQfrhu8Dl6W2nA0Nl4bwMXCuEyEu/t0ro/flfBy5NP87Mup2SfizQZydsGeKYinGI8vgVo0UUuAZ4VAhhQg+B/GWQ9xmBB4UQbnQP9fdSSv8uTu1N6OGT76OHSK7ZBzleAr4uhFiJ3ilx6V7eD3BzWqbvAs8DvYO85x70GcyfpI1lB3p46ETg+0KIBBBCb3m9K8+n3/ea1HvP/wz4AL374yfo3wnon/v2tOwm9HWIr6OvCfxDCPE59H+mbUCQXf5hSilfScfzP0h/nyHgKillgxDix+gTugzoHSdvBLaih4CWSimTw/ieFOME1Z1TodgLQggHEJFSSiHEFcDnpZQjNldXCGEH3kRfXN0tlXUY+1uBlJQyKYRYBNyRXqQdCdn+CDwjpXx9JI6nGBsoj1+h2DuHA7elPXk/+ji9EUNKGRFC/Af6HNZt+3GIWuBfaW89zo5spJFgtTL6Ew/l8SsUCsUkQy3uKhQKxSRDGX6FQqGYZCjDr1AoFJMMZfgVCoVikqEMv0KhUEwy/j+K6erR2qqVSgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "colors = chemi._KELLYS_COLORS\n",
    "for i in range(4):\n",
    "    plt.plot(all_angles[i][12:] + all_angles[i][:12], np.concatenate((all_energies[i][12:], all_energies[i][:12])), color=colors[i])\n",
    "    plt.plot(all_angles[i][12:] + all_angles[i][:12], np.concatenate((all_energies[i][12:], all_energies[i][:12])), '.', color=colors[i])\n",
    "plt.xlabel('Torsion angles (degree)')\n",
    "plt.ylabel('Relative energy (kJ/mol)')\n",
    "plt.savefig('biphenyl_mm_torsion_scan.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get a system\n",
    "integrator = mm.VerletIntegrator(1*unit.femtosecond)\n",
    "platform = mm.Platform.getPlatformByName(\"Reference\")\n",
    "\n",
    "simulation = simulations[0]\n",
    "xyz = off_molecules[0].conformers[0].in_units_of(unit.nanometer)\n",
    "simulation.context.setPositions(xyz)\n",
    "flat_xyz = (xyz / unit.nanometer).flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(387.15522636173904, array([[  512.06140379,  -267.55877391,  -463.59216742],\n",
       "        [  627.46324387,  -336.03510451,  -258.81421554],\n",
       "        [ -274.18330607,   151.41230112,   -53.90962383],\n",
       "        [  -55.53046259,    24.16385982,   216.38021103],\n",
       "        [ -774.42013791,   501.7776067 ,  -255.90490377],\n",
       "        [ -272.2309811 ,    56.23735377,   910.92399018],\n",
       "        [-1039.7011379 ,   538.45777794,   678.17647218],\n",
       "        [-1057.34962256,   570.06114302,   631.76071524],\n",
       "        [ -352.50374849,   189.51027789,   217.09866871],\n",
       "        [  666.61072363,  -359.22475587,  -419.1614252 ],\n",
       "        [  408.62100332,  -219.37034897,  -256.40864825],\n",
       "        [ 1291.38398629,  -682.84559279,  -810.90455818],\n",
       "        [  120.71117643,   -62.97725248,   -73.27896483],\n",
       "        [ 1100.31167363,  -430.52461055,  -592.01922399],\n",
       "        [  937.619622  ,  -654.86094881,  -674.34685279],\n",
       "        [   36.89033338,   -16.82288027,   -79.74005479],\n",
       "        [  -21.64164485,    11.07347168,    26.85930492],\n",
       "        [-1132.0185625 ,   448.12567147,   657.15792121],\n",
       "        [-1023.92348899,   700.70444945,   689.34971333],\n",
       "        [  210.46855524,  -126.32261093,    25.98707822],\n",
       "        [  103.31108869,   -38.01069543,  -220.44832025],\n",
       "        [  -11.9498301 ,     3.02966906,   104.83494617]]))"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fxn_to_minimize(flat_xyz, simulation), jacobian_of_fxn(flat_xyz, simulation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -3.9966931 ,  -7.24618595,  -0.23349645,   0.        ,\n",
       "         0.        ,   0.        ,   0.        ,   0.        ,\n",
       "         0.        ,   0.        ,   0.        ,   0.        ,\n",
       "         4.29775292,   7.06724271,   0.86817156,   0.        ,\n",
       "         0.        ,   0.        ,   0.        ,   0.        ,\n",
       "         0.        ,   0.        ,   0.        ,   0.        ,\n",
       "         7.99896344,  14.15260351,   0.76521164,  -8.30002327,\n",
       "       -13.97366028,  -1.39988675,   0.        ,   0.        ,\n",
       "         0.        ,   0.        ,   0.        ,   0.        ,\n",
       "         0.        ,   0.        ,   0.        ,   0.        ,\n",
       "         0.        ,   0.        ,   0.        ,   0.        ,\n",
       "         0.        ,   0.        ,   0.        ,   0.        ,\n",
       "         0.        ,   0.        ,   0.        ,   0.        ,\n",
       "         0.        ,   0.        ,   0.        ,   0.        ,\n",
       "         0.        ,   0.        ,   0.        ,   0.        ,\n",
       "         0.        ,   0.        ,   0.        ,   0.        ,\n",
       "         0.        ,   0.        ])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "target_dihedral = 0\n",
    "torsion_indices = dih\n",
    "eq_cons = form_eq_constraint(torsion_indices, target_dihedral)\n",
    "eq_cons['jac'](flat_xyz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'type': 'eq',\n",
       " 'fun': functools.partial(<function dihedral_constraint at 0x7f2be828f440>, torsion_indices=[0, 8, 9, 4], target_dihedral=0),\n",
       " 'jac': <function autograd.wrap_util.unary_to_nary.<locals>.nary_operator.<locals>.nary_f>}"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq_cons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimize a scalar function of one or more variables using Sequential\n",
      "Least Squares Programming (SLSQP).\n",
      "\n",
      "Options\n",
      "-------\n",
      "ftol : float\n",
      "    Precision goal for the value of f in the stopping criterion.\n",
      "eps : float\n",
      "    Step size used for numerical approximation of the Jacobian.\n",
      "disp : bool\n",
      "    Set to True to print convergence messages. If False,\n",
      "    `verbosity` is ignored and set to 0.\n",
      "maxiter : int\n",
      "    Maximum number of iterations.\n",
      "finite_diff_rel_step : None or array_like, optional\n",
      "    If `jac in ['2-point', '3-point', 'cs']` the relative step size to\n",
      "    use for numerical approximation of `jac`. The absolute step\n",
      "    size is computed as ``h = rel_step * sign(x0) * max(1, abs(x0))``,\n",
      "    possibly adjusted to fit into the bounds. For ``method='3-point'``\n",
      "    the sign of `h` is ignored. If None (default) then step is selected\n",
      "    automatically.\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import show_options\n",
    "show_options('minimize', method='SLSQP')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Production Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simtk import openmm, unit\n",
    "from dask.delayed import delayed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "forcefield = ForceField('openff-1.0.0.offxml')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Molecule with name 'neutral' and SMILES '[H]c1c(c(c(c(c1c2c(c(nc(c2[H])[H])[H])[H])[H])[H])O[H])[H]',\n",
       " Molecule with name 'cation' and SMILES '[H]c1c(c(c(c(c1c2c(c([n+](c(c2[H])[H])[H])[H])[H])[H])[H])O[H])[H]',\n",
       " Molecule with name 'anion' and SMILES '[H]c1c(c(c(c(c1c2c(c(nc(c2[H])[H])[H])[H])[H])[H])[O-])[H]',\n",
       " Molecule with name 'zwitterion' and SMILES '[H]c1c(c(c(c(c1c2c(c([n+](c(c2[H])[H])[H])[H])[H])[H])[H])[O-])[H]']"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "off_molecules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0, 8, 9, 4], [0, 8, 9, 4], [0, 8, 9, 4], [0, 8, 9, 4]]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dihedrals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_openmm_system(offmol, forcefield):\n",
    "    \n",
    "    top = Topology.from_molecules(offmol)\n",
    "    system = forcefield.create_openmm_system(top)\n",
    "    \n",
    "    integrator = mm.VerletIntegrator(1*unit.femtosecond)\n",
    "    platform = mm.Platform.getPlatformByName(\"Reference\")\n",
    "    simulation = openmm.app.Simulation(top, \n",
    "                                       system, \n",
    "                                       integrator=integrator,\n",
    "                                       platform=platform)\n",
    "    return simulation\n",
    "    \n",
    "def optimize_geometry(index, offmol, forcefield, dihedral):\n",
    "    \n",
    "    simulation = create_openmm_system(offmol, forcefield)\n",
    "    xyz = offmol.conformers[index].in_units_of(unit.nanometer)\n",
    "    simulation.context.setPositions(xyz)\n",
    "\n",
    "    angle = torsions.measure_dihedral_angle(coords=xyz._value, dihedral=dihedral)\n",
    "    flat_xyz = (xyz / unit.nanometer).flatten()\n",
    "\n",
    "    target_dihedral = (angle) *np.pi/180\n",
    "    print(angle)\n",
    "    print(target_dihedral)\n",
    "\n",
    "    # jiggle xyz\n",
    "    minimum = get_energy(xyz, simulation)\n",
    "\n",
    "    print('initial_energy: {}'.format(minimum))\n",
    "    for i in range(5):\n",
    "        flat_xyz = flat_xyz + np.random.normal(0, 0.01, flat_xyz.shape)\n",
    "\n",
    "        new_result = minimize(fxn_to_minimize, flat_xyz, simulation, method='SLSQP', jac=jacobian_of_fxn,\n",
    "           constraints=[form_eq_constraint(dihedral, target_dihedral)],\n",
    "                        options={'ftol': 1e-10, 'disp': True, 'maxiter': 8000})\n",
    "\n",
    "        minimum = min(minimum, new_result.fun)\n",
    "\n",
    "        if new_result.fun == minimum:\n",
    "            result = new_result\n",
    "\n",
    "    print('final_energy: {}'.format(result.fun))\n",
    "    return result\n",
    "\n",
    "def run_torsiondrive(offmol, forcefield, dihedral):\n",
    "    results = []\n",
    "    for j in range(24):\n",
    "        result = delayed(optimize_geometry)(j, offmol, forcefield, dihedral)\n",
    "        results.append(result)\n",
    "    return delayed(results).compute(scheduler=\"processes\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing Molecule with name 'neutral' and SMILES '[H]c1c(c(c(c(c1c2c(c(nc(c2[H])[H])[H])[H])[H])[H])O[H])[H]' with dihedral indices [0, 8, 9, 4]\n",
      "Processing Molecule with name 'cation' and SMILES '[H]c1c(c(c(c(c1c2c(c([n+](c(c2[H])[H])[H])[H])[H])[H])[H])O[H])[H]' with dihedral indices [0, 8, 9, 4]\n",
      "Processing Molecule with name 'anion' and SMILES '[H]c1c(c(c(c(c1c2c(c(nc(c2[H])[H])[H])[H])[H])[H])[O-])[H]' with dihedral indices [0, 8, 9, 4]\n",
      "Processing Molecule with name 'zwitterion' and SMILES '[H]c1c(c(c(c(c1c2c(c([n+](c(c2[H])[H])[H])[H])[H])[H])[H])[O-])[H]' with dihedral indices [0, 8, 9, 4]\n",
      "CPU times: user 1min 37s, sys: 1.56 s, total: 1min 38s\n",
      "Wall time: 35min 6s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "all_energies = []\n",
    "\n",
    "for offmol, dihedral in zip(off_molecules, dihedrals):\n",
    "    print(f\"Processing {offmol} with dihedral indices {dihedral}\")\n",
    "    prod_energies.append(run_torsiondrive(offmol, forcefield, dihedral))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "energies_1 = [i.fun for i in prod_energies[0]]\n",
    "energies_1 = energies_1[12:] + energies_1[:12]\n",
    "energies_2 = [i.fun for i in prod_energies[1]]\n",
    "energies_2 = energies_2[12:] + energies_2[:12]\n",
    "energies_3 = [i.fun for i in prod_energies[2]]\n",
    "energies_3 = energies_3[12:] + energies_3[:12]\n",
    "energies_4 = [i.fun for i in prod_energies[3]]\n",
    "energies_4 = energies_4[12:] + energies_4[:12]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEOCAYAAACn00H/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABKkklEQVR4nO3dd3gc1dXA4d/ZVS+WXGTLvdu4V4oBG1NMh1ADgdBbQqgJoQRCTMhHIJBACC1A6KH3Fjo2xWBcMS5YcpGbLEuy1duuds/3x4zMWpbklbzSSvZ5n2cfaWdmZ86uxzo79945V1QVY4wxprk80Q7AGGNMx2QJxBhjTItYAjHGGNMilkCMMca0iCUQY4wxLRIT7QDaUrdu3XTAgAHRDsMYYzqUBQsWFKpqRv3le1UCGTBgAPPnz492GMYY06GIyLqGllsTljHGmBaxBGKMMaZFLIEYY4xpEUsgxhhjWsQSiDHGmBaxBGKMMaZFLIEYY4xpEUsgxhhjWsQSiDHGmBaxBGKMMaZFLIEYY4xpEUsgxhhjWsQSiDHGmBaxBGKMMaZF2jSBiMg0EXlbRDaJiIrI+fXWi4jMFJFcEakSkVkiMqreNvEi8i8RKRSRCnd/fdryfRhjjGn7K5AUYClwNVDVwPrrgd8BVwL7AvnAxyKSGrLNfcCpwC+AqUAn4F0R8bZe2MYYY+pr0wSiqu+r6h9U9VUgGLpORAS4BrhTVV9T1aXAeUAqcJa7TRpwEfB7Vf1YVRcC5wBjgSPa7p0YY4xpT30gA4FM4KO6BapaBXwBHOgumgTE1ttmA7AiZBtjjDFtoD0lkEz355Z6y7eErMsEAkBhE9vsQEQuFZH5IjK/oKAgUrEaY8xerz0lkDpa77k0sKy+RrdR1UdVdbKqTs7I2GlOeGOMMS3UnhJInvuz/pVEd366KskDvEC3JrYxxhjTBtpTAlmLkyBm1C0QkQSckVZz3EULAH+9bfoAI0K2McYY0wZi2vJgIpICDHGfeoB+IjIe2Kaq60XkPuBmEfkRyAJuAcqB5wFUtURE/gPcLSL5wFbgH8AS4JO2fC/GGLO3a9MEAkwGPg95fpv7eBo4H/gbkAg8CHQG5gJHqmpZyGuuBWqBl9xtPwXOVdVAawdvjDHmJ6K6q/7pPcfkyZN1/vz50Q7DGGM6FBFZoKqT6y8P6wrE7YuYAgzA+dZfACxU1dWRDNIYY0zH0WQCEZGDcMqOnIBzA18JTgmSLkC8iKwBHgUeqdfMZIwxZg/X6CgsEXkLeBlYBxwJpKpqV1Xto6pJwFDgL8DhQJaIzGhsX8YYY/Y8TV2BfAScrqq+hlaq6hpgDfC0WzG3VyvEZ4wxpp1qNIGo6oPh7kRVlwHLIhKRMcaYDqE93UhojDGmA2n0CkREyth1DSoAVLVTxCIyxhjTITTVB3JFm0VhjDGmw2mqD+TptgzEGGNMx9KsUiYichgwEqdpa5mqzmqNoIwxxrR/4d6J3ht4A2dGwFx3cS8RmQ+crKq5jb7YGGPMHincUVj348wEOERV+6pqX5wbCQPuOmOMMXuZcJuwZgDTVXVt3QJVXSMiV+FUwzXGGLOX2d37QIIRicIYY0yHE24C+RS4X0T61i0QkX7AP7ErEGOM2SuFm0CuApKANSKyTkRygNXusqtaKTZjjDHtWFh9IKq6AZjoVtzdBxBguaraNLLGGLOXatZ9IKr6MfBxK8VijDGmAwk7gYjIBOBQoDv1mr5U9foIx2WMMaadC/dGwuuBO3Eml9rCjkUW955J1Y0xxmwX7hXItcCvVfXfrRmMMcaYjiPcUVgebLiuMcaYEOEmkIeBC1ozEGOMMR1LuE1YtwHvi8hi4AfAH7pSVS+McFzGGGPauXATyP8BRwILgc5Yx7kxxuz1wk0glwNnqepLrRmMMcaYjiPcPpAqYFFrBmKMMaZjCTeB3AtcIyLSmsEYY4zpOMJtwpoKTAOOE5Hl7NyJfmKkAzPGGNO+hXsFUgi8DnwG5AFb6z0iQkS8InK7iKwVkWr3519EJCZkGxGRmSKSKyJVIjJLREZFKgZjjDHhafIKREREHW11D8gNwG+A83CGC48FngZqgNvdba4HfgecD6wEbgU+FpHhqlrWRnEaY8xeb1dXIPki8oSInCgiiW0Qz4HAO6r6jqrmqOrbwNvA/uAkNOAa4E5VfU1Vl+Ikm1TgrDaIzxhjjGtXCeQEnCarvwKFIvKmiFwgIt1aKZ6vgENFZB8AERkJHAa8764fCGQCH9W9QFWrgC9wko8xxpg20mQCUdVvVfUPqjoKGI/zB/5CYJOIfCEivxORIRGM5y7gWWC5iPiBZcDTqvqQuz7T/bml3uu2hKzbgYhcKiLzRWR+QUFBBEM1xpi9W7id6Khqtqreo6pTgT7AUzgjs74XkaUiclwE4jkDOBenOWqi+/vlInJR/XDqPZcGltXF/aiqTlbVyRkZGREI0RhjDDRzRsI6qloAPAE8ISJJOGVOfBGI527gHlV90X3+g4j0B24C/oPTnAbO1caGkNd1Z+erEmOMMa0o7CuQxqhqpaq+6U53u7uSgEC9ZQF+inMtThKZUbdSRBJw7lOZE4HjG2OMCdOuhvGWsevCibXAZpyO7dtUtXg34nkHuFFE1uL0f0wAfgs8A6CqKiL3ATeLyI9AFnALUA48vxvHNcYY00y7asK6Iox9eHCalC7AaUo6ezfiuRLnfo+H3H1tBh4D/hyyzd+AROBBnMrAc4Ej7R4QY4xpW6K6e5XZRSRWVf0isi/wnqp2j0xokTd58mSdP39+tMMwxpgORUQWqOrk+svD6gMRkdsbWR4HvOY+XY5zVWCMMWYvEG4n+kUiclXoAhGJxamP1Q9AVStU9bYIx2eMMaadCncY7zHAZyKyVVX/6155vIFzP8hhrRadMcaYdiusBKKq34vIScC7IlKNczd6b+AwVY1YNV5jjDEdR3PuRP8S5w7xF4BeWPIwxpi9WqNXICLydiOrCoEK4Km6CQptQiljjNn7NNWE1djVxYetEYgxxpiOpdEE0oaTSBljjOmAdrsWljHGmL1TowlERD4RkYN3tQMRSReRm0XkysiGZowxpj1rqg/kOeAFEanCmVZ2Pk5tqmqcGlQjgYOBo4E3ceYqN8YYs5doqg/kKRH5L3A68AvgIiCtbjVO6ZIPgQmqurK1AzXGGNO+NHkjoar6ccqkPw8gImk4lXC3uuuMMcbspZo1I6GqlgAlrRSLMcaYDsRGYRljjGkRSyDGGGNaxBKIMcaYFrEEYowxpkXCnZHwPhEZ3drBGGOM6TjCvQLZF/heRL4TkUtFpFNrBmWMMab9CyuBqOpBOHeefw78CcgVkWdE5JDWDM4YY0z71ZwJpVaq6g1AX+BMIAX4SESyReRGEenSWkEaY4xpf1rSiR4LdMIpa+IF1gPnAOtF5KwIxmaMMaYdCzuBiMhkEXkIp6Di34BvgaGqeriqjgJuBu5tnTCNMca0N+GOwvoBmIPTfHU+0F9Vb1bVtSGbPQ9kRDxCY4wx7VK4tbBeBp5Q1U2NbaCqBdh9JcYYs9cIK4Go6u2tHYgxxpiOJawEIiJPNLJKcSaYWgW8pKq5kQrMGGNM+xZuE1YGMBUIAkvdZaMBARYApwB/FpGpqro40kEaY4xpf8Lts/ga+B/QR1Wnqeo0oA/wPvAR0B94D/j77gYkIj1F5GkRKRCRahFZHnrDojhmikiuiFSJyCwRGbW7xzXGGNM84SaQq4E/q2pl3QL39/8DrlVVH3AXMH53ghGRdJxkJcBxwAjgSiA/ZLPrgd+5y/d1130sIqm7c2xjjDHNE24TVgrQE1hRb3mmuw6gtBn7a8z1wGZVPTdk2fahwiIiwDXAnar6mrvsPJwkchbw7908vjHGmDCFewXyBvAfETldRAaISH8ROR34D/C6u81+QNZuxnMSMFdEXhKRfBFZLCJXuIkDYCBO0vqo7gWqWgV8ARzY0A7d4o/zRWR+QUHBboZnjDGmTrgJ5FfAh8BzwGpgjfv7B8Dl7jYrgEt2M55B7v7WAEcB/wTuBH7jrs90f26p97otIet2oKqPqupkVZ2ckWH3ORpjTKTssslJRGKA6cAtOH0Pg3H6KFapakXddhEafeUB5qvqTe7zRSIyFCeBPBCyndYPs4FlxhhjWtEur0BUtRanmSpFVStUdYmqfh+aPCJoM7C83rIVQD/39zz3Z/2rje7sfFVijDGmFYXbhPU9MKQ1A3F9DQyvt2wYsM79fS1OEplRt1JEEnDuUZnTBvEZY4xxhZtAZgJ/F5GTRKSviHQJfUQwnnuBA0TkZhEZ4nbUXwU8CKCqCtwH3Cgip7jT7D4FlOMUczTGGNNGwh12+57783V27Guo63vwRiIYVZ0nIicBdwB/xJlr5I/AQyGb/Q1IxEkqnYG5wJGqWhaJGIwxxoQn3ARyaKtGEUJV3+OnhNXQesW5IprZRiEZY4xpQLjVeGe3diDGGGM6lubMSDhGRB4Qkf+JSE932UkiMqH1wjPGGNNehTsj4ZHAPKA3cBhOHwQ494T8qXVCM8YY056FewVyO/BbVT0Z8IUsn4VTwsQYY8xeJtwEMgqndHt924BIDuM1xhjTQYSbQIpwmq/qmwhsjFw4xhhjOopwE8jzwN0i0gfnvo8Yd5Kne4BnWis4Y4wx7Ve4CeQWnDIi63Dm/1gOfAZ8hTOplDHGmL1MuPeB+IGzReRWYAJO4lmkqtmtGZwxxpj2q1kzCKrqapz5QIwxxuzlwk4gInIGcDhO6fQdmr5U9cQIx2WMMaadCyuBiMjdOHORfw7kYpM3GWPMXi/cK5BzgV+o6qutGYwxxpiOI9xRWB5gcSvGYYwxpoMJN4E8CvyyNQMxxhjTsYTbhJUOnCUiM4AlgD90papeFeG4jDHGtHPhJpCR/NSEtU+9ddahbowxe6FwbyRssxkJjTHGdAxN9oGISNwu1ouIDIpsSMYYYzqCXXWiV4lI97onIvKViIRW5e0OWDkTY4zZC+0qgUi95+OA+F1sY4wxZi8Q9pzoTbBOdGOM2QtFIoEYY4zZC+0qgSg7XmHUf26MMWYvtathvAKsEZG6pJECLAl5bv0fxhizl9pVArmgTaIwxhjT4TSZQFT16bYKxBhjTMdinejGmD1SXrmPhXnl5JX7oh3KHqtdJxAR+YOIqIg8ELJMRGSmiOSKSJWIzBKRUdGM0xjTvuSV+3g7aytzc8t4J3urJZFW0m4TiIgcAFyCU/031PXA74ArgX2BfOBjEUlt2wiNMe1VbrmPgPt7rTrPTeS1ywQiImnAf4GLgKKQ5YIzte6dqvqaqi4FzgNSgbOiEKoxph3qnhS7/XcP0CulybJ+poXaZQLBmcDqVVX9rN7ygUAm8FHdAlWtAr4ADmy78Iwx7VllbRAAj0C3pFgyLYG0irATiIhcLiLLRKSyrgKviNwoIj+PZEAicgkwBPhjA6sz3Z9b6i3fErKu/v4uFZH5IjK/oKAgcoEaY9qtVduqSI71MCYjmcIqPzWBYLRD2iOFlUBE5BrgFpwrg9CbBzcBV0QqGBEZDtwBnK2qTTVa1r8bXhpY5myo+qiqTlbVyRkZGRGK1BjTXlXXBtlQWsOQzokMTE8gqLChpCbaYe2Rwr0C+RVwiar+E6gNWb4QiOQIqClAN2CpiNSKSC1wCHC5+/tWd7v6Vxvd2fmqxBizF1pTXE0QGNolkR7JsSTEeFhbUh3tsPZI4SaQ/sDSBpb7gcTIhcObwBhgfMhjPvCi+3sWkAfMqHuBiCQAU4E5EYzDGNNBZW+rIi3eS7fEGDwiDEiLZ31pDQG1Mn6RFu6c6GuAicC6esuPBZZHKhhVLQaKQ5eJSAWwzR1xhYjcB9wsIj/iJJRbgHLg+UjFYYzpmMp9AXLLfUzOTMEZtAkD0hL4cWsVm8t89OlUfzojszvCTSD3AA+ISBJOf8MUETkH556MC1sruEb8Deeq50GgMzAXOFJVy9o4DmNMO7O6qAqAIV1+ahjp0yker0BOSbUlkAgLK4Go6pMiEoPTwZ0EPIvTgX6Vqr7UivGhqtPrPVdgpvswxpjtVhVV0y0xhs4JP/1pi/UIfTrFs7akhoP66PYrE7P7wh7Gq6qPqWp/nA7rTFXtq6r/ab3QjDEmfCXVteRX+hnaZedu2YFpCZT7Amytqm3glaalwh3Ge6+ITARQ1UJVzW/dsIwxpnlW1TVfdd45gfRPc5qucmw0VkSFewWyPzBfRFa4BQ4HtGJMxhjTLKpKdlE1PVPiSInz7rQ+KdZLj+RY1hZbAomksBKIqh4IDMapT/VLYLWIfCkil4lI59YM0BhjdmVbVS1F1bUM6ZzQ6DYD0xIorKql3BdodBvTPM3pA1mrqn9R1ZE4VXDn4pQbyW2t4IwxJhzZRVUIMDi98QQywF1nzViR09JiirFAPBAHWDo3xkSNqrKqyBmimxi7c/NVnfR4L2nxXnKsGStimlNMcZiI3CYi2cBXwHDgOqBHawVnjDG7sqXCT5kvwNAmmq8ARIQBaQlsKvfhs+KKERHuKKz5wArgBOBhoK+qHqmqz6hqRWsGaIwxTVlVVIVXYGATzVd16oorri+14oqREO6d6B8B56jqitYMxhhjmiPoNl/1T0sgzrvr78N1xRVziqsbHO5rmifcO9H/0NqBGGNMc+WW+aiqDTY5+iqUR4T+neJZW1JNQBWv3ZW+WxpNICJyP3CTqla4vzdKVa+KeGTGGLML2UVVxHqE/mnhJRBwmrFWbqtic7mPPqlWG2t3NHUFMgZntFXd78YY024Egsqa4moGpicQ4wn/SqJPapxTXLG42hLIbmo0gajqoQ39bowx7cH60hp8Ad3l6Kv6Yr0e+qTGk2PFFXdbuKOwbnVLuddfnigit0Y+LGOMadqqoioSYjz0bkGJ9gHpCZT5Amyz4oq7Jdz7QP4EpDSwPMldZ4wxbcYfCJJTXM3g9IQWdYQPcIsr2lS3uyfcBCJAQ/NBTgC2RS4cY4zZtbUlNdTqjhNHNUddccWcErsfZHc0OYxXRMpwEocCa0QkNIl4gQTgkdYLzxhjdrZqWxXJsR56JsfueuNGDEhLYG5uGeW+QIMVfM2u7eo+kCtwrj6eAG4GSkLW+YAcVf2mlWIzxpidVNcG2VBWw5iM5N3qAB+QFs/c3DJySqoZnZEcwQj3Hk0mEFV9GkBE1gJzVNXfJlEZY0wj1hRXE1QanHmwOTonxDjFFUtqLIG0ULh3os+u+11EMnGq8IauXx/huIwxpkHZ26pIj/fSLTHcSkwNqyuu+ENBBb5AMKxSKGZH4Q7j7SQiT4tIFbAJWFvvYYwxra7CFyC33MeQLom7bL5aNW8D79z7BTmLG5+yaEBaPEGFDVZcsUXCTbl/B8YBJwHVwFnA74GNwBmtEpnZI/grFlOZ9xj+isXRDsXsAermPR+6i0KI3766hAcvepnPnpjHv85/gaxv1zW4XWZKHAlesdFYLRTuNeAxwC9U9UsRCQALVPUlEdkMXAa82moRmnbDX7EYf9k8YlP3JTZ5POBM5qOBEoL+LQT9BQT9+dsftVWrqK1YBChsjqHTkEeISz0gqu/BdEx15172tqPolphIekIMqn6C/kL3vHPOP19FHp8+XcmXryW4Nx4IwQA8ee2bnPu3ExkxdeAO+/WIU0crx4ortki4CSQdqEvhJUBXYBXwDfB45MMy7Y2/YjElWRcCPtjswZswBA1WEvTng/p22l686SAx/HT7UC2la64kKfNXJHQ7E4/XOi1NeHzliyjNvpByzaCA4xkX8wJbf3gdrd1K6O1puWu68u5j0yjYlM7QiVtY80NXggEPHk+Q+IQKHr38dSYcPYyTbziM1G4/nX8D0uJZua2KvHIfva02VrOEm0BWA4OA9TgTS50pIt8Bp2A3Eu7xNFhDxaZ7cUZuAwTRYCWxyePwxHYPeWRs/ymeeCfpZF8M6ge8xCQMpTL3Xqq2PEFi9/NIyDgLj7ehAgfGOIL+QirWzwT8rGc6AP283xGXOm37ead047NnSvn8mTWkdE3i0oeOYsjEEpa9fwPrVnSl/z4F9BoWx9dvpPP120F+/HoVJ153OPufPAYRoW+neKe4YkmNJZBmEtWGbjCvt5HItUBAVe8XkcOAd3Eq9XqAq1X1gdYNMzImT56s8+fPj3YYHYq/Ygnl628hUL0G595RBYkjbejj25uxmjL/w1fYmF9Cn+5pTD7qdPwVP1CZ9zD+0i8QbycSM84hIeNsPDGdWvutmA5EVakpeo+KjXegwSo0qHzEI8RRxknDM7afextXbOGFWz4kN6uAfU8cxUnXTyfJLe1evux5Kla8T/KIY0keeSa+kk9ZN+8p3n2kFxuyMhk0PpEzbv853Qd04/1V29hWXcvZozKsuGIDRGSBqk7eaXk4CaSBnfUDJgPZqvpDBOJrE5ZAwqfBGio3P0hV/lN4YjNI6fdn8n3pbCraQu/OPejTbdT2bf2BIBX+IOXVNZSVFlNeXkZ5dTVbiqsoSusHCBKsZVJpDvsePg2A2splVOb9G1/JZ4g3lYSMs0nMOAdPTFqU3rFpL4L+Aso3/BlfyefEJI8jpd/tLM9ay9e1IxnpWcch4w8g4A/w8eNz+fjRuSSnJ/Lzm6cxbGwigZI8aovzqFozj9Ivn4ZgEGLi6H3NayQOmoxqkOqiT/j6v6/zyXP9qPXHcvh53elzxsl8lVvBz0d0o2tiy+9u31NFNIF0VJZAwuOv+IHy9TcTqF5DfNdTSe59HflV8byVVUBQBUHp7N+GX4UqbxK1MTuPiPH6K/HU1uBPSANxBvt5/FUMqSli6kETto+5r638kcq8R/CVfIJ4kknIOIvE7ucSqFm3U4e92bM5Vx3vULHxTjRYQ3Kvq0jI+CXrs5fxfmlnEC8SqGG/5c9S9MV3SFUh3br6SY4rQ6tLm9x3fL9xZF70b2IzBrjHClK4+gPeuOtLVnzbjR4jAqTdeA6TupYxOn62nXf1NDuBiMhvw925qv5jN2JrM5ZAmuZcdTxEVf6T26864jodhD8Q5M15yyiM7QoioEp8RR7pRdkkBqtJCPrQsmrKN5aRv7KYsjw/Kmn0OqgXJUfPIOiNQYIBUgqzKMscg7emnD6VRRw6bRKJsU4NotqqLOeKpPgjkDjQWiAIEh92c5npuAL+fMrX34a/dDYxyRNI7Xc73oQBlNf4eHPeMsqSMp0Ng7UMnv8wPRa8QkL3PnTqN5CY9B5403oSk57pPNIy8RfnsfnhcyDgd85Z8YAGSZ18Mp2Pvpq4zKGAk0gWv/s6b9yTRfoVJ5CY4qP7nHsZPKaAUcfdZeedqyUJJNwbBFVVB+1OcCHHvAmnY344UAN8izOt7tKQbQSnhPylQGdgLvAbVV22q/1bAmmcc9VxC4Hq1cR3PYXk3r/H400lJ28Ls9dupTK+MwSduRM8wVoOzv+M0up9WfJJNllz1xOsDZLeI4UxRwxl3BHDGDihFx6vh4WfzmZj/hb6dO/BiInjmPfC0+R178HWQYfg8VeRWbyFQ6dOolOSU9ygtmoVZTnXE6jOciPzkNTzSpIyL4nSJ2Nak6pSs+1tKjbdhQZ9JPe6moSMswAPS1av5bttSkBiAEFF8ARqSXznM0753S9I6bzTFEU7qFozn6qsOSQOO5DYrn0p/uQhSr58BvVXkzLxRDofcw3xvUY425ZV8cx/PyG432RWX/06lJdz0T2pjDj88tb/EDqADtGEJSIfAi8C83CKOP4ZmAKMVNVt7jY3ALcA5wMrgVuBg4HhqlrW1P4tgexMgz4q8x6kakvdVcdtxHU6mEpfgNlLlpLj6U5SSQ4j180isGY2pb0m0il3EbM/HUleaS+69klj7AwnafQdnYknjKlFfRXlLHjhKTalplIw+HAkGKBb4SYOPXgCXTslukOGLwCc0msp/e8ioctxrfxJmLYW8G2hfMNM/KVfEpM80b3q6M/Wsko+W5pNYXx3umxegP/d+WwuTabL6DSKlpUyfvpUjrh4/5Yds6yQok//TcnsJ9CaCpLHHUuXY68lvu8Y3nnyQzaOG8uWJ7+jZFY2Yw/Zxvn/vAXxNp2o9gYdIoHUJyIpOPednKSq77hXH7nAA6r6f+42iUA+cJ2q/rup/VkC+Ym/YjHVW9/BX/YVQd8m4rucTHKf6xFPCivWb+SbvGr8nngGrn4f2dSFr9/bRmJ1Fr3TNrCppC+pIw7gpN9Pp9fwlo9aCdRUs/jlp1jniSF/2JEoQuctGzhkvzFsWrucrPUF9Ev8mBH95tFp0IPEpkyI8KdgosFXvojq/Kfxlc0BlOSezlVHUIW5K7JYUpVEjK+czO/e4OtnAgS8nQj4A6gqMXExXP7Y6QwY32u3YghUFFH8+WOUfP44wapSkkYfQfWQc3jPO4RAdS1bHv+WmtUFHHvxOg779W14YrtE5s13ULuVQNw/3L8GfgMMBEar6hoRuRFYo6ovRzpg97g9cRLGVFX9SkQG4dyTsp+qzgvZ7j2gUFXPa2Afl+I0d9GvX79J69Y1XNJgb+LcmHUB4DRJJfX+PUndz6Okyscn3y8nP64HaVuW4H/jMxZ/k4F4PPQc2o281VvRYDBi/4nrBGt9LH/jeVZVVbJln2MIeuMg6JyXGggydtuD7NPrE1IH3kN82qEROaaJDl/ZN5SuugwIAkLKgHtI6HwUm7YW83n2JsriupCx6hM2PLmELUX9mHbWBA45dxIFOUWsmr+BIZP7Ruy8AwhUlVIy6z8Uf/Yo25L7Mu/4/zid9arw1hxWvrmeA3+Ww0k3X0dsYv+IHbej2d0Ecg1wPXAXcCcwyk0g5wCXqOq0CMdbd9yXgaHAZFUNiMiBwNdA/9AKwCLyBNBbVY9qan92BeJ0GpZknUPO97ms+zGT/vvk03XoGSzcMJa8zEwk6KfLrOf58YMY+u03jlHTBrPPwQNITk8kZ3Fuq/wn3h5boJbs919innShtOc4p/MzWEv3nGwOH/MotZXLSen7RxK6nR7xY5vWF6wtofjHUwn689wlHmJ6XMv8/P3IDnYhoWwT8W+9zIq5PTjwFwcw7ZyJJKftXsn2sGOrLuer2R+wrOv+4IkBVUbGlbD5le+Z+8YmRh24gbPvupDE9LFtEk9701gCCfdO9F/hJIr3ROQvIcsXAqMaec1uEZF/4PRtHKyqgXqr62e9xqbcNSFUA5Sv/xM53+fy3J1HE/B7SRzYiQGXjSfQpy9d135BevZmBh12ASf/rjfe2B1naRswvlerJI464o1h2AlnU/biv5kXGIF640C8BEs3kdz/MSo3XU/5htsI+PNJyrzcbvjqQIL+QkpWXUqwthCIRTXIRv9UFq2bQk1sCt0XvcK61/MZecJJ3PSnidtvBmwrnoQUMr3VrAjWEsQZtZVd7uGEa4+lS685/O9BeOSS/3LBP7aS3teuguuEm0D6A0sbWO4HIv4VQUTuBc4EDlXVNSGr6r66ZAIbQpZ3B7ZEOo49iQb9lK27EV/xh3w/5yJ6jIHOJ03A328o3sqtjFv5EpNPOBfvqdFv6x253yT0ud9Q0Ht/CvtOpXD8ETz5ny858ZiZJHf5F1V5DxP0F5DS9xZEdm9OCNP6Ar7NlK66hIB/C50GPcT8bwtYGdOL6vQ+pBT+SMIb75A8+mf8+qX9SezUtokjVN9B+zD5+SvZljGO+IotrJ70a95ctY3DThjH2b0688KtX/DgJZ9x4d8L6TnKroIh/CasZcAtqvqGO0/6OLcJ6xrglw1d2rQ4IJF/4iSP6aq6ot66uk70f6nqHe6yBJxO9N9bJ3rDNFhN6drf4i/9gqULr2DB/CBJF//MGRsfrGXS+s/Z75Rzox3mDuqGYEqXvsxeV8rmwUfCkmWMS+vB2P0+ozr/MeI6TSd14N2Ip22aOUzzBarXUbLqYjRQRqfBD/HFnHKyug4BQIK19Jz9PkdedgGJ7aQGVejQ34ryUj7aHKAkYzSTkspJKw7y5NVvERtXw3l/7crgg36911wFN9aEFe58IPcAD4jI2TjNRVNE5E/A/wF3RzDIB4ELgF8ARSKS6T5SwLnhBLgPuFFEThGR0cBTQDnwfKTi2JMEAxWUrr4cX8mXfP3B1Xz0hpBy3lE4/4wOX7fM6AXYiMRBk+ly9FV03u9kjj/qWIavfAXGjOAHbwWvPDQMT9of8JXOpiT7YoK1xdEO1zSgtiqL4uxzUa0mbeiTrNjQmawu7i1j4tzXIb27tJvkAT+dd4mDJtNt7GH8bGRveq6fxYLKFHLi/FzxzDl4Y+J57NpSFr91Bzu3ru9dwkogqvokMBO4A0gCngUuBq5S1ZciGM/lQCrwKbA55HFdyDZ/A/4BPAjMB3oCR+7qHpC9UbC2hNJVl1BdvIj3n7uKRUu99P/jDGJqq5CAD4K1eIK19M1sfwkkVEynDA497TIOyHkV6ZGBf0Y/HvlzkMLS26mtWkFJ1jkEajZFO0wTwl/xAyXZ5yMSQ9rQp/h2hfJ1ZRzJRWuQ2hrn3AvU0qd7j2iH2qTkviM5bvphDF/9BmvoyjfFhVz23GV07R3Ls3+KY9ajf0SDVdEOM2qafR+IiHQDPKqa7z7vq6obdvGydmFvasKq67SsKNrIW49eRGlKMp3PP4Tk4hyO6lJBVUovNuTl0Tczk/7DOsbIElVlzadP80XcWHzxncl/7EsmjEhiv0PvJSY2jqTe1xD05VsdoyjzlX1H2ZorkJgudBr8OJ8sKGRtYg+6rpnFpK69KAmyvTrBxMMPiXa4YQn6qlj47mMs6HMCCcFqZgzuxZs3PEv2vBqmnb6ZY64+kWDVyj323Iv4jYQikgn8EbhQVTtEI/TekkACvjxKV11MUV4Jr9x/FrHjM0g88SC65C3k6KHdSBvasrt424uChf/j44IYSjJGU/HG13hXlXLCRS+T3nUjIFY/K4p8JV9QuvZavPF9SB70KO98s5H8Tr3oseQNpoyfSs/xY6IdYoupKtkfPcmXyZMJxqVweO9kvvnXe8x/r4jBYzfQe2gBA0cVMurYPa+GVov6QEQkXUT+KyIFIpIrIleJ40/AGmA/4MJWitm0QKBmPSXZ55G3xsczd5xByhH9SDzxIHqun8WJ44d0+OQBkDHxGE4cnkHP9Z+TfPJBxB/eh8f/dBSzXxvP1++MZmNWGv6yebvekYmomqIPKV1zFd6EwSQMepJX5zjJo/ecp5k+/cQOnTwARIRhR13IcfEbSCxZz4d5fkZefCgTj/SxeklfvnhtAs/99TBWz/0m2qG2mV2NgbwDmAY8DRwN3AvMAJKBY1R1duuGZ5qjtmo1pasuZu3SNF57+AgG/GoUgX32YdCadzns8OOJ7dwz2iFGTMqA8RzXqTtffPwaWRNPZXCPVBa+0J34/l347uVczsmoZZ/23bWzR6ne+jrl62cSkzwe6f1PXvpmHdUpGfT99BEOu+i3JHXZc+Z5ydz3OE5cvYCPl33Lt3IQnacdQELOOhKH9aBqxRZWzoN9Dot2lG2jySYsEVkHXKSqn7hlRFYB96vqNW0UX0TtyU1YtZXLKVl1KUvn9Oe95w5g2I374u/Zl9FrXufA48/Dm7hnzvgXrK5g4dsPM3/g6agnFoJBNBAk47t/c+I5hxDf+choh7jHq8p/lopNdxGbeiA1Xf/K28u2oDEx9P7sWY649hbiEtvPKKtI8hWu5/PP/8eawSc6E1epooEgWx/6gItv3o+ew8dHO8SIaekw3l7AcgD3hr5q4LHIh2d2R1XhyxSt/CVz3h3Oey8dxPDbDibQPZP91r3OwSdftscmDwBPQjKTTruW3nkLnQVeLxLjwd9pKGU511G97d3oBrgHU1XK1s+kYtNdxKTsS1Hyn3lzZREe9dP/81c46vqZe2zyAIjr1o8Zx59F5pb5TtkdrxeJEdKHp/DABe+zZt630Q6x1e0qgXioq6ntCACVrReOaa7K/GdZ8cnDPHXbDL75dn+G3n4owaQkphV+xMSTfoPExEU7xFYnHi/DEyuQQA1oEDxeKuJ7seKHIylfdxPVW1+Pdoh7HFWlbN0N1Gx9FYCVm9P5YGMNCeV5DFr4JTNu/APeGO8u9tLxeRJTGRlTjAT9oAriZeLoauISAzzy69ks/WRWtENsVbvqAxHgORGpcZ8nAI+JyA5JRFVPbI3gTNNqij7ix0+f5tWXzyX9sIH0278/3poiDg9k0f/YS6MdXpvqO2gf9n3+SrZmjKeqU29yh53I4rxeeDamMZxbUfWT2O2MaIe5R1ANULHhL/iK3merjuCHqrMpSNmftNxFDC4oYv+rr4h2iG2q3+ARTP7v5eT1P4y8wcewPPMQfv5QFW///hWe+N08Tv9DFVPOOCbaYbaKXfWBPBnOTlT1gohF1Ir2pD6Q6q2vU5Yzk6ce/zXJFxyPeL2A0m/hBxx3yUXRDi8q6spQeNO6k/XN+3y/3w34YlIZ4fuYMan3ktz7BhK7nxPtMDs0VT/l626mpuh9SlIu5aOSk53qtcEAozfNZerPTo12iFFRd+6VE8+s2NFUpvVjSnIRX9z+Imt+6MpRlw3kqN+c3GFLn7SoGm9HSQx7m6r8p6nYdDefv/EzUg4dDTHuP2MwQPwe1HHXXImDJpM4yDnHxw87kJRHfs3icZfxY69jKfP3Y7+NN6HqI6nH3plgd5cGayhb+zt8pbOIzbyOL1dNguS6PyFKTecB0QwvqurOvS7Aidnf8cGqpczxjGffm04n6eHn+fDfUFrwAqfdeiYeb7gVpNq/Peed7AVUlYrcf1Gx6W7mzj6dwkGHwJAhzlzlblmSfr4Vu97RXiC2az+GXftf9l/yGIMWPsIm70g+4Sk2b3qdys0P055n4myPNFDp1FQrnYX2mMm7qydQldQFCfi3n3veJTaqHyB96H78bNxgMjd+ySJ/JgMvPpMpJ27km9c38+TVT+Or9u96Jx1Eu57SNtI6chOWapCKjXdSXfg8S5aex6pOBxHoO5B9Vr5IyqpPKO4xni753zPyrJnbv4UbUH8Na+6/hMKqYr6ffjckJTKR+xiR2Z+knld12CaFthSsLaF09a+prVxGVde/8vHGQdR6Ysh87366SxalvSaSlruI7gddxcCj9sy2/pbwF2/h08/eZ+2Ao+nr34Tnm2f5+Im+DBiTyiUPndfmc57sjg45J3qkddQE4rQ7/5GaondZueE3LE89iGByCpNzP2Di8RdTvX7J9hLUljx2psEgWQ/+ntqcD1kw/Z9U9B5Jfz5iSsZW0vtca0mkCXU11QI1aylMupvZxYOJrS4i4eVXOfoPN1DywxwK531Gt30Ps+TRgGB1BV+//yxL+59A15ot9Mh+hnfu60XXPkkcf81RbFm7rdVm+IwkSyB0zASiwRrKcq7DV/I5K7fewNKUA4nxVzC1eilDj/iF/fFrhiUP3EH88gdZOfE6Nk08g1TWMybxe6pihtKncw/6dGuVyTU7rIBvMyWrLibgy2eN/I2FtcNIzV+O7+VvOPWu35LWPSXaIXYIGqhl0ftP8l3mkST7Sxi65Vle/7+u+GtiQSA2zsPlj5/RrpPI7s4HYqLAmcvj1/hKPmdp6Z9Zkn4oSaUbOMK7mWEzzrLk0Uxjr/gDRUNvYui8vzPmgxupIYM5VSewqGw4769PZmPhsmiH2G4EqnOcMvn+EhbV/IOFgX3osmo2Zc9+z+n3XGfJoxnEG8OE4y9meunXVEs8yzMvYvwptYCCgt9Xy8o530U7zBaxBNJOBWuLKV11MTVli5hX9ndWpE6h87pvObRrMv0PPCra4XVYB1x7BRv730ba2i/JXPyyc/MXHgIay4atG6MdXrtQW7WS4uxzqQ7E8EX5PayOHUK3uS+T/0YZ5z70G1K7JkU7xA5HRNjn8LM42rMK8VdRfvgvGHxKOqPPSqXXGKFH34557lkCaYcC/nxKss+nsmI9X5TfT07KWLrNf50p/UfTZ9K+0Q6vwzvipgvJ6jqT1FXfOhNraRDEw8YiobJscbTDiyp/xfeUZF9Aub8nn5bfSUFcLzq//Sibv+7MxY+dT3J6h5i5od3qu/9xHJ9RTVzVNrwnHotvxrGkXnUqC2YV4avYEu3wms0SSDtTU/QRxStOoqSihk8r/0VhUn86v/cfxo48iv6T94l2eHsEEeGk239JxZZUJr13GUPmP0jfJc9SGDuSt7JL2Zb3drRDbHP+isWUrbuVkuwLKagex8e+26mReOKe+DeF64Zz6b/PbFdTz3ZkGSOm0L82F1DweNGYOLRnb5659i5qyhZF/Hj+isWUb7wLf8XiiO/bEkg7UlX4EmU5vyW/vCef+u+h2ptM3JOPMWzSqQw/eEi0w9ujeGO9JI44nJS8ZQxY/ATD5/6DMZ/eSrn25e3cAaxf8wAa9EU7zDbhK19ESdYF1Gx7nTUlU5glNxJbtY2Kv/8XX3AClz50KgnJe35NtbY0qN9Ad1rpACBw4IHkVQ3nrTvup6rw5Yjdp1S99Q1Kss6luuBZSrIujHgSsVFY7YAG/VTmPUTVlsdZkv9zsrqeT1zVVsruf40JJ53MtLMnRjvEPVLO4lxeveZeeiTm4AvEsf+Qhfi7duf74x/CF5PIgYkvMHLopXhiMyJyvJriT6mtWk1cp/3azYx1wdpiSrLPZ9MmYRkXsK3LvqRvXsjqe74jc/R4zrvnBGLjd1Uyz7TEuqwlbNi8mQR/Kcu8falM6w+zvmRg8FUO+vk+pPS5GfG07KqvtnotlZsfwFf8YchSD0k9ryQp85Jm78+G8dI+E0ht9RrKcm6komIjc0v+QH76JFBFamvo9d3nnHilTfjYmnIW55I9bz0FOcV8/94Cjt1vDt1S1/H98Q9R0mkw473PM2nwEcSlTmjR/oO1JdQUvUtVwfMEa9a5S+NJG/afqCcRX+lXlK/7I6u3DGN++h/B44VgLTUPPEdal1Gcc9dxxMTu+RV124Oq/HV8/u0c1vWZTtyWHAYVPMKE/WtIHXgv3rjwJ4IL+DZTmfcwNVvfAk8c8elHUrPtA6AWJLbFUz23qBaWaT2qSnXhS5Rv/DurSw5nacLt+Dslb+/QVW8MXcb1jnaYe7wB43ttH38/YuoAXvhjCiP7rmLKe79h2YE3snjALynJfp9pfVeTlHFaWPtUDeIv+4bqrW/gK/kU1I8ntnvIFjX4Sr+JWgLRQCUVuf+gvOBVVpRcyI/pp4D81Jrde2o/TjnveLwx1sLdVhK79+eY4/uw5NNXmddpNFldb8P/43NMrDmDTgPvIS51vyZfH/Rvo3LL41QXvggoCRlnktTjEjyx3Ujodjr+snnEpu4b8XPOEkgUBP2FlK3/IwX5OSyo/ivb0kaRXrAM7w8rKZp6LEFvDJ5ALUn+jlPqYE8w4eh96Na3M09c/SbrvunBGckvkli8lrXjL6IiZwHTK+6gS7/rEE/D/QGBmk1Ub3uDmq1vEfRvRrxpJHT7OQldTka1mpLsi0FrAMVX8imJ3c/B423b+yn8FT9Qvu4mNmztzCL5D5VpveiyeRFF3Uag7nk3qH9vSx5RIB4v42acQbcl85m9OpvV/S+ieMs4ptT8nm79LyIh45yd7v0KBsqpyn+a6vyn0WA18V1OJKnn5XjjfropMTZ5fKt9WbEmrDZWU/wpRTn/x9L8Y1jV5UxifOV0W/gZc54oJ6BxZI4M0GV0J4qWlTJ++lSOuHj/qMa7NyotrODJa94i5/tcfn76NvzpxSyfciMpwU1M7/Q8PYf8Ea97RaHBanzFn1G99XX85d8CQmzqFBK6nkJc2qE7tGH7KxbjL5uHqp+qvEfxJg4jbfDDeGK7tvp7UvVTmfcoWze+wqKSS8nt4tyUuk/RKt6+u5D0gT4779qRzT9u4vM3XqLs8JPxBqqZVPM3hgzIIKXvTMSbhAZrqC58kcq8x9BAMXHpM0jqeQUxCYNbJR7rAyG6CUQDlZRvupM1q9eyKOYaqpJ70S17Nj8+mkV5WRoHnTmewZP68PR17xLwB/DGern8sdPbdXmDPVmtr5ZX/vwJ3721jAMOS2bo8B9YOP4KxKscFHM3vbp0Ba3FX/4dGijFE9ebhK4nE9/lZ2G1WftKvqB07W/xxmXSacijO3xjjPh7qV5Lac4fWLGxD8tTLiXoiWVY4Vwqlyfz6TMr6ZSRQnlxFRoI2nnXjqxesJF3736CzhdOp6LbUPoXvsmkzP+RkNyH2sqlaO02YlOnkNTramKTRrdqLJZAiF4C8Vd8T97K21lQ+DM2dzuCxOJ1VL7wCeuXpHLwLyYy/dxJpHRx7u7NWZzLqvkbOkSBtT2dqjL72YW8/ffZ9BqSzvFnFDE34yCqUnvRY+tshnpfp++ADBJ7XERsyr6INK/Zx1++kNI1v0E8SXQa8miD3x5353yo62dbn/UKC3yXU5I2km4FSxiTGMs7928kd2UBB5w6hpOun87mrEI779qhRR/8yEt/fJ19r+hF/ugZJFRtoXvZHAbEfsagcaeRmHFWm8RhCYS2TSD+isX4Sr/FX72JpVk1LE+7BPXEkvjl+6x6u4oDTjuAQ86ZaHf2dgA/fp3D079/F2+Mh2N+VcySYec6I5Y0yNDi2UydeirxcS0cblm1kpJVl4IG6DT4YWKTxwBQvq2Sd+79ku/eXAqAJ8bDxfefxIipA3e5T3/FYmqKP6eyeDELN05gbdfTiKspZqJ/NSXrMnnv/m9ISInjjJlHMvpQu7+ovfvsyXm884/ZHHJlHJsnnQYioEFGlH7GwdPOJqYNJqiyBELbJRBf+UK+e/t+NsuBVPYYTWX6QFLWLaLgxaWMnH4w086a2KHmAjCQv3Ybj1/1JpkH11Jz+DHONK6qIILXX8HAmg2MGzaU7pl9mr3vQM16SlZdSrB2K/64v/H1awG+e3MZ/praHbbzeIXxRw1nymljGTy5T4PFNGtKvyT7u3+yyn8chWn74Y9PZ2DBd4wbMJI37llC9tz1jJo+iDNmHmU1rToIVeX1Oz6j2LMa34wdz724mmKGaAHjRo4mPb1Lq8VgCYTWSyCqSmmNn025y9m0YRVbSuMp6zURcP6DJ89+mx5JY5h29mQrB9GBVZVW8/xfnqH21BkEvTFIoJYReQsoS4hnU+cxBGPi6Vq6mlGd4xg+aiIxMeEPcsz5fgUfP/g0K77tiifGy+QTRjFsv368OPMjAv4AHq+HUYcMZuW366guqyFjQGemnDqGfU8cRVKaUlP8Kbk5c1iRP4AN3U9whuVqkIlbXyc2MINXb/+EQG2Qk66fzgGnjrFKzh1MMBDkhb88S/nR05xzL1jL2IofyPPFsqWrc9WaWb6GMb26MmjICDyeyF6V7HEJREQuB34P9ASWAdeo6pdNvaalCWTRolyyNhTTq1M8nVPj2bKtiq1VfsoVauIE7ZyExP00tNPrqyAQm+j8Jw7W0nPFLE4655xmH9e0Px8/+i2LZ39Fl9Gd2La0lMLsOEZNH8SQicn40orJSe5PVUpP4qu3MTiwmXEjx5DetUeD+1JVVs5Zx2dPziN77noSUmKZdMQGJh36BT3H/p6Ebqft1Afiq/Lz/cdZfPPK96xfnkvvfWPpcVBnqgeMpSa5h1tdGKeZI1hL10Wf8s39W+k/tidn//UYMvp1bsNPy0RS9rfreetfL9PZPfeK1ibQY3BXeg+pJXaYl8J+4/AldiGpYgvDYksZO3oCySmpEelX3aMSiIicATwHXA585f68ABipqusbe11LEsiiRbl841fwSt2xt6/zlhWRWrKK1OJskovXEF9VTnKnAZQkDmLFiGnb7+fYt3w9Ew8/pNnv07Q/OYtzeeiSV6j11eLxehi6Xz9yswooLagAIKNfKiOOisc/pAdbe4xBNEhmaRajM9PxemLYlL+FPhkZFK6K5/On5pObVUBa9xQOOWciU04bS3xigNK1v8Vf9hVJvX5LUo+fKhGo1rJt0+esyv6BTb7eFKTvSzA2Ea+/ksS136OrNpPRswvrxxxA0BODJ1hL2f2vMeaggzjikv3t3o4O7pPH5/Le/V+B03pF31GZxCXGUrCuiJL8cmLjahlyqJfEaWMo7zMWT20NqTlLKFm8idhOXkpXlHPKtWe2KInsaQlkLrBEVS8JWZYNvKqqNzX2upYkkJffWUFhzzTEIxAMkrbuG4Yue5yUbauJo4qEoWNJmXAGyWOOJSb1p5pJCz+dzcb8LfTp3sOSxx6m/jc6VSV/7TZWfrOO7G/Xkz1vAzUVPrr29TPo5N5UjpyEP7GzU2UAQJWkrO/wlpTTNSOFjF5pJMTHER8XT0JCInHxcQSKn8BT8xHF/hPIKR1MAKEkthclnZ1ZExMq8uhZuZqBGd3oM3AcP36dxzevLmH1go30Gqvbr5CmnnIUB5w6JoqflomUui8vDQ3zr6n0U7ihiIKcIiehFGZR29tLyT5TCMYmuuWRfAxcupijLjqp2cfeYxKIiMQBlcAvVPWVkOUPAqNVtdG/1i1JIPM/fIUFnSdv/0Y38eMr6ds/jdQJ55A0/DAkxqqUmh0FaoNsWJpH1rfryJq7no3LNzDiqkFU7rO/2zeheAI1BD2xzmiupridpQDJpWvpX5vDkP4j6DloLB7vzq9982+fM/vZhYBzqGOvPNhuCtyDNLc56sN332RN5uTtdc6GblnEEced2Ozj7kkJpBewCThEVb8IWX4rcLaqDq+3/aXApe7T4cDK5hxvUEbM0LiYxE6B2GS8/gp8tVWlawpqs3fvXeygG1AYwf21Nou3BbqndxuckJmZ7gzBVKrz8orziwtXN7Dp9ngH9M8cHUzKiHcGYyieyoKanHV5S5s6TpwnPjk9rtsw3BcV+wqzfMGaioi/oUZi7iD2mnjTO6V069RrQP+68640N2ddcWl5S/bVX1V3KkvdkWth1c980sAyVPVR4NE2iagFRGR+Q5m9vbJ4W1dHixc6XswWb+R0xF61QiAAZNZb3h3oeHNCGmNMB9XhEoiq+oAFwIx6q2YAc9o+ImOM2Tt11CasfwDPish3wNfAr4BewCNRjapl2m3zWiMs3tbV0eKFjhezxRshHa4TvY57I+H1ODcSLgWuDe1UN8YY07o6bAIxxhgTXR2uD8QYY0z7YAnEGGNMi1gCaQMicqmIfC4ixSKiIjKggW1y3HWhjzvrbdNPRN4RkQoRKRSR+90786MVc2cReVZEStzHsyKSHq2Y6x13VgOf54vNjb8ticjlIrJWRKpFZIGITI1WLKFEZGYDn2VeyHpxt8kVkSr3sx/VhvFNE5G3RWSTG9v59dbvMj4RiReRf7nnaIW7v+bX5o9MvE818Hl/G614m2IJpG0kAR8BM3ex3Z9xBgXUPf5St0JEvMB7QCowFfgFcBrw98iHC4QX8/PAROAY4Gj392frVkYh5vqeZMfP87J665uMvy2JUyD0n8AdwAScIen/E5F+0YinASvZ8bMMLbB1PfA74EpgXyAf+FhEUtsothScgTRXA1UNrA8nvvuAU3HO0alAJ+Bd9xxu63gBPmHHz/vYeuvbMt7Gqao92ugBTMa5W35AA+tygOuaeO0xQBDoG7Lsl0A10KmtYwZGuMsPCll2sLtseDRjdo8zC3igifW7jL+Nz425wGP1lmUDf23rWBqIbSawtJF1AmwGbg5ZlgiUAZdFIdZy4PzmxAekAT6cUkh12/R1z92j2jJed9lTwLtNvCZq8dZ/2BVI+3KdiGwVkcUicnO9pp4pwApV3RCy7EMgHpjUplH+FE85O968+TVQARwYsk00Yz7TvcRfJiL31PvGGU78bcL9d56Ec8UX6qO2jqUJg9wml7Ui8qKIDHKXD8SpCrE9dlWtAr6gfcQeTnyTgNh622wAVhC993CwiOSLSJaIPCYi3UPWtZt4O+qNhHui+4FFwFZgP+BOnJP/Ynd9JjuXammsrEtbyAQK1P36A6CqKiL5IfFEM+bngXVALjAK+Cswjp8qGIQTf1vpBnjZ+bPaAhzRxrE0ZC5wPvAjTsmgW4A5bj9C3WfVUOy92yrAJoQTXybOOVm/yOAWovN/6wPgdWAtMACnKfszEZmkqjW0o3gtgbSQiPwFuHkXmx2qqrPC2Z+q/iPk6RIRKQVeEpEbVHVr3WaNvTycY0Q65kaOW7+o5W7FvMOOmxG/OkU06/wgImuAuSIyUVUXNhFDg0U520hYBULbmqr+L/S526G7BjgPqOvcbZexh2hJfFF5D6oaOtjjBxFZgPNl6DicxNKYNo/XEkjL3YczK2JTGp0dMQxz3Z9DcK5K8oCD6m3T2DfXxtxH5GLOA7qLiNR9ixcRATJC4olEzKHuo+Xxz8f51jYUWEh48beVDlUgVFXLRWQZzmf5prs4EwhtqmwvsdeNFmsqvjycc7IbUFBvm6hXt1DVXBHZiPN5QzuK1xJIC6lqIa07p8B49+dm9+c3wC0i0kdVN7rLZgA1OMUldynCMX+DM5pkCj/1I0wBkkOe73bMoXYz/jE4/+lCP89dxd8mVNXnfsucAbwSsmoG8FpbxhIOEUkA9gE+x2lmycOJdV7I+qnA76MVY4hw4lsA+N1tnne36YMz0CLqBVpFpBtOc1vdudt+4m3LHvu99YHz7Wc8cBbOJeax7vMu7vopwLXusoHAz3EmzXorZB9e4AfgM5xhnke42/wrGjG72/zPjekA9z38ALwTrZhDjjsYuBVnBNkAN/YVOFce3nDjb+Nz5AyckTUX4/wh+CdOJ3//dnD+3gMc4p6b+wPvAqV1sQE3uM9PAUYDL+L0PaW2UXwp7rk5Hme20lvd3/uFGx/wsHtuHuGeq58Di0PPl7aI1113j3s+DgCm43zZ2RiteJt8L9E+OfeGB84wSG3gcb67fiJOW3IxzrjwH93XJNXbTz/3P28lTrPWv4D4aMTsbtMFp0mp1H08B6RHK+aQY/YFZrvHqwFW4fxB7lJvu13G38bnyeU4w7nrrtCmRfvcdeOq+4Prc/9ovQaMDFkv7vmyGWeI9myc6aXbKr7pjZyrT4UbH5Dgnptb3XP1HUKGn7dVvDhDjD/EuVfFh9P38VT9WNoy3qYeVkzRGGNMi9h9IMYYY1rEEogxxpgWsQRijDGmRSyBGGOMaRFLIMYYY1rEEogxxpgWsQRiOhx3QqAHoh1HpLiTHS1tg+OcLyKf7WKbB0RkVmvHEikicoWIvB3tOPZWlkDMbmlg5rT6j6da4bCnADe1wn73WG7J+L8At0U7lgh7DJgs7WT2xr2N1cIyu6tnyO/H4/yHDl3W2IxrDRKROFX1NbWNqm5rzj4N4MwEWa2qs6MdCICIxKqqf3f3o6o1IvI8cBXw5e5HZprDrkDMblHVvLoHTimW+svOFJFVIuJzf14S+nr3KuU3IvK6iFQAd4hIrDhzp+eKSI2IbJCQ+eHrN2GJM7f50yJS5M55/UnonNdu0025iBwuIkvdOaQ/F5GBTb03EfmtiCxxt98kIo9LyJzp4e5XRG4SkS3uts+IyJ9EJGcXx75ARJaLMz96lohcKyKekPWXucurRaRARD4Ukaa+EJ4F7NDUIyJecSbaKnIf9+HULwvdRkTkehFZ7X62P4jIL+tts7+ILHRjWSQix7r/rtPd9dPd58eKyHci4gOOCnPfvcWZwKouxvdEZCg7ehs4UUSSmvpMTSuIdp0de+w5D5xvuRry/GScqqFXAMNw5qT2AyeEbKM4dX8uBgbhFOz7HU7p7Wk4tbQOBC4Iec0sQqarBd7CqR82Dafq7tvu6xPd9ee7x/0EZ7KusTiTd324i/dzDXAYTlG7Q4AlwLMh63e5X+BMnPpLF7ufwU1ACZATss1MQqaMBS7Bqdt0mvt5nIBTUfYKd/1koBY4G+iPM1HWtUBME++lGDir3rLr3Vh+jlNd9184NcFmhWzzfzjzoR/txnIWzqyNx7nrU3BKij+PM3HXDGCZ++863d1muvv8B+BI9985I4x9JwFZOLWgxroxPo5THyopJMYknHL4h0f7/8De9oh6APbYcx7snEC+Bp6ot81TwFchz5V61XlxZmf8FJxabQ0cZxZuAsGZI0EJKTyIM2d0CXCx+/x86s117v7x9QGeZry/o3EKHXrC3S9OJdVH6u3nI5pOIOuBc+q95hpgufv7Ke77C6vaLZDuxnloveW57DhXuMf9gz3LfZ6M0wQ5td7r7gPed3+/DNiGm6zdZXUVnKe7z6e7z08N2SacfV+IMy+8hKz34hQQ/Hm9120DLor2/4G97WF9IKY1jQCeqLfsK+DEesvm13v+FPAxkCUiHwHvA/9T1WAjxwji/KEGQFVLROQHYGTIdjWqujLkeS7OvNLpOH98diIih+FcMYzASUpeIA6n1H1umPvdB6dfKNRcnKuRho6ZgVNN+N8i8nDIqhicqrLgfDbrgLUi8iFOQnpdVcsa2idOhVdwroTqjpOG01cV+rkFRWSue3xwPr8E4AMRCa26GotTNRj3/S1VZ57x0PfXkNB/53D2PQnnyqRMREI2IQmnZH+oqpD3adqIJRDT2hoq91x/WcUOK1UXisgAnG/8hwFPA9+LyIwGkojQuNDj1DayrsF+QBHpD7yH88f/VpxvvROBF3CSSHP225yS13Wv+xWNTA6kqmUiMhGnyW4GTpK7Q0T2VdXcBl6y1Y2hczPiCI3lBHae6bGuA7w506iG/juHs28PzhwXZzawr/pJvws7zs5n2oB1opvWtAI4uN6yg4Hlu3qhqpap6iuq+mucuaAPw5net77lOOfxlLoFItIJpy9kl8dpwmScRHGtqn6jqllArxbs50ec/pFQ9Z9vp6pbcObcGKyqq+o/QrarVdXPVPUmnP6BZJxRcA3t04fzWYwMWVaC089yQN0ycb7mh8a2HKfJrn8Dsaxzt1kBjBGR0G//jb6/Zu57Ic6/eWED22xPICIyGOdqZuFORzGtyq5ATGu6G3hFnOlaP8K5ojgbpw2/USLyW5w/botxvo2ehdO5u7H+tqqaLSJv4TT5XIrTWfx/7vbP70bs2TiJ6RoReR3nD+01LdjPP4EnRWQezjDTk3Fm9Stq4jUzgX+JSDFO810sztVPb1X9q4gcj9OE8wXON/FDgVScP+aN+RAned9TL7abRCQLp4P7cpxmrc2w/UrnHuAeN7l8gdNpfgAQVNVHgf/i3F/ymIjcgZNk/+Duv9Erk2bs+zrgLRG5FedKpS/wM5x+pWx3d1OBNSHPTVuJdieMPfacB/U60d1lv8KZEdDv/ryk3noFTqu37BKcb5NlOIlgNnBgyPpZ7DgKqzNOM1cRTlv4J8CokPXnA+X1jjHdPXa3Jt7PVThXA1U4nfo/d18zoDn7xfmDmo8zRe0zwJ3AipD1MwnpRHeX/cL9DKrd9/UVcKa77mCcKUy3urEtJWSUWiPvZR93X6FTEscA9+Ik3WKcUVgPs+MoLMEZPVd3xVCA0wczI2SbA3BGn9W4P091P4P9m/qsw9x3D+BJ9/OrwZnj/Il6n++HwI3RPv/3xofNSGhMGxORN3CG3J7Qxsd9EVimqre38nF+BrwBdFfVwlY+1mic5D5MnWY504asCcuYVuTe3PZr4AOcDvdTcZpgTo1CONfjNKFFlIicB6zBufdmNM5Q3HdaO3m4egHnWvKIDrsCMaYVuZ3L7wATcIaZZgN/U9X/RjWwCBKR6/mp/yQPZ/TaDdr4sGKzh7AEYowxpkVsGK8xxpgWsQRijDGmRSyBGGOMaRFLIMYYY1rEEogxxpgW+X8N7npzh6ftZAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(angles[12:]+angles[:12], np.asarray(energies_1)-min(energies_1), '.', color=colors[0])\n",
    "plt.plot(angles[12:]+angles[:12], np.asarray(energies_1)-min(energies_1), color=colors[0])\n",
    "plt.plot(angles[12:]+angles[:12], np.asarray(energies_2)-min(energies_2), '.', color=colors[1])\n",
    "plt.plot(angles[12:]+angles[:12], np.asarray(energies_2)-min(energies_2), color=colors[1])\n",
    "plt.plot(angles[12:]+angles[:12], np.asarray(energies_3)-min(energies_3), '.', color=colors[2])\n",
    "plt.plot(angles[12:]+angles[:12], np.asarray(energies_3)-min(energies_3), color=colors[2])\n",
    "plt.plot(angles[12:]+angles[:12], np.asarray(energies_4)-min(energies_4), '.', color=colors[3])\n",
    "plt.plot(angles[12:]+angles[:12], np.asarray(energies_4)-min(energies_4), color=colors[3])\n",
    "plt.xticks(fontsize=14)\n",
    "plt.yticks(fontsize=14)\n",
    "plt.xlabel('Torsion angles (degree)', fontsize=14)\n",
    "plt.ylabel('Relative Energy (kJ/mol)', fontsize=14)\n",
    "plt.ylim(0, 110)\n",
    "plt.savefig('biphenyl_mm_min_torsion_scan.pdf', bbox_inches='tight')"
   ]
  }
 ],
 "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}