demo-with-scan.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import Any, Optional\n",
    "\n",
    "import ase\n",
    "import numpy\n",
    "from ase.calculators.calculator import Calculator\n",
    "from sella import Constraints, Internals, Sella\n",
    "\n",
    "\n",
    "def sella_optimize(\n",
    "    atms_obj: ase.Atoms,\n",
    "    calc: Calculator,\n",
    "    order: int = 0,\n",
    "    ints_obj: Optional[Internals] = None,\n",
    "    in_place: bool = False,\n",
    ") -> ase.Atoms:\n",
    "    \"\"\"Optimize a geometry using Sella\n",
    "\n",
    "    :param atms_obj: An ASE Atoms object\n",
    "    :param calc: An ASE Calculator object\n",
    "    :param order: 0 = minimum | 1 = saddle point\n",
    "    :param ints_obj: A Sella Internals object, possibly involving constraints\n",
    "    :param in_place: Modify the atoms object in place?\n",
    "    \"\"\"\n",
    "    if not in_place:\n",
    "        atms_obj = atms_obj.copy()\n",
    "\n",
    "    atms_obj.calc = calc\n",
    "\n",
    "    # Initialize and run the optimization\n",
    "    dyn = Sella(\n",
    "        atms_obj, order=order, internal=(True if ints_obj is None else ints_obj)\n",
    "    )\n",
    "    dyn.run()\n",
    "\n",
    "    return atms_obj\n",
    "\n",
    "\n",
    "def sella_optimize_minimum(\n",
    "    atms_obj: ase.Atoms, calc: Calculator, in_place: bool = False\n",
    ") -> ase.Atoms:\n",
    "    \"\"\"Optimize a minimum-energy structure using Sella\n",
    "\n",
    "    :param atms_obj: An ASE Atoms object\n",
    "    :param calc: An ASE Calculator object\n",
    "    :param in_place: Modify the atoms object in place?\n",
    "    \"\"\"\n",
    "    return sella_optimize(atms_obj, calc, order=0, in_place=in_place)\n",
    "\n",
    "\n",
    "def sella_optimize_ts(\n",
    "    atms_obj: ase.Atoms, calc: Calculator, in_place: bool = False\n",
    ") -> ase.Atoms:\n",
    "    \"\"\"Optimize a TS/saddle-point structure using Sella\n",
    "\n",
    "    :param atms_obj: An ASE Atoms object\n",
    "    :param calc: An ASE Calculator object\n",
    "    :param in_place: Modify the atoms object in place?\n",
    "    \"\"\"\n",
    "    return sella_optimize(atms_obj, calc, order=1, in_place=in_place)\n",
    "\n",
    "\n",
    "def sella_optimize_constrained(\n",
    "    atms_obj: ase.Atoms,\n",
    "    calc: Calculator,\n",
    "    const_coos: Optional[list[tuple[int, ...]]] = None,\n",
    "    in_place: bool = False,\n",
    ") -> ase.Atoms:\n",
    "    \"\"\"Optimize a structure subject to internal coordinate constraints using Sella\n",
    "\n",
    "    :param atms_obj: An ASE Atoms object\n",
    "    :param calc: An ASE Calculator object\n",
    "    :param const_coos: Optionally, constrain a set of coordinates\n",
    "    :param in_place: Modify the atoms object in place?\n",
    "    \"\"\"\n",
    "    if not in_place:\n",
    "        atms_obj = atms_obj.copy()\n",
    "\n",
    "    # Set up constraints\n",
    "    const_obj = Constraints(atms_obj)\n",
    "    for coo in const_coos:\n",
    "        if len(coo) == 2:\n",
    "            const_obj.fix_bond(coo)\n",
    "        elif len(coo) == 3:\n",
    "            const_obj.fix_angle(coo)\n",
    "        elif len(coo) == 4:\n",
    "            const_obj.fix_dihedral(coo)\n",
    "\n",
    "    # Set up internal coordinates\n",
    "    ints_obj = Internals(atms_obj, cons=const_obj)\n",
    "    ints_obj.find_all_bonds()\n",
    "    ints_obj.find_all_angles()\n",
    "    ints_obj.find_all_dihedrals()\n",
    "\n",
    "    # Note: For whatever reason, the Atoms object must be identical to the one passed in\n",
    "    # to Constraints and Internals, so we *must* set in_place=True here\n",
    "    return sella_optimize(atms_obj, calc, ints_obj=ints_obj, in_place=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import automol\n",
    "\n",
    "\n",
    "def scan_relaxed(\n",
    "    geo,\n",
    "    calc: Calculator,\n",
    "    scan_coos: Optional[list[tuple[int, ...]]] = None,\n",
    "    scan_grids: Optional[list[numpy.ndarray]] = None,\n",
    "    const_coos: Optional[list[tuple[int, ...]]] = None,\n",
    "    gra: Optional[Any] = None,\n",
    ") -> tuple[numpy.ndarray, numpy.ndarray]:\n",
    "    \"\"\"Perform a relaxed scan with constraints using Sella\n",
    "\n",
    "    :param geo: An AutoMol molecular geometry data structure\n",
    "    :param calc: An ASE Calculator object\n",
    "    :param scan_coos: The set of coordinates to scan over\n",
    "    :param scan_grids: The grids of coordinate values for the scan\n",
    "    :param const_coos: Optionally, constrain a set of coordinates\n",
    "    :param gra: The molecular graph, for setting coordinate values\n",
    "    :return: The geometries and energies, as arrays over the scan grid\n",
    "    \"\"\"\n",
    "    all_const_coos = list(scan_coos) + list(const_coos)\n",
    "\n",
    "    shape = numpy.shape(scan_grids[0])\n",
    "    geo_grid = numpy.empty(shape, dtype=numpy.object_)\n",
    "    ene_grid = numpy.empty(shape, dtype=numpy.float_)\n",
    "\n",
    "    for idx in numpy.ndindex(shape):\n",
    "        for coo, grid in zip(scan_coos, scan_grids):\n",
    "            assert len(coo) == 2, \"For now, only allowing bond coordinate scans\"\n",
    "            geo = automol.geom.set_distance(geo, coo, grid[idx], gra=gra, angstrom=True)\n",
    "\n",
    "        atms_obj0 = automol.geom.ase_atoms(geo)\n",
    "        atms_obj = sella_optimize_constrained(atms_obj0, calc, const_coos=all_const_coos)\n",
    "        geo = automol.geom.from_ase_atoms(atms_obj)\n",
    "\n",
    "        geo_grid[idx] = geo\n",
    "        ene_grid[idx] = atms_obj.get_total_energy()\n",
    "\n",
    "    return geo_grid, ene_grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import interpolate\n",
    "\n",
    "\n",
    "def scan_saddle_point_index(scan_grids, ene_grid) -> tuple[int, ...]:\n",
    "    \"\"\"Determine the index of the saddle point from the scan data\"\"\"\n",
    "    assert len(scan_grids) == 1, \"2D scans not yet implemented\"\n",
    "\n",
    "    (scan_grid,) = scan_grids\n",
    "\n",
    "    # Underscores mark values sorted by scan coordinate value\n",
    "    idxs_ = numpy.argsort(scan_grid)\n",
    "    ene_grid_ = ene_grid[idxs_]\n",
    "\n",
    "    # Check for delta sign changes to find the roots\n",
    "    deltas = numpy.diff(ene_grid_, prepend=ene_grid_[0])\n",
    "    signs = numpy.sign(deltas)\n",
    "    (root_idxs_,) = numpy.where(numpy.diff(signs))\n",
    "\n",
    "    saddle_idx_ = None\n",
    "    for root_idx_ in root_idxs_:\n",
    "        neg_before = numpy.all(deltas[root_idx_ - 1 : root_idx_ + 1] > 0.0)\n",
    "        pos_after = numpy.all(deltas[root_idx_ + 1 : root_idx_ + 3] < 0.0)\n",
    "        if neg_before and pos_after:\n",
    "            saddle_idx_ = root_idx_\n",
    "\n",
    "    saddle_idxs = (idxs_[saddle_idx_],)\n",
    "    return saddle_idxs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.\n",
      "/home/avcopan/code/other/sella-demo/.pixi/envs/default/lib/python3.12/site-packages/ase/optimize/optimize.py:372: FutureWarning: force_consistent keyword is deprecated and will be ignored.  This will raise an error in future versions of ASE.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:27:34     -525.805458       1.3084       0.0000       0.1000       1.0000\n",
      "Sella   1 12:27:37     -525.903020       0.3646       0.0000       0.1000       1.0946\n",
      "Sella   2 12:27:38     -525.915607       0.2396       0.0000       0.1000       1.3557\n",
      "Sella   3 12:27:39     -525.920654       0.1239       0.0000       0.1000       0.9202\n",
      "Sella   4 12:27:40     -525.921163       0.0491       0.0000       0.1000       1.2178\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/avcopan/code/other/sella-demo/.pixi/envs/default/lib/python3.12/site-packages/ase/optimize/optimize.py:372: FutureWarning: force_consistent keyword is deprecated and will be ignored.  This will raise an error in future versions of ASE.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:27:46     -523.901807       7.5976       0.0000       0.1000       1.0000\n",
      "Sella   1 12:27:48     -524.731295       3.0608       0.0001       0.1150       1.0251\n",
      "Sella   2 12:27:49     -524.991352       1.9133       0.0004       0.1150       1.1487\n",
      "Sella   3 12:27:51     -525.156352       1.3285       0.0004       0.1322       1.0202\n",
      "Sella   4 12:27:52     -525.292944       0.8596       0.0004       0.1521       1.0184\n",
      "Sella   5 12:27:54     -525.418772       0.5524       0.0004       0.1749       1.0167\n",
      "Sella   6 12:27:55     -525.539468       0.3652       0.0002       0.2011       1.0246\n",
      "Sella   7 12:27:57     -525.650621       0.2198       0.0002       0.2011       1.0389\n",
      "Sella   8 12:27:58     -525.725450       0.1934       0.0007       0.2011       1.0558\n",
      "Sella   9 12:27:59     -525.757074       0.1514       0.0007       0.2011       1.0883\n",
      "Sella  10 12:28:01     -525.766113       0.1127       0.0003       0.2011       1.1582\n",
      "Sella  11 12:28:02     -525.768822       0.0860       0.0001       0.2011       1.2713\n",
      "Sella  12 12:28:03     -525.770370       0.0708       0.0001       0.2011       1.4722\n",
      "Sella  13 12:28:04     -525.771853       0.0705       0.0002       0.2011       1.4682\n",
      "Sella  14 12:28:05     -525.773234       0.0850       0.0002       0.2011       1.3123\n",
      "Sella  15 12:28:07     -525.774353       0.0759       0.0001       0.2011       1.0737\n",
      "Sella  16 12:28:08     -525.775070       0.0609       0.0000       0.2011       1.0103\n",
      "Sella  17 12:28:09     -525.775312       0.0364       0.0000       0.2011       1.2595\n",
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:28:10     -525.770433       0.4299       0.0000       0.1000       1.0000\n",
      "Sella   1 12:28:11     -525.775081       0.0860       0.0000       0.1000       1.0950\n",
      "Sella   2 12:28:12     -525.775314       0.0550       0.0000       0.1000       1.4639\n",
      "Sella   3 12:28:13     -525.775479       0.0189       0.0000       0.1000       1.1701\n",
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:28:13     -525.770883       0.4151       0.0000       0.1000       1.0000\n",
      "Sella   1 12:28:14     -525.775116       0.0668       0.0000       0.1000       1.0886\n",
      "Sella   2 12:28:15     -525.775274       0.0414       0.0000       0.1000       1.4286\n",
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:28:16     -525.770142       0.4151       0.0000       0.1000       1.0000\n",
      "Sella   1 12:28:17     -525.774659       0.0832       0.0000       0.1000       1.1130\n",
      "Sella   2 12:28:18     -525.774927       0.0534       0.0000       0.1000       1.4481\n",
      "Sella   3 12:28:19     -525.775107       0.0225       0.0000       0.1000       1.2192\n",
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:28:20     -525.770320       0.4049       0.0000       0.1000       1.0000\n",
      "Sella   1 12:28:21     -525.774471       0.0629       0.0000       0.1000       1.0986\n",
      "Sella   2 12:28:22     -525.774641       0.0382       0.0000       0.1000       1.4297\n",
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:28:23     -525.769405       0.4050       0.0000       0.1000       1.0000\n",
      "Sella   1 12:28:24     -525.773827       0.0766       0.0000       0.1000       1.1225\n",
      "Sella   2 12:28:25     -525.774104       0.0484       0.0000       0.1000       1.4442\n",
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:28:25     -525.768883       0.3942       0.0000       0.1000       1.0000\n",
      "Sella   1 12:28:26     -525.773234       0.0794       0.0000       0.1000       1.1359\n",
      "Sella   2 12:28:27     -525.773563       0.0514       0.0000       0.1000       1.4509\n",
      "Sella   3 12:28:28     -525.773836       0.0362       0.0000       0.1000       1.3300\n",
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:28:29     -525.769249       0.3797       0.0000       0.1000       1.0000\n",
      "Sella   1 12:28:30     -525.773100       0.0553       0.0000       0.1000       1.1169\n",
      "Sella   2 12:28:31     -525.773298       0.0330       0.0000       0.1000       1.4457\n",
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:28:32     -525.768346       0.3827       0.0000       0.1000       1.0000\n",
      "Sella   1 12:28:33     -525.772519       0.0664       0.0000       0.1000       1.1420\n",
      "Sella   2 12:28:34     -525.772823       0.0412       0.0000       0.1000       1.4398\n",
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:28:35     -525.768109       0.3688       0.0000       0.1000       1.0000\n",
      "Sella   1 12:28:36     -525.772154       0.0678       0.0000       0.1000       1.1577\n",
      "Sella   2 12:28:37     -525.772515       0.0440       0.0000       0.1000       1.4426\n",
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:28:38     -525.768148       0.3554       0.0000       0.1000       1.0000\n",
      "Sella   1 12:28:39     -525.772038       0.0660       0.0000       0.1000       1.1705\n",
      "Sella   2 12:28:40     -525.772435       0.0441       0.0000       0.1000       1.4424\n",
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:28:41     -525.768520       0.3417       0.0000       0.1000       1.0000\n",
      "Sella   1 12:28:42     -525.772236       0.0630       0.0000       0.1000       1.1829\n",
      "Sella   2 12:28:43     -525.772661       0.0465       0.0000       0.1000       1.4396\n",
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:28:44     -525.769312       0.3270       0.0000       0.1000       1.0000\n",
      "Sella   1 12:28:45     -525.772833       0.0623       0.0000       0.1000       1.1958\n",
      "Sella   2 12:28:46     -525.773282       0.0489       0.0000       0.1000       1.4347\n",
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:28:47     -525.770627       0.3110       0.0000       0.1000       1.0000\n",
      "Sella   1 12:28:48     -525.773930       0.0650       0.0000       0.1000       1.2097\n",
      "Sella   2 12:28:49     -525.774400       0.0512       0.0000       0.1000       1.4279\n",
      "Sella   3 12:28:50     -525.774902       0.0681       0.0000       0.1000       1.5606\n",
      "Sella   4 12:28:51     -525.775453       0.0631       0.0000       0.1000       1.4696\n",
      "Sella   5 12:28:52     -525.775915       0.0529       0.0000       0.1000       1.4260\n",
      "Sella   6 12:28:53     -525.776281       0.0474       0.0000       0.1000       1.3817\n",
      "Sella   7 12:28:54     -525.776525       0.0268       0.0000       0.1000       1.2128\n",
      "Sella   8 12:28:55     -525.776617       0.0108       0.0000       0.1000       1.1951\n",
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:28:56     -525.774950       0.3144       0.0000       0.1000       1.0000\n",
      "Sella   1 12:28:57     -525.778247       0.0483       0.0000       0.1000       1.1672\n",
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:28:58     -525.776503       0.3314       0.0000       0.1000       1.0000\n",
      "Sella   1 12:28:59     -525.780485       0.0572       0.0000       0.1000       1.1976\n",
      "Sella   2 12:29:00     -525.780743       0.0265       0.0000       0.1000       1.1270\n"
     ]
    }
   ],
   "source": [
    "from tblite.ase import TBLite\n",
    "\n",
    "calc = TBLite(method=\"GFN1-xTB\", verbosity=0)\n",
    "\n",
    "# 0. Define the reaction by SMILES\n",
    "rsmi, psmi = (\"CCCO[O]\", \"[CH2]CCOO\")\n",
    "# rsmi, psmi = (\"CCC[CH2]\", \"CC[CH]C\")\n",
    "rsmis = automol.smiles.split(rsmi)\n",
    "psmis = automol.smiles.split(psmi)\n",
    "\n",
    "# 1. Find minimum energy geometries for the reactants\n",
    "ratms_objs = list(map(automol.geom.ase_atoms, map(automol.smiles.geometry, rsmis)))\n",
    "ratms_objs = [sella_optimize_minimum(o, calc=calc) for o in ratms_objs]\n",
    "rgeos = list(map(automol.geom.from_ase_atoms, ratms_objs))\n",
    "\n",
    "# 2. Map the reaction\n",
    "rxn, *_ = automol.reac.from_smiles(rsmis, psmis, stereo=False)\n",
    "\n",
    "# 3. Generate a TS guess structure\n",
    "rxn = automol.reac.with_structures(rxn, \"geom\", rct_strucs=rgeos)\n",
    "ts_gra = automol.reac.ts_graph(rxn)\n",
    "geo0 = automol.reac.ts_structure(rxn)\n",
    "\n",
    "scan_coos = automol.reac.scan_coordinates(rxn)\n",
    "scan_vals = automol.reac.scan_values(rxn)\n",
    "scan_grids = numpy.meshgrid(*scan_vals)\n",
    "const_coos = automol.reac.constraint_coordinates(rxn)\n",
    "\n",
    "# 4. Do a relaxed scan over the reaction coordinate\n",
    "geo_grid, ene_grid = scan_relaxed(\n",
    "    geo=geo0,\n",
    "    calc=calc,\n",
    "    scan_coos=scan_coos,\n",
    "    scan_grids=scan_grids,\n",
    "    const_coos=const_coos,\n",
    "    gra=ts_gra,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     Step     Time          Energy         fmax         cmax       rtrust          rho\n",
      "Sella   0 12:29:01     -525.772435       3.2565       0.0000       0.1000       1.0000\n",
      "Sella   1 12:29:07     -525.865086       0.6032       0.0000       0.1000       0.8693\n",
      "Sella   2 12:29:08     -525.853138       0.8809       0.0000       0.0650      37.8141\n",
      "Sella   3 12:29:10     -525.878182       0.2114       0.0000       0.0650       0.4964\n",
      "Sella   4 12:29:11     -525.886250       0.1263       0.0000       0.0650       0.6920\n",
      "Sella   5 12:29:12     -525.888936       0.1757       0.0000       0.0650       0.8036\n",
      "Sella   6 12:29:14     -525.886290       0.2807       0.0000       0.0194       5.0839\n",
      "Sella   7 12:29:15     -525.888450       0.1432       0.0000       0.0126       5.2198\n",
      "Sella   8 12:29:16     -525.887876       0.1157       0.0000       0.0126       4.7689\n",
      "Sella   9 12:29:17     -525.888246       0.0577       0.0000       0.0126       0.5997\n",
      "Sella  10 12:29:18     -525.888921       0.0401       0.0000       0.0126       0.9239\n"
     ]
    }
   ],
   "source": [
    "saddle_idxs = scan_saddle_point_index(scan_grids, ene_grid)\n",
    "ts_geo0 = geo_grid[saddle_idxs]\n",
    "\n",
    "ts_atms_obj0 = automol.geom.ase_atoms(ts_geo0)\n",
    "ts_atms_obj = sella_optimize_ts(ts_atms_obj0, calc)\n",
    "ts_geo = automol.geom.from_ase_atoms(ts_atms_obj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10,)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAAGsCAYAAADaNnNsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAB2vklEQVR4nO3dd3hUVcIG8PdOn0lmJj2ZFBJqQkKRIlXF/VSwYF1F1hXUVVbXtsja2GLZXUXwW/uHigqsyroWRHEtCyqiSMcEJUDoKaS3mUkmmXq+PyYZCSlMYJKZSd7f89wn5M65d87hQObNueeeKwkhBIiIiIj6GVmwK0BEREQUDAxBRERE1C8xBBEREVG/xBBERERE/RJDEBEREfVLDEFERETULzEEERERUb/EEERERET9EkMQERER9UsMQURERNQvMQQFWEZGBiRJarM9/PDDnZZ3Op146KGHMHLkSERERCA5ORlz585FaWlpm3Lnn39+u/POnj3b9/o333zT7vXWbceOHQCAlStXdlqmsrLSr/b5W18iIqJQJ/HZYYGVkZGBW2+9FfPmzfPti4yMRGRkZIflzWYzrr32WsybNw+jR49GXV0d5s+fD5fLhZ07d/rKnX/++Rg2bBj++te/+vZptVoYjUYAgMPhQG1tbZtz/+Uvf8GXX36JI0eOQJIkNDU1wWw2tylz8803o7m5Gd98841f7fO3vkRERKFOEewK9EV6vR5JSUl+lTUajVi/fn2bfS+++CImTJiAoqIiDBgwwLdfp9N1el6VStXmNafTibVr1+Luu++GJEkAvKFJq9X6ylRVVeHrr7/GG2+80eZcmzdvxsMPP4wdO3YgLi4OV199NRYtWoSIiIhu1ZeIiCiU8XJYD1i8eDFiY2Nx1lln4YknnoDD4ejW8WazGZIkISoqqs3+VatWIS4uDjk5Obj//vthtVo7PcfatWtRXV2Nm2++udMyb775JnQ6Ha699lrfvp9++gkzZszANddcgx9//BHvvvsuNm3ahLvvvrvb9SUiIgplvBwWYM8++yzGjh2L6OhobN++HQsXLsSVV16J119/3a/jm5ubcc455yArKwtvv/22b/9rr72GgQMHIikpCXv27MHChQsxZMiQdqMyrS699FIAwGeffdbpe+Xk5GDatGlYunSpb9/cuXOh1Wrx6quv+vZt2rQJ06ZNQ2NjIzQajV/1JSIiCnmCTunRRx8VALrcduzY0eGxH3zwgQAgqqurT/k+DodDXHnllWLMmDHCbDZ3WXbnzp0CgNi1a1e714qLi4VMJhMffPBBp8dv3rxZABA7d+5ssz87O1uoVCoRERHh23Q6nQAg9u7de9r1JSIiCjWcE+SHu+++u82dWB3JyMjocP+kSZMAAIcOHUJsbGynxzudTsyaNQtHjx7F119/DYPB0OX7jR07FkqlEgcPHsTYsWPbvLZixQrExsbiiiuu6PT4119/HWeddRbGjRvXZr/H48Htt9+Oe++9t90xJ8736W59iYiIQg1DkB/i4uIQFxd3Wsfm5uYCAEwmU6dlWgPFwYMHsWHDhi7DUqv8/Hw4nc525xVCYMWKFZg7dy6USmWHxzY0NOC9997DokWL2r02duxY5OfnY8iQIQGtLxERUajhxOgA2rJlC5599lnk5eXh6NGjeO+993D77bfjiiuuaDOKkpWVhTVr1gAAXC4Xrr32WuzcuROrVq2C2+1GeXk5ysvLfROqDx8+jL/+9a/YuXMnjh07hs8++wzXXXcdxowZg6lTp7apw9dff42jR4/i1ltv7bSe7777LlwuF37961+3e+2hhx7Cli1bcNdddyEvLw8HDx7E2rVrcc899/hdXyIiorAQ7OtxfcmuXbvExIkThdFoFBqNRmRmZopHH31UNDY2tikHQKxYsUIIIcTRo0c7nWe0YcMGIYQQRUVF4rzzzhMxMTFCpVKJwYMHi3vvvVfU1NS0q8OvfvUrMWXKlC7rOXnyZHHDDTd0+vr27dvFRRddJCIjI0VERIQYNWqUeOKJJ/yuLxERUTjg3WFERETUL/FyGBEREfVLDEFERETUL/HusE54PB6UlpZCr9f7HjtBREREoU0IAavViuTkZMhkXY/1MAR1orS0FGlpacGuBhEREZ2G4uJipKamdlmGIagTer0egPcvkQsBEhERhQeLxYK0tDTf53hXGII60XoJzGAwMAQRERGFGX+msjAEEYULIQBHI9BcDzTVw9NkRmNjA6wNVlgbG2BrtMHtskMpk6CQCShlgFICFHIJCpkEhVwOhVwGuUIOhUIJpUIJhVIJSaYAFGpAoW35qgGUGkAVAagivV+VOoBz44ioj2EIIgo2lwOwlgHWMnjMpWisKYK9thRuSznQWAVFUxXU9hponWbI4fYdJgOgb9l6mgcSmmU6NMv1aFYYYFfoYVcaYVfHwKGOhUcXh4iYJBjjUhBnGgBtdLI3PBERhTCGIKKe5PHAYS6HpfIoGiuL4Kgtgsd8HArrcahtZdDby6F31UIG75ql/gQbp5DDjAhYhA7NUMMpU8EjbxnBkSngEjK4BOD2SHALwC0EhEd4vwoP4PFABg8UcEMOD5RwQym5oIYDajihhhNayQEdmhEpNbfUS0DnaYTO0wg4y/1qug1aWJWxaFLHwx2RCLkxGZqYFBjiB0AblwbJmAroTYCcP4aIKDj404foTDibgPoioO4Y6koPoqr4EGzVRVDbyhHlqkKcqIEKbsQB6OoRvHahQKWIRhliUC5iYFHEoUkdB5c2FohIgNyQCLUhHpFRcYgxRiFOr0GcXoXUCDVUiu4t9yWEgNMtYHe5YXd5vJvT+2eby4Palj/bnW7YnS647DZ47FaIJgtkzfWQ2eshd5ihtNdDba+Fxlnr/Wqvgd5VhzjUtYSoJuicJYCzBGgAUNG+Lh7I0KCMgU2bArcxDfLodEQkDkJk0mBIMYMAYyogk3erfURE/mIIIvJHsxmo2AtU7QOqCiAq98FdWQBF48+jItEt28ncQkIVolEpi0e9MgEN6kTYdSZ49CmQRaVCHZsGfawJcXotMiJVGKtTQSHvuXVMJUmCSiFBpZD1yKU0s82Bwspq1FQUwVJZjKa6UnjMZZDbyqFtrkS0uwbJqEGiVAuV5IbBWQ2Dsxqw7AaK257LBQVq1cmwRQyAK3oQ5AlZiEzNQUz6SMgjY3ug9kTUn/DZYZ2wWCwwGo0wm828O6y/cbuA8t3A8R+A47u8W/WBTotbhRbFIgEliEdzRCp0cQNgTEqHJiYNkfHpiEpIgzFSB5mME4sBwOZwobS+CSW1jaipPI7GymNw1RZBYSmG1laKGGc50qUKpEmVUEuuTs9TCyNK1QPRaBwGhWkEYgedhZRhY6HU9sYsKSIKVd35/GYI6gRDUD8iBFBVABz5Bji6ETi2CbBb2hU7LmJx0JOKAyIVB0UKjkmpiDBlIntQOiYOjsO49GhEqjm4eqYcLg/Kzc0oqbGitvwYmisOQqo9Aq3lMOKajiHFXYIUqbrDYz1CQpnchJrIoXDHZSMi/SwkZ45HZOJg3t1G1E8wBAUAQ1Af53EDRVuAfZ94N8vxNi+bEYEf3EOwWwxGnmcwfvQMRoMiCmPSojBxUCwmDYzBmAHR0Ko4X6W3uT0CVTU1qC38CQ1FP8Jdno+I+gNIth9BnGTu8JgG6FCmGQxb9HAoU0cjfuh4xA08C5JS08u1J6KexhAUAAxBfZAQQPE2YPe/Ifb/B1Jjle+lJqiww52J7z0j8L0nB3tFBlRKBcalR2PiwFhMHBiD0WlR0CgZekKVxyNQerwIpQU7YSvZDWXVPsTbDiLDUwSV5G5X3gk5ypQDUG/IAkyjED14PEzDzoYioqOZXUQULhiCAoAhqA9prIE77x04tq+A1nzIt7teRGC9exw+90zA954RUKi0GJ8Rg4mDYjBxYCxGphi7fecVhZ5aSwMKC/JQf+QHoOInRJn3Y6DrMKKkxg7Ll8mSUBWZBVfiKERmjEVy9mRERif1cq2J6HQxBAUAQ1D4O3r0IOyfP4LBleuhhBMAYBNqfOqeiLWeKdijGoWxGQm+0JOTbOjRu7IodDQ7XDh25CCqDm2Hq2Q3Iur2IrX5IJI7mWtUKcWiTJeJ5riR0AwYi5ThkxCXnNG7lSYivzAEBQBDUBgTAoVfvYro7x6HQbIBAH7yZOBj+UWoTL8co4cMwMSBMRhuMkDOO7aohccjUHK8BGUF29FclAtNzR4kNh7AAFEKmdT+x2Q1onBcNxz2hLMQOWg8Bow8F5HRiUGoORGdiCEoABiCwlTtUTg/ugfKou8AAAcVQ3Fg3KMYMuZ8DE2I5G3q1G01tTU4vm87Go7thLLiR8Q37EeauxjyDoJRqSwJFfoRcCePQ0zmFAwYPhEKtTYItSbqvxiCAoAhKAztfhfik99DcjWhSaiwQn0Dbpz/FAw6fghRYDVazTiWvxXmw9uhKM9DUsNeDBCl7co5hALHVENQFz0a8vSJSB5xLkxpQyDJeNmVqKcwBAUAQ1CYOfQlsGoWINzY4s7GX8Rv8fyd1yAn2RjsmlE/UVVZgaI938N2dCsiq/KQ3rQXMZK1XblKxKBYlw1b4jjoB0/GwJFTYTTyZwxRoDAEBQBDUBip2Au8MR1wWLHGfS4WOG/Hk9eMxq8mDAh2zagf87g9KD6Sj4q9m+Ap3o64ut3IcB2FQvK0KecUchyVZ6DSOBIieRzisqZiUNYoqJXKINWcKLwxBAUAQ1CYsFYAr18AmIvxg5SN2U0P4bIxGXhm1mhIXCGYQkyzzYrCPZthOfg9VGU/ILVxD2JFXbtyZhGBY6qhsMSMgDJtPEzZU5CWPhQy3r1IdEoMQQHAEBQGHDZg5WVA6Q8oU6Ti4oZHEJ+QhI/vmooIPr6CwoEQMJcdQUn+d2g+ug36mt0YYD8ATcuSDicyIwLHVYNREzkUtqhMiPhM6OIHIiZpABKNOsRGqDjxnwjd+/zmJwWFJ48HWPNboPQHNCmMmN24AA6lES//eiwDEIUPSYIxeTCMyYMB3AwAEC4Hyg7lomL/FrhLfkB0/R6kOY/BKDXC6PgRqP0RqAVwxHsKh5CjTMTiEOJgVcTCoY6GRxcLWUQ81IZ46IwxMETFIyY2HrFxidBERgOcmE0EgCGIwtW2V4B9n8AjU2Gu7fcoFEl45uoRGJrIJ4hTeJMUKpiyJsKUNdG3z2lvwqGCXNQeyYWyei/09fsR1VyMaFcVVJIb6VIl0lEJeAA0tWw1HZ/fAwkNUiRscj0cyii4NNHwRCZBbjBBE5MKfXwaIhIHQYodDCjUvdFkoqBhCKLwU3MY+OqvAIDFuAk7PFmYfXYarhmbGuSKEfUMpVqLIaOmAKOmtH3B7QKspXDVFsJacQyNdeWwmyvgslZDslVDYa+D0mmBxt2ASE8DtJIDMggYhBUGlxVwlXoDU/tpSXBDhkp5Emo0A9BgGAxV6lkwZU1C0sAcSDI+Q4/6Bs4J6gTnBIUojwdYeSlQtAU/qsbgCsv9yEoy4KO7pvLhpkRdEELAbG1AVXUF6qorYK2rhK2+Gk5LJaSGcqiaKhHpqEaMpxrpUgUMUlOH52mAFsfVQ2GLGwlNxgSk5JwDg2kwwBsRKERwYnQAMASFqK2vAF88BIdMi/9pegr1KhM+ueccDIyLCHbNiPoEu8uNSnMzaiqKYSvdB09lARQ1+xFl3osM5xFopPaTtuthQGlENuxJY6EfOgVpI86BOjI6CLUnYggKCIagEFR7BFg6BXA14c/OW/C2+yL83w1jcdkoU7BrRtQvNNvtOLz3B9Qc3AapNBfxlj0Y5D4KleRuU84jJJQo0lAdNQpS6njEZ5+HlKFn8TIa9QreHUZ9j8cDfHwP4GrCduRglfsC3DwlgwGIqBdp1GrkjJkMjJns21drtuBI/jZYD26BpiIXqbZ8pEkVGOAuwoCaIqDmP8BuwAotCjXD0RA3BpqBEzFg1HmIief/XwoujgR1giNBIWb7a8Bn96MZalxkfwoxKcPw3h2ToVbwN0uiUCKEQElJEUr3fAtH4XZE1eRhkKMAEZK9XdkiKRllkTlwmcYhetgkDMyZBK2Wz/qjM8PLYQHAEBRCLGXAS+MBRwMecd6Ej5SX4dN7z0VajC7YNSMiPzgcDhTu24nagk2QHd+FJOtPSPMcb1fOLhQ4qhiMuqgRkKWORULmJKQNHQ0FHyFC3cDLYdS3rPsz4GhArmcI3nJfhGW/PosBiCiMqFQqDB09BRj98y3+1roKFP/4LRqObIe2Mg9pTXsRJTUgy10A1BQANauB3UCjUKNAOQQ1xhyI5LEwDpmIgUNHwKhTBbFF1Ff06EhQXV0d7r33XqxduxYAcMUVV+DFF19EVFRUp8cIIfD4449j2bJlqKurw8SJE/F///d/yMnJ8ZWx2+24//778c4776CpqQkXXHABli5ditTUn9eJycjIQGFhYZtzP/TQQ3jqqaf8qjtHgkLEkY3Am1fALSRc4fg7zjn3Aiy8dHiwa0VEASY8HlQWFaAs/zvYi3+AvnYPMuwHoOvgMlqdiESBfCiqDDlwJY2BcchEDB44CGnROj46hELnctgll1yCkpISLFu2DADw29/+FhkZGfjkk086PWbx4sV44oknsHLlSgwbNgx///vf8e2336KgoAB6vXc14N/97nf45JNPsHLlSsTGxuIPf/gDamtrsWvXLsjl3jkiGRkZuPXWWzFv3jzfuSMjIxEZGelX3RmCQoDLAc8rUyGrPoB/ui7CJykL8M5vJ0HJh0gS9QselwtlR/ag+sAWeEp2wVj3E1Lth6CCq13ZEhGHvRiMSn02nImjEDlwAoYMSEFmkh46FS969CchEYL27duH7OxsbN26FRMnepd/37p1KyZPnoz9+/cjMzOz3TFCCCQnJ2P+/Pl46KGHAHhHfRITE7F48WLcfvvtMJvNiI+Px1tvvYXrr78eAFBaWoq0tDR89tlnmDFjBgBvCJo/fz7mz59/WvVnCAoBm54DvnwU1cKAX8pfwL9/fzFMRk6aJOrXXA40FOah6sAWuIp3Ql/zIxLshZCh/UfZMU8i9ol0lGqHwh6XA23aWUjPGIzhyUYkGTSQuMBjnxQSc4K2bNkCo9HoC0AAMGnSJBiNRmzevLnDEHT06FGUl5dj+vTpvn1qtRrTpk3D5s2bcfvtt2PXrl1wOp1tyiQnJ2PEiBHYvHmzLwQB3lGlv/3tb0hLS8N1112HBx54ACpVx9eR7XY77Pafh10tFssZtZ/OkLkErg1PQQFgkesG/O3X5zAAERGgUCFy8AREDp7w875mC9zH81B7aBscRTuhq/4J0fbjyJBVIAMVgGM7UAqgFDBv1eGASMVmWRoaDEOhSMhEbFomBgwahsGmGN5x2s/0WAgqLy9HQkJCu/0JCQkoLy/v9BgASExMbLM/MTHRN7+nvLwcKpUK0dHR7cqceN7f//73GDt2LKKjo7F9+3YsXLgQR48exeuvv97hey9atAiPP/64/w2kHtWw9kFEupuwwzMMKefdjPOGxQe7SkQUqjQGyAefh/jB5/28z1YLlP+IhsJcNBb+AGVVPqIaj8Io2XC2dABn4wBg+QqwADgEuL+WUIZYVCuT4dAlwaVLgFOXCBGZBLkhESp9PDRR8YgwxsMQoYFRq2Rg6gO6HYIee+yxU4aFHTt2AECHQ41CiFMOQZ78uj/HnFzmvvvu8/151KhRiI6OxrXXXovFixcjNja23fELFy7EggULfN9bLBakpaV1+Z7UM+z71yPy8KdwCRk+SLwPT16UFewqEVG40cUAg85H5KDz4ZsJ6mwGag7BUb4X9YU/wVGWD2X9URjtx6GBHamoRqqrGrD86A1HHfAICWZEoExEoEGKgE0WAbtcD4dSD7cyEkKlBzR6yDV6KLRGKHUGqCOM0EQaodNHI9IQDUNULDRqdW/9TVAXuh2C7r77bsyePbvLMhkZGfjxxx9RUVHR7rWqqqp2Iz2tkpKSAHhHe0ymn1cSrays9B2TlJQEh8OBurq6NqNBlZWVmDLlpCcsn2DSpEkAgEOHDnUYgtRqNdT8Rxl8ziZYP/w91ADel1+CP8y9BnLe7UFEgaDUAEkjoEoagYSzZv28XwiIhgpUFxegorAADVXFkNsqoG6qgs5ehUhXDSLcFuhFA2SSQDQaEC01tBwLwNWydfzM2Q7ZhBo2SYtGWSSa5AbYFXo4VEa41VGALhbyiDgo9LHQGhMQEZ0IQ2wS9NGJkCnCe80kIQSanG6YrRY0VRyGQq7AgMwxQatPt0NQXFwc4uLiTllu8uTJMJvN2L59OyZM8F673bZtG8xmc6dhZeDAgUhKSsL69esxZoz3L8XhcGDjxo1YvHgxAGDcuHFQKpVYv349Zs3y/iMuKyvDnj17sGTJkk7rk5ubCwBtwhWFnvz3HkOO4zjKRAwGzVqEBL0m2FUior5OkiDpkxCfnYT47Gmdl3O7gKY6eBqqYLPUwmaphb2hFo6GOria6uFuskA0WwGHFTKHBXJnI5SuRqjdNqiFDTphgxYOAIBOskMHO+I89YAHgBOnDFEeIaFeioRZFoVGRRSaVTFwaGLh0cZBikyAwpAAjTER2hgTjLEmREXHQtkDl+yEELA53DA3OdttltavtmZI1jKoGkqgs5XCYC9DnLMMJncZ0qQKmKQ6AMAu3TkY8OCnAa+jv3psTtDw4cNx8cUXY968eXj11VcBeG+RnzlzZptJ0VlZWVi0aBGuvvpqSJKE+fPn48knn8TQoUMxdOhQPPnkk9DpdLjhhhsAAEajEbfeeiv+8Ic/IDY2FjExMbj//vsxcuRIXHjhhQC8k7K3bt2KX/ziFzAajdixYwfuu+8+XHHFFRgwYEBPNZnO0OG9P2DogdcACcjLeRiXZKUHu0pERD+TK4DIeMgi4xGZBPi34EpbHqcDjQ11sNbXwWatRbOlFs7GGjgb6uCx1ULYaiFrqoPSUQeNsx4RrnoYPBYY4R2FioIVUR4r4CgGHAAaOn8vu1CgQjLAIjOiUe4NTW51FDzaGMgiYqGIjIXGEAtNpBE2KQL1Hi1qnSqY7R6Ym1yob3ajvtmDRlsT7E2NcDQ1wNFsg8fegEhhRTQaECU1IBpWxElmJEl1GCbVIVGqQzzqoZA87St1wgonVugAeXAXvezRxRNWrVqFe++913cn1xVXXIGXXnqpTZmCggKYzWbf9w8++CCamppw5513+hZLXLdunW+NIAB49tlnoVAoMGvWLN9iiStXrvStEaRWq/Huu+/i8ccfh91uR3p6OubNm4cHH3ywJ5tLZ6Ch2QnL6nugktzI00zAjF/OO/VBRERhRqZUQR+dCH10x9NCOmN32GGpqYClphy2ujLYzZVwWSuBxirIbdVQN1dD46xDpKseBmFGJJqgllxIRC0SPbU/jzY1BqARfl6Rc0sK2LQm2CNT4dGnAtEDoIgbAk3iYGgThkCvi8G4IC9TwGeHdYLrBPUeIQTeenUJ5pY/iWao0DxvM6JShga7WkREYcttt8FaWw5rTTls9RWwmyvgtFZBNNYATXVQ2GuhdtRD5WqAxtMIPWyIEI2Qo4PRmxYemQpCoQWUWki6aEi6GEi6GEAbA0TEAXoTYEgG9EneP0cmArLev4MuJNYJIvLXe9/+iEvLXgIkoHrcfKQyABERnRG5Woco0yBEmQb5f5AQgMsOCE/L5gY8bu8lK6UWsiAEmp7GEERB9WNJPaSvHkOczILaiMFIvZSXLImIgkKSvHfQ9SN8CBMFjdnmxNtvvopZsq8BANHXvQjIw/v2TyIiCh8MQRQUQgj87d9f4yH7iwAA+/g7IGVMDXKtiIioP2EIoqB447vDuPLoXxErWdEUmwP1xX8NdpWIiKifYQiiXrersBY16/4X58r3wCXXQDt7JaDgat1ERNS7ODGaelVtowMvvf0+lsnfAwDIL10CxA8Lcq2IiKg/4kgQ9RqPR2DhO5vxiP0fUEpuuLKugDR2brCrRURE/RRDEPWa177eg5sKF2KgrALOyGQornzBe0smERFREPByGPWK7fuOYey3t+JseQGciggoZ78FaKODXS0iIurHGIKox1VVVUD73rUYKTuIJlkktDd/DKSOD3a1iIion+PlMOpRlpoKWF+9FCPFQVgkPaSb1jIAERFRSGAIoh6z6ftvUfXiBRjkOoQaYUD9dR9Ckz4u2NUiIiICwMth1AMqzU34+q2/46qqV6GRnKiRolF+9XvIyZ4Q7KoRERH5MARRwAghsHbTD4j9agFmIw+QgCNRk5F803LkRCcHu3pERERtMARRQBRWN2DtqpdwQ+1LiJWssEOF2il/waCL7uFt8EREFJIYguiMuNwefLjuK6RtfRT3SPmABFRHZiLqxpUwJWUHu3pERESdYgii07avsAx7//0nXG37CErJDQdUsE38PeIueoDPAiMiopDHEETd1uxw4ut3X8DYQy/hl1ItIAHHE89H8vXPISpmYLCrR0RE5BeGIOqWgs1rIfvyEVzqOQpIQJXCBOXMJUg564pgV42IiKhbGILIL43FP+H4e/cj07oVAGCFDsdH3IWsK+8HlJog146IiKj7GIKoaw2VOL7mz0g6/D6GwQOnkGNb3FUY9asnkRWXFOzaERERnTaGIOqYwwbbxuch2/w8UkQTAGCjfBIiL/s7zhl7dpArR0REdOYYgqgtISD2foymTx6CrrkcAJDnGYzd2Q9g1jWzoFXJg1xBIiKiwGAIop/VHUPTxwugPfYVdABKRBzeirgFl86+EzcNiAl27YiIiAKKIYgAlwOe71+AZ+MSaD122IUCr3muhPy8+3D//+RAKedzdomIqO9hCOrvqgpg//dNUNfsgwzA9+4cvJswH/defymGJEQGu3ZEREQ9hiGovxICrp3/hPj8Qag9dlQLA5aIuRh5yW14blIGZDI+74uIiPo2hqD+qKkede/dieijnwIAvnWPxOr0P+PBa6chJUob5MoRERH1Doagfqap6AfY3/4Voh3lcAo5/k92AwZe8xCeOysVEp/2TkRE/QhDUD+Sv+E9DNx4D6LQjGOeRKwe+DhumfVLxESogl01IiKiXscQ1E/sfG8RxuQvhlwS2C6Ngv26FfjDyCHBrhYREVHQMAT1dR439v3zHowvXAVIwPboy5A97w1E6jj3h4iI+jeGoL7MYUPR6zdgeOUGAMBXKb/D/9z6JCQZ1/0hIiJiCOqrbLWofPVKDDD/CLtQ4rMhj+KqG+/m5GciIqIWDEF9UX0R6pddgQTbUZiFDp/kPItfX3c9AxAREdEJGIL6mop8NL5xJaIcVSgVMfjvWf+Hm6+6hAGIiIjoJAxBfcmx72F/axYi3A0o8KRi44RXMO+ycxiAiIiIOsAQ1FcUfA7XuzdB7bFjmycLuVOW4o6LxwW7VkRERCGLIagv2P1veNbcCQXcWO8eh0PnPY/fXTQy2LUiIiIKaQxBYU5sWQrpvwshA7DafS4qf/G/+N3/ZAW7WkRERCGPIShcCQHxzSJIGxcDAJa7Lobror/jd9OGBrliRERE4YEhKBwJAbHuz5C2vAQA+IfzWhhn/BG/PW9wkCtGREQUPhiCwo3HA/HZA5B2vg4AeMw5F+mXLsAtUwcGuWJEREThhSEonHjcEGvvhZT3NjxCwh9dtyJ75j2YOzkj2DUjIiIKOwxB4cLtgvjoDkg/vQ+3kHC/8w6Mu+J3uHFSerBrRkREFJYYgsKB2wXx4W2Q8tfAKeT4vfMunHvVPPxqwoBg14yIiChsMQSFOrcL4sPfQspfA4eQ427X73Hh1b/BrLPTgl0zIiKisMYQFMo8bog1d0DKXw2HkOMu13zMuOY3uHZcarBrRkREFPZkPXnyuro6zJkzB0ajEUajEXPmzEF9fX2Xxwgh8NhjjyE5ORlarRbnn38+8vPz25RZtmwZzj//fBgMBkiS1OE5T+e9Q4rHDfHR7yDteR9OIcc9rntx6bUMQERERIHSoyHohhtuQF5eHr744gt88cUXyMvLw5w5c7o8ZsmSJXjmmWfw0ksvYceOHUhKSsJFF10Eq9XqK2Oz2XDxxRfjj3/8Y0DfO2R4PPB8fBekH9+FS8hwr+seXHrdPFw9hgGIiIgoYEQP2bt3rwAgtm7d6tu3ZcsWAUDs37+/w2M8Ho9ISkoSTz31lG9fc3OzMBqN4pVXXmlXfsOGDQKAqKurO+P3PpnZbBYAhNls9qt8wHg8wv3JAiEeNQjnI1Hizj8+ItbmHe/dOhAREYWp7nx+99hI0JYtW2A0GjFx4kTfvkmTJsFoNGLz5s0dHnP06FGUl5dj+vTpvn1qtRrTpk3r9JhAvbfdbofFYmmzBYPny8ch2/k6PELC/e47cdnsO3D56OSg1IWIiKgv67EQVF5ejoSEhHb7ExISUF5e3ukxAJCYmNhmf2JiYqfHBOq9Fy1a5Js/ZDQakZbW+3dfub99BrLvnwUA/MV9Ky6efQ8uHWnq9XoQERH1B90OQY899hgkSepy27lzJwBAkqR2xwshOtx/opNf9+eYU53jVOdZuHAhzGazbysuLu7W+50p9/bXIf/6cQDAYvcNOP+GB3HxiKRerQMREVF/0u1b5O+++27Mnj27yzIZGRn48ccfUVFR0e61qqqqdiM9rZKSvB/65eXlMJl+HgGprKzs9JjOztPd91ar1VCr1X6/RyC58t6D7LP7AQAvu6/C+BsewwXD/W8vERERdV+3Q1BcXBzi4uJOWW7y5Mkwm83Yvn07JkyYAADYtm0bzGYzpkyZ0uExAwcORFJSEtavX48xY8YAABwOBzZu3IjFixf7XcfTee9gcRz8GrKPfgcZBN7yzEDWDUvwCwYgIiKiHtdjc4KGDx+Oiy++GPPmzcPWrVuxdetWzJs3DzNnzkRmZqavXFZWFtasWQPAewlr/vz5ePLJJ7FmzRrs2bMHN998M3Q6HW644QbfMeXl5cjLy8OhQ4cAAD/99BPy8vJQW1vbrfcONsfx3XD/69dQwIVPPZMx4NcvMgARERH1kh5dJ2jVqlUYOXIkpk+fjunTp2PUqFF466232pQpKCiA2Wz2ff/ggw9i/vz5uPPOOzF+/HgcP34c69atg16v95V55ZVXMGbMGMybNw8AcN5552HMmDFYu3Ztt947mJqrj6Fx+VXQChu2ebIRdcMbmJbJAERERNRbJCGECHYlQpHFYoHRaITZbIbBYAjouZst1ah54RdIcRXhgEhD3fVrMTF7UEDfg4iIqD/qzud3j44EUXtNtkYUvnQlUlxFKBcxaLj2HQYgIiKiIGAI6mW7Vy9GpmMPLEKHqitXYezIkcGuEhERUb/Ep8j3srN/9RdsebkYMeOvw8ixoXWnGhERUX/CENTL5AolJt+zItjVICIi6vd4OYyIiIj6JY4EdaL1prlgPUiViIiIuq/1c9ufm98ZgjphtVoBICgPUiUiIqIzY7VaYTQauyzDdYI64fF4UFpaCr1e3+2Ht4YSi8WCtLQ0FBcXB3y9o1DA9oW/vt5Gti/89fU29rX2CSFgtVqRnJwMmazrWT8cCeqETCZDampqsKsRMAaDoU/84+4M2xf++nob2b7w19fb2Jfad6oRoFacGE1ERET9EkMQERER9UsMQX2cWq3Go48+CrVaHeyq9Ai2L/z19TayfeGvr7exr7evK5wYTURERP0SR4KIiIioX2IIIiIion6JIYiIiIj6JYYgIiIi6pcYgoiIiKhfYggKI99++y0uv/xyJCcnQ5IkfPTRR12W/+abbyBJUrtt//79bcqtXr0a2dnZUKvVyM7Oxpo1a3qwFZ3rifatXLmywzLNzc093Jr2uts+ALDb7fjTn/6E9PR0qNVqDB48GMuXL29TJlT6D+iZNoZzH958880d1j0nJ6dNuVDpw55oXyj1H3B6/0ZXrVqF0aNHQ6fTwWQy4ZZbbkFNTU2bMuHah8Cp2xdqfRhIDEFhpLGxEaNHj8ZLL73UreMKCgpQVlbm24YOHep7bcuWLbj++usxZ84c7N69G3PmzMGsWbOwbdu2QFf/lHqifYB3KfgTXy8rK4NGowlk1f1yOu2bNWsWvvrqK7zxxhsoKCjAO++8g6ysLN/rodR/QM+0EQjfPnz++efb1Lm4uBgxMTG47rrrfGVCqQ97on1A6PQf0P02btq0CXPnzsWtt96K/Px8vP/++9ixYwduu+02X5lw7kN/2geEVh8GlKCwBECsWbOmyzIbNmwQAERdXV2nZWbNmiUuvvjiNvtmzJghZs+eHYBanr5AtW/FihXCaDQGtG6B4E/7Pv/8c2E0GkVNTU2nZUK1/4QIXBvDuQ9PtmbNGiFJkjh27JhvX6j2YaDaF6r9J4R/bXz66afFoEGD2ux74YUXRGpqqu/7cO5Df9oXyn14pjgS1A+MGTMGJpMJF1xwATZs2NDmtS1btmD69Olt9s2YMQObN2/uzSqeka7aBwANDQ1IT09HamoqZs6cidzc3CDUsvvWrl2L8ePHY8mSJUhJScGwYcNw//33o6mpyVcm3PvPnzYC4duHJ3vjjTdw4YUXIj093bcv3PvwRB21Dwjv/psyZQpKSkrw2WefQQiBiooKfPDBB7jssst8ZcK5D/1pHxDefdgVhqA+zGQyYdmyZVi9ejU+/PBDZGZm4oILLsC3337rK1NeXo7ExMQ2xyUmJqK8vLy3q9tt/rQvKysLK1euxNq1a/HOO+9Ao9Fg6tSpOHjwYBBr7p8jR45g06ZN2LNnD9asWYPnnnsOH3zwAe666y5fmXDuP8C/NoZzH56orKwMn3/+ebvLDOHeh606a1+499+UKVOwatUqXH/99VCpVEhKSkJUVBRefPFFX5lw7kN/2hfufdilYA9F0enBaQxVCyHEzJkzxeWXX+77XqlUin/9619tyrz99ttCrVafaRXPSKDadzK32y1Gjx4t7rnnnjOo3Znzp30XXXSR0Gg0or6+3rdv9erVQpIkYbPZhBCh239CBK6NJwunPjzRk08+KWJjY4Xdbm+zP1T7MFDtO1mo9J8Q/rUxPz9fmEwmsWTJErF7927xxRdfiJEjR4rf/OY3vjLh3If+tO9kodSHZ4ojQf3MpEmT2qT3pKSkdr+tVFZWtvutJlyc3L6TyWQynH322WHxG4zJZEJKSgqMRqNv3/DhwyGEQElJCYDw7z9/2niycOrDVkIILF++HHPmzIFKpWrzWrj3IdB1+04Wbv23aNEiTJ06FQ888ABGjRqFGTNmYOnSpVi+fDnKysoAhHcf+tO+k4VbH3aFIaifyc3Nhclk8n0/efJkrF+/vk2ZdevWYcqUKb1dtYA4uX0nE0IgLy+vyzKhYurUqSgtLUVDQ4Nv34EDByCTyZCamgog/PvPnzaeLJz6sNXGjRtx6NAh3Hrrre1eC/c+BLpu38nCrf9sNhtksrYflXK5HIC3LUB496E/7TtZuPVhl4I1BEXdZ7VaRW5ursjNzRUAxDPPPCNyc3NFYWGhEEKIhx9+WMyZM8dX/tlnnxVr1qwRBw4cEHv27BEPP/ywACBWr17tK/P9998LuVwunnrqKbFv3z7x1FNPCYVCIbZu3don2vfYY4+JL774Qhw+fFjk5uaKW265RSgUCrFt27aQb5/VahWpqani2muvFfn5+WLjxo1i6NCh4rbbbvOVCaX+a61zoNsYzn3Y6sYbbxQTJ07s8Jyh1Ic90b5Q6j8hut/GFStWCIVCIZYuXSoOHz4sNm3aJMaPHy8mTJjgKxPOfehP+0KtDwNJEazwFeo8Hg9KS0uh1+shSVKwqwMA+O677zBz5kzf9wsWLAAA/OpXv8Irr7yCwsJCFBUVwWKxAAAsFgsWLFiA0tJSaLVaZGVl4f3338eFF17oKzNixAgsX74cf/vb3/DnP/8ZAwcOxIoVKzB8+HBfmXBuX0VFBW677TZUVFTAYDBg9OjR+Pzzz5GVlRXy7QOANWvW4IEHHsC4ceMQExODq6++Gn/5y19Csv+AnmljuPeh2WzGBx98gMWLF3dY31Dqw55oXyj1H9D9Nl5zzTWorKzE888/jwULFsBoNGLatGl4/PHHQ/L/YU+0L9T68FSEELBarUhOTm43ynUySYhOxrv6uZKSEqSlpQW7GkRERHQaiouLO72s3oojQZ3Q6/UAvH+JBoMhyLUhIiIif1gsFqSlpfk+x7vCENSJ1ktgBoOBIYiIiCjM+DOVhXeHERERUb/EEEQh54eiOpy3ZAMeeH83ahsdwa4OERH1UQxBFFKanW7c/95uFNXa8P6uElzwj2/w/s7iTterICIiOl0MQRRSXvr6EI5UNyJer0ZWkh51Nice+OBH/Oq1rThc1XDqExAREfmJIYhCxr4yC17ZeBgA8LcrR+CTe87Bw5dkQaOUYeuRWlz+4iYU1diCXEsiIuorGIIoJLg9Ag+v/hEuj8DFOUm4eEQSlHIZ7pg2GOvvm4ZskwE2hxvv7yoOdlWJiKiPYAiikLBy8zHsLjFDr1Hg8Stz2ryWFqPD7dMGAQA+zivl/CAiIgoIhiAKuuJaG/73vwUAgD9eOhyJBk27MhdlJ0KnkqOo1oYfiup7uYZERNQXMQRR0D3y8R40Od2YODAG14/v+FElOpUC07MTAQAf5x3vzeoREVEfxRBEQVXdYMeGgioAwJPXjIRM1vkKn1eOSQEA/OfHMjjdnl6pHxER9V0MQRRU3x+qBgBkmwwYHB/ZZdlzh8QhNkKF2kYHNh2s7o3qERFRH8YQREHVGoLOGRp3yrIKuQwzR5kAAB/xkhgREZ0hhiAKGiGEb0TnnCGnDkHAz5fE1uVXoNHu6rG6ERFR38cQREFztLoRpeZmqOQynJ0R49cxY9KikB6rQ5PTjfV7K3q4hkRE1JcxBFHQtF4KG5ceDa1K7tcxkiThytHJAHhJjIiIzgxDEAXNpm7MBzpR6yWx7w5Wo7rBHvB6ERFR/8AQREHhcnuw+XANAP/nA7UaHB+JkSlGuD0Cn/5Y1hPVIyKifoAhiILip+NmWJtdMGgUGJFi7PbxV57lvST26U8MQUREdHoYgigoWucDTRkcB3kXCyR25vzMeADATyVmuD18lhgREXUfQxAFxenOB2o1MC4SOpUcTU43Dlc1BLJqRETUTzAEUa+zOVzYVVgHoPvzgVrJZRJykg0AvKNBRERE3cUQRL1u+9FaON0CKVFapMfqTvs8I1OiAHjnFxEREXUXQxD1utb5QOcOjYMkdX8+UKuRqd6RoD0MQUREdBoYgqjXfdfyqIypp3kprNXIlrvK8kstnBxNRETdxhBEvarKasf+cisAYMrg2DM6FydHExHRmWAIol61+bB3FCgn2YDYSPUZnYuTo4mI6EwwBFGv2nnMe1fY5EFnNgrUqnWhRU6OJiKi7mIICoLXvzuCCktzsKsRFD+2hJXRaVEBOd+oVG8I4uRoIiLqLoagXvZR7nH8/dN9uOT57/BNQWWwq9OrHC4P9pVZAPwcXs4UJ0cTEdHpYgjqZaNSjRhuMqC20YGbV+zA4i/2w+X2BLtaveJAhRUOlwcGjQIDYk5/faATnTg5+ggnRxMRUTeETQhaunQpBg4cCI1Gg3HjxuG7777rtGxZWRluuOEGZGZmQiaTYf78+b1X0VMYFB+JNXdOwY2TBgAAXv7mMGYv24rS+qYg16zntc7bGZUadUbrA53oxMnRP3JyNBERdUNYhKB3330X8+fPx5/+9Cfk5ubi3HPPxSWXXIKioqIOy9vtdsTHx+NPf/oTRo8e3cu1PTWNUo6/XzUSL90wBpFqBXYW1uHSF77DJ7tLIUTfvaTTGoJO56nxXeHkaCIiOh1hEYKeeeYZ3HrrrbjtttswfPhwPPfcc0hLS8PLL7/cYfmMjAw8//zzmDt3LozGwH7gBtLMUcn49N5zMDLFiHqbE/e8k4vb39qFyj46abr1NvZAzQdqxcnRRER0OkI+BDkcDuzatQvTp09vs3/69OnYvHlzwN7HbrfDYrG02XpDemwEVv9uCn5/wVAoZBLW7a3ARc9+i9W7SvrUqJDd5cb+cu/f6cgAjwRxcjQREZ2OkA9B1dXVcLvdSExMbLM/MTER5eXlAXufRYsWwWg0+ra0tLSAnftUVAoZ7rtoGD655xyMSDHA3OTEH97fjbnLt+NAhbXX6tGTCsqtcLoFonVKpEZrA3puTo4mIqLTEfIhqNXJE2mFEAGbXAsACxcuhNls9m3FxcUBO7e/hpsM+OjOqXjw4kyo5DJ8d7Aalzz/HR75eA9qGx29Xp9Aap20PCLFGNB+A05aOZqXxIiIyE8hH4Li4uIgl8vbjfpUVla2Gx06E2q1GgaDoc0WDAq5DHeePwTr7jsPM3IS4fYIvLmlEOc/vQGvf3cEzU53UOp1pvYc75n5QK1aJ0fzDjEiIvJXyIcglUqFcePGYf369W32r1+/HlOmTAlSrXpeRlwEXp0zHv+aNxHDTQZYml34+6f7cN4SbxiyOVzBrmK3tIaTkSlRPXL+1nlBnBxNRET+CvkQBAALFizA66+/juXLl2Pfvn247777UFRUhDvuuAOA91LW3Llz2xyTl5eHvLw8NDQ0oKqqCnl5edi7d28wqn9GpgyOw3/uOQdPXTMSJqMGlVY7/v7pPkx96mu8+NVBmJucwa7iKTU73b65TT01EtR6Xk6OJiIifymCXQF/XH/99aipqcFf//pXlJWVYcSIEfjss8+Qnp4OwLs44slrBo0ZM8b35127duFf//oX0tPTcezYsd6sekDIZRJmTxiAa8am4sMfSvDyxsMorLHhH+sPYOk3hzFzlAmzJ6Rh7IDogM+3CYR9ZRa4PAKxESqYjJoeeY/WydE2h3dy9NBEfY+8DxER9R2S6Ev3YQeQxWKB0WiE2WwO2vygzrjcHnz6Uxle/uYw9pf/fPfY0IRIXH92Gi4bZYLJGNg7sM7Em1uO4ZGP83F+ZjxW3jKhx97nulc2Y8exOjx7/WhcPSa1x96HiIhCV3c+v8NiJIjaUshluPKsFFwxOhk/FNXhne3F+M+PpThY2YC/f7oPf/90H0akGHDh8ERcODwROcmGoI4Qtc4HGhXg9YFOlpVkwI5jdSgo523yRER0agxBYUySJIxLj8G49Bg8cnk21uaVYk3ucfxQVIc9xy3Yc9yC5748iGGJkXh97tkYEBuYh5Z2V+tk5ZGpUT36PsOSvJfACsp7Z6FLIiIKbwxBfYRBo8SNk9Jx46R0VDfY8fX+Sny5twLfHazGgYoGXPfqZqy6bSKGJPTuXJkmx8+TogO9UvTJslpC0IEKjgQREdGphcXdYdQ9cZFqzBqfhmVzx2PjA+djWGIkKix2zHp1a6/fQr63zAyPAOL1aiQa1D36XsNaAt7x+iZYmkP/rjkiIgouhqA+LsGgwbu/nYyRKUbUNjrwq9e2YldhXa+9/4nzgXp6XpJRp/TdfXawjzxuhIiIeg5DUD8QHaHCqnkTcXZGNKzNLsx5Yxu2HanplfdufXL8yB5aH+hkw1pujT/xrjkiIqKOMAT1EwaNEv/8zQScOzQONocb89/NQ4O951edbn2WV0/PB2rlmxfEEERERKfAENSP6FQKLJszHgNidCgzN+N//1vQo+93tLoRh1qe6t5bIYgjQURE5C+GoH5Gq5LjiatHAAD+ueUYdhfX98j7uD0C97+/G0IA5w6NQ4KhZ1aKPlmm7w4xK7gOKBERdYUhqB86d2g8rjorGUIACz/8CS63J+DvseL7o9hVWIdItQKLrhkZ8PN3ZkhCJGQSUGdzospq77X3JSKi8MMQ1E/9eWY2onRK7C2zYMX3xwJ67sNVDXi65VLbny4bjtTo3lukUaOUIyMuAgBQwDvEiIioCwxB/VRcpBp/vGQ4AOCZ9QdQXGsLyHndHoEH3t8Nu8uDc4fGYfbZaQE5b3dkJrauHM0QREREnWMI6seuG5+KiQNj0OR045GP9wRkDs0bm47gh6J6RKoVeOqXo4LyzLLMJIYgIiI6NYagfkySJDx5zUio5DJsKKjC4i8KzigIHaq04n/XHQAA/GXmcKREBedJ9r6RIF4OIyKiLjAE9XOD4yPxyOXZAIBXNh7G3/6z77SC0J7jZvzqtW1wuDyYNiwes8b3/mWwVifeIebx8A4xIiLqGEMQ4cZJ6fj7Vd7b5pd/fxSPfJzfrfCwoaASs17dgiqrHVlJejx9XXAug7VKj42AWiFDs9ODogDNdSIior6HIYgAeIPQkl+OgiQBb20txB/X/ORXEHpnexFu++dO2BxuTB0Si/fumIwEfe+sCdQZuUzC0MRIALwkRkREnWMIIp9ZZ6fhH9eNhkwC/r2jGNe9ugXv7iiCuantE9mFECgot+Jv/9mLhR/+BLdH4JqxKVhx8wQYNMog1b6tYbxDjIiITkER7ApQaLlmbCqUchkWvJeHXYV12FVYh798nI8LshIwPiMGuUV12HqkBtUNDt8x914wFPddODSol8BO1voMMY4EERFRZxiCqJ3LRydjbHo0Pso9jo9yj+NgZQM+31OOz/eU+8polXKcPTAG149Pw2WjTEGsbcdaR4L4IFUiIuoMQxB1KCVKi7t+MQR3nj8Ye8ss+Cj3OA5VNmB0WhSmDI7DWWlRUClC92pqVpIBAHCkuhF2lxtqhTzINSIiolDDEERdkiQJOclG5CT3zlPgAyXRoIZBo4Cl2YUjVY0YbjIEu0pERBRiQvdXeaIzIEmSbzSIk6OJiKgjDEHUZw1L4m3yRETUOYYg6rMyORJERERdYAiiPqv1NvmtR2rw3cGqINeGiIhCDUMQ9VljB0Rj8qBY2Bxu3LxiB97eWhjsKhERUQhhCKI+Sy6TsPI3Z+OasSlwewT+/NEe/PWTvXDzoapERASGIOrj1Ao5/nHdaNw/fRgA7wNif/vmTlibnac4koiI+jqGIOrzJEnC3f8zFC/dMAZqhQxf7a/EFS99j72llmBXjYiIgoghiPqNmaOS8e7tk5Fs1OBodSOuWvo9/rWtCELw8hgRUX/EEET9yllpUfj03nPxP1kJcLg8+OOan/D7f+ehwe4KdtWIiKiXMQRRvxMdocLrc8dj4SVZkMskrN1dipkvfIftR2uDXTUiIupFDEHUL8lkEm6fNhjv/nYSTEYNjtXYMOvVLXjk4z0cFSIi6icYgqhfG58Rg//edx5+NSENAPDmlkLMePZbLq5IRNQPMARRv2fQKLHomlFYddtEpEZrcby+CXPe2I67//UDimttwa4eERH1EIYgohZTh8Thv/PPw81TMiBJwH9+LMMF/9iIJz/bB7ON6woREfU1kuD9wR2yWCwwGo0wm80wGAzBrg71sr2lFjz52T5sOlQNAIjSKXH3L4bghokDoFMpglw7IiLqTHc+vxmCOsEQREIIfHOgCk9+ug8HKxsAeMPQ3MkZuGlyOmIj1UGuIRERnYwhKAAYgqiVy+3B+7tK8MrGwyis8c4RUitkuG58KuZMykBmy9PqiYjIf81ON5qdbkTpVAE9L0NQADAE0cncHoH/5pfj1Y2HsbvE7NufbTLgmrEpuGJ0MhIMmiDWkIgodAghUFRrw95SC/JLLSistaHS0oyqBjuqrHZYm124ZEQSXr5xXEDflyEoABiCqDNCCGw7WosV3x/F1/sr4XR7/wvJJGBUahSMWiV0Kjm0Sjk0KjkcLg8a7S402F1otLvgcHugUyqgVckRoZZDp1IgQa9GRmwE0mN1yIiLQIJeDUmSgtxSIiL/eTwCe0rN+KagCt8fqsbeUgusp1h3beLAGLx7++SA1qM7n9+c4UnUTZIkYdKgWEwaFIu6Rgf+81MZ1vxQgh+K6pFXXB+Q94hQyTFmQDTOzojB2QOjMSYtGlqVPCDnJiIKFLvLja/2VeLLvRX49mAVqhscbV5XyWXITNIj22TAkIRIJBjUSNBrEK9XI16vhkET3BjCkaBOcCSIuquwphF7jlvQ5HSjyeFq+eqBUiEhUq1AhEqBCLUCKoWEJocHtpYyDXYXyuqbcaymEYU1NpTU2eA56X+lQiZhbHo0LhtpwiUjknjZjYiCas9xMz7YVYKP8o6j/oQlRCLVCkwdEotpwxIwNj0Kg+MjoZT37mo8ffJy2NKlS/H000+jrKwMOTk5eO6553Duued2Wn7jxo1YsGAB8vPzkZycjAcffBB33HGH3+/HEETB4nB5cLiqATuP1WL7sTrsOFqLckuz73VJAs5Oj8GlI5Nw1ZiUgE8qJCLqiMPlwUd5x7Hy+2PYW2bx7U8yaHDFWcn4RWYCxqVHQ6UI7hKEfS4Evfvuu5gzZw6WLl2KqVOn4tVXX8Xrr7+OvXv3YsCAAe3KHz16FCNGjMC8efNw++234/vvv8edd96Jd955B7/85S/9ek+GIAoVQggU1zZh3d5yfPpTGXKL6n2vaZQyXD0mBTdNyUBWEv+dElHgNTnc+PeOIrz27RGUmr2/kKnkMlyUk4jrxqXi3KHxkMtCZw5jnwtBEydOxNixY/Hyyy/79g0fPhxXXXUVFi1a1K78Qw89hLVr12Lfvn2+fXfccQd2796NLVu2+PWePRaChACcfBQDnb4ycxP+m1+Bj/JKUVD+829jEzNicfloE9JidDAZtUg0qIP+GxkRha8mhxtvbi3EPzcfQ53NO9cnLlKNm6dk4JdjAzgKrdR5h7gDpE9NjHY4HNi1axcefvjhNvunT5+OzZs3d3jMli1bMH369Db7ZsyYgTfeeANOpxNKpbLdMXa7HXa73fe9xWJpVyYgnDbgyeSeOTf1CyYAN7dsOHFqUHnLRkQUAFoAt7dsvp81LgDftmyB8sdSQBURwBP6L+R/Tayurobb7UZiYmKb/YmJiSgv7/gnfnl5eYflXS4XqqurOzxm0aJFMBqNvi0tLS0wDSAiIqKQFPIjQa1OXjNFCNHlOiodle9of6uFCxdiwYIFvu8tFkvPBCGlzpt6iXqQEAK1jQ6UmptRUmdDUW0TimobUVzThKM1jahusHd4nEwC0mMjkJWkR06yESOSjchJMSBSHTY/Koiom4QQ2Hq0Fm9tOYYNBVUAAJ1Kjt+eNxg3TU6HRtnDy3ModT17/i6E/E+2uLg4yOXydqM+lZWV7UZ7WiUlJXVYXqFQIDY2tsNj1Go11OpeeBaUJAVt2I/6DwlArBqIjQFGDmz/urnJiUOVDThc2YCDlVbsL7diX5kF1Q0O7K12Y291PT7cU+89lwQMjo/EWWlRGJcejfHp0RgcHwlZCE2EpPZcbg/KzM2wNDvR0OyCtdmFRocLqdE65CQbev6DrQ+rabD7lrQorLGhuNbm+7u2NrtgbfnqEQIKmQwKuQS5TIJaIUNshBpxehXiItWIj1QjOUqLgfERGBQXgZQoLRS9eDt5s9ONNbneu70KKqwAAEnS4NqxqXhgRma/WIoj5EOQSqXCuHHjsH79elx99dW+/evXr8eVV17Z4TGTJ0/GJ5980mbfunXrMH78+A7nAxH1N0atEuPSozEuPbrN/kprM/aVWbHnuBk/lZjxY0k9Ss3NOFTZgEOVDfhgV4nv+DEDopCVZMDQhEgMS9RjSEIkF3QMMofLg+8PVeOzn8qwbm8FzE3ODsspZBKyTHqMTo3C2RkxuGRkEtQK9l1Hmp1u5BXXY3dxPXaX1GN3sRnH65v8Pt7h9sDh9v7ZCqC6wYGCio7LKuUS0qJ1iNOrERuhQkyECrERKhi0SkSoFdCp5IhQeb8q5DLIZd6rG/KWKxxOtwcOtwdOt/CtVG+1/xzK6m1OVLc8sqJ1c7g9ALwjP9eOS8VNUzIwOD7yTP7KwkpY3B3Weov8K6+8gsmTJ2PZsmV47bXXkJ+fj/T0dCxcuBDHjx/Hm2++CeDnW+Rvv/12zJs3D1u2bMEdd9zBW+SJTkOV1Y4fS+rxQ1EddhXWIa+4Hs1OT7tykgTERqgRpVMiSqtElE6FKJ0SCXo1Eg0aJBrUSDBoYDJqkKjXcCQpgA5XNWDphsNYt7cc1uafH1OgUshg0Chh0Cig1yigUcpxuKqh3aq+iQY1fnveYPxqQhp0qpD/3bjHmZuc+Hp/Bf67pwIbD1Shyelu87okAclGLQbE6LxbrA4pUVoYtUroNQoYWr7KJQkuj4DbI+DyCDQ53KhptLcJIsW1TThW04ij1Y2wu9r/v+ppqdFa3DQ5A7POToNR2zcGCfrcLfKAd7HEJUuWoKysDCNGjMCzzz6L8847DwBw880349ixY/jmm2985Tdu3Ij77rvPt1jiQw89xMUSiQLA6fZgX5kFu4vrcaDCezntYEUDahodpz64hUohQ1r0zx8igxMiMSQhEkMT9IiLVPG5aX5qtLvw4teH8MamI75n2MXr1bhkRBIuGWHChIEx7dZvEUKg1NyM3cXex7x8srsUZS1rv8REqHDrOQMxd3I69Jq+8YHoL6fbg3X5Ffj3jiJsOVwD1wnLtica1Bg7IBqj06IwOjUKI1IMAf/78XgEyizNKKxpRE2DA7WNDtQ0OlDb6H3QaKPdDZvDhUaHGza7C26PgEcIuIWApyU7KeUSlHKZd1PIEKmWI1KtgF7jDWVGrRLxerX3UpzeezkuJUrb534h6ZMhqLcxBBF1T02DHRUWO+qbHDDbnKhvcqK20YEqqx0VluaWzY5ySzPcJz8X5ATROiWGJuoxItmIESkGjEgxYlBcRK/OlQh1Qgh8+lMZnvh0ny/A/CIzHr87fwjGp0d360PN7nLjwx+O4+VvDqOo1ruGmcmowUs3jMG49JgeqX8oKTc341/bi/Dv7UWotP58w8CwxEjMyEnCjJwk5CQbGMzDCENQADAEEfWM1gm7RbU2FNXacKymEYcrG3Gw0oqiWhs6+omkUcowKiUKY1vmMY0dEIXYyF64kSEE7S214O+f7sXmwzUAvJczHr08BxcOTzijD2qX24P//FiGZ9YfQFGtDQqZhAcvzsS8cwf1uQBQb3Pgy32V+GJPOTYUVPpCeVykGr+akIarxqT0q3kxfQ1DUAAwBBH1vmanG4cqG1BQbsVPx83ILzUjv9QCm8PdrmxqtBbJRq3vqdSJBjWidSpEahQtlwC8m67lwbURLZNJw1WFpRn/WFeA93eVQAjvJcXfTRuM350/OKB3elmbnVj44U/4z49lAIALhyfgf68bHdbPqPN4BA5VNWDrkRr8N78cW4/UthmNnJARgxsnp+PinCSust4HMAQFAEMQUWhwewSOVjfgh6J6/FBYh52FdThU2XBa51IpZIhQyaFVyqFVeTedUgG1Uga1Qg5Ny1e1UgaVXAaV4uevOpXcG6ZaAlWUToUBMboencNkbnKisKYRX++vxLJvj/jC4MxRJjx0cRbSYnpmfRUhBP61vQiPf7IXDpcHSQYNzs+MR2aSHpmJemQm6XtlJE4I4Zss7E/QE0Kg3ubE8fomlNQ1YW+pGbnF9cgrqofV7mpTNitJj+k5SbhspAmZSfoeqT8FB0NQADAEEYWuepsDByoaUGn1zjOqtDSj0mqHucm7Jo6l2YkGuwsNdhdsdrfvNuCeEKGSIy1Gh/RYHRINGsS03NocE6FClFYFlUIGpVzyBSq3ELA7vbcy250eNDpc3kmwDXZUNzhQ3WDH8fomHKtuRJ2t7S3uYwdE4c8zszF2QHQntQms/FIz7lr1A47VtH/eoUGjQHKU1jsiF6WFyahFlK7l7qiWibhalRwSJEiSd+0qALA53L6+abR7215ubkaZpRnlZu/W6HCh2emG3eXxXR7VKGWI0akQpVMhOkIJuUwGp8sDp9u7NTrcKK1v6nDUEAC0SjlGpRrxP1kJmJGThIw4rtfWVzEEBQBDEFHf4XB5YHN4P3ibHG7YHG40Od1oavna+oHb+tUbUty+9VbsLcc32t1otHsXHaxpcKDU3NThHKZAiterMSguAnMmp+OykaZen59jc7iwYX8V9pdbUFBuRUFF53O3QkVcpBopURoMSdBjzIAojBkQhcxEfVhfDiX/MQQFAEMQEZ2K3eVGSV0Timq8k7yrrHbU2hyobbnF2dzkPGEBOw8cLg/kMhnUCu+mUsigUcoRF6lCbIQasZEqxEaqYTJqkBEbgfRYHSJC8JElNocLx+uacLzeu5XWN3lXTG5y+VZNtjQ5YXe5IQQg8POji3Qq75ytCLX38mKUTgWTUYOk1jWkjBro1d41jdRKme8ymNnmvduw1uZAvc0BIfDz7eByCRqlHCajBslRWq6G3c8xBAUAQxAREVH46c7nN8cGiYiIqF9iCCIiIqJ+KfQuNoeI1quEFoslyDUhIiIif7V+bvsz24chqBNWqxUAkJaWFuSaEBERUXdZrVYYjcYuy3BidCc8Hg9KS0uh1+vDesl4i8WCtLQ0FBcX98kJ3mxf+OvrbWT7wl9fb2Nfa58QAlarFcnJyZDJup71w5GgTshkMqSmpga7GgFjMBj6xD/uzrB94a+vt5HtC399vY19qX2nGgFqxYnRRERE1C8xBBEREVG/xBDUx6nVajz66KNQq3v+YYfBwPaFv77eRrYv/PX1Nvb19nWFE6OJiIioX+JIEBEREfVLDEFERETULzEEERERUb/EEERERET9EkNQGPn2229x+eWXIzk5GZIk4aOPPuqy/DfffANJktpt+/fvb1Nu9erVyM7OhlqtRnZ2NtasWdODrehcT7Rv5cqVHZZpbm7u4da01932AYDdbsef/vQnpKenQ61WY/DgwVi+fHmbMqHSf0DPtDGc+/Dmm2/usO45OTltyoVKH/ZE+0Kp/4DT+ze6atUqjB49GjqdDiaTCbfccgtqamralAnXPgRO3b5Q68NAYggKI42NjRg9ejReeumlbh1XUFCAsrIy3zZ06FDfa1u2bMH111+POXPmYPfu3ZgzZw5mzZqFbdu2Bbr6p9QT7QO8q6Ce+HpZWRk0Gk0gq+6X02nfrFmz8NVXX+GNN95AQUEB3nnnHWRlZfleD6X+A3qmjUD49uHzzz/fps7FxcWIiYnBdddd5ysTSn3YE+0DQqf/gO63cdOmTZg7dy5uvfVW5Ofn4/3338eOHTtw2223+cqEcx/60z4gtPowoASFJQBizZo1XZbZsGGDACDq6uo6LTNr1ixx8cUXt9k3Y8YMMXv27ADU8vQFqn0rVqwQRqMxoHULBH/a9/nnnwuj0Shqamo6LROq/SdE4NoYzn14sjVr1ghJksSxY8d8+0K1DwPVvlDtPyH8a+PTTz8tBg0a1GbfCy+8IFJTU33fh3Mf+tO+UO7DM8WRoH5gzJgxMJlMuOCCC7Bhw4Y2r23ZsgXTp09vs2/GjBnYvHlzb1bxjHTVPgBoaGhAeno6UlNTMXPmTOTm5gahlt23du1ajB8/HkuWLEFKSgqGDRuG+++/H01NTb4y4d5//rQRCN8+PNkbb7yBCy+8EOnp6b594d6HJ+qofUB499+UKVNQUlKCzz77DEIIVFRU4IMPPsBll13mKxPOfehP+4Dw7sOuMAT1YSaTCcuWLcPq1avx4YcfIjMzExdccAG+/fZbX5ny8nIkJia2OS4xMRHl5eW9Xd1u86d9WVlZWLlyJdauXYt33nkHGo0GU6dOxcGDB4NYc/8cOXIEmzZtwp49e7BmzRo899xz+OCDD3DXXXf5yoRz/wH+tTGc+/BEZWVl+Pzzz9tdZgj3PmzVWfvCvf+mTJmCVatW4frrr4dKpUJSUhKioqLw4osv+sqEcx/6075w78MuBXsoik4PTmOoWgghZs6cKS6//HLf90qlUvzrX/9qU+btt98WarX6TKt4RgLVvpO53W4xevRocc8995xB7c6cP+276KKLhEajEfX19b59q1evFpIkCZvNJoQI3f4TInBtPFk49eGJnnzySREbGyvsdnub/aHah4Fq38lCpf+E8K+N+fn5wmQyiSVLlojdu3eLL774QowcOVL85je/8ZUJ5z70p30nC6U+PFMcCepnJk2a1Ca9JyUltfttpbKyst1vNeHi5PadTCaT4eyzzw6L32BMJhNSUlJgNBp9+4YPHw4hBEpKSgCEf//508aThVMfthJCYPny5ZgzZw5UKlWb18K9D4Gu23eycOu/RYsWYerUqXjggQcwatQozJgxA0uXLsXy5ctRVlYGILz70J/2nSzc+rArDEH9TG5uLkwmk+/7yZMnY/369W3KrFu3DlOmTOntqgXEye07mRACeXl5XZYJFVOnTkVpaSkaGhp8+w4cOACZTIbU1FQA4d9//rTxZOHUh602btyIQ4cO4dZbb233Wrj3IdB1+04Wbv1ns9kgk7X9qJTL5QC8bQHCuw/9ad/Jwq0PuxSsISjqPqvVKnJzc0Vubq4AIJ555hmRm5srCgsLhRBCPPzww2LOnDm+8s8++6xYs2aNOHDggNizZ494+OGHBQCxevVqX5nvv/9eyOVy8dRTT4l9+/aJp556SigUCrF169Y+0b7HHntMfPHFF+Lw4cMiNzdX3HLLLUKhUIht27aFfPusVqtITU0V1157rcjPzxcbN24UQ4cOFbfddpuvTCj1X2udA93GcO7DVjfeeKOYOHFih+cMpT7sifaFUv8J0f02rlixQigUCrF06VJx+PBhsWnTJjF+/HgxYcIEX5lw7kN/2hdqfRhIDEFhpPWW8JO3m266SQghxE033SSmTZvmK7948WIxePBgodFoRHR0tDjnnHPEp59+2u6877//vsjMzBRKpVJkZWW1CRG9qSfaN3/+fDFgwAChUqlEfHy8mD59uti8eXMvtupn3W2fEELs27dPXHjhhUKr1YrU1FSxYMGCdnNlQqX/hOiZNoZ7H9bX1wutViuWLVvW6XlDpQ97on2h1H9CnF4bX3jhBZGdnS20Wq0wmUzi17/+tSgpKWlTJpz78FTtC7U+DCRJiE7Gu/o5j8eD0tJS6PV6SJIU7OoQERGRH4QQsFqtSE5Obnepr6PCYWfjxo1i5syZwmQydTj73ePxiEcffVSYTCah0WjEtGnTxJ49e7r1HsXFxR2maW7cuHHjxo1b6G/FxcWn/KxXIAy1Lgt+yy234Je//GW715csWYJnnnkGK1euxLBhw/D3v/8dF110EQoKCqDX6/16j9ZyxcXFMBgMAa0/ERER9QyLxYK0tDS/Pu/D/nKYJElYs2YNrrrqKgCAEALJycmYP38+HnroIQDeBzQmJiZi8eLFuP322/06r8VigdFohNlsZggiIiIKE935/O5zt8gfPXoU5eXlbZYwV6vVmDZtWpdLmNvtdlgsljYbBV6j3YU3txyDuckZ7KoQEVE/1+dCUOuCVd1dwnzRokUwGo2+LS0trUfr2V+9/t1RPPJxPp5ZVxDsqhARUT/X50JQq5Pv6BJCdHmX18KFC2E2m31bcXFxT1exX/rpuBkAsOVITZBrQkRE/V1YTozuSlJSEgDviNCJq1meaglztVoNtVrd4/Xr7w5VWgEAByoaYLY5YdQpg1wjIiLqr/rcSNDAgQORlJTUZglzh8OBjRs3hsUS5n1Zs9ONwlqb7/sfiuqCWBsiIurvwnIkqKGhAYcOHfJ9f/ToUeTl5SEmJgYDBgzA/Pnz8eSTT2Lo0KEYOnQonnzySeh0Otxwww1BrDUdrmrAifci7iysxS+yEoJXISIi6tfCMgTt3LkTv/jFL3zfL1iwAABw0003YeXKlXjwwQfR1NSEO++8E3V1dZg4cSLWrVvn9xpB1DMOVjS0+X7nMY4EERFR8IT9OkE9hesEBd7T/92P/9twGBMGxmD70VpolDL89NgMKOV97qosEREFSb9eJ4hCV+tI0MU5SYjSKdHs9CC/lOsxERFRcDAEUa85WOkNQZlJeowbEA0A2HmsNphVIiKifowhiHpFs9ONwppGAMDQhEiMy/CGoF2FnBdERETBwRBEveJIVSM8AjBqlYjXqzE+PQYAsLOwDpyWRkREwcAQRL3iYMsiiUMTIiFJEkalGqGUS6iy2lFc2xTk2hERUX/EEES9onVS9NBE7zIFGqUcI1KMALzrBREREfU2hiDqFSeOBLUan94yOZrzgoiIKAgYgqhXtI4EDUv8ecHKcS3zgnZx0UQiIgoChiDqcXaXG8da7wxL/HkkaFzLSNCBSivMTc6g1I2IiPovhiDqcUervXeGGTQKJOjVvv3xejXSY3UQgg9TJSKi3scQRD3uwAmToiVJavNa62gQL4kREVFvYwiiHneoov2k6FY/rxfEO8SIiKh3MQRRjztw0u3xJxrfsnJ0XnE9nG5Pr9aLiIj6N4Yg6nEd3R7fakh8JIxa78NU9/JhqkRE1IsYgqhHee8MswFoe3t8K5lMwtgBUQC4XhAREfUuhiDqUceqbXB7BPRqBRIN6g7LjM9oWS+I84KIiKgXMQRRj/JdCkuMbHdnWKvWO8R2HuPDVImIqPcwBFGP8k2KTmh/KazV6NQoKGQSKq12lNTxYapERNQ7GIKoRx06YSSoM1qVHDl8mCoREfUyhiDqUV3dHn+i8SdcEiMiIuoNDEHUYyqtzThW7X1m2LAuRoKAn0PQLt4hRkREvYQhiHqE3eXGHW/tgssjkJWkR5JB02X5cS2LJhZU8GGqRETUOxiCKOCEEPjLR3vwQ1E99BoFlv56bKd3hrVK0GswIMb7MNVcPkyViIh6AUMQBdzKzcfw3s4SyCTgpRvGYlB815fCWvGSGBER9SaGIAqo7w9V4++f7gMALLxkOKYNi/f72NZLYpwcTUREvYEhiAKmsKYRd676AW6PwDVjUnDbuQO7dXzrE+X5MFUiIuoNDEEUENZmJ+a9uRPmJidGp0XhyWtGnnIe0MmGJkTCoFGgyenGvjI+TJWIiHoWQxCdMY9H4L53d+NARQMS9GosmzMOGqW82+eRySSM5XpBRETUSxiC6Iw9++UBfLmvAiqFDMvmjkfiKW6H7wonRxMRUW9hCKIz8umPZXjx60MAgEVXj8RZaVFndL5xLfOCdhbW8mGqRJ3YW2rB1iM1/D9CdIYUwa4Aha/8UjPuf383AGDeuQPxy3GpZ3zOs9K8D1OtsHgfppoWozvjcxL1JZsPVeOmFdvhdAuMT4/Gw5dkYXxGTLCrRX2A2ebEtqM1GBQfiYFxEZDLujevMxwxBNFpqW6w47dv7kKT043zhsXj4UuGB+S8WpUcOckG7C4xY1dhHUMQ0QnyS8347Vu74HR7R4B2Ftbh2le24MLhCXhgRhYyk7p+Rh9RZyotzZj92lYcqfI+6kijlCEryYCcZAOykw3INhmQlWSAVtX9+Z6hjCGIus3h8uDOt3/A8fomDIyLwIuzxwT0N4Zx6THYXWLGzsJaXDUmJWDnJQpnxbU23LxiBxrsLkwcGIOnrx2Nlzcewns7S/Dlvkp8tb8S14xJxYLpw5ASpQ12dcNCs9ONgnIr9pSasee4BfmlZhyrbsSkQbG4eWoGJg+K7fZdruGo0tqMX7UEIL1GAZdboMnpRl5xPfKK633lZBIwMC4C2clGZJt+DkhxkergVf4MSaKPX1ReunQpnn76aZSVlSEnJwfPPfcczj333FMeZ7FYYDQaYTabYTAYeqGm4eOPa37Cv7YVQa9WYM1dUzEkwb8Vof312U9luHPVD8hK0uOL+ecF9NxE4aiu0YFfvrIZR6oakZmox3t3TIZRqwQAHKpswD/WFeDzPeUAAJVChrmT0nHXL4YgOkIVzGqHlEa7C/vKLNhz3Iw9pd6vBysb4PZ0/hGYmajHzVMzcNVZKX1uBKRVldWO2cu24HBVI5KNGvz7t5OREq3FsZpG7C21IL/Ugr1lFuwtNaO6wdHhORL06hNGjIzITjYgPUYHWZAup3Xn87tPh6B3330Xc+bMwdKlSzF16lS8+uqreP3117F3714MGDCgy2MZgjr21tZC/OWjPZAkYPlNZ+MXWQkBf49KSzMmPPkVJAnY/eh0GDTKgL8HUahzewSqrHaUmZvw1//sRW5RPZKNGnx451QkGdvfgZlbVIfFX+zH1iO1AAC9WoHbpw3Cb84ZCJ0qfAb9hRBnPPpitjmRX2rGnlIz8lsCz5HqRnT0aRcTocKIFCNGJBswIsWIRIMGH+Uex+ofSmBzuAEASrmE9NgIDIyLwKC4CAyKj8DAOO+8mbhIVdiOFlVZ7fjVa1txqLIBJqMG//7tJKTHRnRavtLa7AtG+8os2FtqwdGajv9eI1RyZJkMbUaMhiXqT2v5lO5iCGoxceJEjB07Fi+//LJv3/Dhw3HVVVdh0aJFXR7bYyFICMBpC9z5etG2I7W47c0dcHkEFlyUiXndXBG6Oy56diNK6pqwbO54nDskrsfepydVN9ix8vtj8AiBLJMew01GDIzVQSHnTZn9ndsjUNXQjHKzHRWWZpSZm1Fhbka5tRnlLX+ubLC3GaUwapR4+7aJXY68CiHw3aFqPLP+IArKvQuOxkWqcef5g3HtuFQog/xvz+0RqG6wo8zsbWe5pRll9c0otzShzGJHeX0T6mwOnJ0Rg3svGOrX3abVDXbsK7dib6kF+8q8oaekrqnDsq0jFsNbRitykvVI1Gs6DDHmJifW5B7Hqm2FnZ4PAPQaBTJaApI3KOmQHhuBjFhdUMKn2yNgbXbC3OSCpcmJ+mYnzDaH988t+8wt2/4yC8oszUgyaPDP30zAgNOYg9lod+FAhRX7y63YV+79eqDcArur/ar/cpmEQXGRGG7SIyvJgCyTHpmJesRERQEBDJIMQQAcDgd0Oh3ef/99XH311b79v//975GXl4eNGze2KW+322G3233fWywWpKWlBT4EORqBJ5MDdz4iIqJw9sdSQNX5CFR3dScE9dlfSaurq+F2u5GYmNhmf2JiIsrLy9uVX7RoEYxGo29LS0vrraoSERFREITPheLTdPIwZ2fXmxcuXIgFCxb4vm8dCQo4pc6benuY0+3B0WobDlRYEB+pwcRBp7eOiMcjMP/dPKzfV4HYCDXev2MSTMaev/OkoMKKq/7ve+hUcmxbeEFYXELyeARe+uYwXv7Gu3jk1MGx+Md1Z8Go885pcnsECmtt2F9mwf5yK/aWWbC/zIKaxo4nG5oMGmSZDBhuMiDbpMfwZANMhrZD9zUNdrz8zWG8u7MYri4meJ4OlUIGrVKOCJUCWpUcOpUcWqXc9+f6ljVFWt9WpZAhLVoHa7MT1mYXmpzugNQjUq2AXq2AXquEQa2EXiuHXqNEldWOrUdrfPMRWt9fLpOgkEuQy2SQS4BCkrXs8/48kMskKFrKyHzfe8vIZBIUJ5SRy71f7U4PdhTWIb/U7FedFTIJCXoNkowtm8H7NVGvRpJRiySjGrER6l5bh6XB7sI/NxdixfdH0Ngyz2XKoFjcd1EmRqR0/puy2yPwwlcH8dbWQr/6UyYB8ZEaJEVpkGRQI8mgRZJRA5NRg0SD92tcpH/tPlLdiJe+Puib8A14L6cMjo9Edsv/i5xkAzKT9IhUh95HWV2jA4W1NhytbkRhTSOOVjfiaLUNRbWNHV4q6ohMAoxaFYxaJQxaJYxahffPGiWMOhWiWvYZtEoYW/YZtAoYNUqoFKH/MxPw/tyUKYO3FAovh3UiXCZGCyFQam5GQbn3g7WgZTtc1eBbS0QmAV/94XwMjOv+cONzXx7Ac18ehEouwzu/nehb0bmneTwCo/+6DtZmFz65+xyMTDX2yvt2h7nJ2TJJ0DsPYXdxPY5Ue9fYuO2cgXj4kqxThjchBCqtdt959pZ5Jx0W1nQ8byxKp/RNNFQr5Pjn5mOw2l0AgF9kxuPXE9PhFgI2hws2hxs2u9v7teX7RocLTY62+3x/trthc7q7vFumI+PTo/HLcam4dKTJd8cS4F1KwdISiCxNTliavfMQTvze0uRq2eeE5aRyzU7/PigmZMTg6rEp7d6/J1RZ7fjuYBU2HqhCmdk7l8LU8iHv/bDXwhSlQVyEOmh3xnSlusGOl74+hFXbCn0/Hy4bZcL90zPb/XywNDtx7zu5+KagCoB3ykZcpBrJLe1MMmqQHOX9c3KUBklGLRL06oDPO9pfbsFPJWYMTdQjK6l3Jtb2JI9HoNTc1BKKGuF0Cxg0Cug1St/XKJ0SRp0SkSpFSP47CnWcE9Ri4sSJGDduHJYuXerbl52djSuvvDJ4E6PPgNnmxP5yCwoqfg47BRVWWJtdHZaPVCsgl0kwNzmx8JIs3D5tcLfe74s9Zbjj7R8AAEt+OQqzzu7dS4Q3Ld+OjQeq8Ojl2bhlas9NwvZHpbUZ+aUW5B/3Bp78UguKatsHFZVChieuGoHrxp/Z35Wl2Yl9Lbemtr7fwQprh6M9I1IM+OMlwzElABPIhRBwuD2+QGSztw1PjQ43mlrCEwD8T1ZCl3eTnAmn2+MLTL7w1Oz0BSWZJGFGThIX1DwNxbU2PLP+AD7KOw4hvCNXsyek4d4LhiJBr8HR6kbc9s8dOFzVCI1ShsW/HIVLRpjCZnSB+jeGoBatt8i/8sormDx5MpYtW4bXXnsN+fn5SE9P7/LYYIYgu8uNQ5UNvqDTOsJTbmnusLyiZYg4M0mPzCTvb0uZSXqkRGnx9tZC/OXjfJydEY3375jidx32l1twzdLNsDncuHlKBh67IidQzfPbi18dxD/WH8Blo0z4vxvG9sp7CiFQUtfkvb22ZfG0PaUWVFntHZZPidJiRIoBOclG5CQbcFZaFGJ7aOEwu8uNgxUNvlGjMnMzLh1pwhWjk/nbIp2WfWUWLPliPza0jPZolXJcNz4VH+Ueh6XZBZNRg2VzxofkSCxRZ7rz+R16F1ID6Prrr0dNTQ3++te/oqysDCNGjMBnn312ygDUWzwegeI6W5vLWAUVVhytbuz0kkRKlBZZSXoMOyHsDIqL7PQ3tAuGJ+IvH+djV2Edahrsfn1A1zY6cNs/d8LmcOOcIXH482WBeSRGd43LaHmi/LG6gKwdcjK3R+BIVYN3LZHjrSMuZlg6GFmTJGBQXARGpHjDTmvoidL13mJ0aoXcu55JihEAJ+7TmRtuMmDFLROw7UgNnvpiP3KL6vHmlkIAwJgBUXh1zjgk6NuvSUTUV/TpkaAz0VMjQXuOm/HWlkLsr7DiYIXVd1nhZEatss2oTlaSHsMS9dCfxsKBM1/8DnuOW7Dk2lGYdYrLNE63B3Pf2I4tR2qQHqvDx3dN7dUP+hPZHC6MfGwd3B6BTQ/9AqnRp3/Zw+5y40B5S+BpGeXZX27pcN6JUi5hWKIeOS2Lp+Uke5+ZExGCky+JAkUIgXV7K7Ds2yPISTbgT5cNh1oR3vNvqH/iSFAIq2l04N2dxb7vVQoZhiZEIjNRf8LlLAMSDeqAjXxcNDwJe45b8OXeilOGoL/9Zy+2HKlBhEqO1+aOD1oAAgCdSoGcZAN+bHmYandDkBACi78owMYDVZ3Op9Gp5L6JxjnJRuSkGDA0Qc+5D9TvSC1zrGbkJAW7KkS9hiGol41INuDe/xmCzCTvrZ0ZvbCC8IXZCXj2ywP47mA1mp3uTu+ueGd7Ed7cUghJAp6bPQbDEoP/ROpx6dH4scSMncfqcOVZ3XuY6t4yC17ZeNj3fZRO2eZS1ogUIzJiI3rtNmUiIgotDEG9LDZSjQXTM3v1PbNNBqREaXG8vgmbDlbjwuzEdmV2HKvFIx/vAQD84aJhuKiDMsEwPj0GK74/hp2Fdd0+ttzsnUg+JCES//zNBCQbO14en4iI+ieO+fcDkiThwuHeB51+ua+i3evH65twx1u74HQLXDbKhLt+MaS3q9ip8S2TowvKLbA2O7t1bGXLHV3pMTqkRGkZgIiIqA2GoH7iomzvdf4v91XCc8LcmCaHG799cydqGh3INhnw9LWjQiosJBo0SI3WwiOA3KL6bh1bafGGoARDz9yyTkRE4Y0hqJ+YOCgGeo0C1Q125JXUA/BOHH7gg93IL7UgNkKF124aH5SnHp/K+HTvaFB3L4lVWr2Xw+J5iy8REXWAIaifUMplOD/Te0ls/V7vJbGl3xzGf34sg1Iu4eUbxyElquefCXY6xmV4H9Wxq7C2W8e1Xg5L0HMkiIiI2mMI6kdaJzuv31uBL/dW4H/XFQAAHr9iBCYM7J1ngp2O1pGg3KJ6uNz+PU8KgG+V53iGICIi6gBDUD8ybVg8FDIJhyobcM87uRACmDMpHTdMHBDsqnVpWKIeerUCNocb+8utfh9XxZEgIiLqAkNQP2LUKjFpUCwAoMnpxqRBMXjk8uwg1+rU5DIJY1rnBR3z75KYEOLnEGTgnCAiImqPIaifuXiE9y6x1Ggtlv56HJQ9vFBjoHR3cnS9zQlHy6Wz+B56oCkREYW30LsViHrU7LPToJBJOG9YPGIigvdIjO5qDUG7/AxBrZOio3VKPgKDiIg6xE+HfkYhl2H2hAFIDtE7wTpz1oAoyGUSyszNOF7fdMryP98ez1EgIiLqGEMQhQWdSoFsk/dpwP7MC/p5UjTnAxERUccYgihsjOvGJTGuEURERKfCEERho/U5YjuP+RGCWh6ZEc9HZhARUScYgihsjE/3Lui4v9yCBrury7Ktc4J4OYyIiDrDEERhI8moQUpU68NUux4NquRq0UREdAoMQRRW/L0kxtWiiYjoVBiCKKz4u14QQxAREZ0KQxCFlXEt84Jyi+o6fZiqzeHyzRniIzOIiKgzDEEUVjKTvA9TbeziYaqtd4bpVHJEqrkoOhERdYwhiMKKXCbhrAFRAIAfOpkczTWCiIjIHwxBFHZab5XvbHI0H5lBRET+YAiisNN6h1hnk6NbL4dxjSAiIuoKQxCFnbPSvA9TPV7fhDJz+4epVjVwjSAiIjo1hiAKOxFqBYab9AA6viTmGwniIzOIiKgLDEEUllrnBXV0SYyPzCAiIn8wBFFYan2i/M7C2navVfGRGURE5AeGIApLrZOj95a2f5gqb5EnIiJ/MARRWDIZtb6HqeYV1fv2O1we1DY6ADAEERFR1xiCKGx1dEmsptE7CqSQSYjWqYJSLyIiCg8MQRS2OlovqPXOsHi9GjKZFJR6ERFReAjLEPTEE09gypQp0Ol0iIqK6rBMUVERLr/8ckRERCAuLg733nsvHA5H71aUelTrSFBuUT3cHgHg5/lAnBRNRESnEpYhyOFw4LrrrsPvfve7Dl93u9247LLL0NjYiE2bNuHf//43Vq9ejT/84Q+9XFPqSVlJBkSqFWiwu7C/3ALgxNvjGYKIiKhrYfmI7ccffxwAsHLlyg5fX7duHfbu3Yvi4mIkJycDAP7xj3/g5ptvxhNPPAGDwdBbVaUeJJdJGDMgCt8drMauwjrkJBtPuBzGNYKIiKhrYTkSdCpbtmzBiBEjfAEIAGbMmAG73Y5du3Z1eIzdbofFYmmzUejzTY5uWTm69ZEZHAkiIqJT6ZMhqLy8HImJiW32RUdHQ6VSoby8vMNjFi1aBKPR6NvS0tJ6o6p0hk5eOZqPzCAiIn+FTAh67LHHIElSl9vOnTv9Pp8ktb8zSAjR4X4AWLhwIcxms28rLi4+7bZQ7zlrQBRkEnwPU63iIzOIiMhPITMn6O6778bs2bO7LJORkeHXuZKSkrBt27Y2++rq6uB0OtuNELVSq9VQqzl6EG4i1QoMNxmQX2rBzmN1vDuMiIj8FjIhKC4uDnFxcQE51+TJk/HEE0+grKwMJpMJgHeytFqtxrhx4wLyHhQ6xqdHt4SgWt9zwzgniIiITiVkLod1R1FREfLy8lBUVAS32428vDzk5eWhoaEBADB9+nRkZ2djzpw5yM3NxVdffYX7778f8+bN451hfdC4DO+8oC/3VcLVsl5QXCRDEBERdS1kRoK645FHHsE///lP3/djxowBAGzYsAHnn38+5HI5Pv30U9x5552YOnUqtFotbrjhBvzv//5vsKpMPWh8yx1ix+ubAAAxESqoFGGZ74mIqBdJQggR7EqEIovFAqPRCLPZzNGjMDBl0VcoNXsnRWcl6fHF/POCXCMiIgqG7nx+89dl6hNaL4kBnBRNRET+YQiiPqH1khjAEERERP5hCKI+YdwJIYhrBBERkT8YgqhPyErSI0IlB8Db44mIyD8MQdQnKOQyTB7sXWdqcEJkkGtDREThICxvkSfqyOJfjsTesnRMHRyYRTeJiKhvYwiiPiM2Uo1zh8YHuxpERBQmeDmMiIiI+iWOBHWidQ1Ji8US5JoQERGRv1o/t/1ZC5ohqBNWqxUAkJaWFuSaEBERUXdZrVYYjcYuy/CxGZ3weDwoLS2FXq+HJEnBrs5ps1gsSEtLQ3FxcZ98/AfbF/76ehvZvvDX19vY19onhIDVakVycjJksq5n/XAkqBMymQypqanBrkbAGAyGPvGPuzNsX/jr621k+8JfX29jX2rfqUaAWnFiNBEREfVLDEFERETULzEE9XFqtRqPPvoo1Oq++SgJti/89fU2sn3hr6+3sa+3ryucGE1ERET9EkeCiIiIqF9iCCIiIqJ+iSGIiIiI+iWGICIiIuqXGILCyLfffovLL78cycnJkCQJH330UZflv/nmG0iS1G7bv39/m3KrV69GdnY21Go1srOzsWbNmh5sRed6on0rV67ssExzc3MPt6a97rYPAOx2O/70pz8hPT0darUagwcPxvLly9uUCZX+A3qmjeHchzfffHOHdc/JyWlTLlT6sCfaF0r9B5zev9FVq1Zh9OjR0Ol0MJlMuOWWW1BTU9OmTLj2IXDq9oVaHwYSQ1AYaWxsxOjRo/HSSy9167iCggKUlZX5tqFDh/pe27JlC66//nrMmTMHu3fvxpw5czBr1ixs27Yt0NU/pZ5oH+BdBfXE18vKyqDRaAJZdb+cTvtmzZqFr776Cm+88QYKCgrwzjvvICsry/d6KPUf0DNtBMK3D59//vk2dS4uLkZMTAyuu+46X5lQ6sOeaB8QOv0HdL+NmzZtwty5c3HrrbciPz8f77//Pnbs2IHbbrvNVyac+9Cf9gGh1YcBJSgsARBr1qzpssyGDRsEAFFXV9dpmVmzZomLL764zb4ZM2aI2bNnB6CWpy9Q7VuxYoUwGo0BrVsg+NO+zz//XBiNRlFTU9NpmVDtPyEC18Zw7sOTrVmzRkiSJI4dO+bbF6p9GKj2hWr/CeFfG59++mkxaNCgNvteeOEFkZqa6vs+nPvQn/aFch+eKY4E9QNjxoyByWTCBRdcgA0bNrR5bcuWLZg+fXqbfTNmzMDmzZt7s4pnpKv2AUBDQwPS09ORmpqKmTNnIjc3Nwi17L61a9di/PjxWLJkCVJSUjBs2DDcf//9aGpq8pUJ9/7zp41A+Pbhyd544w1ceOGFSE9P9+0L9z48UUftA8K7/6ZMmYKSkhJ89tlnEEKgoqICH3zwAS677DJfmXDuQ3/aB4R3H3aFIagPM5lMWLZsGVavXo0PP/wQmZmZuOCCC/Dtt9/6ypSXlyMxMbHNcYmJiSgvL+/t6nabP+3LysrCypUrsXbtWrzzzjvQaDSYOnUqDh48GMSa++fIkSPYtGkT9uzZgzVr1uC5557DBx98gLvuustXJpz7D/CvjeHchycqKyvD559/3u4yQ7j3YavO2hfu/TdlyhSsWrUK119/PVQqFZKSkhAVFYUXX3zRVyac+9Cf9oV7H3Yp2ENRdHpwGkPVQggxc+ZMcfnll/u+VyqV4l//+lebMm+//bZQq9VnWsUzEqj2ncztdovRo0eLe+655wxqd+b8ad9FF10kNBqNqK+v9+1bvXq1kCRJ2Gw2IUTo9p8QgWvjycKpD0/05JNPitjYWGG329vsD9U+DFT7ThYq/SeEf23Mz88XJpNJLFmyROzevVt88cUXYuTIkeI3v/mNr0w496E/7TtZKPXhmeJIUD8zadKkNuk9KSmp3W8rlZWV7X6rCRcnt+9kMpkMZ599dlj8BmMymZCSkgKj0ejbN3z4cAghUFJSAiD8+8+fNp4snPqwlRACy5cvx5w5c6BSqdq8Fu59CHTdvpOFW/8tWrQIU6dOxQMPPIBRo0ZhxowZWLp0KZYvX46ysjIA4d2H/rTvZOHWh11hCOpncnNzYTKZfN9PnjwZ69evb1Nm3bp1mDJlSm9XLSBObt/JhBDIy8vrskyomDp1KkpLS9HQ0ODbd+DAAchkMqSmpgII//7zp40nC6c+bLVx40YcOnQIt956a7vXwr0Pga7bd7Jw6z+bzQaZrO1HpVwuB+BtCxDefehP+04Wbn3YpWANQVH3Wa1WkZubK3JzcwUA8cwzz4jc3FxRWFgohBDi4YcfFnPmzPGVf/bZZ8WaNWvEgQMHxJ49e8TDDz8sAIjVq1f7ynz//fdCLpeLp556Suzbt0889dRTQqFQiK1bt/aJ9j322GPiiy++EIcPHxa5ubnilltuEQqFQmzbti3k22e1WkVqaqq49tprRX5+vti4caMYOnSouO2223xlQqn/Wusc6DaGcx+2uvHGG8XEiRM7PGco9WFPtC+U+k+I7rdxxYoVQqFQiKVLl4rDhw+LTZs2ifHjx4sJEyb4yoRzH/rTvlDrw0BiCAojrbeEn7zddNNNQgghbrrpJjFt2jRf+cWLF4vBgwcLjUYjoqOjxTnnnCM+/fTTdud9//33RWZmplAqlSIrK6tNiOhNPdG++fPniwEDBgiVSiXi4+PF9OnTxebNm3uxVT/rbvuEEGLfvn3iwgsvFFqtVqSmpooFCxa0mysTKv0nRM+0Mdz7sL6+Xmi1WrFs2bJOzxsqfdgT7Qul/hPi9Nr4wgsviOzsbKHVaoXJZBK//vWvRUlJSZsy4dyHp2pfqPVhIElCdDLeRURERNSHcU4QERER9UsMQURERNQvMQQRERFRv8QQRERERP0SQxARERH1SwxBRERE1C8xBBEREVG/xBBERERE/RJDEBEREfVLDEFERETULzEEERERUb/EEERERET90v8D8xIYfh6wajwAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot\n",
    "\n",
    "print(scan_saddle_point_index(scan_grids, ene_grid))\n",
    "\n",
    "scan_grid, = scan_grids\n",
    "idxs_ = numpy.argsort(scan_grid)\n",
    "x_points = scan_grid[idxs_]\n",
    "ene_grid_ = ene_grid[idxs_]\n",
    "\n",
    "y = interpolate.Akima1DInterpolator(x_points, ene_grid_)\n",
    "x = numpy.linspace(numpy.min(*scan_grids), numpy.max(*scan_grids), 100)\n",
    "\n",
    "fig, ax = pyplot.subplots(3)\n",
    "ax[0].plot(*scan_grids, ene_grid)\n",
    "ax[0].plot(x, y(x))\n",
    "ax[1].plot(x, y(x, nu=1))\n",
    "ax[1].plot(x, numpy.zeros_like(x))\n",
    "ax[2].plot(x, y(x, nu=2))\n",
    "ax[2].plot(x, numpy.zeros_like(x))\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4f48986bec3341f28f91036deeaebcce",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(Image(value=b\"<?xml version='1.0' encoding='iso-8859-1...\", format='svg+xml', height='300', wid…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Starting geometry:\n"
     ]
    },
    {
     "data": {
      "application/3dmoljs_load.v0": "<div id=\"3dmolviewer_17177777597189531\"  style=\"position: relative; width: 400px; height: 400px;\">\n        <p id=\"3dmolwarning_17177777597189531\" style=\"background-color:#ffcccc;color:black\">3Dmol.js failed to load for some reason.  Please check your browser console for error messages.<br></p>\n        </div>\n<script>\n\nvar loadScriptAsync = function(uri){\n  return new Promise((resolve, reject) => {\n    //this is to ignore the existence of requirejs amd\n    var savedexports, savedmodule;\n    if (typeof exports !== 'undefined') savedexports = exports;\n    else exports = {}\n    if (typeof module !== 'undefined') savedmodule = module;\n    else module = {}\n\n    var tag = document.createElement('script');\n    tag.src = uri;\n    tag.async = true;\n    tag.onload = () => {\n        exports = savedexports;\n        module = savedmodule;\n        resolve();\n    };\n  var firstScriptTag = document.getElementsByTagName('script')[0];\n  firstScriptTag.parentNode.insertBefore(tag, firstScriptTag);\n});\n};\n\nif(typeof $3Dmolpromise === 'undefined') {\n$3Dmolpromise = null;\n  $3Dmolpromise = loadScriptAsync('https://cdnjs.cloudflare.com/ajax/libs/3Dmol/2.1.0/3Dmol-min.js');\n}\n\nvar viewer_17177777597189531 = null;\nvar warn = document.getElementById(\"3dmolwarning_17177777597189531\");\nif(warn) {\n    warn.parentNode.removeChild(warn);\n}\n$3Dmolpromise.then(function() {\nviewer_17177777597189531 = $3Dmol.createViewer(document.getElementById(\"3dmolviewer_17177777597189531\"),{backgroundColor:\"white\"});\nviewer_17177777597189531.zoomTo();\n\tviewer_17177777597189531.addModel(\"\\n     RDKit          3D\\n\\n 12 11  0  0  0  0  0  0  0  0999 V2000\\n   -2.1718   -0.9476   -0.2462 C   0  0  0  0  0  0  0  0  0  0  0  0\\n   -0.7518    0.2605    1.9213 C   0  0  0  0  0  0  0  0  0  0  0  0\\n    1.6113    1.5221    0.8229 C   0  0  0  0  0  0  0  0  0  0  0  0\\n    3.5557   -0.3490    0.1467 O   0  0  0  0  0  0  0  0  0  0  0  0\\n    2.3085   -1.9870   -1.6159 O   0  0  0  0  0  1  0  0  0  0  0  0\\n   -1.0376   -2.5986   -1.0336 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -1.9172    0.2229   -1.9822 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -4.0997   -0.8711   -0.0452 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -1.4137    0.9972    3.4696 H   0  0  0  0  0  0  0  0  0  0  0  0\\n    0.1522   -1.5057    2.8302 H   0  0  0  0  0  0  0  0  0  0  0  0\\n    1.2232    2.3351   -1.0605 H   0  0  0  0  0  0  0  0  0  0  0  0\\n    2.5660    2.9075    1.6786 H   0  0  0  0  0  0  0  0  0  0  0  0\\n  1  2  1  0\\n  1  6  1  0\\n  1  7  1  0\\n  1  8  1  0\\n  2  3  1  0\\n  9  2  1  0\\n  2 10  1  0\\n  3  4  1  0\\n  3 11  1  0\\n  3 12  1  0\\n  4  5  1  0\\nM  RAD  1   5   2\\nM  END\\n\",\"sdf\");\n\tviewer_17177777597189531.setStyle({\"stick\": {}, \"sphere\": {\"radius\": 0.3}});\n\tviewer_17177777597189531.setBackgroundColor(\"0xeeeeee\");\n\tviewer_17177777597189531.zoomTo();\nviewer_17177777597189531.render();\n});\n</script>",
      "text/html": [
       "<div id=\"3dmolviewer_17177777597189531\"  style=\"position: relative; width: 400px; height: 400px;\">\n",
       "        <p id=\"3dmolwarning_17177777597189531\" style=\"background-color:#ffcccc;color:black\">3Dmol.js failed to load for some reason.  Please check your browser console for error messages.<br></p>\n",
       "        </div>\n",
       "<script>\n",
       "\n",
       "var loadScriptAsync = function(uri){\n",
       "  return new Promise((resolve, reject) => {\n",
       "    //this is to ignore the existence of requirejs amd\n",
       "    var savedexports, savedmodule;\n",
       "    if (typeof exports !== 'undefined') savedexports = exports;\n",
       "    else exports = {}\n",
       "    if (typeof module !== 'undefined') savedmodule = module;\n",
       "    else module = {}\n",
       "\n",
       "    var tag = document.createElement('script');\n",
       "    tag.src = uri;\n",
       "    tag.async = true;\n",
       "    tag.onload = () => {\n",
       "        exports = savedexports;\n",
       "        module = savedmodule;\n",
       "        resolve();\n",
       "    };\n",
       "  var firstScriptTag = document.getElementsByTagName('script')[0];\n",
       "  firstScriptTag.parentNode.insertBefore(tag, firstScriptTag);\n",
       "});\n",
       "};\n",
       "\n",
       "if(typeof $3Dmolpromise === 'undefined') {\n",
       "$3Dmolpromise = null;\n",
       "  $3Dmolpromise = loadScriptAsync('https://cdnjs.cloudflare.com/ajax/libs/3Dmol/2.1.0/3Dmol-min.js');\n",
       "}\n",
       "\n",
       "var viewer_17177777597189531 = null;\n",
       "var warn = document.getElementById(\"3dmolwarning_17177777597189531\");\n",
       "if(warn) {\n",
       "    warn.parentNode.removeChild(warn);\n",
       "}\n",
       "$3Dmolpromise.then(function() {\n",
       "viewer_17177777597189531 = $3Dmol.createViewer(document.getElementById(\"3dmolviewer_17177777597189531\"),{backgroundColor:\"white\"});\n",
       "viewer_17177777597189531.zoomTo();\n",
       "\tviewer_17177777597189531.addModel(\"\\n     RDKit          3D\\n\\n 12 11  0  0  0  0  0  0  0  0999 V2000\\n   -2.1718   -0.9476   -0.2462 C   0  0  0  0  0  0  0  0  0  0  0  0\\n   -0.7518    0.2605    1.9213 C   0  0  0  0  0  0  0  0  0  0  0  0\\n    1.6113    1.5221    0.8229 C   0  0  0  0  0  0  0  0  0  0  0  0\\n    3.5557   -0.3490    0.1467 O   0  0  0  0  0  0  0  0  0  0  0  0\\n    2.3085   -1.9870   -1.6159 O   0  0  0  0  0  1  0  0  0  0  0  0\\n   -1.0376   -2.5986   -1.0336 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -1.9172    0.2229   -1.9822 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -4.0997   -0.8711   -0.0452 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -1.4137    0.9972    3.4696 H   0  0  0  0  0  0  0  0  0  0  0  0\\n    0.1522   -1.5057    2.8302 H   0  0  0  0  0  0  0  0  0  0  0  0\\n    1.2232    2.3351   -1.0605 H   0  0  0  0  0  0  0  0  0  0  0  0\\n    2.5660    2.9075    1.6786 H   0  0  0  0  0  0  0  0  0  0  0  0\\n  1  2  1  0\\n  1  6  1  0\\n  1  7  1  0\\n  1  8  1  0\\n  2  3  1  0\\n  9  2  1  0\\n  2 10  1  0\\n  3  4  1  0\\n  3 11  1  0\\n  3 12  1  0\\n  4  5  1  0\\nM  RAD  1   5   2\\nM  END\\n\",\"sdf\");\n",
       "\tviewer_17177777597189531.setStyle({\"stick\": {}, \"sphere\": {\"radius\": 0.3}});\n",
       "\tviewer_17177777597189531.setBackgroundColor(\"0xeeeeee\");\n",
       "\tviewer_17177777597189531.zoomTo();\n",
       "viewer_17177777597189531.render();\n",
       "});\n",
       "</script>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Scan coordinates:\n",
      " - (4, 5) distance: 1.826187177856278\n",
      "Constraint coordinates:\n",
      " - (3, 4) distance: 1.434203250708657\n",
      "\n",
      "Scan saddle geometry:\n"
     ]
    },
    {
     "data": {
      "application/3dmoljs_load.v0": "<div id=\"3dmolviewer_17177777597948208\"  style=\"position: relative; width: 400px; height: 400px;\">\n        <p id=\"3dmolwarning_17177777597948208\" style=\"background-color:#ffcccc;color:black\">3Dmol.js failed to load for some reason.  Please check your browser console for error messages.<br></p>\n        </div>\n<script>\n\nvar loadScriptAsync = function(uri){\n  return new Promise((resolve, reject) => {\n    //this is to ignore the existence of requirejs amd\n    var savedexports, savedmodule;\n    if (typeof exports !== 'undefined') savedexports = exports;\n    else exports = {}\n    if (typeof module !== 'undefined') savedmodule = module;\n    else module = {}\n\n    var tag = document.createElement('script');\n    tag.src = uri;\n    tag.async = true;\n    tag.onload = () => {\n        exports = savedexports;\n        module = savedmodule;\n        resolve();\n    };\n  var firstScriptTag = document.getElementsByTagName('script')[0];\n  firstScriptTag.parentNode.insertBefore(tag, firstScriptTag);\n});\n};\n\nif(typeof $3Dmolpromise === 'undefined') {\n$3Dmolpromise = null;\n  $3Dmolpromise = loadScriptAsync('https://cdnjs.cloudflare.com/ajax/libs/3Dmol/2.1.0/3Dmol-min.js');\n}\n\nvar viewer_17177777597948208 = null;\nvar warn = document.getElementById(\"3dmolwarning_17177777597948208\");\nif(warn) {\n    warn.parentNode.removeChild(warn);\n}\n$3Dmolpromise.then(function() {\nviewer_17177777597948208 = $3Dmol.createViewer(document.getElementById(\"3dmolviewer_17177777597948208\"),{backgroundColor:\"white\"});\nviewer_17177777597948208.zoomTo();\n\tviewer_17177777597948208.addModel(\"\\n     RDKit          3D\\n\\n 12 11  0  0  0  0  0  0  0  0999 V2000\\n   -2.0695   -1.0057   -0.3097 C   0  0  0  0  0  0  0  0  0  0  0  0\\n   -0.7059    0.4053    1.7629 C   0  0  0  0  0  0  0  0  0  0  0  0\\n    1.9298    1.2954    0.8213 C   0  0  0  0  0  0  0  0  0  0  0  0\\n    3.2852   -0.8290    0.0683 O   0  0  0  0  0  0  0  0  0  0  0  0\\n    2.3531   -1.7164   -2.3168 O   0  0  0  0  0  1  0  0  0  0  0  0\\n   -0.5385   -2.0001   -1.4173 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -3.0196    0.2565   -1.6270 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -3.3532   -2.4427    0.4045 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -1.7470    2.0896    2.3514 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -0.4207   -0.8171    3.4023 H   0  0  0  0  0  0  0  0  0  0  0  0\\n    1.7066    2.6301   -0.7481 H   0  0  0  0  0  0  0  0  0  0  0  0\\n    2.9929    2.1618    2.3674 H   0  0  0  0  0  0  0  0  0  0  0  0\\n  1  2  1  0\\n  1  6  1  0\\n  1  7  1  0\\n  1  8  1  0\\n  2  3  1  0\\n  9  2  1  0\\n  2 10  1  0\\n  3  4  1  0\\n  3 11  1  0\\n  3 12  1  0\\n  4  5  1  0\\nM  RAD  1   5   2\\nM  END\\n\",\"sdf\");\n\tviewer_17177777597948208.setStyle({\"stick\": {}, \"sphere\": {\"radius\": 0.3}});\n\tviewer_17177777597948208.setBackgroundColor(\"0xeeeeee\");\n\tviewer_17177777597948208.zoomTo();\nviewer_17177777597948208.render();\n});\n</script>",
      "text/html": [
       "<div id=\"3dmolviewer_17177777597948208\"  style=\"position: relative; width: 400px; height: 400px;\">\n",
       "        <p id=\"3dmolwarning_17177777597948208\" style=\"background-color:#ffcccc;color:black\">3Dmol.js failed to load for some reason.  Please check your browser console for error messages.<br></p>\n",
       "        </div>\n",
       "<script>\n",
       "\n",
       "var loadScriptAsync = function(uri){\n",
       "  return new Promise((resolve, reject) => {\n",
       "    //this is to ignore the existence of requirejs amd\n",
       "    var savedexports, savedmodule;\n",
       "    if (typeof exports !== 'undefined') savedexports = exports;\n",
       "    else exports = {}\n",
       "    if (typeof module !== 'undefined') savedmodule = module;\n",
       "    else module = {}\n",
       "\n",
       "    var tag = document.createElement('script');\n",
       "    tag.src = uri;\n",
       "    tag.async = true;\n",
       "    tag.onload = () => {\n",
       "        exports = savedexports;\n",
       "        module = savedmodule;\n",
       "        resolve();\n",
       "    };\n",
       "  var firstScriptTag = document.getElementsByTagName('script')[0];\n",
       "  firstScriptTag.parentNode.insertBefore(tag, firstScriptTag);\n",
       "});\n",
       "};\n",
       "\n",
       "if(typeof $3Dmolpromise === 'undefined') {\n",
       "$3Dmolpromise = null;\n",
       "  $3Dmolpromise = loadScriptAsync('https://cdnjs.cloudflare.com/ajax/libs/3Dmol/2.1.0/3Dmol-min.js');\n",
       "}\n",
       "\n",
       "var viewer_17177777597948208 = null;\n",
       "var warn = document.getElementById(\"3dmolwarning_17177777597948208\");\n",
       "if(warn) {\n",
       "    warn.parentNode.removeChild(warn);\n",
       "}\n",
       "$3Dmolpromise.then(function() {\n",
       "viewer_17177777597948208 = $3Dmol.createViewer(document.getElementById(\"3dmolviewer_17177777597948208\"),{backgroundColor:\"white\"});\n",
       "viewer_17177777597948208.zoomTo();\n",
       "\tviewer_17177777597948208.addModel(\"\\n     RDKit          3D\\n\\n 12 11  0  0  0  0  0  0  0  0999 V2000\\n   -2.0695   -1.0057   -0.3097 C   0  0  0  0  0  0  0  0  0  0  0  0\\n   -0.7059    0.4053    1.7629 C   0  0  0  0  0  0  0  0  0  0  0  0\\n    1.9298    1.2954    0.8213 C   0  0  0  0  0  0  0  0  0  0  0  0\\n    3.2852   -0.8290    0.0683 O   0  0  0  0  0  0  0  0  0  0  0  0\\n    2.3531   -1.7164   -2.3168 O   0  0  0  0  0  1  0  0  0  0  0  0\\n   -0.5385   -2.0001   -1.4173 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -3.0196    0.2565   -1.6270 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -3.3532   -2.4427    0.4045 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -1.7470    2.0896    2.3514 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -0.4207   -0.8171    3.4023 H   0  0  0  0  0  0  0  0  0  0  0  0\\n    1.7066    2.6301   -0.7481 H   0  0  0  0  0  0  0  0  0  0  0  0\\n    2.9929    2.1618    2.3674 H   0  0  0  0  0  0  0  0  0  0  0  0\\n  1  2  1  0\\n  1  6  1  0\\n  1  7  1  0\\n  1  8  1  0\\n  2  3  1  0\\n  9  2  1  0\\n  2 10  1  0\\n  3  4  1  0\\n  3 11  1  0\\n  3 12  1  0\\n  4  5  1  0\\nM  RAD  1   5   2\\nM  END\\n\",\"sdf\");\n",
       "\tviewer_17177777597948208.setStyle({\"stick\": {}, \"sphere\": {\"radius\": 0.3}});\n",
       "\tviewer_17177777597948208.setBackgroundColor(\"0xeeeeee\");\n",
       "\tviewer_17177777597948208.zoomTo();\n",
       "viewer_17177777597948208.render();\n",
       "});\n",
       "</script>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Scan coordinates:\n",
      " - (4, 5) distance: 1.609520654478949\n",
      "Constraint coordinates:\n",
      " - (3, 4) distance: 1.4342001464642709\n",
      "Energy: -525.7724354244447\n",
      "*** ***\n",
      "\n",
      "\n",
      "TS geometry:\n"
     ]
    },
    {
     "data": {
      "application/3dmoljs_load.v0": "<div id=\"3dmolviewer_17177777598378289\"  style=\"position: relative; width: 400px; height: 400px;\">\n        <p id=\"3dmolwarning_17177777598378289\" style=\"background-color:#ffcccc;color:black\">3Dmol.js failed to load for some reason.  Please check your browser console for error messages.<br></p>\n        </div>\n<script>\n\nvar loadScriptAsync = function(uri){\n  return new Promise((resolve, reject) => {\n    //this is to ignore the existence of requirejs amd\n    var savedexports, savedmodule;\n    if (typeof exports !== 'undefined') savedexports = exports;\n    else exports = {}\n    if (typeof module !== 'undefined') savedmodule = module;\n    else module = {}\n\n    var tag = document.createElement('script');\n    tag.src = uri;\n    tag.async = true;\n    tag.onload = () => {\n        exports = savedexports;\n        module = savedmodule;\n        resolve();\n    };\n  var firstScriptTag = document.getElementsByTagName('script')[0];\n  firstScriptTag.parentNode.insertBefore(tag, firstScriptTag);\n});\n};\n\nif(typeof $3Dmolpromise === 'undefined') {\n$3Dmolpromise = null;\n  $3Dmolpromise = loadScriptAsync('https://cdnjs.cloudflare.com/ajax/libs/3Dmol/2.1.0/3Dmol-min.js');\n}\n\nvar viewer_17177777598378289 = null;\nvar warn = document.getElementById(\"3dmolwarning_17177777598378289\");\nif(warn) {\n    warn.parentNode.removeChild(warn);\n}\n$3Dmolpromise.then(function() {\nviewer_17177777598378289 = $3Dmol.createViewer(document.getElementById(\"3dmolviewer_17177777598378289\"),{backgroundColor:\"white\"});\nviewer_17177777598378289.zoomTo();\n\tviewer_17177777598378289.addModel(\"\\n     RDKit          3D\\n\\n 12 11  0  0  0  0  0  0  0  0999 V2000\\n   -2.0640   -0.9866   -0.2923 C   0  0  0  0  0  0  0  0  0  0  0  0\\n   -0.7299    0.4503    1.7712 C   0  0  0  0  0  0  0  0  0  0  0  0\\n    1.9205    1.2808    0.8229 C   0  0  0  0  0  0  0  0  0  0  0  0\\n    3.1745   -0.9088   -0.0015 O   0  0  0  0  0  0  0  0  0  0  0  0\\n    2.1195   -1.7502   -2.2024 O   0  0  0  0  0  1  0  0  0  0  0  0\\n   -0.2829   -1.7960   -1.4706 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -3.0953    0.2058   -1.6094 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -3.1794   -2.5834    0.3515 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -1.7687    2.1483    2.3260 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -0.4634   -0.7521    3.4282 H   0  0  0  0  0  0  0  0  0  0  0  0\\n    1.7382    2.6456   -0.7267 H   0  0  0  0  0  0  0  0  0  0  0  0\\n    3.0441    2.0739    2.3622 H   0  0  0  0  0  0  0  0  0  0  0  0\\n  1  2  1  0\\n  1  6  1  0\\n  1  7  1  0\\n  1  8  1  0\\n  2  3  1  0\\n  9  2  1  0\\n  2 10  1  0\\n  3  4  1  0\\n  3 11  1  0\\n  3 12  1  0\\n  4  5  1  0\\nM  RAD  1   5   2\\nM  END\\n\",\"sdf\");\n\tviewer_17177777598378289.setStyle({\"stick\": {}, \"sphere\": {\"radius\": 0.3}});\n\tviewer_17177777598378289.setBackgroundColor(\"0xeeeeee\");\n\tviewer_17177777598378289.zoomTo();\nviewer_17177777598378289.render();\n});\n</script>",
      "text/html": [
       "<div id=\"3dmolviewer_17177777598378289\"  style=\"position: relative; width: 400px; height: 400px;\">\n",
       "        <p id=\"3dmolwarning_17177777598378289\" style=\"background-color:#ffcccc;color:black\">3Dmol.js failed to load for some reason.  Please check your browser console for error messages.<br></p>\n",
       "        </div>\n",
       "<script>\n",
       "\n",
       "var loadScriptAsync = function(uri){\n",
       "  return new Promise((resolve, reject) => {\n",
       "    //this is to ignore the existence of requirejs amd\n",
       "    var savedexports, savedmodule;\n",
       "    if (typeof exports !== 'undefined') savedexports = exports;\n",
       "    else exports = {}\n",
       "    if (typeof module !== 'undefined') savedmodule = module;\n",
       "    else module = {}\n",
       "\n",
       "    var tag = document.createElement('script');\n",
       "    tag.src = uri;\n",
       "    tag.async = true;\n",
       "    tag.onload = () => {\n",
       "        exports = savedexports;\n",
       "        module = savedmodule;\n",
       "        resolve();\n",
       "    };\n",
       "  var firstScriptTag = document.getElementsByTagName('script')[0];\n",
       "  firstScriptTag.parentNode.insertBefore(tag, firstScriptTag);\n",
       "});\n",
       "};\n",
       "\n",
       "if(typeof $3Dmolpromise === 'undefined') {\n",
       "$3Dmolpromise = null;\n",
       "  $3Dmolpromise = loadScriptAsync('https://cdnjs.cloudflare.com/ajax/libs/3Dmol/2.1.0/3Dmol-min.js');\n",
       "}\n",
       "\n",
       "var viewer_17177777598378289 = null;\n",
       "var warn = document.getElementById(\"3dmolwarning_17177777598378289\");\n",
       "if(warn) {\n",
       "    warn.parentNode.removeChild(warn);\n",
       "}\n",
       "$3Dmolpromise.then(function() {\n",
       "viewer_17177777598378289 = $3Dmol.createViewer(document.getElementById(\"3dmolviewer_17177777598378289\"),{backgroundColor:\"white\"});\n",
       "viewer_17177777598378289.zoomTo();\n",
       "\tviewer_17177777598378289.addModel(\"\\n     RDKit          3D\\n\\n 12 11  0  0  0  0  0  0  0  0999 V2000\\n   -2.0640   -0.9866   -0.2923 C   0  0  0  0  0  0  0  0  0  0  0  0\\n   -0.7299    0.4503    1.7712 C   0  0  0  0  0  0  0  0  0  0  0  0\\n    1.9205    1.2808    0.8229 C   0  0  0  0  0  0  0  0  0  0  0  0\\n    3.1745   -0.9088   -0.0015 O   0  0  0  0  0  0  0  0  0  0  0  0\\n    2.1195   -1.7502   -2.2024 O   0  0  0  0  0  1  0  0  0  0  0  0\\n   -0.2829   -1.7960   -1.4706 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -3.0953    0.2058   -1.6094 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -3.1794   -2.5834    0.3515 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -1.7687    2.1483    2.3260 H   0  0  0  0  0  0  0  0  0  0  0  0\\n   -0.4634   -0.7521    3.4282 H   0  0  0  0  0  0  0  0  0  0  0  0\\n    1.7382    2.6456   -0.7267 H   0  0  0  0  0  0  0  0  0  0  0  0\\n    3.0441    2.0739    2.3622 H   0  0  0  0  0  0  0  0  0  0  0  0\\n  1  2  1  0\\n  1  6  1  0\\n  1  7  1  0\\n  1  8  1  0\\n  2  3  1  0\\n  9  2  1  0\\n  2 10  1  0\\n  3  4  1  0\\n  3 11  1  0\\n  3 12  1  0\\n  4  5  1  0\\nM  RAD  1   5   2\\nM  END\\n\",\"sdf\");\n",
       "\tviewer_17177777598378289.setStyle({\"stick\": {}, \"sphere\": {\"radius\": 0.3}});\n",
       "\tviewer_17177777598378289.setBackgroundColor(\"0xeeeeee\");\n",
       "\tviewer_17177777598378289.zoomTo();\n",
       "viewer_17177777598378289.render();\n",
       "});\n",
       "</script>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Scan coordinates:\n",
      " - (4, 5) distance: 1.329220053597468\n",
      "Constraint coordinates:\n",
      " - (3, 4) distance: 1.3661323494287925\n"
     ]
    }
   ],
   "source": [
    "# Visualize the results\n",
    "def log_geometry_info(geo, scan_coos, const_coos, gra=None):\n",
    "    \"\"\"Print and display information about a geometry\n",
    "\n",
    "    :param geo: An automol geometry data structure\n",
    "    :param coo: The reaction/scan coordinate\n",
    "    :param const_coos: Any constrained coordinates\n",
    "    \"\"\"\n",
    "    # Display the geometry\n",
    "    automol.geom.display(geo, gra=None)\n",
    "\n",
    "    # Print distances for the scan coordinates\n",
    "    print(\"Scan coordinates:\")\n",
    "    for coo in scan_coos:\n",
    "        dist = automol.geom.distance(geo, *coo, angstrom=True)\n",
    "        print(f\" - {coo} distance: {dist}\")\n",
    "\n",
    "    # Print distances for the constraint coordinates\n",
    "    print(\"Constraint coordinates:\")\n",
    "    for coo in const_coos:\n",
    "        dist = automol.geom.distance(geo, *coo, angstrom=True)\n",
    "        print(f\" - {coo} distance: {dist}\")\n",
    "\n",
    "\n",
    "automol.graph.display(ts_gra, label=True, exp=True)\n",
    "\n",
    "print(\"\\nStarting geometry:\")\n",
    "log_geometry_info(geo0, scan_coos, const_coos, gra=ts_gra)\n",
    "\n",
    "print(\"\\nScan saddle geometry:\")\n",
    "ene = ene_grid[saddle_idxs]\n",
    "log_geometry_info(ts_geo0, scan_coos, const_coos, gra=ts_gra)\n",
    "print(\"Energy:\", ene)\n",
    "print(\"*** ***\\n\")\n",
    "\n",
    "print(\"\\nTS geometry:\")\n",
    "log_geometry_info(ts_geo, scan_coos, const_coos, gra=ts_gra)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}