mol_gen_test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This notebook is licensed under the GNU General Public License, Version 2 (GPL-2.0).\n",
    "#See the `LICENSE_GPL` file in the root directory of this project for the full text.\n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from molecular_generation_utils import *\n",
    "from invert_CM import *\n",
    "import torch\n",
    "from Model import Multi_VAE\n",
    "from torch.distributions.categorical import Categorical\n",
    "from tqdm import tqdm\n",
    "import copy\n",
    "\n",
    "reproduce_paper = True\n",
    "\n",
    "if reproduce_paper:\n",
    "    paper_path = '_paper'\n",
    "else:\n",
    "    paper_path = ''\n",
    "\n",
    "\n",
    "properties = torch.load('./data{}/data/properties_total.pt'.format(paper_path))\n",
    "p_means = torch.load('./data{}/data/properties_means.pt'.format(paper_path))\n",
    "p_stds = torch.load('./data{}/data/properties_stds.pt'.format(paper_path))\n",
    "norm_props = (properties - p_means)/p_stds\n",
    "\n",
    "properties_list =  ['eAT', 'eMBD', 'eXX', 'mPOL', 'eNN', 'eNE', 'eEE', 'eKIN', 'DIP', 'HLgap', 'HOMO_0', 'LUMO_0', 'HOMO_1', 'LUMO_1', 'HOMO_2', 'LUMO_2', 'dimension']\n",
    "p_arr = np.array(properties_list)\n",
    "\n",
    "PATH = './models_saved/masked/epoch=2597-step=145487.ckpt'\n",
    "\n",
    "modello = Multi_VAE.load_from_checkpoint(\n",
    "    PATH,\n",
    "    map_location=torch.device('cpu'),\n",
    "    structures_dim = len(torch.load('./data/data_val/CMs.pt')[0,:]),\n",
    "    properties_dim = len(torch.load('./data/data_val/properties.pt')[0,:]),\n",
    "    latent_size = 21,\n",
    "    extra_dim = 32 - len(torch.load('./data/data_val/properties.pt')[0,:]),\n",
    "    initial_lr = 1e-3,\n",
    "    properties_means = p_means,\n",
    "    properties_stds = p_stds,\n",
    "    beta = 4.,\n",
    "    alpha = 1.,\n",
    "    decay = .01\n",
    ")\n",
    "\n",
    "# if you want the non-masked model uncomment this\n",
    "\n",
    "# if reproduce_paper == True:\n",
    "#     PATH='./special/VAE_reduced_21'\n",
    "#     modello.VAE.load_state_dict(torch.load(PATH,map_location=torch.device('cpu')))\n",
    "#     PATH='./special/prop_ecoder_reduced_21'\n",
    "#     modello.property_encoder.load_state_dict(torch.load(PATH,map_location=torch.device('cpu')))\n",
    "\n",
    "modello.eval()\n",
    "modello.freeze()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit the dataset property distribution with 91 (optional number) Gaussians"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 1/1 [02:18<00:00, 138.54s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "using 91 components\n"
     ]
    }
   ],
   "source": [
    "gm, labels = props_fit_Gaussian_mix(\n",
    "    norm_props, \n",
    "    min_components = 91,\n",
    "    max_components = 92\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# recall that here the target value is the normalized one\n",
    "\n",
    "generated = start_generation(\n",
    "    modello,\n",
    "    {\n",
    "        'mPOL': 2.,\n",
    "        'eMBD': 2.\n",
    "    },\n",
    "    p_arr,\n",
    "    11,\n",
    "    int(5e3),\n",
    "    gm.means_,\n",
    "    gm.covariances_,\n",
    "    cm_diff = 5,\n",
    "    deltaz = 6,\n",
    "    check_new_comp = True\n",
    "    verbose = True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the coulomb matrices in the dataset\n",
    "CMs = torch.load('./data{}/data/CMs_total.pt'.format(paper_path))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [],
   "source": [
    "#find the closest mmatrix in the dataset\n",
    "\n",
    "where_closest = []\n",
    "for cm in generated:\n",
    "    where_closest.append(torch.argmin(torch.norm(CMs - cm.view(1,-1), dim = 1)))\n",
    "    #print(torch.norm(CMs - cm.view(1,-1), dim = 1).min())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "index 0 is out of bounds for dimension 0 with size 0",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-32-4010e6ec26b6>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0mn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m \u001b[0mpos\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcomp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_cartesian\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgenerated\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtolist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     13\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcomp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m \u001b[0matom\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mAtoms\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcomp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpos\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: index 0 is out of bounds for dimension 0 with size 0"
     ]
    }
   ],
   "source": [
    "#some visualization and optimization\n",
    "\n",
    "import rmsd\n",
    "from ase.visualize import view\n",
    "from ase.optimize import BFGS, FIRE\n",
    "from ase import Atoms\n",
    "from ase.io import write, read\n",
    "from openbabel import openbabel as ob\n",
    "from ase.constraints import FixAtoms\n",
    "\n",
    "# here choose the calculator you prefer, or do not optimize by commenting out the opt steps in the code\n",
    "#from spooks.spookynet import SpookyNetCalculator, spookynet\n",
    "#calc=SpookyNetCalculator(load_from=\"./spooks/qm7x_known_unknown_0.pth\", charge=0, magmom=0)\n",
    "\n",
    "n = 0\n",
    "pos, comp = get_cartesian(generated[n,:].tolist())\n",
    "print(len(comp))\n",
    "atom = Atoms(comp, pos)\n",
    "write('./temp.xyz', atom)\n",
    "obConversion = ob.OBConversion()\n",
    "obConversion.SetInAndOutFormats(\"xyz\", \"mol2\")\n",
    "mol = ob.OBMol()\n",
    "obConversion.ReadFile(mol, \"temp.xyz\")  \n",
    "obConversion.WriteFile(mol, 'temp.mol2')\n",
    "obConversion = ob.OBConversion()\n",
    "obConversion.SetInAndOutFormats(\"mol2\", \"xyz\")\n",
    "mol = ob.OBMol()\n",
    "obConversion.ReadFile(mol, \"temp.mol2\")  \n",
    "mol.ConnectTheDots()\n",
    "mol.PerceiveBondOrders()\n",
    "mol.AddHydrogens()\n",
    "obConversion.WriteFile(mol, './temp.xyz')\n",
    "atom = read('./temp.xyz')\n",
    "print(comp, len(comp))\n",
    "c = FixAtoms(indices=[atomo.index for atomo in atom if atomo.symbol != 'H'])\n",
    "atom.set_constraint(c)\n",
    "atom.set_calculator(calc)\n",
    "#try:\n",
    "opt = FIRE(atom)\n",
    "opt.run(steps = 5)\n",
    "# atom.set_constraint()\n",
    "# opt.run(steps = 10)\n",
    "#except:\n",
    "#    pass\n",
    "view(atom, viewer = 'x3d')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<html>\n",
       "\n",
       " <head>\n",
       "\n",
       "  <title>ASE atomic visualization</title>\n",
       "\n",
       "  <link rel=\"stylesheet\" type=\"text/css\"\n",
       "\n",
       "   href=\"https://www.x3dom.org/x3dom/release/x3dom.css\">\n",
       "\n",
       "  </link>\n",
       "\n",
       "  <script type=\"text/javascript\"\n",
       "\n",
       "   src=\"https://www.x3dom.org/x3dom/release/x3dom.js\">\n",
       "\n",
       "  </script>\n",
       "\n",
       " </head>\n",
       "\n",
       " <body>\n",
       "\n",
       "  <X3D>\n",
       "\n",
       "   <Scene>\n",
       "\n",
       "    <Transform translation=\"0.00 0.00 0.00\">\n",
       "\n",
       "     <Shape>\n",
       "\n",
       "      <Appearance>\n",
       "\n",
       "       <Material diffuseColor=\"0.188 0.314 0.973\" specularColor=\"0.5 0.5 0.5\">\n",
       "\n",
       "       </Material>\n",
       "\n",
       "      </Appearance>\n",
       "\n",
       "      <Sphere radius=\"0.71\">\n",
       "\n",
       "      </Sphere>\n",
       "\n",
       "     </Shape>\n",
       "\n",
       "    </Transform>\n",
       "\n",
       "    <Transform translation=\"2.79 0.05 1.02\">\n",
       "\n",
       "     <Shape>\n",
       "\n",
       "      <Appearance>\n",
       "\n",
       "       <Material diffuseColor=\"0.188 0.314 0.973\" specularColor=\"0.5 0.5 0.5\">\n",
       "\n",
       "       </Material>\n",
       "\n",
       "      </Appearance>\n",
       "\n",
       "      <Sphere radius=\"0.71\">\n",
       "\n",
       "      </Sphere>\n",
       "\n",
       "     </Shape>\n",
       "\n",
       "    </Transform>\n",
       "\n",
       "    <Transform translation=\"0.74 -1.25 0.23\">\n",
       "\n",
       "     <Shape>\n",
       "\n",
       "      <Appearance>\n",
       "\n",
       "       <Material diffuseColor=\"0.565 0.565 0.565\" specularColor=\"0.5 0.5 0.5\">\n",
       "\n",
       "       </Material>\n",
       "\n",
       "      </Appearance>\n",
       "\n",
       "      <Sphere radius=\"0.76\">\n",
       "\n",
       "      </Sphere>\n",
       "\n",
       "     </Shape>\n",
       "\n",
       "    </Transform>\n",
       "\n",
       "    <Transform translation=\"2.28 -1.08 0.19\">\n",
       "\n",
       "     <Shape>\n",
       "\n",
       "      <Appearance>\n",
       "\n",
       "       <Material diffuseColor=\"0.565 0.565 0.565\" specularColor=\"0.5 0.5 0.5\">\n",
       "\n",
       "       </Material>\n",
       "\n",
       "      </Appearance>\n",
       "\n",
       "      <Sphere radius=\"0.76\">\n",
       "\n",
       "      </Sphere>\n",
       "\n",
       "     </Shape>\n",
       "\n",
       "    </Transform>\n",
       "\n",
       "    <Transform translation=\"2.92 -0.46 -1.07\">\n",
       "\n",
       "     <Shape>\n",
       "\n",
       "      <Appearance>\n",
       "\n",
       "       <Material diffuseColor=\"0.565 0.565 0.565\" specularColor=\"0.5 0.5 0.5\">\n",
       "\n",
       "       </Material>\n",
       "\n",
       "      </Appearance>\n",
       "\n",
       "      <Sphere radius=\"0.76\">\n",
       "\n",
       "      </Sphere>\n",
       "\n",
       "     </Shape>\n",
       "\n",
       "    </Transform>\n",
       "\n",
       "    <Transform translation=\"3.57 0.58 -0.14\">\n",
       "\n",
       "     <Shape>\n",
       "\n",
       "      <Appearance>\n",
       "\n",
       "       <Material diffuseColor=\"0.565 0.565 0.565\" specularColor=\"0.5 0.5 0.5\">\n",
       "\n",
       "       </Material>\n",
       "\n",
       "      </Appearance>\n",
       "\n",
       "      <Sphere radius=\"0.76\">\n",
       "\n",
       "      </Sphere>\n",
       "\n",
       "     </Shape>\n",
       "\n",
       "    </Transform>\n",
       "\n",
       "    <Transform translation=\"5.08 0.49 0.03\">\n",
       "\n",
       "     <Shape>\n",
       "\n",
       "      <Appearance>\n",
       "\n",
       "       <Material diffuseColor=\"0.565 0.565 0.565\" specularColor=\"0.5 0.5 0.5\">\n",
       "\n",
       "       </Material>\n",
       "\n",
       "      </Appearance>\n",
       "\n",
       "      <Sphere radius=\"0.76\">\n",
       "\n",
       "      </Sphere>\n",
       "\n",
       "     </Shape>\n",
       "\n",
       "    </Transform>\n",
       "\n",
       "   </Scene>\n",
       "\n",
       "  </X3D>\n",
       "\n",
       " </body>\n",
       "\n",
       "</html>\n",
       "\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# visualize the closest molecules in the dataset by coulomb matrix\n",
    "\n",
    "import rmsd\n",
    "from ase.visualize import view\n",
    "from ase import Atoms\n",
    "\n",
    "pos, comp = get_cartesian(CMs[where_closest[n],:].tolist())\n",
    "print(len(comp))\n",
    "atom = Atoms(comp, pos)\n",
    "view(atom, viewer = 'x3d')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[['opt-new_mol_1                  89.152023       -0.349659      101.345923        6.679409        0.258495']\n",
      " ['opt-new_mol_2                  83.942242       -0.275254      103.457106        6.612292        0.490576']\n",
      " ['opt-new_mol_3                  81.672425       -0.320542       96.750483        6.537104        0.130669']\n",
      " ['opt-new_mol_6                  87.693071       -0.300381      100.507797        6.651634        0.531671']]\n",
      "[['opt-new_mol_1-pp               97.006998       -0.433016      106.544639        6.457454        0.514589']\n",
      " ['opt-new_mol_5-pp               91.010867       -0.405348      105.686245        7.151042        0.180508']]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAikAAAFLCAYAAAATRxzSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABmtElEQVR4nO29e5wcZZ3v/366p5PMJUlPNmEuPQQIkEAYE4xRkEsAlQ2Kcl+DCILuklXC1YCyiMKPPbuuN0CRVRMFZInKOSxyRA4KKARRWAwh6CQwBJIAc2EkZCowM7n0dD+/P7qr09NTVX2r6q7q+b5fr3n1dHVdnqqurvrU96q01giCIAiCIPiNULUHIAiCIAiCYIWIFEEQBEEQfImIFEEQBEEQfImIFEEQBEEQfImIFEEQBEEQfImIFEEQBEEQfImIFEEQBEEQfIlvRIpS6hyl1G1KqT8opd5RSmml1D028x6qlPqyUur3Sqk3lFJ7lVIDSqn/q5Q6Kc92LlRKPauUGlJK7VRKPaGU+rg3eyUIgiAIQqkovxRzU0ptABYCQ0APcBiwRmt9vsW8vwCWAZuAp4AdwDzgNCAMXKG1/p7Fct8GVqbXfx8wCTgXmAFcprX+vus7JgiCIAhCSfhJpJxESjy8ApwAPI69SLkIeEFr/XzO9BOARwENHKi17s/67Bjgj8CrwPu11oPp6QcCzwGNwGFa621u75sgCIIgCMXjG3eP1vpxrfVmXYBq0lrflStQ0tPXAk+QspAck/Px59Ov/2YKlPQy24DbgcnAZ0sbvSAIgiAIbuMbkeIi8fTraM70D6Vff2OxzMM58wiCIAiCUGVqSqQopQ4APgyMAE9mTW8EYsBQtgsoi83p17meD1IQBEEQhIKoq/YA3EIpNRlYQ8pt86Vslw4wPf2602Zxc3rUm9EJgiAIglAsNSFSlFJh4L+AY4F7gW97sI3lwHKAxsbG9x122GFub0IQBKHm2bVrl+1n9fX1FRxJ8cTjcUZHcyMJrPH7vhTDc889t11rPasa2w68SEkLlHuAfwD+N3C+RfCtaSmZjjXmdMNuO1rrVcAqgMWLF+t169aVOmRBEIQJS3d3N/F4fNz0SCTCvHnzqjCi0qmlfbGjt7eXv//7v59Zre0HOiZFKRUBfk6q1snPgPO01uNkrtZ6GOgFmpRSbRarOjT9+rJXYxUEQRCgpaUFpdSYaUopWlpaqjSi0vHjvhiGQXd3N11dXXR3d2MYRsnr6u3tZXBwMP+MHhJYS4pSahIpy8npwN3AZ7XWSYdFfg9cAJwC3Jnz2Uez5hEEQRA8IhqNAjAwMEA8HicSidDS0pKZHiRK2ZfcG39zczOxWMyV8RiGQW9vL6YzIR6P09vbO2asxVBtgQIBFSnpINn7gY8BPwGW5xEoAD8kJVK+opR6IKeY2wpgD+PFiyAIgpAHwzCKulFHo9FAihIritkXK8uE+d4NoTIwMEButIPWmoGBgcAeb9+IFKXUGcAZ6bet6dcPKqXuSv+/XWt9dfr/H5ISKNtJuXG+lmtyA57QWj9hvtFa/0kpdTPwReAvSimzLP4y9pXF3+beHgmCINQ+bj+91zJ2lonBwUFXRIpVfIzT9CDgG5ECHAlcmDNtTvoP4DXAFCkHpV9nAl9zWOcT2W+01iuVUn8lZTlZDiSB9cC3tNa/LnXggiAIE5VafHovlGItSF4TiURsA3lLobm5ueouH9+IFK31jcCNBc57YhnbuQu4q9TlBUEQhH3U4tN7IeS6buLxOD09PUD1LEgtLS1jrFpQXiCvW7Ey5RDo7B5BEAShutg9pZf69B4EDMOwtTD09fXZLtfc3FzU9GKJRqPEYrHMsY9EIsRisbJEUywWY9OmTc+5MsAS8I0lRRAEQQgebj+9BwEnIZJM7svhsHIHAZ5l90BtBSWDiBRBEAShDGoppbhQsoWIHXYBxbFYzBdulKAgIkUQBEEoi1p7ei+HcDgMTOyAYjeRmBRBEARBKAJTiFjR1pYqaj5RA4rdRiwpgiAIglAEbW1t4+JwIBVfYlpJ3E4HriS5sTTNzc0zqjUWsaQIgiAIQhFYubdyA2D92NenEAzDoKenJyOw4vE4ra2tB1ZrPGJJEQRBEDJ42VsmCBSy/4ZhjGvcZxgGjY2NGfES1IBiq8wlZVHSvVKISBEEQRAA73vL+J1C998uKLanp4eenp4xgsSvosSuWm4hmUuVRNw9giAIAuDcW2YiUOj+5wt+NdONc60tfsFMj8526fh1vGJJEQRBEDyl1lxIdkGx2fg53djJEqSUGvdZNRGRIgiCIHhGEFxIpuujUKyq7Frhh3RjK4FY7Lh0FVWLuHsEQRAEwJveMn53IeW6PqzI3f/cHjl2VDvd2Ekg2qG1JhwOj+n/8+abb27zaoz5EEuKIAiCAOyzbNSSayYfVq6PbOz23wyK7erqsl222unGToLEya2TSCQ4/PDDs9ezw/XBFYiIFEEQBCHDROst42RB6ezszLu8XXxKKBTyZTyKSSwWo6enx/KzXAuQFHMTBEEQahIvXEhuYueScSp9n42dtaS9vb3kMVWCaDRKR0dHQQXnZs2aVTXVKpYUQRAEwTP87kKyC4JNJBIYhsHw8HDVxm5Xy6RQmpubLV0+DQ0NdHd3E4/HCYfDaK1JJpNEIhGampoYGBgYU++lrq5ukou7VRQiUgRBEISSKSS92M8upGg0Sl9fn2V8hpU7JDczyapCqzm9HHePGdBrjsusZWKOuRCsBGJDQwO7du3KrDeRSKCUoqOjA8Bym8lkcrTkHSkTESmCIAhCSQQhvTgbO0FVbJXVwcHBzP7ZLVtu5Va7Wia5tVfyicRcgdjd3W273kQiYflZNRGRIgiCIBRFthvCiuybuF9wElSFFGerNHbjyZ5eikgsZL25hEKhqmkFCZwVBEEQCqaQuiJ+xKlei1XH4kKxC7AtNPDWDruA3uzppdSgKaV2y+jo6N6iF3IJESmCIAg1hGEYdHd309XVRXd3t+v9WPLVFQkiZnG2YjCPbVtbm2WGTFtbW1ljshJOVpk3bqw3H2+99VZvWRstAxEpgiAINUIlGscVYkHxS3pxMRQb5JodyJpdfTYSiRCLxcqukZJb1dbt9RZKOByWYm6CIAhBptxUUbfo7++3DHzs7+93bTz54jf8lF6cjV06bragsts3UyjkfmYGnM6bN8+T79usamtHIftULolEgrlz577HtRUWiYgUQZjg+OUGGxRyj1dTUxOGYZSVKuoWiUSiqOnFsGXLFkZGRiw/U0q58pTvJYXUa7GqmWK6WOyqs1YzNqfUGjTFNFMEpE6KIAjVwY1aDBMJq+Nl9SRrlSoaZJwESpCEbb56LeY+WIl2u2ymajcRLKUGTZCCnkWkCMIEptBaDEKKYoJGq3EjCIVClvU5QqHywg/tBArAvHnzylq337BzsThZWbzEC0unH1Ou7ZDAWUGYwJRSM2EiU8xxqcYTtl2/mHL6yOQLujUtb7WCXXaUV4Gs+cbiRSC0XeZQR0cHnZ2dVbcOZSOWFEGYwOQLFBTGUswTaFNTk8ejGY+Tu6IQrJ7a88UvOBUP82u8k12VVsMwxsSexOPxzHvTwlLJ8Rdi6SzlGOc7T+z6GVUDESmCMIGplgk7qBRz8R4aGvJ8PHY3qEJupLk36tyeLuZTeyH7alVh1m/xTk5Vcs3jYFcEraenpypjzmfpLOcYO50nuSKmmsXcRKQIwgSm3CfviUbu8XLCa5dZOTcoq3LqVnEnxTxJG4YxZrt+infKPVZWOFVprRb5LJ1eHuNsEfPyyy//tayVlYGIFEGY4FTahB10so+X2e7eCq9dZuXcoLy4IecKJD/FO1nVjymVQro+u0U+S6efjrFXiEgRBEEoESf3j9cus0rdoEzrWj5LRK5AqmS8k1NchmEYrtSJgcp3fc5n6bQ7xuFwOCOgg24dFZEiCIJQIuaFv7+/P3MjDIVCtLe3e35TqIQIMJ/ao9EofX19ea0R2eOpVLxTPrdXX1+fa9txaujnlTXFydJpd4wTiUTmfKx2LFC5iEgRBEEog2q5y8oRAXbl1LPJfgLv7e21rL9itYxJpeKdnNxeQEHjLgS7irPVxOoYJxIJy+PR09NDT09P4CwrIlIEQZiwWJnv/dp7JpdyRIBVOfVsIpFIpkibkwUhGyuBVAkB5+T2KrT8u1LKF+m2pZB7jLu6uhznD5plRUSKIAgTDsMw6Ovrs3zK9jLGwG3KEQGxWIzGxsa81phCbvTVfDp3issoJD7HDYHip67PhdTyCVJVaREpgiBMKApNRw2CSMnGyioEzpahQqwx+W54HR0dVb3Z2QX1FhosW65ACYVCvjpXmpqaCrJ8BSUDSESKIAgThtyKorWCnUCB/JahfNaYcDjseMPPrsjqJU6F66r5nboV8+IWhRYRDEpVaREpglAj+LUEeTWwOhZQe31mTPI9OZdjGSrE0tDX1+fJuWZXJbYScRV2zRpz8dvNvlAXV1CqSotIEYSAYxjGmBRYCF5wnJvY9V8Jh8MFm/b9FGPgFYWK2kJu1F5YE/K55dyMq8iNS1FKZZoymscoFAqhtfZ1C4lCRHjQHmBEpAhCgHG6kAcpOM5N7OpiFBqjUMnsnkpav5xcQn4UtVapxbnE43G6u7vL2k52I0Wr7yH7ePjZWun0/SqlPO/Y7BUiUgQhwOS7kAclOM5Nyn2qr6RAcasBX766J5FIJK9LyErUGoZR1DjcpNBz124+pRTRaJSdO3fanhPZxeoKOeZ+biHh9P0GVaCAiBRBCDT5LuR+85dXGz/Vwyik906hT+5OdU8KKdxmkns+FVpnJBQKFTRfMRSSSuu0rHmsSrl5V7I/TzGUaskJqkABESmCEGicLuR+85d7iV2AZTbhcJhp06Y53rQaGhq8GJ4l+Xrv2MXWgPVNp7GxEcMwxsVMNDY2FixSckVtoSLBjN9wE6d+QfmCWs1CdKVQ6f48hWJneRseHi44oyeIuC9/BUGoGC0tLSilxk03azcE+QmqUMybeb4baltbGzt37nScZ9euXRVzcdhZuczpdrE1dtPtLDPFpOfmitpqWuKi0SixWCwzhkgkQkdHBx0dHY7WMLO5XldXV954ld7e3nHfd75U7mph9/0ODg46nvuhUKiqbrtyEUuKIASYSvVH8SOFWE9MTCGXL17FvKn39/fT1tbm6XHM13vHbqzJZBLDMDI9ddy6eYZCoXH7W0j3Y8CzAO3cGJB8dW6smus5YX7fAwMDvv/dlOr6SiaTvguKLgYRKYIQcPwczOcVhVSNzSa74VwhJBIJzy/s5QjM3t7eggVaodi5bAo5xpUI0C5EkJUac5QdtOxXyonRCXKmn4gUQfARfk5x9BOlVBiNx+N5q6dmU4kLu5PAdBqr1tp1YdDT08Pw8PCYuIv+/v6ClvXaLVRIk8NybuKw7/u2CzSudu2cQq1adgQ10883MSlKqXOUUrcppf6glHpHKaWVUvfYzBtRSl2hlLpTKbVBKbU3Pf8/FbCdC5VSzyqlhpRSO5VSTyilPu7+HglCcZjWAfNiYj7dBdmf7AWlPvFGIhHa2tosY3jsqOaFvdixusHg4OCY41uIoKtEgHYhVrCmpqayxVI8HicWi40TJH7I7rGK0SmGoGb6+cmScj2wEBgCeoDDHOZtBG5N/z8AvAnsn28DSqlvAyvT618NTALOBR5USl2mtf5+qYMXBCcKsZAUkpI60bA6bqXGYGQf8+wqok5xKtUOHIXSrEblUEwJ/UpZ+woRi2acTm6GExTuBjK/71gsVnVRYkWu5a2rq6vgZROJROYYBQk/iZSrSImHV4ATgMcd5h0BPgZs0Fr3K6VuBG5wWrlS6hhSAuVV4P1a68H09G8BzwHfVkr9Wmu9rcz9EIQxFFq0K19Kaq2RT7g5HbdiaW5uHlNBtJBgzGqkcFsdk46OjrLM/OVgJ+JCoRDz58+v2DgKceWYmS7Nzc0MDQ1l3Hta6zH7EIlEaGpqskzXDlLKvpOFNRQKZYKITYIaQOsbd4/W+nGt9WZdwC9Ra71Xa/2w1rowh2mKz6df/80UKOl1bQNuByYDny1mzIJQCE4WEhOnm29QzbROFOLasjtuxWCmrTo9FUej0XHm/WqkcNsdE0g92WcXTPOieJoVdsG0XtRFccIu1d6KwcFBmpqa6OjoIJlMjhEophAxLSXZrpMgpeyb54oTTq0yit3W3Llz31PUQi7iJ0uK13wo/fobi88eBr6ansfRIiMIxZLPQpKv50aQnu4KpRDXVrkWpM7OzoLmszr+1bBaOB2TlpaWMZ950dDPCr+kuJvby22kacfg4KDlbyr7HMu2qJkWrJ6enqrto2EY9PX1Zb7bcDicSYMvNtXc6fww+x0V8n2aYqiurm5SUTvjIhNCpCilGoEYMGRjfdmcfp1buVEJEwU7U3UoFMrrU66vrw/M010xFCLc7Cikcmyh2GWNVCMWyOmYFNJszyu8SHG36tydfVO2o9DMLCdyj7ObPZRKxcrlaKbB79ixg5GREVe3Z2Wts9rXap53JhNCpADT06925SbN6VG7FSillgPLAWbPnu3awITaxy51sJCnYbcvTn7BSbiZT3l2aK2JxWKOjeMKLW/vZPqudCyQ3TEpN7W2WNxOtc2Os3FKq85Xm8au0m6x5LpP/RCwbnceaq09vwY47asf4uEmikgpG631KmAVwOLFi/3RoUzwnNwnvlAoRHt7e1EXL3PebFPuRMIqGNRJuOU7RmbFVaf5zPL2+b4np4twpWOB7I5JpQWKm1ktuW6KfJYQpxumW7+dXPepHwLWK7EtJ7FrN73SAtkK3wTOeoxpKZlu87k53fB+KEJQMM3AuRHyPT09RdcuiUajJV9kg1onxTAMNm3aNKavjlum9HxpuYUGCDoJkUrFAhmGQXd3Nz09PRWvi2LS0dFBZ2enqwKlkAJsVljdFN38DQwMDIxZX74eSpWg1G2Zy4XD4bzzzZs3r+h99UM83ISwpGith5VSvUBMKdVmEZdyaPr15QoPTfAxTv7YSpuCgxaX4lS2XmtNf3+/5xkqhTwB2lkvslOWvST3ONWSpa3YLBITqxtmqeuyIlco5+uh5Ca5VsWmpqZMunQpZHd7duprZK6/2H2NRqMFVx32ioliSQH4ffr1FIvPPpozjyA4XjhqzRTsNvkC7hKJhOf7VcjTqVUVz3wpy27ih8BEsO4GXC6lfL92N0y3z5VsS5vVOeBFOrJVinm+DsZO5MYORaPRvJaSUvbVjWDlcpgQlpQ0PwQuAL6ilHogq5jbgcAKYA9wZ/WGJ/iJfDUISjHP2vUEMbErnBXEOimFXHi99HcX8yTsRfZKofhFgHoRKFrs92tm9wDjUmS9OFfi8XgmbqkS54DbgnTnzp3jxHQhlpJi97Xa1x/fiBSl1BnAGem3renXDyql7kr/v11rfXXW/Neyr3T+kenXzyqljkv//5TW+sfm/FrrPymlbga+CPxFKXUfqbL4y4AZwGVSbVYwyedLN2sNRCKRMdH3DQ0NxONxyxoE5gUle93ZgYpWLpKg1kkp5Kbi5Q06KIW5/BCYaOL2OIppiGfWtLFLB66vr/fkOFUy1djt8Vs90HhR1yZdo6dqfkjfiBRSQuPCnGlz0n8ArwFXZ312Cqny+dkck/4z+XH2h1rrlUqpv5KynCwHksB64Fta61+XM3hh4mGKkWyyBUs8Hs/4iLOFip0rwS+Fs9ygpaWl4j1nTCKRSGCOWTWPUy5uPTHnphznlqV32m5fX59lOrBXabiVTDUupgN3oVgVZcu2lPT29tLT05M5x0rJ3opGo7z55puvuTrwIvCNSNFa3wjcWMT8J5a4nbuAu0pZVpgYdHd3u7q+vr6+gi+C1XQ91Arm07cfG8TlEo1GPSnWVSxuWexyLSHZqftTpkyx3M+mpqbMstUIHK6UJcuL2KPs+BZTjJiCZXh4eJxF2Hxf7G9jcHBwhzsjLh7fiBRB8IpCOhBnUwmzbK2Re4xHR0erOp5SL8ZukH0sshu9mece7LOW+QE3LXZ2cRfJZNJWiA0NDWWWrQb50nfdolLXAVOk24miYrpc+wERKUJNk5uWZ+WCEcrDKo7AD1TjYuyUUmyee0opX2T0gPvF20r57uPxeN72EF5Sie/CzcypQjtC1woTKQVZmIDYldKuZBxApZ7UqoVf0mj9QCHHwk/HanBwkO7ubtduokE81722cBTSsbhQOjs7x9RGmQiISBFqGqcLkN2Fw82UO6VUJq2yVvGL5cQPBPFYmO6BcoVKbnXmatHR0VHUb9jrFFu3RLxbxQ/d7s3kNSJShAmLXZqxU/noXMwLh1kELPsC6VVRKL9R7ToKdlTjYuzHY9HR0ZF3nkLbCNhRSAl8pRTNzc2eVxru6emxFEtWLQcqkeLvlnDNHn++80wpRUdHx7jfgNvuvUogMSlCzWAVIFtq2l9uuWm7QDTTUpOdUTIRzLG5jeOqSUNDA5MnT7atP1NJ/JRSnE2+QoJQ3s20EIGjtWZoaKgi/YmSySRKKZRSJJNJy6Blr1P83f6NZF/HnGrQ5KYjB02U5CIiRagJci8I2QGy5WJVv2R0dNSynkN/f39NWk5yM1a89OOHQiG01nlN5EqpMZaqci/GQ0N7aWqaZPt57jlmJYSi0ahl6mc1GRgYyAhnp3GVYwUqVOBU0h2mtaauro758+ePmV7u77OQbEEvRHz291NLNZXyISJFCCyGYdDf31+WH7yrq6ugp+7c+iV22Qh+8Mm7QfaFOBevAw2TySQdHR22jdi8uCD397/LwQd/j1dfvZy2tqnjPre66QwODlqKllgsRmNjo2/SjM0xmGPzorKxF4XK3MDt429XERfGih8vRGru9zNRaiqJSBECiZMLpliqWVPDjzh1U60Ulb4A//KXL7Fr1ygPPPASX/jC+8d9XshNJ/s8MsdvVgStJrkWEi+ewv2UsZSN2zFCVkGwlahaG8SsKbcQkSIEglwTayKRcPXCWGxNDTuXh9dBgW7hZLK2S9uuNZ57ro8330wVEvvBD/4MwH/+5zpmz54OQGtrE+97X3tR68w9j6otUOwsJG6JQPM88mPBQi+CYu2+T6+/50QiUdE+Q35CRIrge6xMrKHEMM3Gb5i09zX2TjqAndFTSIYbKzam9vZ2S2tDe3txN7VqYHU8s/t7+IGuri7P/exf/OIjPPnkazQ27nva3rp1kGXL7mN4OM6SJQewdu1FRa+3q6uLjo4OotGo5/E7Trh9/HKFbVNTE4ZhVMyKUsyxLHff7US8XSG1XItNIYHKxaK1tvydBjFjpxhEpAi+J9fE2jC8ngO2XgI6SVjvIqHqae3/Fq8d9J+MNC6qyJiCHLgWlOJrdv5+t3j00Qv40pceZdWq59i1K1XGf3g4Tn19HVdeeRTf+MbJmXmLven09PQwPDxcVQuDm1lmVsK20sHBhR7LcDhc1r47xZ1YZdUopWhqahrX7K9Sx6fW3dUiUoSqky9rIvvJJZQY5oCtlxBODmemhfUu0HDA1kvoPvz3JMMNRY+hlJoaQQ1cq7YLohi89PdPmhTm1ltP4U9/eoM//3mfi6uzcz9uueWUMfOa52MxN55qZvhUIhbDr5QbwOsUd2KKHyeLUraoscPtQOOg9eMpBhEpQkWwyhaJRCJEIpFxjccGBwfZs2cPc+bMGVcFc7rxG9A2T1Q6yfSdv2Fwxlm24zAvKn6oqVEtCun94Se8HKth7Oa55/qZMqWOpqYIQ0NxnnuuH8PYTTQ6Zcy8ZnbMpk2bfBmDYVLJWIxaJF/cSe7DSXd3t6WosSNov79qE4woPyHQmObT3B9mPB637Yw6MjJCb2/vOP/rpL2vpSwnFoT1LiJ7XnccS0tLC7FYjM7OzszfRBIokDoGlSio5RZeVHE1DIPu7m5+/OPHUQquu+5o3nzzam666SSUgkceedV2OT9bFEKhkCdVjv1YSdeOcoPX7fbVbrqT4Mj9nZkCMkjHs9qIJUXwnFJNxVbm8r2TDiCh6i2FSkLVE58823Z9zc3NlXXPxN+F1+6FdzfD1EPhgGUQGV+DoxLkBgJGo9FM3RE/44VVIDvmYPHiWdx//0c4+ODpvPvuO1xzzTF84hNzaW6eYrmsX90eXsdE+bWSrhXlBq/bVXNNJBIYhjHuGDtZRrTWYwJ+TdHi9vH0ugVEc3PzDE834ICIFMFz3LwR7oyeQmv/t8DqPqFC7Jy+L5bA9PtWJaj1b0/BEx9LuaYSwxBuhPVfhBP/H+x3nKubssq6yC58ZuUz91NFVNiXuREKhVBKefq9ZQuNmTOnMHPmlDGxL4cdNtN2Wb+JukgkUrE2DEopXwq0bBoaGso+X8zlcwtFJpNJy0BupxL15nK564jFYq5lAHntrjYMg9bW1gM820AeRKQIrmIVBOumDzYZbuS1g/6TA7ZegkITSo5AuBGNoveQH5EMN1Q/0yb+bkqgjL67b1oiHej7xMfgzD6INLmyqXxZF34UJNnklravBMXWunCqvltt3LQyOQWw+9WClIvpJi73ph2NRhkYGBgX3GoVyG3+X6hlRGtdVi2izs7OkpcthYGBAZRSVQsNEZEiuIZd6fCGhoaiL/BKKaLRqOUNdqRxEa+/72kOCj0L774CUw9BzV7GbJdu/GXz2r2Owb28fi8c/I+ubCooNw87/NYl2mx3EA6HaWtrA3CtsrHbuOm+tPvtQuo78qNAs8OtTBcvC7f5OfA6G6tYwkojIkVwBaeGWnbBsXZkW0IaGxvp6+vL/KjNm0fq4vyeMkftEe9u3mc5ySUxnBJWLlHtC0g5mAXPKkm+1FCTRCLh2xiMUChEe3u7q8fO7rdr3vBrOSPFzoLkVLjNjb5hdjilJ1eyorVfOp2LSBHKxq2T2cr0H8haJFMPTcWgWAmVcCNMPcS1TQX15lHxIOY0frjolkuln8INw2B0dLSi23TCfIhxQ0Q6WZCcCrd5aV1ra2sb82BWafzm3pQUZKEs3BIo4XDYd6b/kjlgGdi5cFUIZi9zbVNBSSc2Uy4jkQgdHR1VSfvOrbkTZAYGBiq2Lb+5u+bNm0c0GrXNaCkm08XJghSNRjMWFUidu7FYjKGhIU+PRzQatRUoXgsXu3IR1UQsKULB5GaRWBViK4ZKB4BVjMjUVBZPbnaPCqWmuxg7E41G2bFjR1nfQyWoVAaKHebFt1Zw+yZil2lSzd5DVmQLEKsqwG5nulhZcr12A3Z3dzu6fKzSoN3CjzFuIlKEgrDKIinnQlnzxYz2Oy6VxfP6vZngXmYvc1WgGIZRVbNwMbiRcVEOfrz4loPbvx+7G76f3GMNDQ00NjZmeuSY6ergfZ0YNzBTuPOVxM93XfWyn5WfLCgmIlKEDLk3vewg1VIu8g0NDZZP+F4U6PIlkaaCs3jsgvfsurHmika/U+kmaH4J+vMKL34/Ztl/2Gc19RN79uwZc85ni/NSm1HaCTE7l1E51ri6ujrmzZtHd3d3WQG3XvWz8qs7VETKBKMYl00ikcj8KItV2Nk3WfvsHAHsg/dyK8LG43F6enrYsWMH8Xg8MALFpFJN0GpFoNi5WrwOOvarAM53Yy/l5l2My8gwjLLOK/O37Ia1wguLh1Ptlv7+/q2ub7BARKRMAOwu2oW4bMwffjFZJNk/8kBm51QYuwuf3fH2e/xJJbGyQNWCQIHxfV9MhoaGXFl/rYi5bJyuUXZWyWwLUr7lysF00bnRAdltd59hGI5u48HBwR2ubrAIRKTUOG5ciOLxOB0dHZbpePX19WNumhOto3Ah2F0chfJwSh8NMmYdFLsATTeeooMoUAoJ4rW6eVvVNCnUPeSmVSkej2eKBZaDF+7yfK69+fPnv8/VDRaBiJQax636JeYPWW62xWEVcOxl4JufcbsJWtBusk40NDQwZ86cMdPsnt7deIoO2rEzhZuTYLC6eTuJDLM8vdM1zS8F/UxLtlfXXT8GzJqISAk4Tv023ML8gYvrpnisAo7Ni6N5LGvJRWGHWNjssTs2dsXEghR0vhvN+ojmLaWZpRWL4oopFFfXRyk1psKuKSrsmlEW455JJpMZ60zuA4SfUta9fjD0c1FIESkBxsnc3djY6Lvo/IlEPnN6MpnMRNO/8847FRpV5fC6Bk6+TISgCD+n4xR06+WrYc0P6pNoYK+CSVpz/2TNF3aFODjhLFRM107uPud7UCrXPaO1pqenx1cVV028tMK6VcHXC0SkBBinaomGYfguOn+iUKi/v6+vD611zX1PXtfAyVeYraGhIRMMadbU8CMNDQ155/HKeum1iNtNSqDsydIie9P//6A+yb8NhZjsYFEJh8PMnz+/6O26VQ/Hr+dMuenHdvFx0WjUtyJFyuLXKKX8UN0oMy0U7u9PJpOBFyi5GSheuiN6e3vp6uqip6fH8bjt3r2b7u5uurq6fNVzJpdq3ghjsZinv+v1EY3dN6TTnztR6rHxq7goFicLW6n7mFvy3rTMmFbJjo6OktbrNWJJEYB9fSnA2zLTQm0Ri8Uq4o4oJhslO87AzyKw2jfUxsbGcbV43OItpTOWk1z2qtTnXuDn2IpCMcWjUxfmYnD67eRaZsyquH5CREqAKbbfRiQSoampadwy2U+/+WoGCPuY6KnFkUjE82DqIKbKFko1W0N4XbBtllZM0tZCZZJOfe4FVsHG4E5tkkqQ/VDoRuB0Ib+feDzu2wJ+ICIl0NhZPhobG21P7mg0mgmqnag3VzdwSi0OStBmOZTj1ilU3NWyQAFvSttDYce3v7/f0xvSorji/sk2qcLpz/PR3d1d9LXJnDe3p5XfBYpSalwXeDcCpwv5/YRCIV/3thKR4iNKeTJ3snzYrUtSicvHLrV4YGCAlpaWmry5mk+j5QjbYurG1OIxzMaL36Dd8R0eHvbMtWPFFBRf2BXKye5JCZQv7HIOmjUpJpsl99rp1xuuFU6/p0pcq7Pdo35ERIpPcLvolwgRb7G72Mfj8ZpL/XbT2mYn7np6ejLZBRIHVTp2x7cagu/ghOLfhkLj6qQUIlBMCslmsbp2BoVIJMK8efOqPQxfIyLFJzg9mYvY8B9OAXpBukgWQlNTk2vnYCHHptYtKOBdPIrfzr3JKD5YgGvHiXz75GdXRT68+L5qzU0qIsUnTJQbXtCptQtAIbjZvbjQ7IvBwcGaje1xM007181RSH+boJFP0AX5Gum2WHXr+uSnLB8RKT7BrXQzwX3c6oIq2GdfWGEVGB50zB40blimguzmKIampibHz4OcdpxMJunq6nLNpZrvtxIOhx3rM4XDYQ4//PDMez8UQxSR4hNqoU9HLeLn1LwgkpuxkA8va3lUEjfFiUmQ3RzFsHPnzsw5YHUztxO+QRAvZtZRpRqPJhIJx3TsRCLBpk2bMi0J/HAMRaT4hKD36ahVJsqNAFJl2kdGRsZNd7Myaa452u4i2NzcXBMCMRwO09bW5snvuNo3j0rh1ATQfN2xY8e4czdox6dSMYj50rGzj7UfjqGIFB8hGTn+wjAMX/xIK0EoFGLOnDmudNW2S6W38pfbHd/BwcEx87rRTbcahEIhz37TtRh/UghWN3MrcR1Eyr3e1GIcl4gUQciiluNP7CwlAO3t7UD5FYedUulLvXiW00232rh5HhmGQX9/v+8Lk1WC7OOaryN2kCg3BjEWi9WEezQbESmCkKYW3At2hMNhS0uJVaXLcnBKpS+FcrvpVptSbzq51qimpqYJ2dncLn7CPK6GYfi2e28pJBIJVwNpawERKYJA7V3sslFK0dbWBnjfm8ntVPpCuumWW4fDK0oNfLeyRtWaCb8QzBu1VUJBU1MTL774Ys1ZlZxib5yoZQuwiBRhQpH9Yw6FQiilauJCl20RqWbjQ7ezAarVTbdYcutKlBMwO5GCte3I7jUG1IxVKRQKMX/+fLq6uvLOW2ggbS1bgEFEijCByP0x10rQoZ96MxVTB6UQqtVNt1i01q6VOC9V5AU1uDiXfOdzd3d3YG/IyWSS7u7ugoOezXPB6cGj1kVtwSJFKTVJa73Xi0Eopc4BTgCOBBYCU4E1WuvzHZY5BrgeOBqoBzYDdwC3aa0tH42VUh8HrgbeC4SBjcB/aq1/6trOCL4kyO6cfBH7fur9UWwdlHy40U23UpS7v+aNqBSCHFxskk/k1YpLIx6Po1Rh30koFMrULclePtsVFPTjkY9QEfO+opT6R6VUMcsUyvXApaRESm++mZVSpwNPAkuAXwLfByYBtwC/sFnmUuBBoBO4B1gNtAN3KaW+XfYeCL4lqAIlEonQ0dFBLBazDcD0Y0XiaDTKvHnzXBmb2U13sk5ZTiD1OlkX3k23UpSzv729vfT09JR0w8kOLjYtTnsV7FGp4OI9tlE9/iFf/I5pBa2VG7LWmnA4nDlnTNez1XxWFpfsYHQ/XgPcpBh3z9ukbuxfUkp9TWt9r9VMSql/AKZqre8oYt1XAT3AK6QsKo/bzaiUmpYeRwI4UWu9Lj39q8DvgXOUUudqrX+RtcyBwLeBHcBirfW29PSbgD8DK5VS/621frqIMQs+JuhPXZ2dnWPeB6UisRfH3Y1uul5TzndhGEZZgbFBDi42yeeirEWXRiKRIJFIZGoR5bp0ksmkY7yc+Rtramqq6cDqgkWK1vq9SqmzgBuAnymlrgWu11o/lDPrPwMnkXK9FLrujCgpwAx2DjALuNsUKOl17FZKXQ/8DvgCYy0qnwMmA98wBUp6mUGl1L8DPwE+D4hICSi5AbFa68Be1MLh8LhpQahI7GUAnxvddL2k2DRuN8VcUIKLnRgaGnL8PKgPG4VgFi7M/U3nC641LSi1LFCg+MBZBfwVOBRYAPxKKfUX4A1Slo2DSblTXnZzkDl8KP36G4vPngRGgGOUUpO11nsKWObhnHmEgGAYBn19fePMoUEOiM1OF87F7xWJ+/r6AisMy6VYgTIRg4udiMfjmWZ2VgLcDz1kvCY31iTfPjc1NdHd3V2p4VWNYgJnLycV82Ge8UlSMS0L0386/dlrpKwpXmFGVo0TQlrrUaXUVuAIYA7wYgHL9CulhoEOpVSD1ro26ivXOH6NM3Gq6loIbhZWcwurUvkw9gmuoaEh0OKwHIrtbeS268Kt4GLz+6uW0DRdG4lEgtdff53e3l5CoVBmbE7de2uJ119/nf7+/rz7/Pbbb7tyPEzvhXms/UYxlpQvAHHgn4Bfaq2HIRMjchNwCSlryvVa67VuDzSL6enXnTafm9OjRS7TmJ7P8g6jlFoOLAeYPXt2gUMVysEp7c6PAgVg8uTJzJgxw9LCY9LQ0MCuXbvGxZcEQaCAtXm5VnqnFEuhvY28jJEyg4vHZvekBEqhwcXJZJJQKEQsFmPy5MkFZ59UkonUq6i+vh6A0dFRTy1IWmv27NlDb29v5hzwG8WIlAOA32qt78meqLV+B7hSKfVTUnEgdymltmutf+viOKuO1noVsApg8eLFtS/nq4xp+jWJx+P09PTQ09Pjyx+SiXkDnz9/PmBthbAKkvNbfIlJrfu7S6WYmihWQs9t3AgujsViTJkyxcNRlsdEESjZArGurs5TkaKUYsqUKcRiMd544w3PtlMOxYiUHU7za62fV0qdBGwCvg54JVJMa8h0m8/N6UbOMjPTn73tsIydpUXwmGKeNP1wsXJ6qhscHMw8XduVofd7fImwj9xqssVk8pSbuVMM5QQXa62ZPHmyyyMSSiH3ISz3/POCyZMn+9aVVoxIeQz4pFJqntbaMlpHa92nlHoE+Jgro7OmG1gMzAWey/5AKVUHHASMAltylpmZXubpnGXaSLl6eiQepXIEPRvHD0LJS/wa8+M1zc3NmS6ypoULSs+qKrU4WzXwo4unVlBKFWwVSSQSjI6OUleXuj17bU0xx+dXihEp3wQ+BTyulLpUa32/zXwHlD8sR34PfBo4Bfh5zmdLgAbgyazMHnOZY9PL5KYZfzRrHqEC1Gp5+lphogoUwDa+pFSrV61npBRLOBwmFApNmOMSCoXGWagK2Xdznrq6uoxYGR0dDdSDnFsU7NzXWm8CVgKtwP9RSj2tlFqplHq/UiqmlHqPUuqHwPvZl1XjBfcB24FzlVKLzYlKqSnA/0q//UHOMncCe4BL04XdzGWagevSb3/o1YCFsfi9MJMbFRwNwyh/IBWit7eXrq6uzN9EFShuYhgGL77o5WUwmCQSiQkjUCD1ADY6Opp5X1dXRyQSKchyEY/HM8vW1dUxZcoUyxpKtU5REYha6+8D55OK6ziKlHXlGeB1YANwcXrWW4tZr1LqDKXUXUqpu4Br05M/aE7LLlufDtS9mFTvnSeUUj9WSn0zvf0PkhIxY6rhaq23AtcAM4B1SqnblVK3AH8hVdvlO1JttnJU6yJV6A/cjfH19vYGQqhUIqgzKJQakG0YBt3d3XR1ddHd3U1vby+9vb010V07KPT29rJixQoOOeQQpk+fzsEHH8wll1wyRnAvX76choYGPv/5z49b/vrrr6ehoYGzzjorM23p0qU0NDSM+3vf+95nO44bbriBQw89lJ0794U3xuNxPvWpT7FgwQJGRkYygqO+vj6vWMm9Fk2aNCnvsag1iv5Vaq1/BswGLiCVzfMKsJeUpWI9cL7Wek2Rqz0SuDD9tzQ9bU7WtHNyxvAAqfL5TwJnA5eRSo/+InCutnhM11rfBpxGqqngZ0ilE78JXKS1vrrI8QplUOleE6FQiI6ODg4//HA6Ozvp6OhwJUNIKUVHR4fl/mT31vAzIlD2kUwmixaWuT1l4vE4g4ODvrYU1hrbtm3juOOOY9OmTaxevZquri7uuOMOXnzxRY4//nhee+21zLwdHR3cf//9DA8PZ6aNjo6yZs0a9t9//zHr/fnPf86WLVsyfy+99BJTp04dI2Ry+cpXvsKMGTNYuXJlZtoDDzzA/fffz+rVq1FKjbGsFHIdGh0dZffu3ezatYvdu3dXPLsxHA6zadOm5/LP6Q3FVpwFUiXogTXpv7LRWt8I3FjkMn+kyABdrfWDpJoMClXEqg+NF3R0dIyLJSin2mdDQwPxeHxcAKWde8RvZu0tW7ZM2HomhTIwMFBU/InfXZcTgSuvvJJQKMRDDz1EQ0MDAPvvvz8PPfQQ73nPe7jyyiv55S9/CaR6YvX39/Pf//3ffOYznwHg4YcfZsqUKRx33HG8/fa+5M8ZM2aM2c4vfvELRkZGuPDCC23HMmnSJFavXs2SJUv4xCc+wdFHH83ll1/OV77yFRYsWACkrgt1dXWMjo4WZG3Lvo5UI8HArgJ2pShJpAhCqZhZPV7+0BoaGpgzZ47lZ8Vsu9BaGHblq/3UnTS33ftEJV9BMDthaVfXxm9CtNq8++677NixI5OdMmPGDKZOnerZ9nbs2MGjjz7KDTfckBEoJg0NDSxfvpybbrppjMXwwgsv5O67786IlLvvvpsLLriAbdu2OW7rzjvv5OSTT6ajo8NxvgULFnDddddx+eWXs2DBAubMmcMXv/jFcfNlW1T8SjgcrnqpBBEpgidYXdSBilhQ4vE4hmEwPDw8rpBaoTeVYmph+L1DsVlNshqERkZI5tw8qkUoFGL+/PmOgi1XWFr1h8ouLCjs49133+Wtt97K/A5GR0d56623ADwTKq+88gpaa9uHicMOOwytNa+88kpm2rJly7juuut45ZVXaGpq4tFHH+Xmm2/mX//1X223s3nzZv7whz9w77332s6TzcqVK7n77rt54oknWL9+/bh4uF27dhW0nmri1EeskohIEVzDriCb2TjLrIdSLs3NzTQ2Nma2FQ6Hx5hNzZtILoXGXxRbC8NPHYqtxGG14k7q3nqLuR/9KC8//DCjs2ZVZQzZtLe3Z16tzo9cYTmRU7FLYceOHeN+31prduzY4ak1pRCyA06bm5s57bTTuPvuu5k+fTpLliwZF4+Sy5133klraysf/ehHM9P++Mc/csYZZ2Te33bbbZx77rkArF27ltdff526ujqeeeYZDj30UHd3yGOyr2GGYTB37tz3VGssIlIEV8gX66G1dvS/dnZ2jiuFb0V2rxRTBHR3d5eVSeFG35xqV5B1euKvFtN+9ztCe/Yw7fe/Z8eyZVUbR65oNF+zj1c4HKatrW3Md9jX11fhkQYbO/eFl26NQw45BKUUL730Eqeffvq4z1966SXq6uo48MADx0z/zGc+w8UXX0xTUxNf/epXHbexd+9e1qxZw0UXXZSpWQKwaNEinnnmmcz7/fbbD4CdO3fy+c9/nssvv5z99tuPL33pS5x00kl53UR+IVeg9Pb2UldXV7W0IhEpgiuUE2dimtidAmqdrBPlxAX4uW9OoZQTDOw2UzZupC4dfDgjbRqf8YtfsDdtNh79u79j9xFHVGw8dnFFhYhKieEpDjMY1Gq6V8yYMYOTTz6ZVatWcdlll42JSxkZGeFHP/oRp512GtOnj+2ictJJJzFp0iTefvttPvGJTzhu48EHH2T79u1cdNFFY6bX19dz8MEHj5v/mmuuYfr06Xz1q18lEonwq1/9iksuuYRf/epXpe9oGqVUJq5Ka+1JyXzT8g3Q399f9euKiBShJHLdCoUIBavy99kmdju3iTmtp6eHgYGBMaKinFoknZ2dJS/rJ/yUYdL6rW/R9NxzJNJdXAEivb3sf/XVhHftYuh972PbXXdVZCx+iguaCMyYMWNMTAqkvoPcLBm3ufXWWznxxBM59dRTueGGGzjkkEPYsmULN954I5FIhG9/+9vjllFK8eyzzxbUs+iOO+7gpJNO4qCDDso7lgcffJB7772XJ554IrPeVatWcdRRR3HHHXfwuc99rrSdZF8zwGy86pKstXbs5F5J/NtOVvAtVrUh8qGUor29nVgslrGcRCKRcW6WaDTKvHnz6OzszDwB527LLJRmjqMU/JR5Uy5+yjB5bfVqtp9/PkprwungwPCuXSit2X7++by2enVFxmF1bgneMnXqVGbNmjWm58ysWbM8j0c58MADeeqppzj88MP5x3/8Rw477DBOOeUUQqEQzzzzDK2trbbjnTZtmuO6t27dyhNPPMFnP/vZvOPYvn07l112GV/60pd473vfm5k+Z84c/vVf/5V/+Zd/4fXXXy9u57LQWrNr1y727t2bmeblb98PAgVA+eUJLEgsXrxYr1u3rtrDqBqFxI5kU45LxW5bpsiwG0dHR8e47B4TN2JQKoFhGPT392fibUKhEO3t7ePGXez3UQnmfOpTNHR1Zd6PdHay5ee5rbbcw6omTrl0ZY2/1kkkEsydO7fk5evr69m7d69vquz+4Ac/4F/+5V+45557+PjHP17t4QSCl19+2bYq97Jly9i4cWNVuhCKu0coGqcboun6cSvWw25b+W7KZsxBLBazrXHhZ6ziTJLJ5JhA2OzsJj8Reucd6jdtIjl5Msn6ekK7dlG/aROhd94hmefJtVSKLcKWj1ItdAC70ayPaN5SmllasSiumIJ/u8yWi1na3S8CBeALX/gC++23Hy+++CIf/vCHqc9yPwrBQkSKUDROxcsKKX7m1rbAWqzkunKqnXlTCk5xJn19fWNie/x0cwBo+tOfQCkGLr2Uty+4gJl3303Ld79L05/+xDunnOLJNt20JBmGUXLa9qthzQ/qk2hgr4JJWnP/ZM0XdoU4OFGbQsXLwNhyOPvss6s9hECQr8BhtfHn2SX4mkoWL8u3LT8XUXMin3XH6abr5wsKwPDixWy+/372pqv+bv/sZ3nnhBNIeGRFgfJijEoJArdiNymBsidLi+xN//+D+iT/NhRicg1ZVJRS1NXV+VakCPnJrTHlR+TsEoqmksXLCtlW0F05ZjDw8PAwQ0NDvosvKZbEzJkkZs4cM22vTZsCtyhVmFp9F6WyPqKxi/DT6c8/GK8NkRIOh8d15A3CDU8YSxC+LxEpQklU0oXitK1aceVoraUjcQnYBRMXipvp228pnbGc5LJXpT6vFRKJRCBKuwvBR0SKIFSYoFtK/EB25eFycPO7mKUVk7S1UJmkU58LglAcIlIEQfA9bgdlm3EobrIorrh/srW1RKU/FwShOESkCILge9wIhu7t7fXUpTYFxRd2hXKye1IC5Qu7aitoVqhNIpEIyWTSV7EqIlIEwWOys0dCISnybIVZYM+qAF85hdrsOnN7xcEJxb8NhcbVSRGBIvidSZMmMW/ePN8VMRSRIggekps94vf04UpiVfjPLMDnBtVqvDgZVTNZPMLEwbw2lZOG7wUiUgTBQ/zU/K9a5HZq9bItgdcuHUGoVUwXj1M3+mogtmdB8BA/PZFUg46OjrxNJd1CBIrgJ5YvX85ZZ51V7WEUjNY64x71i0ABsaQIQlnk3hjdSo2tBSKRSEaMVKKWjQgUwS38Xpjuscce47TTTnOcZ/Xq1Xz6058ueJ1KKV9ZUExEpAhCiVg9uQ8ODjI4OEgkEqGpqalKI6s+ue0JRMwJQWLSpEm+LlZ37LHHsmXLlsz7JUuWcPbZZ3PFFVdkps2YMaOodWb3A/MT4u4RhBJxenKPx+MT9sk+16VjJ+bK6TTc29tLV1dX5q+cdQnuMZJM8ptdO7nj3bf4za6djFQgUPypp57ihBNOYNasWbS2tnL88cezceNGAB555BE+8pGP0N7eTiwW47TTTuOll17KLLt06VIuv/xyrr32WmKxGLNnz+b2229nz549XHnllbS1tTF37lx+9rOfjdnm0qVLueyyy7j66qtpb2+nvb2d6667zjEwXmvNzTffzBFHHMGMGTN4//vfz89//vOC9iOX+vp6WltbaW1tpb6+nt7eXo4++ujMtNbW1nFtC4KKWFIEoQgqndIaNDo7O+nt7aWnp4eenh7HeQcHB0uyptiJHqG6bNy7i6/t7EHrVLPFKShWD/2Nm6Z3cMSkek+2OTo6yic/+UkuvPBC7rzzTuLxOBs2bCAcDgMwMjLCpZdeSmdnJ7t27eIb3/gG55xzDuvXr8/cxO+9914uu+wy1q5dy0MPPcQ111zDo48+ysknn8xTTz3FPffcwyWXXMJJJ51EW1tbZtv33nsv559/Po8//jhdXV2sWLGC1tZWLr/8csux3njjjTzwwAPccsstzJ07l//5n/9hxYoVRKNRTj75ZMf9cGLDhg1orVm0aJELR9R/iEgRBAdya5z41STqB0KhkCvBq1u2bGFkZCTzvqGhgTlZDQpFkDgTDoeZNm0ahmFU7FwdSSb52s4edmVtbzcaNHxtZw/3/N3B1HtQI+idd97BMAw+9rGPZc6R7MrEZ5xxxpj5f/SjH9HS0sK6des45phjADj88MO5/vrrAbj88sv5zne+Q11dHStWrADguuuu4+abb+aZZ57hzDPPzKyrtbWV73znOyilmDdvHps3b+Z73/uepUgZHh7mtttu48EHH+TYY48F4MADD2TdunX86Ec/4qijjnLcDyeef/55Zs6cyf777z/us2XLlvGHP/yBE088cZw1KCiIu0cQsjAMg+7ubrq6uti4cSM9PT0Zq0kymRSB4kB7e7vrAgVST8Om/90wjLLWPxFoa2tjaGiooufqk3vexW5zWqc+94IZM2Zw/vnnc9ppp3HmmWfyve99jzfeeCPz+ZYtW7jooos44ogjaGlp4cADDySZTI6Zp7OzM/O/UopZs2ZxxBFHZKZFIhGam5v529/+NmbbH/jAB1BqXz2co446ir6+Pt55551x43zxxRfZvXs3p59+OrNmzcr8rV69mq1bt+bdDyc2bNjAwoULLT9bsWIFq1evLmg9fkVEijChyRYlmzZtore3NyNKRJDYEw6HM9VzI5FISVVhm5ubx03LFSi5093ut1OLZAvrStE3ujdlObFgN5q+xF7Ptr1q1SrWrl3Lcccdx0MPPcTChQt59NFHATj77LPZvn07t912G2vXruXpp5+mrq6OvXv3jcdMjzdRSo2bBuUVYjSXve+++3jmmWcyf8899xwPPvhg3v1wYsOGDbz3ve+1/GzJkiVMnTq15HH7AXH3CBOWXNeEVIMtjOwnz1IpJbvHb+W6hX20101iCspSqExB0R72NohzwYIFLFiwgJUrV3L66aezZs0aFi1aRHd3N7feeisnnHACkHKNjI6OurLNP//5z2itM9aUZ599lra2NqZNmzZu3sMPP5zJkyfz+uuvc+KJJxa1HyeffLLt/ENDQ2zevJkjjzyy3N3xLSJShAmJYRgS2+ABzc3NlsdVUo5rmyWTp7J66G9YGVOUSn3uBdu2beMnP/kJp556Ku3t7WzdupWuri4uvvhimpubmTlzJnfeeScdHR309fVx3XXXUVfnzm2vv7+fa665huXLl7Nx40ZuvfVWvvzlL1vOO3XqVK644gquu+46tNYcd9xxDA0N8eyzzxIKhfjwhz9sux9OvPDCC2itRaQIQq2RL/NEsMbKDJ6NKURKrYnit74hQmE0hELcNL1jXHaPUnDT9A5PgmYhlYq7efNmPv3pT/P222+z3377ce6557Jy5UpCoRB33303V199NYsXL+bggw/m61//Ouedd54r2162bBmJRIITTjgBpRQXXnghl112me38N9xwAy0tLXz3u9/liiuuYNq0aSxYsICrrrrKcT+c2LBhA9OnT+eggw5yZZ/8iBK/e/EsXrxYr1u3rtrDEMpAXAfF42XPHZPu7m4RKVUgkUgwd+7cstezK5nkyT3v0pfYS3t4EksmT/VMoFSTpUuXMn/+fG655ZZqDyUvTz75JD/84Q/zZve8/PLLtinPy5YtY+PGjVXpmimWFEGoAn4vuw0pC8jQ0NC4TsVekR20LAST+lCIpfXTqz0MIc2pp57KX//6V4aHhznkkENYs2YNRx11VLWHVRQiUgShCvhZoFQjfkSaAwqC+zz00EPVHkLZiEgRJiR2AZ61QigUIplMZl4LpZRU4nIwDIP+/n5fizZByOW3v/1ttYcwYRCRIkw4DMNgaGjI8+10dHQAVLSMvpVbplArRXbX4kog8SeC4C/MuDPwT3KBiBRhQmEYxrh25EoplFIl1Ulpbm4eV348N8DUfPUiWLeQWJFYLJZXpOR2LXaL7LYCQYjDEYSJjNaagYEBWlpa6OjoyPx2R0dHvavGlwcRKcKEYmBgYFwl2UL68SilbIVIY2Nj5sdsJxrcLuceiUQK7u1RyLq8CIrNteCIQBEE/xOPx+nt7SUWi2WuMS+//PJfqzUeESnChKIU90Ku2yb3ph6NRvPe4N0s516K1cMpBqepqclVgWIYBn19fVLB1yXMQGbDMHxjghf8S7FxaFaYFpVKun/tEJEiTCjsioWFw+FxDQTt3Da5ZLs07KwSpYijUCiEUmqMBaJUq4dVkTUTc5obGT1W7jQhP7mWuuzpjY2NwL7zT4SK4IRbDwd+iRcTkSJMKFpaWixjUtra2gB7a4kduTdl01QKY0VNMZVU3XTlZOMUmzI4OFiQSMkVZE1NTWNqqUin6OJpbm6msbHRUnz46YlWmFhkV5dubm6eUa1xiEgRJhTmxd7JdVMMdjEuuTcWK3FkhVcBrG5gJciyRY9fnryCRHbKt52FJPu4ShdooVKMjo6yadMmkskkbW1tVau7LyJFmHAUEkNSKHY35tzpVuKoqamJd955J+POCYVCtLe3+/ap2UqQCeVRiLUt+4lWhKDgNqFQiHA4nOnmbFJIQkElqL2mCoJQQewa7llNj0ajzJs3j87OTlpaWjAMY0y8STKZZMeOHZ6NFVKuhWKmZyM3SPfJzvpqaWkZd6PItqy5nSEmCKFQiMmTJ7vWGdoLRKQIQhnku7HYYWeVGBkZycS0eEEsFhsnSAopg+/lmCYy2e6b4eHhcedELBZjeHiYrq4uCZi14LHHHqOhocHxb82aNdUepi9RSjF58uTMez9YTazwr3wShACQL8bFDierRKFBrKUSi8Uc158dHBtKd7CVdGJvMM8Du6rAIkycOfbYY9myZUvm/ZIlSzj77LO54oorMtNmzKhazKev0VozOjqasaLkPmz5BREpglAmpcS4FJPtUwmyhUk2Ik68xXQL1lofqaGhvTQ1TfJ8O/X19dTX1wOwc+dOent7Ofroo2ltbfV827WAKVJGR0erPRRbxN0jCFXATxk8W7Zsoaenx1eiaSLg50yucujvH2L27O/T3+99f6xsNmzYgNaaRYsWVXS7QcZ08YyOjvrW3SMiRRCqQDQapaGhwfKzQoJYy8UwDLq7u+nq6mJkZMTz7U1kwuEwsC+LAlIWFLNQYK3F+/zqV5vZvXuUBx/cXNHtPv/888ycOZP9999/zPRly5bR3t7OeeedV9HxBAHTxeNXgQLi7hEEz8mNNzADVefMmWP7mZdIVdjKYIpN8/tNJpPjvl/DMGrC1bN+/ZsMDAwDsHr18wCsWvU8++8/DYCWlkYWLfLWBbNhwwYWLlw4bvqKFSv4zGc+IwG0FpgxZ34msCJFpSTgP6X/jgAU8CLwY2CV1nqcM10p9XHgauC9QBjYCPyn1vqnlRq3MLGwCojMLkOfL4jVC6TeSeWw++6zm1LWAtde+zhPPfUGDQ37Uu+3bdvJBRf8ipGROMcdtz+PPPIpT8ewYcMGTjvttHHTlyxZwpNPPunptoNKIpFg165d1R6GI/6XUfbcA6wCDgR+TkqcNAA/AO7KnVkpdSnwINCZXnY10A7cpZT6dkVGLEw4nMrQVxrDMNi0aVPN3Bj9jtN339vbW1Pfw69//UlWrHgfWmtGRlL7NTISR2vNpZe+j1//+pOebn9oaIjNmzdz5JFHerodofIEUqQopc4EzgO2AkdorS/WWl8BHAn8GrhAKXVW1vwHAt8GdgCLtdYrtNZXAQuAV4GVSqkPVnYvBKFymB10JVunMuSLK6o1S9akSWG+9a0PM3/+zDHT58+fyTe/+WEmTQp7uv0XXngBrbWIlBokkCIFODP9+h2t9XZzotZ6L/DV9NtLs+b/HDAZ+L7WelvW/IPAv6ffft6z0QpClTADZKXeRuUwg2InGoaxm+efH2DKlDB/93f1TJkS5vnnBzCM3Z5ve8OGDUyfPp2DDqpaixnBI4IqUswIrC0Wn5nTjldKmYn6H0q//sZi/odz5hEE1yinDH0pGIbBiy++SFdXV6ZKaS25FYLAvHnzJmQJ+8ce24pS8LWvHc+2bSv46lePRyn43e+2eb7tFStW0N/f79uCZELpBDVw1rSeWMnmOenXuvT/LwFm3/uXc2fWWvcrpYaBDqVUg9Za8jEF1zCfqCuRwdPd3S2CpIoopTLfa19fX5VHU3mOP34269Z9jnnz/g6Aq676AB/72MFEo1OqOq5TTz2Vv/71rwwPD3PIIYewZs0ajjrqqKqOyQuUUtTV1dXcNSCoIuUh4FPAF5VSv9Ba7wBQSkWA/y9rPvNxdXr6dafN+nYCjen5LEWKUmo5sBxg9uzZZQ1emFh4lcFjVyVW8IaGhgbmzEk9AxmGQV9fXybGJxwO09bWlql7MhFjf1paGmlpaRwzzRQs1eShhx6q9hA8RymVKXNfawRVpPwCuABYCmxSSv1fYDfwEaANeB2YDbh2pdBaryKVTcTixYtrK+pNCBxS66TyZItBsw2CKVQSiQQ9PT0S+yNUBfM6UIvXg0CKFK11Qin1CeCLwPnAhaREyhPA2cB96Vn/ln7dCcwkZSl522KV+SwtguAbzEwdobLE43G6urqqPQxBmFAEUqQAaK3jwDfSfxmUUlOAQ4HtWuut6cndpETKXODpnPnbSLl6eiQeRfAzuS4GoTwaGhqYMWOGCD5B8DFBze5x4lxgEqkCbya/T7+eYjH/R3PmEQTfYbp3qiFQzE69tYQZX1Js92pBCBqhUCjQWU+BFSlKqWkW044EvgUMAv+R9dGdwB7g0nRhN3P+ZuC69NsfejVWQSiXapWyD4VCzJs3L/+MAaK5uTkTAGu+F4RaJZlMUlcXWKdJcN09wKNKqV1AF/AucDhwKrAL+ITWOpMDqLXeqpS6BvgesE4pdS+wFzgH6CBVFO7p3A0Igl+oVgZPMpmsqTiMSjRwFAS/EeQMwCCLlPtIuXbOB+qBXlLZN1/XWo9zMmutb1NKbSPVYPAzpKxIm4DrpcGg4HcikUigLzTVJjt9OJda6EIsCLVKYEWK1vpbpFw7xSzzIKkmg4IQKFpaWspOOTZrKUw07KwnhmHQ399fhREJgvuEQqGaDKoPrEgRhImEGeBpFm+LRCK0tLQA2GanZBcYK4QXX3yRRCLhxnB9gZNrR+rMCLVEKBRi8uTJ7N27t6Z+wyAiRRACQzQatRUcdtVPi2HatGk15fpwij3p6+sTgSLUBOFwmEmTJuWfMaCISBGEgOMkXgrFMIwJ0xTPMIyaNIsLE5datKCYiEgRBKGgFOeGhgZGRoJT77C7u5umpiZ27tw5xsokFhShEixfvpzt27dz//33e7qdWhUnJoGtkyIIgnsUkjlUikCpZiG4eDzO4ODgGKtJIpEQK4oQeB577DEaGhoc/9asWVPtYbqCWFIEQXA1xTkSiYwpANfd3V30uiORCE1NTTUVIyMIbnHssceyZcuWzPslS5Zw9tlnc8UVV2SmzZgxoxpDcx2xpAiCQEtLi2uls+PxON3d3ZkYF6t1K6Xo6Oigo6PDcrumFaQWS/IL7vLUU09xwgknMGvWLFpbWzn++OPZuHEjAI888ggf+chHaG9vJxaLcdppp/HSSy9lll26dCmXX3451157LbFYjNmzZ3P77bezZ88errzyStra2pg7dy4/+9nPxmxz6dKlXHbZZVx99dW0t7fT3t7Odddd52il01pz8803c8QRRzBjxgze//738/Of7+ve4rQfSilCoX236/r6elpbW2ltbaW+vp7e3l6OPvrozLTW1taaCaYVkSIIAtFolFgsVpIoCIfD46bF43F6e3sxDGPcuiORCLFYLBPwG4vFLNdhrkcIEPF3CW+7k7qu6wlvuxPi73q6udHRUT75yU9yzDHH8D//8z+sXbuWSy+9NHM+jYyMcOmll/Lkk0/ym9/8hmnTpnHOOeewd+/ezDruvfdempqaWLt2LStXruSaa67hk5/8JIceeihPPfUUn/70p7nkkkvG1dS59957SSaTPP7449x2223ccccdfP/737cd64033shPf/pTbrnlFtavX8/VV1/NZZddxsMPP2y7H5MnT6a+vp4pU6YwefJky9/Jhg0b0FqzaNEil46qv1ASRFY8ixcv1uvWrav2MATBM4qpI9LR0ZGp35JLruun3G0J3pBIJJg7d25Z6wht/xOTnj4DtEYlhtHhRlCKvR98gOTMY9wZaA47duygo6OD3/72txx//PF55x8eHqalpYVHHnmEY445hqVLl7Jnzx6eeOIJIGXtOOCAA/jABz7AfffdB6SE8syZM7nrrrs488wzgZQlpb+/nxdeeCFjCfyP//gPfvzjH/PKK68AYwNnh4eH2X///XnwwQc59thjM+O55ppr2Lx5M3fccce4/bAqvmg17dZbb+Xmm2/m9ddfHzO9p6eHf/zHf+Stt96irq6Oa6+9lrPOOsv22Lz88su2DwvLli1j48aNVelSKDEpgiCMw0xptisUl42TwCjEElKt5omCi8TfZdLTZ6BGhzKTVGIYgElPn8Huj26BuibXNztjxgzOP/98TjvtNE488UROOukkzjzzTPbff38AtmzZwk033cSf//xntm/fTjKZJJlM8sYbb2TW0dnZuW/MSjFr1iyOOOKIzLRIJEJzczN/+9vfxmz7Ax/4wBhX5VFHHcVNN93EO++8w7RpY/vfvvjii+zevZvTTz99zDLxeJwDDjgg736YWP1ONmzYwMKFC8dND4fDfPOb32ThwoW8+eabHHvssSxdupTGxkbHY+o3xN0jCIIl0WjUNmYkm3IFhrh0gk+49z6wOw+0Jtxzn2fbXrVqFWvXruW4447joYceYuHChTz66KMAnH322Wzfvp3bbruNtWvX8vTTT1NXVzfG3ZPr4lRKWbo9y8kKM5e97777eOaZZzJ/zz33HA8++GDe/XBiw4YNvPe97x03va2tLSNeWltbmTlzZiAD0UWkCIJgixkzkh20Vyy9vb2On0twbPBRQ69mLCfjPksMo4a2WH7mFgsWLGDlypUZd8maNWt4++236e7u5pprruFDH/oQhx12GO+++y6jo6OubPPPf/7zGIH+7LPP0tbWNs6KAnD44YczefJkXn/9dQ4++OAxf7Nnz3bcDyeGhobYvHkzRx55pON869evJ5FI0NHRUdxO+gBx9wiC4IgZ4GoYRkHun1wGBwcdS9S70TxRqC666WB0uNFSqOhwI7rJugN1uWzbto2f/OQnnHrqqbS3t7N161a6urq4+OKLaW5uZubMmdx55510dHTQ19fHddddR12dO7e9/v5+rrnmGpYvX87GjRu59dZb+fKXv2w579SpU7niiiu47rrr0Fpz3HHHMTQ0xLPPPksoFOLDH/6w7X448cILL6C1dhQpO3bs4OKLL+b2228vZ3erhogUQRAKwoxTcVtQ5DZPFIJHInYOkb9a36BRikTHOZ5st76+ns2bN/PpT3+at99+m/32249zzz2XlStXEgqFuPvuu7n66qtZvHgxBx98MF//+tc577zzXNn2smXLSCQSnHDCCSiluPDCC7nsssts57/hhhtoaWnhu9/9LldccQXTpk1jwYIFXHXVVY774cSGDRuYPn06Bx10kOXne/bsYdmyZaxcuZKjjz66rP2tFpLdUwKS3SNMZAzDKFpQZAcnWq0nEokQiUQCVXa/lghqdk+1WLp0KfPnz+eWW24pe11WGTtuoLXmoosu4tBDD+X666/PO79k9wiCUBOY7p+urq6C5m9ubh43LTftOB6PixUl4CRnHsPuj24h3HMfamgLumlOyoLiQVZPLaG1JhQKud6u4emnn+a+++6js7MzE5z7k5/8xPKBwc+ISBEEoSQKKaXf3NxsGY8iacc1Sl0TiQMvqvYoAocX/aSOOeYYhoetg5mDhIgUQRBKwirgVSmVqSYLKYuJ2bsnEonQ0tJCNBoVq4kQaH77299WewgTBhEpgiCURG7Aa7YIAWuXjpmO7IV5WxCE2kNEiiAIJWPGp1hh5dLRWjMwMCACRRCEgpBiboIgeIKdS0dcPYIgFIqIFEEQPEEqyQqCUC4iUgRB8ISWlpZqD0EQhIAjIkUQBE+wi1URBEEoFBEpgiB4gmEY1R6CIAgBR0SKIAiuU2ozQkEQhGxEpAiC4DoDAwPVHoIgCDWAiBRBEFxH0owFQXADESmCILiOpB97g1KKjo4OOb5pHnvsMRoaGhz/1qxZU+1hCmUgFWcFQXCdlpYWV2JSwuEw06ZNY2hoKFN6v6mpCcMwJmSDQq21xPpkceyxx7Jly5bM+yVLlnD22WdzxRVXZKbNmDGjGkMTXEJEiiAIrhONRunv7yeRSDjO19zcPEaAZPf+ycfg4KALIxU8YWgImpo830x9fT319fUA7Ny5k97eXo4++mhaW1s937ZQGUSkCILgCfkESjgcJhaLFb1ewzAkvdnP9Pcz5Ygj2L1xI7S1VWyzGzZsQGvNokWLKrZNwXskJkUQBE9wiptQStFW4g3MqnGh4B/Cv/oVavduwg8+WNHtPv/888ycOZP9998/M62np4elS5eyaNEiPvCBD3D//fdXdExC+YglRRAET2hpaaG3t3ecoAiHw7S1tZVckdYvmUMNDQ2MjIxUexi+QK1fj0qnndetXp16XbUKnRYMuqUF7bGFY8OGDSxcuHDMtHA4zDe/+U0WLlzIm2++ybHHHsvSpUtpbGz0dCyCe4hIEQTBE0wRMjAwUFLMiR2RSCSvUFFKeW5tEYGyj8i11xJ+6il0Q0Nmmtq2jUkXXIAaGSFx3HHsfeQRT8ewYcMGTjvttDHT2traMha71tZWZs6cyeDgoIiUACHuHkEQPCMajTJv3jw6OzuZN2+eK/18WlpaUEpZfhaJROjo6CAWi2XcTZFIhObmZknb9ZC9v/418RUrQGtUWrypkRHQmvill7L317/2dPtDQ0Ns3ryZI4880nae9evXk0gk6Ojo8HQsgruIJUUQhEBRqIXGShAZhmHpghLKZNIkRr/1LULPPEP4uecyk5Pz5zP6zW96vvkXXngBrbWtSNmxYwcXX3wxt99+u+djEdxFRIogCIEjGo2WZJUZHh4WgeIVhkHo+efRU6ZAYyMMDxN6/nkwDPC4I/aGDRuYPn06Bx100LjP9uzZw7Jly1i5ciVHH320p+MQ3EfcPYIgTAh6e3ultoqHhB97DJQi/rWvsXvbNuJf/SooRfh3v/N82ytWrKC/v3+cG1BrzfLlyznhhBM477zzPB+H4D5iSREEYUIgAsVbEscfT3LdOvS8ean3V11F8mMfQ3tsRXHi6aef5r777qOzs5MH0ynRP/nJT+js7KzamITiEJEiCMKEp7Ozk66urmoPI9i0tKBbWsZMMgVLtTjmmGMYHh6u6hiE8hB3jyAIgiAIvkREiiAIE4Lm5mbbzwzDcPwcsE17FgTBO0SkCIIwIXDqEzQwMEAsFrMVKkopYrEYDVnFygRB8B6JSREEYcJjVrCNxWLEYjEMwxhXhwWkyqwgVBoRKYIgTBjsSurnVqPNrcNiGAY9PT1eD6+qaK3FpTVB8XPtIBEpgiBMGKyaHiqlMpYSSAmSvr4+kslk5nM/X8TdQCnFnj17mDJlSrWHIlSBPXv2+FagSkyKIAgThmg0Oq6vTywWy1hNTIuJKVDA30+ZbtLb28vu3bsnzP4KqXN79+7d9Pb2VnsotoglRRCECYVTSf2BgYHKDsYnhEIhkskkb7zxhoiUCYZpQQmF/GmzCLRIUUqdClwBzAf+DugHngNu1lo/bTH/McD1wNFAPbAZuAO4TWudqNS4BUHwJ1bxKhMFv96khIlNYM9KpdQ3gF8Di4DfAN8F1gOnA39USp2fM//pwJPAEuCXwPeBScAtwC8qN3JBEPyIYRjVHoIgCDkE0pKilGoFrgYGgAVa679lfXYS8HvgJuCe9LRpwGogAZyotV6Xnv7V9LznKKXO1VqLWBGECYhhGL72ywvCRCWolpQDSI39f7IFCoDW+nHgXWBW1uRz0u9/YQqU9Ly7Sbl/AL7g6YgFQfAtAwMDEoshCD4kqCJlM7AX+IBSamb2B0qpJcBU4LGsyR9Kv/7GYl1PAiPAMUqpyR6MVRAEn+NFLEokEhlXf0UQhOIIpEjRWu8Avgy0AJuUUquUUl9XSv1v4BHgUeCfsxYxW3G+bLGuUWArKdfXHE8HLgiCL/FCTMTj8QkdiCsIbhDImBQArfWtSqltpLJzLs766BXgrhw30PT0606b1ZnTo3bbU0otB5an3+5RSklfd8EtZgLbqz2IWqO5uXnGrFmzYnV1dZNGR0f3vvXWW72Dg4M7rOaNxWKzp0+fPsvqsyAyODiYt2GiIBTK1q1bq7btwIoUpdSXgH8HvkcqU+dN4DDg68AapdSRWusvubU9rfUqYFV62+u01ovdWrcwsZHzSXAbpdS6vr4+OacEV1BKrcs/lzcE0t2jlDoR+AbwK631F7XWW7TWI1rr9cCZQC+wUillum9MS8n0cSsbO93wZsSCIAiCIBRLIEUK8PH06+O5H2itR4BnSe3be9OTu9Ovc3PnV0rVAQcBo8AW10cqCIIgCEJJBFWkmFk4dj5kc/re9Ovv06+nWMy7BGgA/qS13lPg9lcVOJ8gFIKcT4LbyDkluEnVzicVxNoASqlPAveSKub2Pq11b9ZnHwUeAvYAHVrrt9PF3F4FpgHHZhVzm0JKwHwQ+JQUcxMEQRAE/xBUkRICfgt8hFThtl+SCpw9nJQrSAFXaq2/m7XMGcB9wG5SZfB3AKeRSk++D/ikDuLBEARBEIQaJZAiBUApFQFWAOeSajDYQEp4PAt8T2v9iMUyxwJfIWU5mUIqXfmO9PzSYFAQBEEQfERQY1LQWse11rdqrY/WWk/TWtdprffTWn/cSqCkl/mj1vpjWutmrXW91vo9WutbnASKUuoYpdT/U0rtUErtUkr9RSl1pVIqXMx4lVLa4e+ZYvdfCCZunU8W670+63z6iFvjFfyNG+eTUmqaUupWpdQflFJ9SqndSqm/KaWeTa+r0ct9EPyFS+dUTCl1mVLqYaXUNqXUHqXU20qpR5VSZxU1nqBaUipBunPyf5NyEd1LylLzCdIuIq31PxSxLg28Btxl8XGP1vrHZQ9Y8DVunk85610EPEMqDqsJOFlr/ZjzUkLQcet8UkodCGwiZYXeDLxFqizDh0jVntoEfFBr/Y7LuyD4DBfPqf8gVRV+K7CWVDjGAcBZpBJfbtFaf7GgdYlIsSYdbPsKqR9r2cG2aZGyVmt9ojcjFvyM2+dT1nqnAM+RqvHzKnABIlJqHjfPp/QTckhrPa6Gv1LqHuDTwJe11t90cRcEn+HyOXUW8LbWem3O9MNJPVBNAxZrrZ/Lt67AunsqgHROFtzEq/Pp66Tq/FwEJMscoxAcXDuftNYJK4GS5v+kXw8tdaBCYHDznLo/V6Ckp79IykIDcGIh6wpsWfwKUHDn5CLqq0SVUp8DWklVwX1Oay3xKBMD188npdSHgCuAq7TWm5VS7oxUCAJeXJ+s+ET69S9lrEMIBpU6p0xBPFrIzCJS7HHsnKyU2gocQapz8osFrnMh8JPsCUqpF4ALtNZ/LWOsgv9x9XxSSk0nFd/0B1L9q4SJhevXp3T1bfOJeQZwPHAkqcreq8scr+B/vLjnjSHtUjob0IBlgksuIlLsKbtzcg43kwpIeplUUNJhpAKLzgF+n26I2OuwvBBs3D6fbiN1IzlR6vtMSNw+nyB1P7ghZ9p/AZekTf5CbePFOZVBpUy9PwZagP9Mu37yUtMxKenUJ6fU39y/e7wai9Z6pdb6T1rr7VrrIa31unSk9H8DM4Grvdq24A5+OZ+UUmeTCpD9ktZa+k0FFL+cTyZa691aa0XqvtBBKs7pI8C6dAaQ4HP8dk7l8B3gH0hZfwvK7IHat6S8SspqUSh9Wf9XqnPyD0mZv5aUuR7Be6p+PimlZpA6Z34H/KCIsQj+o+rnkxVpy1wv8FOlVDfwNPB99jV2FfyLL88ppdQ3gatIxbacWkxMS02LFK31h8tYvBtYTKpz8pg0KeVu5+S30q9SMMnn+OR8mk3K8vZhIGkTLPtoevpVWutbyxiz4CE+OZ8c0Vo/o5QyKDATQ6gufjynlFK3AFeSim36uNZ6pJjla9rdUyZud0624+j0q5jtaxu3zqe3SQVfW/1tTs/zcPp9V5ljFvxLRa5PSqmppGpaFJSJIQQaV88pleJ2UgLlUVIWlKIECgBaa/mz+CP1w3yLVBXPxVnTpwB/IhWdfG7OMg2kAmJn50xfAEQstrEA2J5e13nV3mf5C8b55LCNu9Lr+Ui191f+vP1z+fr0HmCKxTYmAT9Nr2tNtfdZ/gJ1TilSGWEa+H9W51ehf1Jx1gFVZOdkpdSJpExaYyrLKqXuIlVv4A/AG6ROgsNIKdYwqS/zn7V8GTWNW+eTw/rvAi5EKs5OCFy8Pt0KfBb4I6nWHQbQDvw9qZpO3cBJWut+T3dIqDounlM3ADcCu4Bbgb0Wm9ugtX4g35hqOialXLTWDyilTiDVOfls9nVO/iKpzsmFiooHSKnUBaQK5kwhZbZ/GFittf6Vy0MXfIiL55MguHk+/R9SPZ8+mP6bCrxDqmfPd0ilixZvphcCh4vn1EHp13rgX2zm+Smpe6MjYkkRBEEQBMGXSOCsIAiCIAi+RESKIAiCIAi+RESKIAiCIAi+RESKIAiCIAi+RESKIAiCIAi+RESKIAiCIAi+RESKIAiCIAi+RESKIAiCIAi+RESKIAiCIAi+RESKIAiuo5T6rFJKF/D3eLXHKgiCf5HePYIgeMGi9OvzwF8c5nukAmMRBCGgiEgRBMELTJFyi9b6v6o6EkEQAou4ewRBcBWlVAhYmH7bVc2xCIIQbESkCILgNvOARiABvFjlsQiCEGBEpAiC4Damq+dVrfXuqo5EEIRAIyJFEAS3MUWKuHoEQSgLESmCILiNKVLOypN+PM9cQCl1Q3rau1YrVEo9kf78qZzpN+SsM66U2qaU+l9KKeUw37BSaqNS6ttKqRZPjoIgCGUj2T2CILjNkenXh4DtNvPsBV7Oer8QeBeYqpTaX2v9hvmBUupU4INAHHghZz0LScW9fC79vhG4EPgK8CpwZ9Z83cBF6fdTgfcBVwDnK6WO0VpvKXgPBUGoCEprXe0xCIJQIyilDgZeSb+dpbW2Eym5y70K/Ak4D/io1vqR9PQQsIGUODkf+LzW+kc5y/1Ra/2ZrGmTgV3A97TWV2bN94zW+tM5252bXvdvtNZnFr3DgiB4irh7BEFwE9PV80YRAqUJOAh4lpT14/Csj88HDgQeTL9/wWK5v+assh1QwOY886G1fpmUxecUpZRYlgXBZ4hIEQTBTUyR8lwRyywgJSpeIBVsezhkLCI3Ad8iJTySjBUa5nIvKqXqlFJTlFKLgP9Kz3dnznzjREqaV4EpQHMRYxYEoQKISBEEwU1MkbK+iGXMwm9/IUukAJeQEg83p+fZorUetljuQVLxKrtIiaM4sERrPZIzn51IiaRfh4oYsyAIFUBEiiAIbvLe9GsxImUBKfeQQUpIHK6UmkYq+PWmtDBZwPig2QWkAnPfn/47BricVJDt13Pm26m1ft1m+wcDvVrrXUWMWRCECiA+WEEQXEEptT8wK/32n5RS/+Aw+//SWpsBtgvZ14SwK72ObwKDwGqlVBiYDzyQs46FwF+01uuypj2tlHo/cJFS6hKdygxYiI0VRSk1BTjeYt2CIPgAESmCILjFoqz/z3CYL0nK4kG6lsl7gO+lP9sM7AH+GThXax1XSs0n5fbJdFPOWu4Oi/XvAiYBIaVUMj2fXZPDL5CKRbndaccEQagOIlIEQXAFrfX/JRWgWgxzgCbSAkRrPaqU+jEpkfG/0/MsSL++YLHcGAuJUqoB+Hvgz1rrRDoletx86XlPBf4dWKW1LibQVxCECiEiRRCEapIdNAuA1vpSi3l2aq23WSwXVkodTSq+7kBSxdlaSaUuZ88XSs83GZgNnEnK2vNfpK06giD4DxEpgiBUk4XAbsZWn81lAeMtIab4+GH6NQn8DXgC+Cet9V9z5vt++vVdoJ9U4bgTtNZ/KHnkgiB4jlScFQRBEATBl0gKsiAIgiAIvkREiiAIgiAIvkREiiAIgiAIvkREiiAIgiAIvkREiiAIgiAIvkREiiAIgiAIvkREiiAIgiAIvkREiiAIgiAIvkREiiAIgiAIvuT/B6oD3VL9JU3PAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 582.492x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "font = {'family' : \"Times New Roman\",\n",
    "        'weight' : 'normal',\n",
    "        'size'   : 20}\n",
    "\n",
    "matplotlib.rc('font', **font)\n",
    "#plt.rcParams[\"font.family\"] = \"Times New Roman\"\n",
    "\n",
    "gold = (1+(5**0.5))/2\n",
    "plt.figure(figsize=(5*gold,5))\n",
    "\n",
    "plt.scatter(properties[:,1], properties[:,3], c = 'lightgrey', label = 'QM7-X')\n",
    "index = [0,1,2,4]\n",
    "mols = pd.read_csv('/work/projects/tcp/lmedrano/alessio/vae/mol-new/code/fhi/extract/xyzfiles2/12/data-12.dat', header = None).values[index]\n",
    "print(mols)\n",
    "for i in range(0, len(mols)):\n",
    "    if i == 0:\n",
    "        plt.scatter(float(mols[i][0].split()[2]), float(mols[i][0].split()[3]), c = 'turquoise', s = 50, label = r'samples $T_1$')\n",
    "    else:\n",
    "        plt.scatter(float(mols[i][0].split()[2]), float(mols[i][0].split()[3]), c = 'turquoise', s = 50)\n",
    "        \n",
    "plt.scatter(-1*p_stds[1] + p_means[1], 2*p_stds[3] + p_means[3], c = 'navy', marker = '*', s = 70, label = r'$T_1$')\n",
    "\n",
    "mols = pd.read_csv('/work/projects/tcp/lmedrano/alessio/vae/mol-new/code/fhi/extract/xyzfiles/set2/data-set2.dat', header = None).values\n",
    "print(mols)\n",
    "for i in range(0, len(mols)):\n",
    "    if i == 0:\n",
    "        plt.scatter(float(mols[i][0].split()[2]), float(mols[i][0].split()[3]), c = 'orange', s = 50, label = r'samples $T_2$')\n",
    "    else:\n",
    "        plt.scatter(float(mols[i][0].split()[2]), float(mols[i][0].split()[3]), c = 'orange', s = 50)\n",
    "                \n",
    "plt.scatter(-1.5*p_stds[1] + p_means[1], 1*p_stds[3] + p_means[3], c = 'red', marker = '*', s = 70, label = r'$T_2$')\n",
    "plt.legend(loc = 'lower right', prop = {'size': 14})\n",
    "plt.xlim(-.5, -.2)\n",
    "plt.ylim(80, 120)\n",
    "plt.locator_params(axis='x', nbins=4)\n",
    "plt.locator_params(axis='y', nbins=4)\n",
    "plt.xlabel(r'$E_{MBD}$', fontdict={'size':22})\n",
    "plt.ylabel(r'$\\alpha$', fontdict={'size':22})\n",
    "plt.tight_layout()\n",
    "plt.savefig('./new_chem.pdf', dpi = 100)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.13"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}