{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example 7: Rotational barrier prediction in KRAS inhibitors\n", "\n", "## Objective\n", "The objective of this task is to use the `molli` workflow to transition from `.CDXML` molecular drawings to full-fledged 3D structures with subsequent CREST//GFN2-XTB conformer generation, ORCA equilibrium structure minimization, as well as ORCA relaxed surface scan.\n", "This workflow is useful to assess semi-quantitatively the possibility of hindered rotation in a number of molecules (see e. g. [*J. Med. Chem.* **2020**, *63*, 1, 52](https://doi.org/10.1021/acs.jmedchem.9b01180))\n", "\n", "## Prerequisites\n", "- `molli` installation\n", "- `pandas` installation (it is NOT a dependency of molli and must be installed separately)\n", "- `OpenBabel` installation (it is NOT a dependency of molli and must be installed separately)\n", " Note: if you experience problems with openbabel refusing to optimize the structures, do not forget to set the `BABEL_DATADIR` environment variable to point to a valid location with openbabel data files. With conda installations, it seems to be consistently pointing to `$CONDA_PREFIX/share/openbabel`\n", "- `barriers.csv` file that has the experimental rotational barriers listed\n", "- `cladosporin.cdxml` file that contains the structural drawings\n", "\n", "## Methods\n", "\n", "### Conformer generation\n", "Rotational barrier profiles are generated using b97-3c method, as reported by Grimme et al." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import molli as ml\n", "from molli.external import openbabel\n", "from math import degrees\n", "\n", "import pandas as pd\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from scipy.signal import find_peaks\n", "\n", "# This configures on-screen preview of molecules\n", "ml.visual.configure(bgcolor=\"black\", theme=\"dark\", width=500)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 1. `CDXML` Parsing\n", "\n", "In this step, all of the molecules are read from the CDXML file and converted into their 3D versions.\n", "This can be done from both the command line interface, but for the sake of exercise, this workflow will show the script that utilizes the key function of `ml.Molecule.join(...)`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To make sure that the parsing step was successful, we will perform a basic inspection of the resulting `.mlib` file" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Dihedral 117.02 for Molecule(name='23', formula='C28 H29 Cl1 F1 N5 O2')\n", "Dihedral 113.62 for Molecule(name='25', formula='C27 H27 Cl1 F1 N5 O2')\n", "Dihedral 134.54 for Molecule(name='28', formula='C24 H24 Cl1 F1 N6 O2 S1')\n", "Dihedral 144.75 for Molecule(name='27', formula='C24 H25 Cl1 F1 N7 O2')\n", "Dihedral 131.87 for Molecule(name='18', formula='C27 H27 Cl1 F1 N5 O2')\n", "Dihedral 128.08 for Molecule(name='29', formula='C25 H27 Cl1 F1 N7 O2')\n", "Dihedral 111.87 for Molecule(name='24', formula='C28 H29 Cl1 F1 N5 O2')\n", "Dihedral 139.62 for Molecule(name='26', formula='C26 H26 Cl1 F1 N6 O2')\n", "Dihedral 126.61 for Molecule(name='22', formula='C27 H25 Cl1 F1 N5 O2')\n" ] } ], "source": [ "cdxf = ml.CDXMLFile(\"rot_barrier.cdxml\")\n", "core = cdxf[\"core\"]\n", "core_ap = core.get_attachment_points()[0]\n", "names = set(cdxf.keys())\n", "names.remove(\"core\")\n", "core.add_implicit_hydrogens()\n", "\n", "# This is the destination \n", "lib = ml.MoleculeLibrary(\"kras_inhibitors.mlib\", overwrite=True, readonly=False)\n", "\n", "with lib.writing():\n", " for n in names:\n", " sub = cdxf[n]\n", " sub_ap = sub.get_attachment_points()[0]\n", " \n", " mol = ml.Molecule.join(\n", " core,\n", " sub,\n", " core_ap,\n", " sub_ap,\n", " name=n,\n", " )\n", "\n", " mol.add_implicit_hydrogens()\n", "\n", " openbabel.obabel_optimize(\n", " mol,\n", " ff=\"MMFF94\",\n", " coord_displace=0.01,\n", " inplace=True,\n", " )\n", "\n", " # just a little extra: calculate the preliminary dihedral angle)\n", " d = mol.dihedral(\"1\", \"2\", \"3\", \"4\")\n", " print(f\"Dihedral {degrees(d):0.2f} for {mol}\")\n", "\n", " lib[n] = mol" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "application/3dmoljs_load.v0": "
\n

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n jupyter labextension install jupyterlab_3dmol

\n
\n", "text/html": [ "
\n", "

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n", " jupyter labextension install jupyterlab_3dmol

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/3dmoljs_load.v0": "", "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/3dmoljs_load.v0": "", "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%mlib_view kras_inhibitors.mlib 22" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will now import the molecule library into our python code. the success of this step will be demonstrated by rendering the parsing result of an arbitrary molecule. Give it a try and substitute the name with any names that were obtained in the inspect stage!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 2. XTB relaxed surface scan.\n", "\n", "In this step we will be attempting to generate the guess structure of the transition state of the hindered rotation by performing the *Relaxed Surface Scan* of the rotation." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Submitting jobs: 0it [00:00, ?it/s]\n", "Waiting for jobs: 0it [00:00, ?it/s]\n", "Finalizing the calculations: 100%|███████████████| 9/9 [00:00<00:00, 180.65it/s]\n" ] } ], "source": [ "# Jupyter does not play nice with the Thread parallelism, so we placed the workflow in a separate file.\n", "!python scan.py " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "9 2 43 44 70.75379386832473\n", "Minima: [], maxima: [26] 255.47593799748057\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAGwCAYAAABPSaTdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABkd0lEQVR4nO3dd3hT5d8G8PskbdK996C7ILTsWZW9FRUcoKigiIiCIiAKqIAKKKKgIG5BFAT9CY4XBdkIyCpUyi7QQksX3TttkvP+URqorKRNejLuz3XlojnnJLkzaL59zjMEURRFEBEREVkpmdQBiIiIiEyJxQ4RERFZNRY7REREZNVY7BAREZFVY7FDREREVo3FDhEREVk1FjtERERk1eykDmAOtFotMjMz4erqCkEQpI5DREREehBFEaWlpQgKCoJMdvP2GxY7ADIzMxEaGip1DCIiImqA9PR0hISE3HQ/ix0Arq6uAGpfLDc3N4nTEBERkT5KSkoQGhqq+x6/GRY7gO7UlZubG4sdIiIiC3O7LijsoExERERWjcUOERERWTUWO0RERGTVWOwQERGRVWOxQ0RERFaNxQ4RERFZNRY7REREZNVY7BAREZFVY7FDREREVo3FDhEREVk1FjtERGRetmwBWras/ZfICFjsEBGR+RBFYMYM4OTJ2n9FUepEZAVY7BARkfn46y/g4MHanw8erL1O1EgsdoiIyDyIIvDGG9DKar+atDIZ1DNmsnWHGo3FDhERmYcrrToyrRYAINNqYXc4ETuXfA+tlgUPNRyLHSIikp4oQvv669AI9b+W1IIMHvPfxojP/8GZnFKJwpGlY7FDRETS++svyA4dglzU1ttsJ2rRJjsFDtu3YPBHf2PBxlOorNZIFJIsFYsdIiKSliii8tXp0EK48W5BhjmH1kKt0WLZjnPov3gndpzObeKQZMkkLXbmz5+PTp06wdXVFX5+fnjggQdw+vTpeseMHj0agiDUu3Tt2rXeMSqVChMnToSPjw+cnZ1x3333ISMjoymfChERNZC6sgpV59Mgw4375QiiFhGVBfhyRDyC3B2QXlCJ0csP4oXVh5FbUtXEackSCaIoXTf3gQMHYsSIEejUqRPUajVmzpyJ5ORknDhxAs7OzgBqi52cnBwsX75cdzuFQgEvLy/d9fHjx+P333/HihUr4O3tjSlTpqCgoACJiYmQy+W3zVFSUgJ3d3cUFxfDzc3N+E+UiIhu6rOd5/Dt2r/RTFOGTx/vAC9n5fUH+fkBISEoV6mxaPMZLN+bBo1WhKvSDq8MbI6RXcIgl924ZYisl77f35IWO/91+fJl+Pn5YefOnejevTuA2mKnqKgIv/zyyw1vU1xcDF9fX3z33XcYPnw4ACAzMxOhoaH4448/MGDAgOtuo1KpoFKpdNdLSkoQGhrKYoeIqImdv1yGQR/9DZVai/cfao2HO4bqdbvjmcWYsf4Y/k0vAgC0CXHH3KHxiAt2N2FaMjf6Fjtm1WenuLgYAOq12gDAjh074Ofnh9jYWIwdOxa5uVfP1SYmJqKmpgb9+/fXbQsKCkJcXBz27t17w8eZP38+3N3ddZfQUP3+cxERkfFotSJe+zkZKrUWd8f44KEOIXrftlWQO9aNT8Db97eCq9IO/2YUY+iyPUjhiC26AbMpdkRRxOTJk3HXXXchLi5Ot33QoEFYtWoVtm3bhg8++AAHDx5E7969dS0z2dnZUCgU8PT0rHd//v7+yM7OvuFjTZ8+HcXFxbpLenq66Z4YERHd0Kr9F3AgrQBOCjnmDY2HIBh2GkouE/BEt3BsndIDccFuqNGI2HnmsonSkiWzkzpAnQkTJuDo0aPYvXt3ve11p6YAIC4uDh07dkRYWBg2bNiAYcOG3fT+RFG86X8cpVIJpfIG54SJiKhJXCqqxLt/ngIAvDqwBUK9nBp8X35uDujd3A/HLpXgbG6ZsSKSFTGLlp2JEyfit99+w/bt2xEScutmzMDAQISFhSElJQUAEBAQgOrqahQWFtY7Ljc3F/7+/ibLTEREDSOKImauT0Z5tQYdwzzxRNewRt9ntL8rACCFxQ7dgKTFjiiKmDBhAtatW4dt27YhIiLitrfJz89Heno6AgMDAQAdOnSAvb09Nm/erDsmKysLx44dQ0JCgsmyExFRw6w/cgk7Tl+Gwk6G9x5qDZkRRlHF+LkAAFJySmFG427ITEh6GuuFF17A6tWr8euvv8LV1VXXx8bd3R2Ojo4oKyvD7Nmz8eCDDyIwMBBpaWmYMWMGfHx8MHToUN2xY8aMwZQpU+Dt7Q0vLy9MnToV8fHx6Nu3r5RPj4iI/uNyqQpv/d8JAMBLfWIQ5etilPuN8HGGTABKqtS4XKqCn5uDUe6XrIOkxc6nn34KAOjZs2e97cuXL8fo0aMhl8uRnJyMlStXoqioCIGBgejVqxfWrl0LV1dX3fGLFi2CnZ0dHnnkEVRWVqJPnz5YsWKFXnPsEBFR05n923EUVdSgVZAbnu0eabT7dbCXI8zbGal55TibW8Zih+qRtNi5XVOjo6MjNm3adNv7cXBwwJIlS7BkyRJjRSMiIiPbeCwbG5KzIJcJeO/B1rCXG7cnRbSfC1LzypGSW4aEaB+j3jdZNrPooExERNatuKIGb/x6DAAwrnukSSb/i67rt5PLuXaoPhY7RERkcnP/OIHLpSpE+jrjxT4xJnmMq52UOSKL6mOxQ0REJrU7JQ8/HsqAIAALHmwNB3vT9KeM8avty8m5dui/WOwQEZHJlKvUeG3dUQDAqG7h6BjudZtbNFyUX+0C0vnl1cgvU93maLIlLHaIiMhklu04i4zCSgR7OOKVAc1N+lhOCjuEeDoCYOsO1cdih4iITCK/TIXle9IAAG8OaQlnpekHAOv67bDYoWuw2CEiIpP4fNd5VFRrEB/sjv4tm2b5nhh/9tuh67HYISIio8strcLKf9IAAJP7xRq8onlD1Q0/Z7FD12KxQ0RERvfZjvOoqtGiXTMP9Gzu22SPG8O5dugGWOwQEZFR5ZRU4fv9FwA0basOcLVlJ6dEheLKmiZ7XDJvLHaIiMioPtl+FtVqLTqFe+KuJl62wdXBHgFX1sXiqSyqw2KHiIiM5lJRJdYcSAcAvNzErTp1Yvzr+u3wVBbVYrFDRERG88n2s6jWaNEt0hsJUdIsxhnNZSPoP1jsEBGRUaQXVODHg1dbdaRSt2wE59qhOix2iIjIKJZsS4FaK+LuGB90jjDdshC3c/U0FosdqsVih4iIGi0trxw/H74EQNpWHQCI9q0tdi4VVaJMpZY0C5kHFjtERNRoH29LgUYroldzX7Rv5ilpFk9nBXxclACAc2zdIbDYISKiRjp3uQy/HDGPVp06MZxJma7BYoeIiBrloy0p0IpA3zv80TrEQ+o4AK7222EnZQJY7BARUSOcySnF70czAQCT+sZInOaqqy07nGuHWOwQEVEjfLQlBaIIDGwVgLhgd6nj6ET5sWWHrmKxQ0REDXIyqwQbkrMgCObTV6dO3Vw7FwsqUFWjkTgNSY3FDhERNciizWcAAPfEB6J5gKvEaerzcVHAw8keoljbgZpsG4sdIiIyWHJGMf46kQNBMK++OnUEQeCILNJhsUNERAZbvKW2Vef+NkGI9jOvVp06dbm4Rhax2CEiIoMkpRdh66lcyGUCXuprXn11rhWj66TMEVm2jsUOEREZ5MMrfXWGtgtGhI+zxGlujmtkUR0WO0REpLdDaQXYdeYy5DIBL/Y2v74616obkZWWX4FqtVbiNCQlFjtERKS3RVf66jzcIQTNvJ0kTnNr/m5KuCrtoNGKSMsvlzoOSYjFDhER6WX/+XzsOZsPe7mACb2jpY5zW4IgXJ1ckJ2UbRqLHSIi0ssnO84BAB7uGIoQT/Nu1anDTsoEsNghIiI9HLtUjF1nLkMmAM91j5I6jt64ICgBLHaIiEgPn+2sbdW5t3WQ2ffVuVZdJ+WzPI1l01jsEBHRLaXlleOP5CwAwPieltOqAwDRV05jnc8rg1rDEVm2isUOERHd0hd/n4dWBHo198UdgW5SxzFIsIcjHO3lqNGIuFBQIXUckoidIQeLooidO3fi77//RlpaGioqKuDr64t27dqhb9++CA0NNVVOIiKSQG5JFf53KAMAML6n+Y/A+i+ZTEC0nwuSLxXjbG4ZonxdpI5EEtCrZaeyshLz5s1DaGgoBg0ahA0bNqCoqAhyuRxnz57FrFmzEBERgcGDB2Pfvn2mzkxERE3k6z2pqNZo0SHME53CPaWO0yBcEJT0atmJjY1Fly5d8Nlnn2HAgAGwt7e/7pgLFy5g9erVGD58OF5//XWMHTvW6GGJiKjpFFfWYNW+iwCA8T2iIAiCxIkaJrpuRFYOh5/bKr2KnT///BNxcXG3PCYsLAzTp0/HlClTcOHCBaOEIyIi6Xy/7wLKVGrE+rugdws/qeM0WN2ILA4/t116nca6XaFzLYVCgZgY814vhYiIbq2qRoPle1IBAM/1iIJMZpmtOsDVEVlnc8ug0YoSpyEp6NWyc/ToUb3vsHXr1g0OQ0RE5uGnxAzklVUj2MMRQ9oESR2nUUI9HaGwk0Gl1uJSYaVFzRNExqFXsdO2bVsIggBRvHFFXLdPEARoNBqjBiQioqal1mjxxa7aSQSf7R4Je7llz1JiJ5ch0scZp7JLkZJbymLHBulV7KSmppo6BxERmYkNyVlIL6iEl7MCj3S0jilFYvxdrxQ7Zehzh7/UcaiJ6VXshIWFmToHERGZAVEU8emVBT+fSgiHo0IucSLjiOHq5zbNoEkF65w7dw6LFy/GyZMnIQgC7rjjDrz00kuIirKsacSJiKi+Hacv41R2KZwVcjzZLVzqOEZzda4dDj+3RQafiN20aRNatmyJAwcOoHXr1oiLi8P+/fvRqlUrbN682RQZiYioidS16jzWpRncna6fU81S1a1+fja37Kb9T8l6Gdyy89prr+Hll1/Gu+++e932V199Ff369TNaOCIiajqH0gpwIK0A9nIBY+6KlDqOUYV5O8NOJqC8WoOs4ioEeThKHYmakMEtOydPnsSYMWOu2/7000/jxIkTRglFRERNr65V58H2IQhwd5A4jXHZy2WI8HEGwMkFbZHBxY6vry+SkpKu256UlAQ/P8udYZOIyJadyi7B1lO5EITa4ebWKIbLRtgsg09jjR07Fs8++yzOnz+PhIQECIKA3bt347333sOUKVNMkZGIiEzs853nAQCD4gIQaaUrg0f7ckFQW2VwsfPGG2/A1dUVH3zwAaZPnw4ACAoKwuzZs/Hiiy8aPSAREZlWekEFfvs3E0Dt0hDWKtqfa2TZKoOLHUEQ8PLLL+Pll19GaWltU6Crq6vRgxERUdP46u/z0GhF3BXtg9YhHlLHMZmrc+2U6mb9J9vQqDnAXV1dWegQEVmwvDIV1hxMBwCM72m9rToAEOHjDJkAlFSpcblUJXUcakIGFzv5+fl44YUX0LJlS/j4+MDLy6vehYiILMeKPWlQqbVoHeKOhChvqeOYlIO9HGHeHJFliww+jfX444/j3LlzGDNmDPz9/RvVDDh//nysW7cOp06dgqOjIxISEvDee++hefPmumNEUcScOXPwxRdfoLCwEF26dMEnn3yCVq1a6Y5RqVSYOnUqfvjhB1RWVqJPnz5YtmwZQkJCGpyNiMjalVbVYOU/aQCA53tG2cRpnWg/F6TmleNsbhnujPaROg41EYOLnd27d2P37t1o06ZNox98586deOGFF9CpUyeo1WrMnDkT/fv3x4kTJ+DsXFt9L1iwAB9++CFWrFiB2NhYvPPOO+jXrx9Onz6tO4U2adIk/P7771izZg28vb0xZcoU3HvvvUhMTIRcbh3ruhARGdsPBy6ipEqNSF9n9G8ZIHWcJhHj54LNJ3KQwmUjbIrBxU6LFi1QWVlplAffuHFjvevLly+Hn58fEhMT0b17d4iiiMWLF2PmzJkYNmwYAODbb7+Fv78/Vq9ejXHjxqG4uBhff/01vvvuO/Tt2xcA8P333yM0NBRbtmzBgAEDjJKViMiaqNQafPV3KgDgue5RkMmsv1UHuHauHZ7GsiUG99lZtmwZZs6ciZ07dyI/Px8lJSX1Lo1RXFwMALq+P6mpqcjOzkb//v11xyiVSvTo0QN79+4FACQmJqKmpqbeMUFBQYiLi9Md818qlcqouYmILM1PhzKQW6qCv5sS97cLkjpOk4nxqz0jwLl2bIvBLTseHh4oLi5G7969622vG8an0WgaFEQURUyePBl33XUX4uLiAADZ2dkAAH9//3rH+vv748KFC7pjFAoFPD09rzum7vb/NX/+fMyZM6dBOYmILF21Wotl288CqJ1XR2lnO6f7o3xdIAhAfnk18stU8HZRSh2JmoDBxc7IkSOhUCiwevXqRndQvtaECRNw9OhR7N69+7p9/30MfeZHuNUx06dPx+TJk3XXS0pKEBoa2oDURESW56fEdGQWV8HXVYlHOzeTOk6TclTIEezhiIzCSpzNLWOxYyMMLnaOHTuGI0eO1Bsx1VgTJ07Eb7/9hl27dtUbQRUQUNthLjs7G4GBgbrtubm5utaegIAAVFdXo7CwsF7rTm5uLhISEm74eEqlEkolP+BEZHtqW3VqF/wc3yMKDva206pTJ8bPBRmFlUjJLUOXSOsebk+1DO6z07FjR6SnpxvlwUVRxIQJE7Bu3Tps27YNERER9fZHREQgICAAmzdv1m2rrq7Gzp07dYVMhw4dYG9vX++YrKwsHDt27KbFDhGRrfpfYgYuFVXC11WJx7rYVqtOnRh/9tuxNQa37EycOBEvvfQSXnnlFcTHx8Pe3r7e/tatW+t9Xy+88AJWr16NX3/9Fa6urro+Nu7u7nB0dIQgCJg0aRLmzZuHmJgYxMTEYN68eXBycsJjjz2mO3bMmDGYMmUKvL294eXlhalTpyI+Pl43OouIiGpbdT65pq+OLbbqALVz7QDg8HMbYnCxM3z4cADA008/rdsmCEKDOih/+umnAICePXvW2758+XKMHj0aADBt2jRUVlbi+eef100q+Ndff9VbpmLRokWws7PDI488optUcMWKFZxjh4joGj8frm3V8XFRYqSNtuoA166RxZYdWyGIoigacoO6UVA3ExYW1qhAUigpKYG7uzuKi4vh5uYmdRwiIqOrVmvR+4MdyCisxOv33IFn7o6UOpJkSqtqED/7LwDAv7P6w93R/ja3IHOl7/e33i07M2bMwAMPPIDOnTsbJSARETWddYczkFFY16pjeX+UGpOrgz0C3R2QVVyFs7ll6BDmefsbkUXTu4NyVlYW7r33XgQGBuLZZ5/Fhg0boFJx1VgiInNXo9Fiqa6vTiQcFTzFX9dv5yz77dgEvYud5cuXIycnBz/++CM8PDwwZcoU+Pj4YNiwYVixYgXy8vJMmZOIiBroaquOwuZbderUzaTMfju2waCh54Ig4O6778aCBQtw6tQpHDhwAF27dsWXX36J4OBgdO/eHQsXLsSlS5dMlZeIiAxwbavOuO5RbNW5QrdGFoef2wSD59m51h133IFp06Zhz549SE9Px6hRo/D333/jhx9+MFY+IiJqhPWHLyG9oBLezgqM7Gq7I7D+6+ppLBY7tsDgoec34+fnhzFjxmDMmDHGuksiImqEeq06PSLhpDDar3yLF+1bW+xcKqpEmUoNFyVfG2um97s7bNiw29+ZnR0CAgLQr18/DBkypFHBiIiocX45cgkXCyrg7azA413ZV+dans4K+LgokVemwunsUo7IsnJ6n8Zyd3e/7cXR0REpKSkYPnw43nzzTVPmJiKiW1Bf06rzbHe26txIu2YeAIBDaQXSBiGT0/vTv3z5cr3vdMOGDRg/fjzeeuutBoUiIqLG+SUpExfyK+DlrMAT3diqcyNdIryw+UQODqQWYFyPKKnjkAkZ3EF569atN923dOlSAMCdd96Jjh07NjwVERE1mFqjxZJtKQDYqnMrXSJqVzw/kFYAjdagxQTIwhhc7Dz44IM4ePDgddsXL16MGTNmAAA8PDywbt26xqcjIiKD/Xptqw776tzUHYGucFHaobRKjZNZJVLHIRMyuNhZtGgRBg8ejBMnTui2LVy4ELNmzcKGDRuMGo6IiAxzbavO2Lsj4cxRRjdlJ5ehY3htx+QDqey3Y80M/l/w1FNPIT8/H/3798fu3buxdu1azJs3D3/++ScSEhJMkZGIiPT027+ZSMuvgKeTPZ5kX53b6hzhhR2nL2N/aj6evitC6jhkIg0q+adOnYr8/Hx07NgRGo0Gf/31F7p06WLsbEREZIDaVp3aEVhju7NVRx+6fjupBRBFEYIgSJyITEGv/wkff/zxddsCAwPh5OSE7t27Y//+/di/fz8A4MUXXzRuQiIi0svvRzORmlcODyd7PNktXOo4FiE+2B0O9jIUVtQgJbcMsf6uUkciE9Cr2Fm0aNENt8vlcuzZswd79uwBULt2FosdIqKmp9GKWLL1SqvO3ZGcEVhPCjsZOoR5Ys/ZfOxPLWCxY6X0+t+Qmppq6hxERNQIv/+bifNXWnVGJYRLHceidA73ri12zudz9JqVatRCoEREJD2NVsTH14zAYquOYbpEegG42m+HrI9exc67776L8vJyve5w//79HIJORNSEfj6cgfOXy+HuyBFYDdE21AMKuQy5pSqk5VdIHYdMQK9i58SJEwgLC8P48ePx559/4vLly7p9arUaR48exbJly5CQkIARI0bAzc3NZIGJiOiqymoNPvjrNABgQq9ouDrYS5zI8jjYy9E21AMAcCA1X9owZBJ6FTsrV67Etm3boNVqMXLkSAQEBEChUMDV1RVKpRLt2rXDN998g9GjR+PUqVO4++67TZ2biIgAfLMnFTklKgR7OHINrEboHFF7Kmv/eU4uaI30PrHbunVrfP755/jss89w9OhRpKWlobKyEj4+Pmjbti18fHxMmZOIiP4jv0yFT3ecAwBMG9gcDvZyiRNZri6RXli6HdjPmZStksG92ARBQJs2bdCmTRtT5CEiIj0t2XYWZSo14oLdMKR1kNRxLFr7Zp6QywRcKqpERmEFQjydpI5ERsTRWEREFigtrxzf77sAAJgx6A7IZJz5tzGclXaID3YHwHWyrBGLHSIiC/T+ptNQa0X0bO6LhGh2IzCGLuy3Y7VY7BARWZgjFwuxITkLggC8NqiF1HGshm6+nTQWO9aGxQ4RkQURRRHz/jgJAHiofQhaBHCqD2PpEOYFQQBS88qRW1IldRwyIoOLnRUrVqCigpMuERFJYfOJHBxMK4SDvQyT+8dKHcequDvao2VgbfG4j/12rIrBxc706dMREBCAMWPGYO/evabIREREN6DWaPHuxlMAgDF3RSDQ3VHiRNanS4Q3AE4uaG0MLnYyMjLw/fffo7CwEL169UKLFi3w3nvvITs72xT5iIjoirWH0nH+cjm8nBUY1yNK6jhWiZMLWieDix25XI777rsP69atQ3p6Op599lmsWrUKzZo1w3333Ydff/0VWq3WFFmJiGxWuUqNRZtrF/t8sXc03LgshEnUFTspuWXIL1NJnIaMpVEdlP38/HDnnXeiW7dukMlkSE5OxujRoxEVFYUdO3YYKSIREX3593nklakQ5u2Ex7pwWQhT8XJWINbfBQBwkKOyrEaDip2cnBwsXLgQrVq1Qs+ePVFSUoL/+7//Q2pqKjIzMzFs2DCMGjXK2FmJiGxSbmkVvth1HgAwbUALKOw4kNaU6vrtcOkI62Hw/5ghQ4YgNDQUK1aswNixY3Hp0iX88MMP6Nu3LwDA0dERU6ZMQXp6utHDEhHZosVbUlBRrUHbUA8Mjg+QOo7VY78d62Pw2lh+fn7YuXMnunXrdtNjAgMDkZqa2qhgREQEnM0txdqDtX88zhh8BwSBy0KYWt1MyiezS1BcWQN3R/aPsnQGt+x8/fXXtyx0gNrFQsPCeE6ZiKix3tt4GhqtiH4t/XUtDmRafm4OiPBxhigCh9hvxyoY3LLz8ccf33C7IAhwcHBAdHQ0unfvDrlc3uhwRES27EBqATafyIFcJuDVgVwWoil1ifBCal45DqQWoM8d/lLHoUYyuNhZtGgRLl++jIqKCnh6ekIURRQVFcHJyQkuLi7Izc1FZGQktm/fjtDQUFNkJiKyetcuCzG8Uyii/VwkTmRbOkd4Yc3BdM6kbCUMPo01b948dOrUCSkpKcjPz0dBQQHOnDmDLl264KOPPsLFixcREBCAl19+2RR5iYhswp/HspGUXgQnhRyT+sZIHcfmdImsHZF17FIxylVqidNQYxlc7Lz++utYtGgRoqKuzt4ZHR2NhQsXYvr06QgJCcGCBQuwZ88eowYlIrIV1Wot3ruyLMTYuyPh5+ogcSLbE+zhiGAPR2i0IhIvFEodhxrJ4GInKysLavX1Va5ardYtGREUFITS0tLGpyMiskGr91/AhfwK+Lgo8Wz3SKnj2KwukbUdwg/wVJbFM7jY6dWrF8aNG4cjR47oth05cgTjx49H7969AQDJycmIiIgwXkoiIhtRXFmDj7edBQC83C8GzkqDu1aSkdQNQd/PRUEtXoOGnnt5eaFDhw5QKpVQKpXo2LEjvLy88PXXXwMAXFxc8MEHHxg9LBGRtVu85QwKyqsR7eeC4R05yENKdTMp/5tejKoajcRpqDEM+pNBFEWoVCr8+uuvSE9Px+nTpyGKIlq0aIHmzZvrjuvVq5fRgxIRWbvT2aVY+c8FAMCsIS1hJ+eyEFIK83aCn6sSuaUqHLlYhG5R3lJHogYyuNiJiYnB8ePH0bx583oFDhERNZwoipjz+3FotCIGtPLH3TG+UkeyeYIgoEukN37/NxP7U/NZ7Fgwg/5skMlkiImJQX4+z18SERnTxmPZ2HsuH0o7GV6/p6XUceiKulmr2UnZshncRrpgwQK88sorOHbsmCnyEBHZnMpqDd7ZUDuB4LgeUQj1cpI4EdXpeqXYOXyxENVqrcRpqKEM7ub/+OOPo6KiAm3atIFCoYCjo2O9/QUFrH6JiAzx2c5zuFRUiSB3B4zvEXX7G1CTifZzgZezAgXl1Ui+VIQOYVyfzBIZXOwsXrzYBDGIiGxTekEFPtt5DgAw856WcFRwXUFzIggCOod7YePxbOw7X8Bix0IZXOyMGjXKFDmIiGzS3A0noVJr0S3SG4PjA6SOQzfQJbK22DmQWoAXONjYIjVoXOO5c+fw+uuv49FHH0Vubi4AYOPGjTh+/LhRwxERWbPdKXnYeDwbcpmAWfe1hCAIUkeiG6jrpHworQBqDfvtWCKDi52dO3ciPj4e+/fvx7p161BWVgYAOHr0KGbNmmX0gERE1qhGo8Xs32v/QHyiaxhaBLhJnIhupkWAG9wc7FBercGJrBKp41ADGFzsvPbaa3jnnXewefNmKBQK3fZevXrhn3/+MWo4IiJrtfKfCzibWwYvZwVe7hsrdRy6BblMQKfwK0tHnOcgHEtkcLGTnJyMoUOHXrfd19eX8+8QEekhr0yFxZvPAABeGdAc7k72Eiei26lbFHQ/59uxSAYXOx4eHsjKyrpu+5EjRxAcHGyUUERE1mzBxlMoVakRF+yGR7j+lUXofGWdrINpBdBqRYnTkKEMLnYee+wxvPrqq8jOzoYgCNBqtdizZw+mTp2KJ5980qD72rVrF4YMGYKgoCAIgoBffvml3v7Ro0dDEIR6l65du9Y7RqVSYeLEifDx8YGzszPuu+8+ZGRkGPq0iIiaRFJ6EX48VPs7as59rSCXsVOyJYgLcoOTQo7iyhqczimVOg4ZyOBiZ+7cuWjWrBmCg4NRVlaGli1bonv37khISMDrr79u0H2Vl5ejTZs2WLp06U2PGThwILKysnSXP/74o97+SZMmYf369VizZg12796NsrIy3HvvvdBouEItEZkXrVbE7N9qOyUPaxfMOVssiJ1chg5hngCA/efZZcPSGDzPjr29PVatWoW33noLR44cgVarRbt27RATE2Pwgw8aNAiDBg265TFKpRIBATeee6K4uBhff/01vvvuO/Tt2xcA8P333yM0NBRbtmzBgAEDDM5ERGQqPx/OQFJ6EZwVcrw2qIXUcchAXSO98XdKHg6kFWD0nRFSxyEDGFzs1ImKikJUlOmnNd+xYwf8/Pzg4eGBHj16YO7cufDz8wMAJCYmoqamBv3799cdHxQUhLi4OOzdu/emxY5KpYJKpdJdLynhUEIiMq2Sqhq8t/E0AODFPjHwc3OQOBEZ6tpFQUVR5LxIFsTgYkej0WDFihXYunUrcnNzodXWn2Bp27ZtRgs3aNAgPPzwwwgLC0NqaireeOMN9O7dG4mJiVAqlcjOzoZCoYCnp2e92/n7+yM7O/um9zt//nzMmTPHaDmJiG7n4y0pyCtTIdLHGU+xVcAitQ5xh71cQF5ZNS4VVSLEkwu2WgqDi52XXnoJK1aswD333IO4uDiTVrbDhw/X/RwXF4eOHTsiLCwMGzZswLBhw256u9tV3NOnT8fkyZN110tKShAayhERRGQaZ3NLsWJvGgDgjSEtobBr0OT1JDGlnRxBHo64kF+BS4UsdiyJwcXOmjVr8OOPP2Lw4MGmyHNLgYGBCAsLQ0pKCgAgICAA1dXVKCwsrNe6k5ubi4SEhJvej1KphFKpNHleIiJRFDHn9xNQa0X0vcMPvZr7SR2JGiHIvbbYySyulDoKGcDgPy8UCgWio6NNkeW28vPzkZ6ejsDAQABAhw4dYG9vj82bN+uOycrKwrFjx25Z7BARNZW/TuTg75Q8KOQyvH5PS6njUCMFezoCAC4VstixJAYXO1OmTMFHH30EUWz8pEplZWVISkpCUlISACA1NRVJSUm4ePEiysrKMHXqVPzzzz9IS0vDjh07MGTIEPj4+OhmcHZ3d8eYMWMwZcoUbN26FUeOHMHjjz+O+Ph43egsIiKpVFSr8fb/nQAAPHN3BMJ9nCVORI0V5HGl2ClisWNJDD6NtXv3bmzfvh1//vknWrVqBXv7+tOcr1u3Tu/7OnToEHr16qW7XtePZtSoUfj000+RnJyMlStXoqioCIGBgejVqxfWrl0LV1dX3W0WLVoEOzs7PPLII6isrESfPn2wYsUKyOVyQ58aEZFRfbQlBRmFlQhyd8ALvaRpESfjCtEVO1USJyFDGFzseHh43HBtrIbo2bPnLVuINm3adNv7cHBwwJIlS7BkyRKjZCIiMobjmcX4ancqAOCt++PgrGzwTB9kRq6exqqQOAkZwuD/fcuXLzdFDiIiq6HRipixLhkarYhBcQHo29Jf6khkJHWnsTKLqjjXjgVp0PhHtVqNLVu24PPPP0dpae0aIZmZmSgrKzNqOCIiS/TdP2n4N6MYrko7zL6vldRxyIgC3Wsng6ys0aCwokbiNKQvg1t2Lly4gIEDB+LixYtQqVTo168fXF1dsWDBAlRVVeGzzz4zRU4iIouQWVSJ9zfVzpQ8bVAL+HOmZKviYC+Hr6sSl0tVuFRYCS9nhdSRSA8Gt+y89NJL6NixIwoLC+Ho6KjbPnToUGzdutWo4YiILM2s346jvFqD9s08MLJzM6njkAlwRJbladBorD179kChqF/NhoWF4dKlS0YLRkRkaTYey8bmEzmwkwmYP6w1ZDL257BGIR6O+De9iMWOBTG4ZUer1UKj0Vy3PSMjo96QcCIiW1JaVYNZvx0DAIzrEYnmAfx9aK2CPGpPTWay2LEYBhc7/fr1w+LFi3XXBUFAWVkZZs2aJckSEkRE5uD9TaeRU6JCuLcTJvaOkToOmVCwB2dRtjQGn8ZatGgRevXqhZYtW6KqqgqPPfYYUlJS4OPjgx9++MEUGYmIzNrhi4X4bt8FAMDcofFwsOekptYs+MoCoDyNZTkMLnaCgoKQlJSENWvWIDExEVqtFmPGjMHIkSPrdVgmIrIFNRotZqxLhigCw9oH485oH6kjkYnxNJbladCUno6Ojnjqqafw1FNPGTsPEZFF+ervVJzKLoWnkz0X+rQRIR61LTv55dWorNbAUcGWPHPXoEkFiYgIuJBfjsVbzgAAZt7TknOu2Ag3Rzs4XylwMovZumMJWOwQETWAKIp4/ZdjUKm1SIjyxoPtg6WORE1EEIRr1shisWMJWOwQETXAr0mZ+DslDwo7GeYOjecaSTYmmBMLWhQWO0REBiosr8bb/3cCAPBi72hE+DhLnIia2tUFQVnsWIIGFTtFRUX46quvMH36dBQUFAAADh8+zBmUicgmzPvjJPLLqxHr74Jnu0dJHYckwNNYlsXg0VhHjx5F37594e7ujrS0NIwdOxZeXl5Yv349Lly4gJUrV5oiJxGRWfjnXD5+SswAAMwbGg+FHRvIbRFPY1kWg/+XTp48GaNHj0ZKSgocHK6u5jto0CDs2rXLqOGIiMxJVY0GM9cnAwBGdmmGjuFeEiciqbDYsSwGFzsHDx7EuHHjrtseHByM7Oxso4QiIjJHy7afxfm8cvi6KjFtYAup45CE6k5jZRdXQaMVJU5Dt2NwsePg4ICSkpLrtp8+fRq+vr5GCUVEZG6SM4qxbMc5AMDsIa3g7mgvcSKSkp+rA+xkAtRaEbmlVVLHodswuNi5//778dZbb6GmpgZA7XwDFy9exGuvvYYHH3zQ6AGJiKRWVaPByz8mQa0VcU98IAbHB0gdiSQmlwkIcK/tysFOyubP4GJn4cKFuHz5Mvz8/FBZWYkePXogOjoarq6umDt3rikyEhFJauGm0zibWwYfFyXefiCOc+oQAPbbsSQGj8Zyc3PD7t27sW3bNhw+fBharRbt27dH3759TZGPiEhS+87n4+s9qQCA9x6M55IQpMNix3IYXOykpaUhPDwcvXv3Ru/evU2RiYjILJRW1WDqT/9CFIERnULR5w5/qSORGeFcO5bD4NNYkZGRuOuuu/D555/rJhQkIrJG7/zfSWQUViLE0xGv38sVzak+zqJsOQwudg4dOoRu3brhnXfeQVBQEO6//3789NNPUKlUpshHRCSJrSdzsPZQOgQBWPhwG7goDW4IJyvH01iWw+Bip3379nj//fdx8eJF/Pnnn/Dz88O4cePg5+eHp59+2hQZiYiaVEF5NV79uXbywGfuikDXSG+JE5E5uvY0lihyrh1z1uB5zgVBQK9evfDll19iy5YtiIyMxLfffmvMbERETU4URbz+SzLyylSI8XPBlP7NpY5EZirIvbbYKa/WoKRSLXEaupUGFzvp6elYsGAB2rZti06dOsHZ2RlLly41ZjYioib327+Z+CM5G3YyAR8+0hYO9nKpI5GZclTI4X1ldF5GUYXEaehWDD4J/cUXX2DVqlXYs2cPmjdvjpEjR+KXX35BeHi4CeIRETWdrOJKvPHLMQDAi31iEB/iLnEiMndBHo7IL69GZlEVWgXx82KuDC523n77bYwYMQIfffQR2rZta4JIRERNTxRFTPvfUZRUqdEmxB3P94ySOhJZgGAPRyRfKsalQrbsmDODi52LFy9y9lAisjrf77+Iv1PyoLST4YNH2sJO3uCz/GRD6jopZxZzfSxzplexc/ToUcTFxUEmkyE5OfmWx7Zu3doowYiImkpqXjnmbTgJAHhtUAtE+7lInIgsRd1cO5xY0LzpVey0bdsW2dnZ8PPzQ9u2bSEIQr1hdnXXBUGARqMxWVgiImPTaEVM+TEJlTUadIv0xqhu4VJHIgtSN9dOBufaMWt6FTupqanw9fXV/UxEZC0+33UOhy8WwVVph4WPtIFMxtP0pL9gzqJsEfQqdsLCwnQ/X7hwAQkJCbCzq39TtVqNvXv31juWiMicncgswaLNZwAAs+5rpfviItJXXZ+dy6UqVNVoOFWBmTK4B16vXr1uuCZWcXExevXqZZRQRESmplJrMPnHJNRoRPRr6Y8H2wdLHYkskKeTPRyvFDhZ7KRstgwudur65vxXfn4+nJ2djRKKiMjUPtx8BqeyS+HtrMD8YfEcZUoNIggCgjwcAPBUljnTe+j5sGHDANS+saNHj4ZSqdTt02g0OHr0KBISEoyfkIjIyHaeuYzPd54HAMwbFg8fF+VtbkF0c8GeTjh3uZwjssyY3sWOu3vtzJCiKMLV1RWOjlfPbSsUCnTt2hVjx441fkIiIiPKKanC5LVJAIAnuoZhQKsAaQORxQu+0rLD1c/Nl97FzvLlywEA4eHhmDp1Kk9ZEZHF0WhFvLTmCPLLq3FHoBtm3nOH1JHICtR1bGexY74MnkF51qxZpshBRGRyS7alYN/5Ajgp5PjksXYcOUNGUTcii6exzJfBxQ4A/O9//8OPP/6Iixcvorq6ut6+w4cPGyUYEZEx7T2Xh4+2pgAA5g2NR6QvZ0km4whyr1sygsWOuTJ4NNbHH3+Mp556Cn5+fjhy5Ag6d+4Mb29vnD9/HoMGDTJFRiKiRskrU2HSmiSIIvBIxxA80I7DzMl46lp2soqqoNWKtzmapGBwsbNs2TJ88cUXWLp0KRQKBaZNm4bNmzfjxRdfRHFxsSkyEhE1mFYrYvKP/yK3VIUYPxfMvq+V1JHIygS4OUAmANUaLfLKVFLHoRswuNi5ePGiboi5o6MjSktLAQBPPPEEfvjhB+OmIyJqpM93nceuM5fhYC/DJyPbw0nRoLP3RDdlJ5chwK12RBbXyDJPBhc7AQEByM/PB1C7jMS+ffsA1K6Zde3ioEREUku8UICFf50GAMwe0gqx/q4SJyJrxU7K5s3gYqd37974/fffAQBjxozByy+/jH79+mH48OEYOnSo0QMSETVEUUU1Jq4+Ao1WxH1tgjC8U6jUkciKBXFBULNmcHvuF198Aa1WCwB47rnn4OXlhd27d2PIkCF47rnnjB6QiMhQoihi6k9HkVlchXBvJ8wdGsflIMikONeOeTO42JHJZJDJrjYIPfLII3jkkUeMGoqIqDGW70nDlpM5UMhlWPpYe7g62Esdiaxc3WkstuyYJ72KnaNHj+p9h61bt25wGCKixjqaUYT5f54EAMy85w7EBbtLnIhsQd1prAz22TFLehU7bdu2hSAIt+2ALAgCNBqNUYIRERmqpKoGE1YfQY1GxIBW/niyW5jUkchGhPA0llnTq9hJTU01dQ4iokYRRRHT1yXjYkEFgj0cseDBNuynQ02mrmWntEqNkqoauPHUqVnRq9gJC+NfR0Rk3lYfuIgNR7NgJxOw5LF2cHfilw01HWelHTyc7FFUUYPMokq4BfDzZ04MHnoOAN999x3uvPNOBAUF4cKFCwCAxYsX49dffzVqOCIifZzMKsGc308AAKYNbI72zTwlTkS2KJjDz82WwcXOp59+ismTJ2Pw4MEoKirS9dHx8PDA4sWLjZ2PiOiWiitq8Nz3iahWa9GzuS+euStS6khko+pOZXFiQfNjcLGzZMkSfPnll5g5cybkcrlue8eOHZGcnGzQfe3atQtDhgxBUFAQBEHAL7/8Um+/KIqYPXs2goKC4OjoiJ49e+L48eP1jlGpVJg4cSJ8fHzg7OyM++67DxkZGYY+LSKyQFqtiElrj+BCfm0/nQ8faQuZjP10SBp1LTtcMsL8GFzspKamol27dtdtVyqVKC8vN+i+ysvL0aZNGyxduvSG+xcsWIAPP/wQS5cuxcGDBxEQEIB+/frp1uMCgEmTJmH9+vVYs2YNdu/ejbKyMtx7770cFUZkAxZvTcH205ehtJPh8yc6wMtZIXUksmFXT2NVSZyE/svgSQUjIiKQlJR0XaflP//8Ey1btjTovgYNGoRBgwbdcJ8oili8eDFmzpyJYcOGAQC+/fZb+Pv7Y/Xq1Rg3bhyKi4vx9ddf47vvvkPfvn0BAN9//z1CQ0OxZcsWDBgw4Ib3rVKpoFJdXZm2pKTEoNxEJL0tJ3Lw8dYUAMDcofGcT4ckd3V9rAqJk9B/Gdyy88orr+CFF17A2rVrIYoiDhw4gLlz52LGjBl45ZVXjBYsNTUV2dnZ6N+/v26bUqlEjx49sHfvXgBAYmIiampq6h0TFBSEuLg43TE3Mn/+fLi7u+suoaFcM4fIkpy/XIaX1yYBAJ7sFoaHOoRIG4gIbNkxZwa37Dz11FNQq9WYNm0aKioq8NhjjyE4OBgfffQRRowYYbRg2dnZAAB/f/962/39/XUjwLKzs6FQKODp6XndMXW3v5Hp06dj8uTJuuslJSUseIgsRLlKjXHfJaJUpUbHME+8fo9hLcpEplLXQTmntArVai0Udg0a8EwmYHCxAwBjx47F2LFjkZeXB61WCz8/PwDApUuXEBwcbNSA/50UTBTF204UdrtjlEollEqlUfIRUdMRRRGv/O9fpOSWwc9ViWUj2/MLhcyGj4sCSjsZVGotsour0MzbSepIdEWjfkv4+PjAz88P2dnZmDhxIqKjo42VCwEBAQBwXQtNbm6urrUnICAA1dXVKCwsvOkxRGQ9Pt91Hn8kZ8NeLuDTx9vDz81B6khEOoIgcPVzM6V3sVNUVISRI0fC19cXQUFB+Pjjj6HVavHmm28iMjIS+/btwzfffGO0YBEREQgICMDmzZt126qrq7Fz504kJCQAADp06AB7e/t6x2RlZeHYsWO6Y4jIOuxOycOCjacAAG/e2xIdwrwkTkR0vSAWO2ZJ79NYM2bMwK5duzBq1Chs3LgRL7/8MjZu3Iiqqir8+eef6NGjh8EPXlZWhrNnz+qup6amIikpCV5eXmjWrBkmTZqEefPmISYmBjExMZg3bx6cnJzw2GOPAQDc3d0xZswYTJkyBd7e3vDy8sLUqVMRHx+vG51FRJYvvaACE384DK0IPNQhBI935RI2ZJ44i7J50rvY2bBhA5YvX46+ffvi+eefR3R0NGJjYxs1a/KhQ4fQq1cv3fW6TsOjRo3CihUrMG3aNFRWVuL5559HYWEhunTpgr/++guurq662yxatAh2dnZ45JFHUFlZiT59+mDFihX1JjwkIstVVaPB+FWJKKyoQVywG955II4LfJLZ4izK5kkQRVHU50B7e3tcuHABQUFBAAAnJyccOHAAcXFxJg3YFEpKSuDu7o7i4mK4ublJHYeIrhBFEVN/OoqfD2fAy1mB3ybciRBPdvok8/W/xAxM/elf3BXtg++f6SJ1HKun7/e33n12tFot7O2vruIql8vh7OzcuJRERLfw/b4L+PlwBmQCsOTRdix0yOzxNJZ50vs0liiKGD16tG7IdlVVFZ577rnrCp5169YZNyER2aRDaQW6lcxfHdgCd0b7SJyI6PauHY2lz1Qp1DT0LnZGjRpV7/rjjz9u9DBERACQU1KF8asOQ60VcU98IJ7tzpXMyTIEuDtAEACVWov88mr4uHBON3Ogd7GzfPlyU+YgIgIAVKu1eH7VYVwuVSHW3wULHmrNv47JYijsZPBzVSKnRIVLhZUsdswEpx4lIrMhiiJmrE9G4oVCuCrt8PkTHeGsbNBE70SS4cSC5ofFDhGZjUWbz+B/ibUdkj9+tB0ifDgIgixP8JWO9OykbD5Y7BCRWVi9/yI+3lY7yeg7D8SjVws/iRMRNUyQR+0yJhmca8dssNghIsltPZmD139JBgBM7B2Nx7o0kzgRUcOF8DSW2WGxQ0SSSkovwoTVR3RLQUzuFyt1JKJGCeJcO2bH4J5/+fn58Pb2BgCkp6fjyy+/RGVlJe677z7cfffdRg9IRNYrLa8cY1YcRGWNBt1jfTF/WDxHXpHFC/Zky4650btlJzk5GeHh4fDz80OLFi2QlJSETp06YdGiRfjiiy/Qq1cv/PLLLyaMSkTWJL9MhdHLDyC/vBpxwW5YNrI97OVsbCbLVzcaq6iiBuUqtcRpCDCg2Jk2bRri4+Oxc+dO9OzZE/feey8GDx6M4uJiFBYWYty4cXj33XdNmZWIrERltQZPf3sIafkVCPF0xDejO8GFQ8zJSrg62MPVofbzzFNZ5kHvYufgwYOYO3cu7rrrLixcuBCZmZl4/vnnIZPJIJPJMHHiRJw6dcqUWYnICqg1Wkz84TD+TS+Ch5M9VjzVGX6uDlLHIjKqutadDBY7ZkHvYqegoAABAQEAABcXFzg7O8PLy0u339PTE6WlpcZPSERWQxRFvPnbcWw5mQulnQxfPdkR0X4uUsciMjouCGpeDDpB/t+Og+xISESGWLbjHFbvvwhBAD4a0Q4dw71ufyMiC6TrpMy5dsyCQSfJb7XquUqlMn46IrIaPydm4P1NpwEAs4e0wsC4AIkTEZkOW3bMi97FzpNPPlmvJedGq54/+eSTxklFRFbl75TLePXnowCAcT0iMSohXNpARCYWxIkFzYrexc6KFStMGIOIrNXxzGKM//4w1FoR97UJwqsDWkgdicjkeBrLvOjdZ0culyM3N9eUWYjIyqQXVOCp5QdRplKjW6Q33n+4NWQy9vUj61e3ZER2SRXUGq3EaUjvYkcURVPmICIrk11chZFf7UduqQotAlzx+ZMdoLSTSx2LqEn4uCihkMugFWsLHpIWpyslIqPLL1Nh5Ff7cLGgAmHeTvj26c5wc7CXOhZRk5HJBAReWf08s4jFjtQMGo21adMmuLu73/KY++67r1GBiMiyFVfW4ImvD+Dc5XIEujtg1TNd4O/GSQPJ9gS5O+JCfgUuFVUA4DQLUjKo2Bk1atQt9wuCAI1G06hARGS5ylVqjF5+ACeySuDjosCqZ7ogxNNJ6lhEkmAnZfNh0Gms7OxsaLXam15Y6BDZrqoaDZ759hCOXCyCu6M9vhvTBZG+nB2ZbFewbvg5T2NJTe9ih7MlE9HNVKu1eH7VYfxzPh/OCjm+fboz7gh0kzoWkaSCOdeO2eBoLCJqFI1WxMtrk7DtVC4c7GX4ZnQntA31kDoWkeTqTmNxFmXp6V3sjBo1Co6Ojrc8pqamptGBiMhyaLUiXv35KDYkZ8FeLuDzJzqiS6S31LGIzIJuFuXCSjYYSEzvYkej0dzyzTp06BDatWtnlFBEZP5EUcSc34/jf4kZkMsELHm0HXrE+kodi8hsBLrXjkKsrNGgsIKNAVLSu9g5duwYWrZsiU2bNtXbXlNTgxkzZiAhIQF33XWX0QMSkXlasOk0vv3nAgQBWPhwawyMC5Q6EpFZcbCXw9e1dvFsnsqSlt7FzoEDB/D0009jyJAhGDduHMrKynDo0CG0bdsWq1evxoYNG/DZZ5+ZMisRmYlPtp/FpzvOAQDeeSAOQ9uFSJyIyDzVncrK4PBzSeld7NjZ2eGtt97CP//8gz179iA2NhYJCQm48847kZycjH79+pkyJxGZieV7UvH+ptMAgJmD78DILmESJyIyX3VrZLFlR1oGLxehVCphb2+P4uJiKBQK3HnnnXB1dTVFNiIyMz8eTMec308AAF7qE4Ox3SMlTkRk3oKuLBnB4efSMmjo+fz589GxY0e0bdsWmZmZWLBgASZMmID777+fK6ITWbkfD6bj1XVHAQDP3BWBSX1jJE5EZP6CPTiLsjnQu9jp1q0blixZgp9++gnLly+Hu7s7nn/+efz7778oKipCy5YtsXbtWlNmJSKJLN+Timk/H4UoAo93bYaZ99zBiUaJ9BB8ZbmUzGIWO1LSu9gJDw/HsWPHMGTIkHrbIyMjsWPHDsycORNjxowxekAiktYn28/qTl09c1cE3r4/joUOkZ50p7HYsiMpvYudNWvWwMur/qqte/bsgUqlgiAIePnll3HkyBGjByQiaYiiiPc3ndJ1Rn6xTwxbdIgMFOpV27KTX16NkirOtSMVgzsoX2vQoEG4dOmS7npMDM/hE1mD2gkDT+CT7bXDy6cPaoHJ/WJZ6BAZyM3BHv5utXPtnM0tkziN7WpUscPpr4msj0Yr4rWfk7FibxoA4O37W2FcjyhpQxFZsBi/2hHLLHak06hih4isS41Gi0lrk7D2UDpkArDw4TZ4olu41LGILFq0nwsAFjtSsmvMjT///HP4+/sbKwsRSaiqRoMJq49gy8kc2MkEfDSiHe5pzSUgiBqrrthJySmVOIntalSx89hjjxkrBxFJqKJajXHfJeLvlDwo7GT47PH26N2Cf8gQGUNMXbHDlh3JGHQaKysrC99//z3++OMPVFdX19tXXl6Ot956y6jhiMj0SqtqMOqbA/g7JQ9OCjlWjO7EQofIiGL8a/vsZBRWoqJaLXEa26R3sXPw4EG0bNkSL7zwAh566CHExcXh+PHjuv1lZWWYM2eOSUISkWkUlldj5Ff7cTCtEK4OdvhuTGckRPtIHYvIqng5K+DtrAAAnL9cLnEa26R3sTNjxgwMGzYMhYWFyMnJQb9+/dCjRw/OrUNkoXJLqzDii304mlEMTyd7/DC2KzqEed3+hkRksCjdqSz225GC3n12EhMT8cknn0Amk8HV1RWffPIJwsLC0KdPH2zatAnNmjUzZU4iMqJzl8swZsVBpOVXwM9ViVXPdNE1tROR8cX4ueBAagFScthvRwoGdVCuqqqqd33atGmQyWTo378/vvnmG6MGIyLT2H4qFy/+cASlKjWCPRyx6pkuCPdxljoWkVVjJ2Vp6V3sxMXFYe/evWjdunW97VOnToUoinj00UeNHo6IjEcURXy28zwWbDoFUQQ6hXti2cgO8HVVSh2NyOrVtZyeY7EjCb377Dz55JPYs2fPDfe98soreOutt3gqi8hMVVZr8OKaJLy3sbbQebRzM6x6pisLHaImUjfXTlp+OVRqjcRpbI8gcs0HlJSUwN3dHcXFxXBzc5M6DpFRXSqqxLMrD+F4ZgnsZAJm39cKj3cNkzoWkU0RRRGt5/yF0io1Nk66Gy0C+F1jDPp+f3O5CCIrdiC1APct2Y3jmSXwclZg1TNdWOgQSUAQhKv9dthJucnp3WenV69et13xWBAEbN26tdGhiKjxvt93AbN/Ow61VkTLQDd88WQHhHg6SR2LyGbF+Lni8MUirpElAb2LnbZt2950X0lJCX744QeoVCpjZCKiRqhWazH79+NYvf8iAOCe1oF4/6HWcFI0anUYImokLggqHb1/+y1atOi6bWq1Gp988gnmzp2L4OBgvP3220YNR0SGyStTYfz3iTiYVghBAKb2b47ne0bdtlWWiEwv2p8TC0qlwX/qrVq1Cm+++SYqKysxe/ZsPPvss7Cz41+ORFI5dqkYz648hMziKrgo7fDRiLbocwfXuCIyF3V9dlLzyqHWaGEnZ7fZpmJwdbJx40a89tprSE1NxdSpUzF58mQ4O3NCMiIp/Zp0Ca/+fBRVNVpE+Djjyyc7INqPMyITmZMgd0c4KeSoqNbgQkEFonxdpI5kM/QuKw8cOIBevXph6NCh6NWrF86dO4c33njDpIXO7NmzIQhCvUtAQIBuvyiKmD17NoKCguDo6IiePXvWW5yUyNoVlFdjwurDeGlNEqpqtOgR64tfXriThQ6RGZLJBF2BwxFZTUvvlp2uXbvC0dER48ePR3h4OFavXn3D41588UWjhQOAVq1aYcuWLbrrcrlc9/OCBQvw4YcfYsWKFYiNjcU777yDfv364fTp03B15S97sm5/Hc/GjPXHkFemglwm4IVe0XipTwzkMvbPITJXMX4uSL5UjLO5pQACbns8GYfexU6zZs0gCALWr19/02MEQTB6sWNnZ1evNaeOKIpYvHgxZs6ciWHDhgEAvv32W/j7+2P16tUYN27cTe9TpVLVGzlWUlJi1MxEplRcUYM5vx/HuiOXANT+8vzgkTZoHeIhbTAiuq2rnZTZstOU9C520tLSTBjj5lJSUhAUFASlUokuXbpg3rx5iIyMRGpqKrKzs9G/f3/dsUqlEj169MDevXtvWezMnz8fc+bMaYr4REa1/XQuXvv5KHJKVJAJwLPdozCpbwwc7OW3vzERSS7myilmDj9vWmbdFbxLly5YuXIlNm3ahC+//BLZ2dlISEhAfn4+srOzAQD+/vVHm/j7++v23cz06dNRXFysu6Snp5vsORAZQ0lVDV7931E8tfwgckpUiPRxxk/PJeC1QS1Y6BBZkGvn2tFobX61piZj0GmsI0eOwNvbGwCwdOlSPPnkkyZdS2rQoEG6n+Pj49GtWzdERUXh22+/RdeuXQHguvlDRFG87ZwiSqUSSiUXQCTLsDslD9P+9y8yi6sgCMDTd0bglQHNWeQQWaBQT0co7GRQqbW4VFiJZt6c1bwp6N2yk5GRAY3m6kqtM2bMQF5enklC3YyzszPi4+ORkpKi68fz31ac3Nzc61p7iCxRuUqNmeuT8fjX+5FZXIVmXk5Y+2w3vHFvSxY6RBbKTi5DpE/tKGZOLth0GnwaS4rF0lUqFU6ePInAwEBEREQgICAAmzdv1u2vrq7Gzp07kZCQ0OTZiIxp3/l8DPxoF1ZdWfLhyW5h2DjpbnSO8JI4GRE1Vow/++00NbOe8njq1KkYMmQImjVrhtzcXLzzzjsoKSnBqFGjIAgCJk2ahHnz5iEmJgYxMTGYN28enJyc8Nhjj0kdnahBSqtqsHDTaXz7zwUAQLCHI95/qDUSon0kTkZExhLtyxFZTc2gYuerr76Ci0vtm6RWq7FixQr4+NT/JWzMoecZGRl49NFHkZeXB19fX3Tt2hX79u1DWFgYAGDatGmorKzE888/j8LCQnTp0gV//fUX59ghi7TpeDZm/Xoc2SVVAIBHOzfDjMEt4OpgL3EyIjKmGA4/b3KCqOf5qPDw8Nt2/BUEAefPnzdKsKZUUlICd3d3FBcXm7TDNdGNZBVXYtavx/HXiRwAQJi3E955IA53x/hKnIyITCElpxT9Fu2Ci9IOybP7c6HeRtD3+9vs59khslYarYjv/knDwr/OoEylhp1MwLPdI/FiH86bQ2TNwrydYScTUKZSI7ukCoHujlJHsnpm3WeHyFqdyCzB9PXJ+De9CADQvpkH5g2LR4sAtiwSWTuFnQxh3k44d7kcKTllLHaagN7FTmVlJbZu3Yp7770XQO3EfNcuuSCXy/H222/DwcHB+CmJrERltQaLt5zBV7tTodGKcFXaYdqgFhjZuRlkXNOKyGbE+LnWFju5Zegey1PWpqZ3sbNy5Ur83//9n67YWbp0KVq1agVHx9qK9NSpUwgKCsLLL79smqREFm7nmct4/ZdkpBdUAgAGxwdg1pBW8HfjHwhEtibG3wUbj+PKgqBkanoXO6tWrbqukFm9ejUiIyMBAN9//z0++eQTFjtE/3G5VIW3/+8Efvs3EwAQ5O6At+6PQ9+WnPySyFZdu2wEmZ7exc6ZM2cQGxuru+7g4ACZ7OqchJ07d8YLL7xg3HREFkyjFbH2YDre/fMkSqrUkAnA6IQITOkfC2clu8sR2bK6YudMTpleyxxR4+j9G7e4uBh2dlcPv3z5cr39Wq22Xh8eIlt25GIhZv12HEczigEArYLc8O6w1ogPcZc4GRGZgyhfFwgCUFxZg7yyavi6cr1GU9K72AkJCcGxY8fQvHnzG+4/evQoQkJCjBaMyBLll6mwYONprD2UDgBwVdrh5X6xeLJbGOzkDV6dhYisjIO9HM28nHAhvwIpuaUsdkxM79++gwcPxptvvomqqqrr9lVWVmLOnDm45557jBqOyFJotCJW/pOGXgt36AqdYe2DsXVqDzx9VwQLHSK6TsyVU1nn2G/H5PRu2ZkxYwZ+/PFHNG/eHBMmTEBsbCwEQcCpU6ewdOlSqNVqzJgxw5RZicxS4oUCvPHLcZzIKgEA3BHohrfvb4WO4Vy0k4huLsrPBVtO5nLZiCagd7Hj7++PvXv3Yvz48Xjttdd0q54LgoB+/fph2bJl8Pfn6BKyHZdLVXj3z1P4+XAGAMDNwQ5TBzTHY52bsSWHiG4rxq92HceUHBY7pmbQkJCIiAhs3LgRBQUFOHv2LAAgOjoaXl78C5Zsh1qjxcp/LmDR5jMoVakBAI90DMG0gS3g48Lz7kSkn7rTWGzZMb0GjX/18vJC586djZ2FyOztP5+PWb8dx6ns2onA4oLd8Nb9cWjfzFPiZERkaaKuFDt5ZSoUVVTDw0khcSLrxck+iPRwNrcUH24+gz+SswEAHk72eGVAc4zo1AxyLvNARA3gorRDkLsDMourcDa3jP38TIjFDtEtZBRWYPGWFKw7nAGtCAgCMKJTM0wb0ByezvwrjIgaJ9rfFZnFVUhhsWNSLHaIbuByqQqfbD+LVfsvoEZT2xm/X0t/TOkfy5XJichoYvxcsOvMZXZSNjEWO0TXKK6owRd/n8M3u9NQWaMBACREeeOVAc3Rjv1yiMjI6jopn73MYseUWOwQAaioVmP5njR8vvMcSqpqR1i1CfXAtAHNcWe0j8TpiMha6RYEzeHq56bEYodsmkqtwQ/7L2Lp9nPIK6td2y3W3wVT+zdHv5b+XJyPiEyqrtjJLK5CaVUNXB3sJU5knVjskE1Sa7RYf+QSFm9JwaWiSgBAMy8nTO4XiyFtgjjCioiahIeTAr6uSlwuVeHc5XK0DfWQOpJVYrFDNqW0qgZrD6Zj+Z40XZHj56rEi31i8EjHUCjsOPMxETWtGD8XXC5V4WxuGYsdE2GxQzbhUlEllu9OxZqD6Si7Muuxt7MC43pE4slu4XCwl0uckIhsVbSfC/aey0dKLvvtmAqLHbJqRzOK8OXfqfgjOQsabe0Q8mg/FzxzVwQeaBfMIoeIJKcbkcXh5ybDYoesjlYrYuupXHz593kcSC3QbU+I8sbYuyPRI9YXMvbJISIzEV23ICjXyDIZFjtkNSqrNfjf4Qx8szsVqXnlAAA7mYAhbYLwzN0RaBXkLnFCIqLrxfjXtuykF1agqkbDFmcTYLFDFi+jsAJrD6bj+30XUFhRAwBwc7DDY13CMDohHAHuDhInJCK6OW9nBTyc7FFUUYNzl8v4h5kJsNghi1RVo8HGY9n4KTEde8/lQ6ztjoNQL0c8fWcEHukYCmclP95EZP4EQUCMnwsOphXibC6LHVPgtwFZDFEUkZRehJ8SM/D7v5kovTLTMVDbH+fxrmEY0CqAc+QQkcWJ9nPFwbRCrpFlIix2yOzlllbhlyOX8NOhjHod+II9HPFQhxA81CEEoV5OEiYkImoc3YgsdlI2CRY7ZJaq1VpsO5WL/yWmY/vpy7ph4w72MgyKC8TDHULQNdKbo6qIyCrULRvBuXZMg8UOmQ2NVsThi4X4MzkbvyZdQn55tW5fu2YeeLhDKO5tEwg3rh1DRFambkRWWn4FqtVazuZuZCx2SFLVai32nsvDpuM52HwiR7cYJwD4uioxrH0wHu4QopuHgojIGgW4OcBFaYcylRpp+eWI9efvPGNisUNNrlylxo7Tl7HpeDa2n8pFqepqR2NXBzv0aeGHIW2C0CPWF3Zy/nVDRNZPEARE+7kgKb0IZ3PLWOwYGYsdahIF5dXYciIHm45n4++zeahWa3X7fF2V6N/SHwNaBaBrpDebb4nIJtUVOyk5ZUC81GmsC4sdMglRFHEmpwy7z+Zh84lsHEgtwJU+xgCAcG8nDGgVgP6tAtAu1IMdjYnI5sWwk7LJsNgho0kvqMCes3nYcy4f/5zLQ15Zdb39LQPdMDAuAANaBSDW3wWCwAKHiKhOXSdlDj83PhY71GCXS1X453w+9p7Nw55zeUgvqKy339Fejk4RXuge44MBrQI4Fw4R0S3EXBmIcT6vHGqNln0WjYjFDumtpKoGB84XYM+5POw9m4/TOfWbWu1kAtqGeiAh2gd3RnmjbTMPKO24oB0RkT6CPRzhYC9DVY0W6YWViPBxljqS1WCxQzek1mhxOqcUSelFSLpYVDtC4HKZbg2qOncEuuHOKG/cGe2DThFecOF6VEREDSKTCYjydcHxzBKk5JSy2DEifjMRACC7uApJ6YU4crEIR9KLkJxRjMoazXXHhXs7ISHaBwlR3ugW6Q1vF6UEaYmIrFOM35ViJ7cM/VtJncZ6sNixQXllKpzJLkXypWIkpRfhyMUiZJdUXXecq9IOrUPd0S7UE21DPdC2mQd8WNwQEZlMzJX5dc6xk7JRsdixYiVVNUjJKcXp7DKcySnF6exSnMkprbcMQx2ZADQPcEPbUA+0C/VAu2YeiPJ14ZBwIqImFOVbN/ycxY4xsdixAhXVapy/XF5b0OSU4kx2bWGTWXx9aw0ACALQzMsJLQJc0a5ZbatNfLA7nNnfhohIUtcOP9dqRf7BaST8drMQVTUaXMivQGpeOdLyy5GWV677OadEddPbBbg5IDbAFc39XRDr74rmAa6I9nOBk4JvPRGRuQnzcoK9XEBljQaXiio5ZYeR8BvPDNRotCiurNFd8suqa4uZK0VNWl45skqqrhsJdS1PJ3tdMaP7188V7k5cIZyIyFLYyWWI9HHB6ZxSnL1cxmLHSFjsmND+8/k4mVWC4kp1vWKmpO7fqtp/K6qvH/V0I64OdojwcUa4tzPCfZwR4eOEcG9nRPg4w8NJYeJnQ0RETSHa70qxk1OGXs39pI5jFVjsmNC6w5ew9lC63se7OtjB3dEeXs4KNPNy+k9h4wxPJ3susUBEZOWiuUaW0bHYMaG2zTxQplLDzdEe7re4uDnawdXBHnJ2RCMisnl1nZT/OZ+PE5klaBnkJnEiyyeI4q16gtiGkpISuLu7o7i4GG5u/FAREZF0Mgor0PuDnahWawEA98QHYlLfGN0cPHSVvt/fXGWMiIjIjIR4OuHPl+7GkDZBEARgQ3IW+i/ehUlrjiA1r1zqeBaJLTtgyw4REZmnU9klWLT5DDYdzwEAyGUCHmwfjIm9YzhSC/p/f7PYAYsdIiIyb8cuFePDzWew7VQuAMBeLuCRjqGY0Dsage6OEqeTDosdA7DYISIiS3D4YiE+/OsMdp/NAwAo7GQY2aUZxveMgp+rg8Tpmh6LHQOw2CEiIkuy/3w+PvjrDA6kFQAAHOxlGNUtHIPiAxHq6QgvZ4VNTFVic8XOsmXL8P777yMrKwutWrXC4sWLcffdd+t1WxY7RERkaURRxO6zefjgrzNISi+qt89JIUeIpyNCPZ1q//VyQsg1P7s7Wsfs+jZV7KxduxZPPPEEli1bhjvvvBOff/45vvrqK5w4cQLNmjW77e1Z7BARkaUSRRHbT+fim91pSMktveV6iXXcHOx0xY+3ixIeTvbwdLKHh6MC7k728HC0h4eTAp5O9nB3sofSTt4Ez8RwNlXsdOnSBe3bt8enn36q23bHHXfggQcewPz58297exY7RERkLapqNMgsqkR6YSUyCiuQXlCJ9MIKZBRWIqOgAvnl1Qbfp6O9HB5OtQWQu6MdXJR2cFbawUlhB2eFHM5KOzgr5XBS1O5zUshr/1Ve3e/lrICDvXGLJn2/vy1+BuXq6mokJibitddeq7e9f//+2Lt37w1vo1KpoFJdrXxLSkpMmpGIiKipONjLEenrgkhflxvur6hWI6OwEukFFbhUVInC8hoUVVajuKIGhRXVKKqsQXFFDYoqa1BUUQ2tCFTWaFBZrEFWcVWDc825rxVGJYQ3+PaNYfHFTl5eHjQaDfz9/ett9/f3R3Z29g1vM3/+fMyZM6cp4hEREZkVJ4UdYv1dEavHjMxarYiyajWKrhRERVeKoAqVGuXVGpSr1CivVqNCdfXncpVGt61MpUZFde2xzkrpSg6LL3bq/LfXuSiKN+2JPn36dEyePFl3vaSkBKGhoSbNR0REZGlkMgFuDvZwc7BHMzRuEkMpe81YfLHj4+MDuVx+XStObm7uda09dZRKJZRKZVPEIyIiIlzfKNGULH5tLIVCgQ4dOmDz5s31tm/evBkJCQkSpSIiIiJzYfEtOwAwefJkPPHEE+jYsSO6deuGL774AhcvXsRzzz0ndTQiIiKSmFUUO8OHD0d+fj7eeustZGVlIS4uDn/88QfCwsKkjkZEREQSs4p5dhqL8+wQERFZHn2/vy2+zw4RERHRrbDYISIiIqvGYoeIiIisGosdIiIismosdoiIiMiqsdghIiIiq8Zih4iIiKwaix0iIiKyaix2iIiIyKpZxXIRjVU3iXRJSYnESYiIiEhfdd/bt1sMgsUOgNLSUgBAaGioxEmIiIjIUKWlpXB3d7/pfq6NBUCr1SIzMxOurq4QBMFo91tSUoLQ0FCkp6fb7Jpbtv4a2PrzB/ga2PrzB/ga8Pmb7vmLoojS0lIEBQVBJrt5zxy27ACQyWQICQkx2f27ubnZ5Af8Wrb+Gtj68wf4Gtj68wf4GvD5m+b536pFpw47KBMREZFVY7FDREREVo3FjgkplUrMmjULSqVS6iiSsfXXwNafP8DXwNafP8DXgM9f+ufPDspERERk1diyQ0RERFaNxQ4RERFZNRY7REREZNVY7BAREZFVY7FjQsuWLUNERAQcHBzQoUMH/P3331JHahKzZ8+GIAj1LgEBAVLHMqldu3ZhyJAhCAoKgiAI+OWXX+rtF0URs2fPRlBQEBwdHdGzZ08cP35cmrAmcLvnP3r06Os+E127dpUmrAnMnz8fnTp1gqurK/z8/PDAAw/g9OnT9Y6x9s+APq+BNX8OPv30U7Ru3Vo3cV63bt3w559/6vZb+/sP3P41kPL9Z7FjImvXrsWkSZMwc+ZMHDlyBHfffTcGDRqEixcvSh2tSbRq1QpZWVm6S3JystSRTKq8vBxt2rTB0qVLb7h/wYIF+PDDD7F06VIcPHgQAQEB6Nevn25dNkt3u+cPAAMHDqz3mfjjjz+aMKFp7dy5Ey+88AL27duHzZs3Q61Wo3///igvL9cdY+2fAX1eA8B6PwchISF49913cejQIRw6dAi9e/fG/fffrytorP39B27/GgASvv8imUTnzp3F5557rt62Fi1aiK+99ppEiZrOrFmzxDZt2kgdQzIAxPXr1+uua7VaMSAgQHz33Xd126qqqkR3d3fxs88+kyChaf33+YuiKI4aNUq8//77JckjhdzcXBGAuHPnTlEUbe8zIIrXvwaiaHufA09PT/Grr76yyfe/Tt1rIIrSvv9s2TGB6upqJCYmon///vW29+/fH3v37pUoVdNKSUlBUFAQIiIiMGLECJw/f17qSJJJTU1FdnZ2vc+DUqlEjx49bObzAAA7duyAn58fYmNjMXbsWOTm5kodyWSKi4sBAF5eXgBs8zPw39egji18DjQaDdasWYPy8nJ069bNJt///74GdaR6/7kQqAnk5eVBo9HA39+/3nZ/f39kZ2dLlKrpdOnSBStXrkRsbCxycnLwzjvvICEhAcePH4e3t7fU8Zpc3Xt+o8/DhQsXpIjU5AYNGoSHH34YYWFhSE1NxRtvvIHevXsjMTHR6maVFUURkydPxl133YW4uDgAtvcZuNFrAFj/5yA5ORndunVDVVUVXFxcsH79erRs2VJX0NjC+3+z1wCQ9v1nsWNCgiDUuy6K4nXbrNGgQYN0P8fHx6Nbt26IiorCt99+i8mTJ0uYTFq2+nkAgOHDh+t+jouLQ8eOHREWFoYNGzZg2LBhEiYzvgkTJuDo0aPYvXv3dfts5TNws9fA2j8HzZs3R1JSEoqKivDzzz9j1KhR2Llzp26/Lbz/N3sNWrZsKen7z9NYJuDj4wO5XH5dK05ubu51lb0tcHZ2Rnx8PFJSUqSOIom6kWj8PFwVGBiIsLAwq/tMTJw4Eb/99hu2b9+OkJAQ3XZb+gzc7DW4EWv7HCgUCkRHR6Njx46YP38+2rRpg48++sim3v+bvQY30pTvP4sdE1AoFOjQoQM2b95cb/vmzZuRkJAgUSrpqFQqnDx5EoGBgVJHkURERAQCAgLqfR6qq6uxc+dOm/w8AEB+fj7S09Ot5jMhiiImTJiAdevWYdu2bYiIiKi33xY+A7d7DW7E2j4H/yWKIlQqlU28/zdT9xrcSJO+/5J0i7YBa9asEe3t7cWvv/5aPHHihDhp0iTR2dlZTEtLkzqayU2ZMkXcsWOHeP78eXHfvn3ivffeK7q6ulr1cy8tLRWPHDkiHjlyRAQgfvjhh+KRI0fECxcuiKIoiu+++67o7u4urlu3TkxOThYfffRRMTAwUCwpKZE4uXHc6vmXlpaKU6ZMEffu3SumpqaK27dvF7t16yYGBwdbzfMfP3686O7uLu7YsUPMysrSXSoqKnTHWPtn4HavgbV/DqZPny7u2rVLTE1NFY8ePSrOmDFDlMlk4l9//SWKovW//6J469dA6vefxY4JffLJJ2JYWJioUCjE9u3b1xuCac2GDx8uBgYGivb29mJQUJA4bNgw8fjx41LHMqnt27eLAK67jBo1ShTF2qHHs2bNEgMCAkSlUil2795dTE5Olja0Ed3q+VdUVIj9+/cXfX19RXt7e7FZs2biqFGjxIsXL0od22hu9NwBiMuXL9cdY+2fgdu9Btb+OXj66ad1v+99fX3FPn366AodUbT+918Ub/0aSP3+C6IoiqZvPyIiIiKSBvvsEBERkVVjsUNERERWjcUOERERWTUWO0RERGTVWOwQERGRVWOxQ0RERFaNxQ4RERFZNRY7REREZNVY7BAREZFVY7FDRE1i9OjREAQBgiDAzs4OzZo1w/jx41FYWFjvuCNHjuDee++Fn58fHBwcEB4ejuHDhyMvL093zM8//4wuXbrA3d0drq6uaNWqFaZMmXLLx9++fTt69eoFLy8vODk5ISYmBqNGjYJarQYArFixAh4eHkZ/3kQkPRY7RNRkBg4ciKysLKSlpeGrr77C77//jueff163Pzc3F3379oWPjw82bdqEkydP4ptvvkFgYCAqKioAAFu2bMGIESPw0EMP4cCBA0hMTMTcuXNRXV1908c9fvw4Bg0ahE6dOmHXrl1ITk7GkiVLYG9vD61Wa/LnTUQSa5IVuIjI5o0aNUq8//77622bPHmy6OXlpbu+fv160c7OTqypqbnp/bz00ktiz549DXrsRYsWieHh4Tfdf6OFTGfNmiWKoiiqVCrxlVdeEYOCgkQnJyexc+fO4vbt23W3Xb58ueju7i6uX79ejImJEZVKpdi3b1+rWeCSyBqwZYeIJHH+/Hls3LgR9vb2um0BAQFQq9VYv349xJusURwQEIDjx4/j2LFjej9WQEAAsrKysGvXrhvuT0hIwOLFi+Hm5oasrCxkZWVh6tSpAICnnnoKe/bswZo1a3D06FE8/PDDGDhwIFJSUnS3r6iowNy5c/Htt99iz549KCkpwYgRI/TOR0QmJnW1RUS2YdSoUaJcLhednZ1FBwcHXQvKhx9+WO+4GTNmiHZ2dqKXl5c4cOBAccGCBWJ2drZuf1lZmTh48GARgBgWFiYOHz5c/Prrr8WqqqqbPrZarRZHjx4tAhADAgLEBx54QFyyZIlYXFysO6auheZaZ8+eFQVBEC9dulRve58+fcTp06frbgdA3Ldvn27/yZMnRQDi/v37DX6diMj42LJDRE2mV69eSEpKwv79+zFx4kQMGDAAEydOrHfM3LlzkZ2djc8++wwtW7bEZ599hhYtWiA5ORkA4OzsjA0bNuDs2bN4/fXX4eLigilTpqBz5866fj3/JZfLsXz5cmRkZGDBggUICgrC3Llz0apVK2RlZd007+HDhyGKImJjY+Hi4qK77Ny5E+fOndMdZ2dnh44dO+qut2jRAh4eHjh58mRjXi4iMhIWO0TUZJydnREdHY3WrVvj448/hkqlwpw5c647ztvbGw8//DA++OADnDx5EkFBQVi4cGG9Y6KiovDMM8/gq6++wuHDh3HixAmsXbv2lo8fHByMJ554Ap988glOnDiBqqoqfPbZZzc9XqvVQi6XIzExEUlJSbrLyZMn8dFHH9U7VhCE625/o21E1PRY7BCRZGbNmoWFCxciMzPzpscoFApERUWhvLz8pseEh4fDycnplsf8l6enJwIDA3W3USgU0Gg09Y5p164dNBoNcnNzER0dXe8SEBCgO06tVuPQoUO666dPn0ZRURFatGihdx4iMh07qQMQke3q2bMnWrVqhXnz5mHp0qX4v//7P6xZswYjRoxAbGwsRFHE77//jj/++APLly8HAMyePRsVFRUYPHgwwsLCUFRUhI8//hg1NTXo16/fDR/n888/R1JSEoYOHYqoqChUVVVh5cqVOH78OJYsWQKgtmAqKyvD1q1b0aZNGzg5OSE2NhYjR47Ek08+iQ8++ADt2rVDXl4etm3bhvj4eAwePBgAYG9vj4kTJ+Ljjz+Gvb09JkyYgK5du6Jz585N80IS0a1J3WmIiGzDjYaei6Iorlq1SlQoFOLFixfFc+fOiWPHjhVjY2NFR0dH0cPDQ+zUqZO4fPly3fHbtm0TH3zwQTE0NFRUKBSiv7+/OHDgQPHvv/++6WMfPnxYfPzxx8WIiAhRqVSK3t7eYvfu3cXffvut3nHPPfec6O3tXW/oeXV1tfjmm2+K4eHhor29vRgQECAOHTpUPHr0qCiKVzs2//zzz2JkZKSoUCjE3r17i2lpaY1+zYjIOARRvMn4TiIiuq0VK1Zg0qRJKCoqkjoKEd0E++wQERGRVWOxQ0RERFaNp7GIiIjIqrFlh4iIiKwaix0iIiKyaix2iIiIyKqx2CEiIiKrxmKHiIiIrBqLHSIiIrJqLHaIiIjIqrHYISIiIqv2/9aptF49B0fKAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "scan_lib = ml.ConformerLibrary(\"kras_inhibitors_scan.clib\")\n", "\n", "with scan_lib.reading():\n", " ens = scan_lib[\"25\"]\n", "\n", "# nrg = [\n", "# x[\"ORCA/SCF_Energy\"] + x[\"ORCA/VdW_Correction\"]\n", "# for x in ens.attrib[\"conformer_attrib\"]\n", "# ]\n", "\n", "nrg = ens.attrib['XTB/Conformer_Energies']\n", "\n", "nrg_rel = 2625.5 * (nrg - np.min(nrg))\n", "\n", "idx = ens.get_atom_indices(\"2\", \"1\", \"3\", \"4\")\n", "print(*idx, np.degrees(ens[28].dihedral(*idx)))\n", "\n", "minima, _ = find_peaks(-nrg_rel, distance=10)\n", "maxima, _ = find_peaks(nrg_rel, distance=10, height=10)\n", "\n", "plt.plot(nrg_rel)\n", "plt.plot(minima, nrg_rel[minima], \"ob\")\n", "plt.plot(maxima, nrg_rel[maxima], \"^r\")\n", "plt.xlabel(\"RSS Step\")\n", "plt.ylabel(\"GFN2-XTB Relative energy (kJ/mol)\")\n", "print(f\"Minima: {minima}, maxima: {maxima}\", np.max(nrg_rel[maxima]))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Approx. barrier 18: 125.8 kJ/mol 30.1 kcal/mol\n", "Approx. barrier 22: 127.5 kJ/mol 30.5 kcal/mol\n", "Approx. barrier 23: 162.3 kJ/mol 38.8 kcal/mol\n", "Approx. barrier 24: 203.4 kJ/mol 48.6 kcal/mol\n", "Approx. barrier 25: 246.9 kJ/mol 59.0 kcal/mol\n", "Approx. barrier 26: 81.1 kJ/mol 19.4 kcal/mol\n", "Approx. barrier 27: 57.9 kJ/mol 13.8 kcal/mol\n", "Approx. barrier 28: 38.1 kJ/mol 9.1 kcal/mol\n", "Approx. barrier 29: 69.3 kJ/mol 16.6 kcal/mol\n" ] } ], "source": [ "with scan_lib.reading():\n", " for k in sorted(scan_lib.keys()):\n", " ens = scan_lib[k]\n", " # nrg = [\n", " # x[\"ORCA/SCF_Energy\"] + x[\"ORCA/VdW_Correction\"]\n", " # for x in ens.attrib[\"conformer_attrib\"]\n", " # ]\n", " \n", " nrg = ens.attrib['XTB/Conformer_Energies']\n", " nrg_rel = 2625.5 * (nrg - np.min(nrg))\n", "\n", " eq = 0\n", " ts_min = np.max(nrg_rel)\n", " barrier = ts_min - nrg_rel[eq]\n", " print(f\"Approx. barrier {ens.name}: {barrier:6.1f} kJ/mol {barrier/4.184:6.1f} kcal/mol\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have the approximate structures, we will write out equilibrium structure guesses in one library (`rb_eq_g0.mlib`), and our transition state structures --- into a separate one (`ts_eq_g0.mlib`)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "eq_lib = ml.MoleculeLibrary(\"rb_eq_g0.mlib\", readonly=False)\n", "ts_lib = ml.MoleculeLibrary(\"rb_ts_g0.mlib\", readonly=False)\n", "\n", "with scan_lib.reading(), eq_lib.writing(), ts_lib.writing():\n", " for key, ens in scan_lib.items():\n", " nrg = ens.attrib['XTB/Conformer_Energies']\n", " nrg_rel = 2625.5 * (nrg - np.min(nrg))\n", "\n", " eq1 = int(np.argmin(nrg_rel[:6]))\n", " eq2 = int(np.argmin(nrg_rel[-6:]) + len(nrg_rel) - 6)\n", " ts = int(np.argmax(nrg_rel))\n", "\n", " eq_lib[key + \"_eq1\"] = ens[eq1]\n", " eq_lib[key + \"_eq2\"] = ens[eq2]\n", "\n", " ts_lib[key] = ens[ts]" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "application/3dmoljs_load.v0": "
\n

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n jupyter labextension install jupyterlab_3dmol

\n
\n", "text/html": [ "
\n", "

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n", " jupyter labextension install jupyterlab_3dmol

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/3dmoljs_load.v0": "", "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/3dmoljs_load.v0": "", "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/3dmoljs_load.v0": "
\n

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n jupyter labextension install jupyterlab_3dmol

\n
\n", "text/html": [ "
\n", "

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n", " jupyter labextension install jupyterlab_3dmol

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/3dmoljs_load.v0": "", "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/3dmoljs_load.v0": "", "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# This will show both\n", "%mlib_view rb_eq_g0.mlib 26_eq1\n", "%mlib_view rb_ts_g0.mlib 26" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we need to perform geometry optimizations with ORCA. For transition states we will be using the eigenvector following." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "================================================================================\n", "Optimizing the equilibrium structures\n", "Submitting jobs: 0it [00:00, ?it/s]\n", "Waiting for jobs: 0it [00:00, ?it/s]\n", "Finalizing the calculations: 100%|██████████████| 18/18 [00:00<00:00, 37.60it/s]\n", "================================================================================\n", "Optimizing the transition structures\n", "Submitting jobs: 0it [00:00, ?it/s]\n", "Waiting for jobs: 0it [00:00, ?it/s]\n", "Finalizing the calculations: 100%|████████████████| 9/9 [00:00<00:00, 12.20it/s]\n" ] } ], "source": [ "!python rb_w1_opt.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 3: Extracting preliminary barrier and TS geometry\n", "This is the final step in which we will analyze the obtained outputs" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Barrier for 23: 3.7, 4.0, freq[6]=-74.12\n", "Barrier for 25: 34.9, 35.5, freq[6]=-39.53\n", "Barrier for 28: 16.7, 16.5, freq[6]=-33.68\n", "Barrier for 27: 17.2, 17.7, freq[6]=-35.69\n", "Barrier for 18: 25.9, 25.5, freq[6]=-33.52\n", "Barrier for 29: 24.5, 25.1, freq[6]=-30.06\n", "Barrier for 24: 47.6, 48.1, freq[6]=-24.92\n", "Barrier for 26: 22.1, 21.4, freq[6]=-27.45\n" ] }, { "data": { "application/3dmoljs_load.v0": "
\n

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n jupyter labextension install jupyterlab_3dmol

\n
\n", "text/html": [ "
\n", "

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n", " jupyter labextension install jupyterlab_3dmol

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/3dmoljs_load.v0": "", "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "Molecule(name='26', formula='C26 H26 Cl1 F1 N6 O2')" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# We will analyze the tra\n", "import pandas as pd\n", "\n", "ts_lib = ml.MoleculeLibrary(\"rb_ts_g1.mlib\")\n", "eq_lib = ml.MoleculeLibrary(\"rb_eq_g1.mlib\")\n", "\n", "with ts_lib.reading(), eq_lib.reading():\n", " for key in ts_lib.keys():\n", " if key == \"22\":\n", " continue\n", " ts = ts_lib[key]\n", " freq = ts.attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\"Vibrational frequencies\"]\n", " \n", " G_ts = ts.attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\"Gibbs Energy (Hartree)\"]\n", " G_eq1 = eq_lib[key + \"_eq1\"].attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\"Gibbs Energy (Hartree)\"]\n", " G_eq2 = eq_lib[key + \"_eq2\"].attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\"Gibbs Energy (Hartree)\"]\n", "\n", " dG_act1 = 2625.5*(G_ts - G_eq1)\n", " dG_act2 = 2625.5*(G_ts - G_eq2) \n", " print(f\"Barrier for {key}: {dG_act1/4.184:5.1f}, {dG_act2/4.184:5.1f}, {freq[6]=:5.2f}\")\n", " \n", " ts = ts_lib[\"26\"]\n", "\n", "ts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 4. Replace the substituents on the transition state\n", "\n", "In the previous step we were able to identify a problem: some of the transition states did not converge to what we needed. This can happen, especially considering that XTB geometries can substantially differ from the DFT ones. In order to fix this problem we will use molli's ability to manipulate structures in a streamlined fashion. \n", "\n", "First, we will remove the substituent from atom \"2\" in the direction of atom \"3\" (and write the corresponding file)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "ts.remove_substituent(\"2\", \"3\", ap_label=\"ap1\")\n", "with open(\"ts_core.mol2\", \"wt\") as f:\n", " ts.dump_mol2(f)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "# Create new library\n", "from math import radians\n", "new_ts_lib = ml.MoleculeLibrary(\"rb_ts_g2.mlib\", readonly=False, overwrite=True)\n", "with ts_lib.reading(), new_ts_lib.writing():\n", " for k, m in ts_lib.items():\n", " m.remove_substituent(\"3\", \"2\", ap_label=\"ap2\")\n", " fused = ml.Molecule.join(ts, m, \"ap1\", \"ap2\", optimize_rotation=False,)\n", " fused.rotate_dihedral((\"1\", \"2\", \"3\", \"4\"), radians(180))\n", " fused.name = k\n", " new_ts_lib[k] = fused\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "##" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimizing the transition structures\n", "Submitting jobs: 0it [00:00, ?it/s]\n", "Waiting for jobs: 0it [00:00, ?it/s]\n", "Finalizing the calculations: 100%|████████████████| 9/9 [00:00<00:00, 16.42it/s]\n" ] } ], "source": [ "!python rb_w2_opt.py" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "application/3dmoljs_load.v0": "
\n

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n jupyter labextension install jupyterlab_3dmol

\n
\n", "text/html": [ "
\n", "

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n", " jupyter labextension install jupyterlab_3dmol

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/3dmoljs_load.v0": "", "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/3dmoljs_load.v0": "", "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# In this line one can make sure that all of the transition states look solid now!\n", "%mlib_view rb_eq_g1.mlib 23_eq1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 5. Compile the data in a presentable form." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
barrier (R->TS)barrier (S->TS)TS low freqTS N imag freqEQ1 N imag freqEQ2 N imag freq
18108.3106.7-32.2100
22103.2102.6-36.5100
23132.0133.3-48.3110
24150.0151.9-30.0100
25146.1148.4-35.0100
2692.689.4-27.8100
2773.776.0-38.6100
2869.969.2-35.0100
29101.3103.8-25.5100
\n", "
" ], "text/plain": [ " barrier (R->TS) barrier (S->TS) TS low freq TS N imag freq \\\n", "18 108.3 106.7 -32.2 1 \n", "22 103.2 102.6 -36.5 1 \n", "23 132.0 133.3 -48.3 1 \n", "24 150.0 151.9 -30.0 1 \n", "25 146.1 148.4 -35.0 1 \n", "26 92.6 89.4 -27.8 1 \n", "27 73.7 76.0 -38.6 1 \n", "28 69.9 69.2 -35.0 1 \n", "29 101.3 103.8 -25.5 1 \n", "\n", " EQ1 N imag freq EQ2 N imag freq \n", "18 0 0 \n", "22 0 0 \n", "23 1 0 \n", "24 0 0 \n", "25 0 0 \n", "26 0 0 \n", "27 0 0 \n", "28 0 0 \n", "29 0 0 " ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ts_lib = ml.MoleculeLibrary(\"rb_ts_g3.mlib\")\n", "eq_lib = ml.MoleculeLibrary(\"rb_eq_g1.mlib\")\n", "\n", "with ts_lib.reading(), eq_lib.reading():\n", " data = {}\n", " for key in sorted(ts_lib.keys()):\n", " ts = ts_lib[key]\n", " ts_freq = ts.attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\"Vibrational frequencies\"]\n", " G_ts = ts.attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\"Gibbs Energy (Hartree)\"]\n", "\n", " eq1 = eq_lib[key + \"_eq1\"]\n", " G_eq1 = eq1.attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\n", " \"Gibbs Energy (Hartree)\"\n", " ]\n", " eq1_freq = eq1.attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\"Vibrational frequencies\"]\n", "\n", " eq2 = eq_lib[key + \"_eq2\"]\n", " G_eq2 = eq2.attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\n", " \"Gibbs Energy (Hartree)\"\n", " ]\n", " eq2_freq = eq2.attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\"Vibrational frequencies\"]\n", "\n", " dG_act1 = 2625.5 * (G_ts - G_eq1) # This is the (R) -> TS barrier\n", " dG_act2 = 2625.5 * (G_ts - G_eq2) # This is the (S) -> TS barrier\n", "\n", " nim_eq1 = np.count_nonzero(eq1_freq < 0)\n", " nim_eq2 = np.count_nonzero(eq2_freq < 0)\n", " nim_ts = np.count_nonzero(ts_freq < 0)\n", "\n", " data[key] = {\n", " # \"barrier (expt)\": \n", " \"barrier (R->TS)\": dG_act1, \n", " \"barrier (S->TS)\": dG_act2,\n", " \"TS low freq\": ts_freq[6],\n", " \"TS N imag freq\": nim_ts,\n", " \"EQ1 N imag freq\": nim_eq1,\n", " \"EQ2 N imag freq\": nim_eq2 \n", " }\n", " # f\"{key}: Barrier={dG_act1:5.1f} ({dG_act2:5.1f}) kJ/mol, TS low freq: {ts_freq[6]:5.2f}, TS #im freq: {nim_ts}, EQ1 #im freq: {nim_eq1}, EQ2 #im freq: {nim_eq2}\"\n", "\n", " ts = ts_lib[\"26\"]\n", "\n", "pd.DataFrame.from_dict(data=data, orient=\"index\").round(1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Step 6. Resolving the imaginary frequency problem\n", "\n", "As can be seen from the table above, the project is nearly complete now. The only remaining problem is that for compound **23** we have accidentally converged to the transition state (as indicated by a single imaginary frequency). Upon closer inspection one can notice that this is the tert-butyl rotation. Here we have a couple of options for dealing with this, but we are going to use the similar approach as before: we will create a new collection and we will repeat the calculation." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "application/3dmoljs_load.v0": "
\n

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n jupyter labextension install jupyterlab_3dmol

\n
\n", "text/html": [ "
\n", "

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n", " jupyter labextension install jupyterlab_3dmol

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/3dmoljs_load.v0": "", "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "Molecule(name='23', formula='C28 H29 Cl1 F1 N5 O2')" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "with eq_lib.reading():\n", " mol_fail = eq_lib[\"23_eq1\"]\n", "\n", "# This will rotate the tert-butyl into the correct orientation.\n", "mol_fail.rotate_dihedral((50,49,44,43), 0)\n", "\n", "mol_fail" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "eq_fix = ml.MoleculeLibrary(\"rb_eq_fix.mlib\", readonly=False, overwrite=True)\n", "with eq_fix.writing():\n", " eq_fix[\"23_eq1\"] = mol_fail" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "================================================================================\n", "Optimizing the equilibrium structures (again)\n", "Submitting jobs: 0it [00:00, ?it/s]\n", "Waiting for jobs: 0it [00:00, ?it/s]\n", "Finalizing the calculations: 100%|████████████████| 1/1 [00:00<00:00, 36.99it/s]\n" ] } ], "source": [ "!python rb_w3_opt.py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that the final optimization is done, we will reassemble the library." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "application/3dmoljs_load.v0": "
\n

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n jupyter labextension install jupyterlab_3dmol

\n
\n", "text/html": [ "
\n", "

You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
\n", " jupyter labextension install jupyterlab_3dmol

\n", "
\n", "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/3dmoljs_load.v0": "", "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/3dmoljs_load.v0": "", "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "eq_fix_g1 = ml.MoleculeLibrary(\"rb_eq_fix_g1.mlib\")\n", "rb_eq_g1 = ml.MoleculeLibrary(\"rb_eq_g1.mlib\")\n", "rb_eq_g2 = ml.MoleculeLibrary(\"rb_eq_g2.mlib\", readonly=False, overwrite=True)\n", "\n", "with eq_fix_g1.reading(), rb_eq_g1.reading(), rb_eq_g2.writing():\n", " for k in rb_eq_g1.keys():\n", " if k in eq_fix_g1:\n", " rb_eq_g2[k] = eq_fix_g1[k]\n", " else:\n", " rb_eq_g2[k] = rb_eq_g1[k]\n", "\n", "%mlib_view rb_eq_g2.mlib 23_eq1\n" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
barrier (R->TS)barrier (S->TS)TS low freqTS N imag freqEQ1 N imag freqEQ2 N imag freq
18108.3106.7-32.2100
22103.2102.6-36.5100
23141.0133.3-48.3100
24150.0151.9-30.0100
25146.1148.4-35.0100
2692.689.4-27.8100
2773.776.0-38.6100
2869.969.2-35.0100
29101.3103.8-25.5100
\n", "
" ], "text/plain": [ " barrier (R->TS) barrier (S->TS) TS low freq TS N imag freq \\\n", "18 108.3 106.7 -32.2 1 \n", "22 103.2 102.6 -36.5 1 \n", "23 141.0 133.3 -48.3 1 \n", "24 150.0 151.9 -30.0 1 \n", "25 146.1 148.4 -35.0 1 \n", "26 92.6 89.4 -27.8 1 \n", "27 73.7 76.0 -38.6 1 \n", "28 69.9 69.2 -35.0 1 \n", "29 101.3 103.8 -25.5 1 \n", "\n", " EQ1 N imag freq EQ2 N imag freq \n", "18 0 0 \n", "22 0 0 \n", "23 0 0 \n", "24 0 0 \n", "25 0 0 \n", "26 0 0 \n", "27 0 0 \n", "28 0 0 \n", "29 0 0 " ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ts_lib = ml.MoleculeLibrary(\"rb_ts_g3.mlib\")\n", "eq_lib = ml.MoleculeLibrary(\"rb_eq_g2.mlib\")\n", "\n", "with ts_lib.reading(), eq_lib.reading():\n", " data = {}\n", " for key in sorted(ts_lib.keys()):\n", " ts = ts_lib[key]\n", " ts_freq = ts.attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\"Vibrational frequencies\"]\n", " G_ts = ts.attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\"Gibbs Energy (Hartree)\"]\n", "\n", " eq1 = eq_lib[key + \"_eq1\"]\n", " G_eq1 = eq1.attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\n", " \"Gibbs Energy (Hartree)\"\n", " ]\n", " eq1_freq = eq1.attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\"Vibrational frequencies\"]\n", "\n", " eq2 = eq_lib[key + \"_eq2\"]\n", " G_eq2 = eq2.attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\n", " \"Gibbs Energy (Hartree)\"\n", " ]\n", " eq2_freq = eq2.attrib[\"ORCA/THERMOCHEMISTRY_Energies\"][\"Vibrational frequencies\"]\n", "\n", " dG_act1 = 2625.5 * (G_ts - G_eq1) # This is the (R) -> TS barrier\n", " dG_act2 = 2625.5 * (G_ts - G_eq2) # This is the (S) -> TS barrier\n", "\n", " nim_eq1 = np.count_nonzero(eq1_freq < 0)\n", " nim_eq2 = np.count_nonzero(eq2_freq < 0)\n", " nim_ts = np.count_nonzero(ts_freq < 0)\n", "\n", " data[key] = {\n", " # \"barrier (expt)\": \n", " \"barrier (R->TS)\": dG_act1, \n", " \"barrier (S->TS)\": dG_act2,\n", " \"TS low freq\": ts_freq[6],\n", " \"TS N imag freq\": nim_ts,\n", " \"EQ1 N imag freq\": nim_eq1,\n", " \"EQ2 N imag freq\": nim_eq2 \n", " }\n", " # f\"{key}: Barrier={dG_act1:5.1f} ({dG_act2:5.1f}) kJ/mol, TS low freq: {ts_freq[6]:5.2f}, TS #im freq: {nim_ts}, EQ1 #im freq: {nim_eq1}, EQ2 #im freq: {nim_eq2}\"\n", "\n", " ts = ts_lib[\"26\"]\n", "\n", "pd.DataFrame.from_dict(data=data, orient=\"index\").round(1)" ] } ], "metadata": { "kernelspec": { "display_name": "dev-blake", "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.11.6" } }, "nbformat": 4, "nbformat_minor": 2 }