From b0653bba212e8e2db62d466f622d4515481f3ccf Mon Sep 17 00:00:00 2001 From: Anton Zamyatin Date: Tue, 19 Aug 2025 08:50:32 +0200 Subject: [PATCH 1/6] Add deep learning barrier heights tutorial (ChemTorch) --- .../deep_learning_barrier_heights.ipynb | 1231 +++++++++++++++++ docs/index.md | 1 + 2 files changed, 1232 insertions(+) create mode 100644 docs/contents/notebooks/deep_learning_barrier_heights.ipynb diff --git a/docs/contents/notebooks/deep_learning_barrier_heights.ipynb b/docs/contents/notebooks/deep_learning_barrier_heights.ipynb new file mode 100644 index 00000000..d6ff5626 --- /dev/null +++ b/docs/contents/notebooks/deep_learning_barrier_heights.ipynb @@ -0,0 +1,1231 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "5f7bd60d", + "metadata": {}, + "source": [ + "# Deep Learning Of Reaction Barrier Heights\n", + "This two-part tutorial showcases how ORCA integrates into downstream deep learning workflows by serving as a data source training and evaluation.\n", + "1. First, we will show how to calculate the barrier height of a chemical reaction using ORCA with the ORCA Python inferface (OPI). \n", + "2. Second, we use the ChemTorch framework to train and evaluate a graph neural network (GNN) on a curated subset of the popular [RGD1 dataset](https://www.nature.com/articles/s41597-023-02043-z) which contains precomputed barrier heights.\n", + "\n", + "The term *barrier height*, also called *activation energy*, refers to the energy difference between reactant and transition state (TS)." + ] + }, + { + "cell_type": "markdown", + "id": "ab7189fb", + "metadata": {}, + "source": [ + "## Part I: Obtaining Barrier Heights via ORCA/OPI\n", + "The first part of this tutorial consists of six steps:\n", + "1. Import Dependencies\n", + "2. Define Base Directory\n", + "3. Pick a Reaction\n", + "4. Setup and Run Calculations\n", + "5. Visualize the TS and IRC\n", + "6. Compute the Barrier Height" + ] + }, + { + "cell_type": "markdown", + "id": "1a2f1185", + "metadata": {}, + "source": [ + "### Step 1: Import Dependencies\n", + "We start by importing the modules needed for:\n", + "\n", + "- Interfacing with ORCA input/output\n", + "- Numerical calculations and data handling\n", + "- Plotting results\n", + "\n", + "You might first need to install the additional required packages via\n", + "```bash\n", + "pip install py3Dmol rdkit matplotlib pandas\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "5edf5999", + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "import shutil\n", + "\n", + "# > OPI imports for performing ORCA calculations and reading output\n", + "from opi.core import Calculator\n", + "from opi.output.core import Output\n", + "from opi.input.structures.structure import Structure\n", + "from opi.input.arbitrary_string import ArbitraryStringPos\n", + "from opi.input.simple_keywords import SimpleKeyword, BasisSet, Scf, Dft, Task, Neb, DispersionCorrection, Wft\n", + "\n", + "# > for visualization of molecules\n", + "from rdkit import Chem\n", + "from rdkit.Chem import Draw\n", + "import py3Dmol\n", + "import matplotlib.pyplot as plt\n", + "\n", + "# > Pandas for data handling\n", + "import pandas as pd" + ] + }, + { + "cell_type": "markdown", + "id": "c7d99c35", + "metadata": {}, + "source": [ + "### Step 2: Define Base Directory\n", + "We will create a `RUN` folder which will serve as the base directory for all ORCA calculations." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "f24e7bb2", + "metadata": {}, + "outputs": [], + "source": [ + "# > Calculation is performed in `RUN`\n", + "base_dir = Path(\"RUN\")\n", + "base_dir.mkdir(exist_ok=True)" + ] + }, + { + "cell_type": "markdown", + "id": "34efb8cb", + "metadata": {}, + "source": [ + "### Step 3: Pick a Reaction\n", + "For this example, we will remove a datapoint from the RGD1 dataset and calculate its barrier height from scatch.\n", + "We chose the tautomerization of acetone.\n", + "\n", + "We will use the reactant and product geometries provided in the RGD1 dataset as a starting point for our calculations.\n", + "Please refer to the [original paper](https://www.nature.com/articles/s41597-023-02043-z) for a detailed explenation of how these structures were obtained." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "bcc52329", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABAoAAAHqCAYAAACN2KwJAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAdJlJREFUeJzt3Xd8VFX+//H3pPdCElqAhA7SizRpoQgooCyCChZA17rqd+3rWtddXV3brr0gIiogNnqV0KUjvZdQAwmk18nM+f0xv1yICRBCQgqv5+PBQ3Pnzr1nJpM5n3nPuefYjDFGAAAAAAAAktzKuwEAAAAAAKDiICgAAAAAAAAWggIAAAAAAGAhKAAAAAAAABaCAgAAAAAAYCEoAAAAAAAAFoICAAAAAABgISgAAAAAAAAWggIAAAAAAGAhKABwXi+//LJsNpuWLFlS3k05L5vNpt69e5d3MwAAqDAOHTokm82mMWPGlHdTKpzK8NyMGTNGNptNhw4dKu+m4CpGUIASyX+TPfefp6enIiMjNXLkSK1fv768m3hB5fnhcsmSJbLZbHr55ZdLdD+bzaaOHTued7+5c+da+/EBGgCAslNUPeTl5aW6detq1KhR2rJlS3k38YqJjo5WdHR0ie6X/9xt27atyH0cDociIyOt/fgADZQ9j/JuACq3hg0b6o477pAkZWRkaMOGDZo2bZp++eUXLVq0SD179iznFlY9Hh4e2rBhg7Zs2aLWrVsXun38+PHy8PBQXl5eObQOAICrz7n1UHp6ulavXq3Jkyfrp59+0q+//qrrrruunFtYsbm5ub67/PLLL/XOO+8Uun3u3Lk6fvw49Q1wBTGiAJelUaNGevnll/Xyyy/rP//5jxYvXqzXX39ddrtdL7zwQnk3r0oaMGCA3Nzc9OWXXxa6LTExUTNnztQNN9xQDi0DAODqdG499NZbb2nFihX6+9//rpycHP39738v7+ZVeJ6enurbt6+++eYb2e32Qrd/+eWXCg4OVrdu3cqhdcDViaAApe6ee+6RJG3YsKHQbbm5uXrnnXfUvn17+fv7KzAwUD169NCMGTMK7btnzx49/fTTat++vcLCwuTj46MmTZro2WefVXp6epHnTktL0yuvvKLWrVvLz89PwcHBateunV544QXZ7XZr+L4kLV26tMBQwa+++kqSlJKSojfeeEO9evVS7dq15eXlpdq1a+uuu+7S/v37C53z3Ov4v/vuO7Vt21a+vr6qVauWHnvsMWVlZRXYNyYmRpL0yiuvFDh/cYfR1alTR/3799e3336r3NzcArd98803ys3N1bhx4857/8TERP3f//2f6tevL29vb1WvXl0jR44873C/89myZYtuu+021apVS15eXoqKitIjjzyi06dPF7n/5s2bNXr0aNWpU0fe3t6qVauWBg4cqJkzZ1r7XGhOhK+++qrA7+liLuW1BgBAaXvkkUckSevWrbO25V8WeOzYMd11112qWbOm3NzcCvR7EyZMUOfOnRUQEKCAgAB17tz5vH2fw+HQG2+8oUaNGsnHx0eNGjXS66+/LqfTWeT+F7os8XyXDuTm5urdd9/Vtddeq8DAQAUEBOiaa67R448/rqSkJOvyi7i4OMXFxRWobS7lMstx48YpISGhQF0gSQkJCZo1a5Zuv/12+fr6nvf+l/K8nU9aWppeeukltWjRQr6+vgoJCdGAAQO0YsWK8+5/obpTuvicCJd6qeiyZcs0ZMgQhYeHy9vbW40bN9bzzz+vzMzMS3qswMVw6QHKjIdHwZdXTk6OBg4cqCVLlqht27a65557ZLfbNXv2bN100016//339Ze//MXa/6efftL48eMVExOj3r17y+l0avXq1XrjjTe0dOlSLVu2TJ6entb+p06dUq9evbRr1y61bdtWDz74oJxOp3bt2qU33nhDTzzxhKKjo/XSSy/plVdeUVRUVIE37bZt20qSdu7cqRdffFExMTEaNmyY/P39tWvXLn333XeaPXu2Nm7cqKioqEKP94MPPtC8efN00003qU+fPpo3b57+97//KTExUd9++60kqXfv3jp06JAmTpyoXr16FegYQkJCiv3cjhs3TvPnz9fMmTM1fPhwa/uXX36pFi1aqHPnzkXeLyEhQV27dtX+/fvVu3dv3XbbbTp48KB++OEHzZ49W/Pnz1f37t0vev4ZM2Zo5MiRcnNz00033aS6detqx44d+uCDDzR//nytWbNGoaGh1v4//vijRo0aJWOMhgwZoqZNm+rUqVNas2aNxo8fryFDhhT7sRfHpb7WAAAoK/lfUOQ7ffq0unbtqmrVqum2225Tdna2goKCJEmPPvqo3n//fUVGRlpfvPz4448aO3asNm3apP/+978FjnXffffpyy+/VP369fXwww8rOztb77zzjlatWlUqbc/KylL//v21cuVKNW7cWGPHjpW3t7f27t2rTz/9VHfddZdVW7333nuSpP/7v/+z7n8pH4CHDRum0NBQTZgwQX/605+s7ZMmTZLdbte4cePOO1r1Up+3opw5c0Y9e/bU9u3bdd111+mBBx5Qamqqpk+frpiYGE2bNk0333yztX9x6s5Lqe2K4+OPP9bDDz+skJAQDRkyRNWrV9f69ev1r3/9S7GxsYqNjZWXl1epnhNXMQOUwMGDB40kM2DAgEK3vfbaa0aSufHGGwtsf+6554wk88ILLxin02ltT01NNR07djReXl7m2LFj1vajR4+anJycQsd/5ZVXjCTzzTffFNg+fPhwI8k899xzhe4THx9v7Ha79bMk06tXryIfW3Jysjl9+nSh7YsXLzZubm7m3nvvLbD9pZdeMpJMcHCw2bVrl7U9MzPTNGnSxLi5uRV4XLGxsUaSeemll4o8//nk3+/+++83OTk5JiwszNxwww3W7WvXrjWSzNtvv21OnDhR5GMcO3askWT+9re/Fdg+e/ZsI8k0atTIOByOQo8tNjbW2paYmGiCgoJMZGSkOXToUIHjTJ482Ugyf/nLX6xt8fHxxt/f3/j7+5uNGzcWelxHjhy54PnyTZgwwUgyEyZMKLC9qMd5qa81AABK4kL10IsvvmgkmZiYGGubJCPJjB071uTl5RXYf+nSpUaSad68uUlOTra2nzlzxjRp0sRIMsuWLbO259cFbdq0Menp6db2o0ePmvDwcCPJ3H333QXOcaH6JyoqykRFRRXY9sQTTxhJ5s477yzU3uTkZJOWlnbB+xdHVFSU8fb2NsYY85e//MV4eHiYEydOWLe3aNHCtGrVyhhjzIABA4wkc/DgQev2S33e8n9nf3xuRo0aZSSZzz//vMD2kydPmrp165qIiAiTlZVlbS9u3Xm+8+Ur6ndy9913F3qc27dvNx4eHqZNmzYmMTGxwP6vv/66kWTeeuutIs8BlASXHuCy7Nu3z7om76mnnlKfPn303HPPqUaNGvrPf/5j7ed0OvXxxx+rYcOG1pD7fIGBgXrxxReVm5urn376ydoeGRlZZCqa/03wokWLrG3x8fH66aef1LBhwyKHudWoUaPQCIfzCQ4OVrVq1Qptj4mJUYsWLQqc91yPPfaYmjZtav3s6+ur22+/XU6ns8jLMC6Hl5eXRo8erfnz5+v48eOSXKMJPD09deeddxZ5n9zcXE2ePFlhYWF6/vnnC9x2ww03qH///tq3b59Wrlx5wXN//fXXSk1N1euvv15oZMVtt92m9u3ba8qUKda2iRMnKiMjQ0888YTatWtX6Hh16tQp1mMurpK81gAAuBx/rId69uypf/zjH/Lx8dG//vWvAvt6eXnpzTfflLu7e4HtEydOlOS6DC84ONjaHhoaqpdeekmSCgyl//rrryVJL774ovz9/a3tkZGReuyxxy77MeXl5emzzz5TcHCw/vvf/xZqb3BwsAICAi77POcaN26c8vLyrOdizZo12r59+wUvqbzU560oiYmJmjp1qvr06aN77723wG3Vq1fXU089pYSEBKsGLM26s7g+/fRT5eXl6f3331dYWFiB255++mlFRERo8uTJpXpOXN249ACXZf/+/XrllVcKbKtZs6aWL1+uRo0aWdt2796tpKQk1a5du9D+kmtIvCTt2rXL2maM0YQJE/TVV19p27ZtSklJKXDNXf4HZElav369jDGKiYkpcDlCSS1ZskTvvfee1qxZo8TExAIz7J5vSFeHDh0Kbcv/EJycnHzZbfqjcePG6X//+58mTpyov/71r5oyZYoGDx6siIgIxcfHF9p/165dys7OVkxMjPz8/ArdHhMTo4ULF+r3339Xjx49znve1atXS3J13kXN2ZCdna3ExEQlJiYqPDxca9eulSRdf/31JX2ol6QkrzUAAC7HufWQp6enatSooVGjRunZZ59Vq1atCuxbv359hYeHFzrGpk2bJBU9XD9/fqPff//d2rZ582ZJKrLPvlA/Xly7du1SWlqa+vXrV+BywrLUrl07tW3bVhMmTNAzzzyjL7/8Ul5eXtaKEkW51OetKOvWrZPD4VBOTk6RH/z37t0ryfWcDB48uNTrzuLIr7/mz5+vX3/9tdDtnp6e1DYoVQQFuCwDBgzQvHnzJLk+gE2cOFHPPPOMhg4dqrVr11pJ85kzZyRJ27dv1/bt2897vIyMDOv/H330UX3wwQeqW7euhg4dqlq1asnb21uSayLAnJwca9+UlBRJrhT9ck2bNk233nqrAgICNGDAAEVHR8vPz8+aSC8uLq7I++VfX3iu/DTZ4XBcdrv+qE2bNmrfvr0mTJigevXqKTk5+YKJe2pqqiRXyl2UWrVqFdjvfPJ/lx9++OEF98vIyFB4eHip/m6KoySvNQAALse59dDFnK8fTk1NlZubmyIiIoq8j81mK9BHp6SkyM3NrcjQ4XznuBRXuv/ON27cOD366KNatGiRpkyZYk3cdz6X+rwVJb92WLly5QVHVubXDuXx3OS38Y8jVICyQlCAUhMREaEnn3xSKSkp+uc//6nnn3/emtgm/0P08OHD9cMPP1z0WKdOndKHH36o1q1b67fffivwDXh8fHyhb4rzJ4s5duzYZT+Ol19+WT4+PtqwYYMaN25c4LZzh9RXBPfcc48efvhhPfPMM6pdu7YGDRp03n3zfwcnT54s8vb8UQhFBR5FHWfr1q1q2bLlRdt47u+mqNmUz5W/jnJRayTnd8oXc6mvNQAArqQ/Tm6YLygoSE6nUwkJCapevXqB206dOiVjTIE+Ojg4WE6nU4mJiYU+JJ+vr7fZbEX2sZKrnz136H5p1laXYvTo0Xrqqac0ZswYpaamWpMTns+lPm/nO4YkPfHEE3rrrbcu2sZLeW5Ko7Y5t42pqakKDAws9v2AkmKOApS65557TrVr19ZHH31kLfnXvHlzBQUFaf369UWuj/tHBw4ckDFG/fr1KzRMfvny5YX279ixo9zc3BQbG1us47u5uZ33W/79+/erefPmhUKCEydO6MCBAxc99sXkX+NXGqMMRo0aJR8fH2uZpT9eP3iuZs2aycfHR+vWrStyCZ38pZnyV384n/wVFX777bditbFTp06SpAULFlx03/yhjUV1vPlDCy/mUl9rAABUBPnz+BS1RHBRfXSbNm0kFV0XFbVNcvWzRfWxhw4dKnSZZNOmTRUUFKR169YpKSnpou13d3cvldqmWrVquvnmm3Xs2DFFRkZqwIABF9z/Up+3olx77bWy2WzFrm0upe68UKhQ3NpGOlt/5V+CAJQ1ggKUOl9fXz3zzDOy2+169dVXJbmG4D/44IOKi4vTk08+WeSb6rZt23Tq1ClJsibJW7VqVYF5CY4ePaq//e1vhe5bo0YNDR8+vMg5EyRXonxuklutWjUdPXq0yPZHRUVp3759BdL47OxsPfjgg6XywTN/osQjR45c9rFCQkI0f/58/fzzz/rrX/96wX29vLx0++23KzExUa+//nqB2+bNm6f58+erUaNGuu666y54nLFjxyowMFB///vfixzan5mZWaATu/vuuxUQEKC33367yGsEz+04r732WkmuCZrO/b3/9ttv1hKTF3OprzUAACqCu+++W5Lr8so/XmKQX9vk7yPJmrz4H//4R4HL6Y4dO3be5QCvvfZaHTp0SEuXLrW25ebm6vHHHy+0r4eHh+6//36lpKToscceKxQCpKSkKD093fq5WrVqSkxMVHZ2drEf8/n8+9//1s8//6xffvnF+kb+fC71eStKzZo1NXLkSK1atUr/+c9/ZIwptM+aNWusL1oupe4MCgpS06ZNtWLFCu3bt8/aJy0trcia9nweeugheXh46JFHHtHhw4cL3Z6cnHxJwQNwMVx6gDJx33336Y033tDXX3+t5557zpqBfuPGjfrf//6n2bNnq2fPnqpevbqOHTumrVu3avPmzfrtt99UvXp11apVS8OHD9ePP/6ojh07qm/fvjp58qRmzZqlvn37FjmJ3kcffaRt27bpX//6l+bMmaM+ffrIGKM9e/ZowYIFOnnypJXq9unTR99//71uvvlmtWvXTu7u7ho6dKhat26tRx55RI888ojatWunW265RXl5eVq4cKGMMWrTpo01eVBJNWvWTLVr19aUKVPk7e2tOnXqyGaz6ZFHHikw5K+4evbsWex933jjDS1dulT//Oc/tWrVKnXu3FmHDh3StGnT5OfnpwkTJly0Q86fVXfEiBFq06aNBg4cqGbNmiknJ8cqPrp162Zdq1m9enV9/fXXuu2229SpUycNHTpUTZs2VWJiotasWaPo6Gj98ssvkqQuXbrouuuu0+LFi9W1a1f17NlTcXFxmj59uoYMGaKff/65WI/zUl5rAABUBD179tQjjzyi999/Xy1bttTw4cNljNGPP/6oo0eP6tFHHy3Q58fExGjs2LGaMGGCWrVqpWHDhiknJ0dTp05Vly5dNGvWrELnePzxx7VgwQLdcMMNuv322+Xn56eFCxcqJCTEmqvoXP/4xz+0evVqTZo0SatXr9agQYPk7e2tAwcOaN68eVqxYoX1bX2fPn20fv16DRo0SD169JCXl5d69ux5SXVKvujo6IterljS5+18PvroI+3evVtPP/20Jk2apK5duyokJERHjhzR+vXrtXfvXp04ccIa6XopdecTTzyh++67T127dtWIESPkdDo1d+5c6wuS4mjZsqU++ugjPfjgg2ratKluuOEGNWzYUGlpaTpw4ICWLl2qMWPG6JNPPin2MYELKp9VGVHZXWjd4Hzvv/++tfZuvry8PPPpp5+a6667zgQFBRlvb29Tr149M3DgQPPxxx8XWAc4LS3NPPHEEyY6Otp4e3ubxo0bm1dffdXk5uaedx3glJQU88ILL5hmzZoZb29vExwcbNq2bWtefPFFk5uba+134sQJM3LkSBMeHm7c3NyMJDNhwgRjjDFOp9N88sknpkWLFsbHx8fUrFnT3HPPPebUqVOmV69e5o9/Ni+99JKRZGJjYwu1Z8KECQWOnW/16tWmV69eJjAw0FpT+dy1couSv17y/ffff8H98h/f+Z6jhIQE8+ijj5qoqCjj6elpwsPDzS233GK2bt1aaN8LPbZdu3aZe+65x0RFRRkvLy8TGhpqWrVqZR599FGzdu3aQvtv2rTJjBw50tSoUcN4enqaWrVqmUGDBplZs2YV2C8xMdHcddddplq1asbX19d06dLFzJ8//7zP5fke56W81gAAKIni1EPnOl+fda4vv/zSXHvttcbPz8/4+fmZa6+91nz55ZdF7puXl2def/1106BBA+Pl5WUaNGhgXnvtNbNv3z4jydx9992F7jNt2jTTqlUr4+XlZWrWrGkeeeQRk5aWZqKiokxUVFSh/bOzs81bb71l2rZta3x9fU1AQIC55pprzBNPPGGSkpKs/dLS0syf//xnU6tWLePu7m4kmZdeeumiz0lUVJTx9va+6H7GGDNgwIDz1kzFfd7yf2dFPTeZmZnmzTffNB06dDD+/v7G19fX1K9f39x8883m66+/Nna7vcD+xa07jTHmww8/NI0bNzaenp6mXr161j5FvSbuvvvu8z7OtWvXmttuu83Url3bquPat29vnn32WbNz586LPodAcdmMKWJsDQAAAAAAuCoxRwEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAAtBAQAAAAAAsBAUAAAAAAAAC0EBAAAAAACwEBQAAAAAAAALQQEAAAAAALAQFAAAAAAAAItHeTcAuKo5ndLx49KWLdLGjdLevVJSkuuf3S5VqyaFhkp160pt20odOkhRUZKXV3m33NX2o0eljz6SjCl4m5ub1KaNdNtt5dM2AABQ8VTmuscYVxtff13Kyipc+zRuLN14o1SrVtH3X7tWWrNGOnHC9XNwsDRihNSgQcnac/SoNHOmdOTI2eNde63Up0/Jjgf8AUEBUB4cDmn5cumnn6RNm6RTp6S0NFfHk5fn+ud0Sp6ekoeH5O0tBQS4Os/GjaWBA6Wbb3Z1CuUlJ0d67TXp558L3+buLg0bRlAAAACqRt1jt0vffCN9+qnr8fxRz55St27nDwp27ZJ++UXavt31c+3aUvfuJQ8KEhJcQcHGjWeP5+5OUIBSQ1AAXEnGuJLkTz+VVqyQduyQEhNdHWRRcnPP/v+pU5LNJu3bJ+3e7UqlH3lEat78yrT9XFlZ0rJl0vTprnb9kZublJJy5dsFAAAqjqpS9+TlSYcPu0ZR5o8I+KOkpPM/LknKzJROn5ZOnnT97OlZ8PFeKru94PG8vKSMjJIfD/gDggLgSjFGio+X3ntPmjLFNfQuP5H28HCl5FFRUo0akp+f68N2ZqZ05owUF+fqWHNzXR/Af/9dOnTI1SE9+qjUooWrM70SHA7X4xg/3vVfyZX4e3qeHToIAACublWl7jFGSk2Vvv/e1Q7J1daoKFdokJ19ZdoBXGEEBcCVYIyrI5k6VfrkE9dwO8k1RCwy0pWOt2jh+m+9eq7hdh4ero7p5EnXcLVt21zD1Q4fdg37T0qSvv5a8vWV/v53KTzc1XGVtaQkackSae7cs9v69XN14BkZBAUAAFztqlLdk5npmlNh0qSzQUfLllLTplJsLEEBqiyCAuBKyMtzdXj//vfZztLDw3Xd3eDB0i23uCbscXc//zF27nQN9f/5Z2nrVtfw/5wc6cMPpY4dpeHDXZ1nWSbsOTmuYYPjx0vp6a5tYWGuoYCvvXbl0n0AAFBxVZW6x+FwTRY4ZYorvJBc8yfcc49rtMTKlWV3bqCcsTwiUNacTik5WXr11bPXkUmuIWvPPis984zUqdOFO0vJlbr/3/9J//iH1K7d2RTd4ZCef96VuF/o2rjLZYxrKODcuWc7Rjc314SFnTtXjBmJAQBA+aoqdY/kGuGwerX07bdnt7VqJd11l2sUBFCFERQAZS0z0zUBz8yZZ7d5e7uW1xk61PWNfHH5+Ei9e7vuGxFxdvvhw9J337mG5f1xuZ7SkpUlLV4sffGF62ebTapeXXrxRdf8BAAAAFWl7nE6pVWrXBMx5o+i9PR0jZIIDGQUJao8ggKgLBnjmoznq68Kbr/nHum666SQkEs/pre369q4F14ouP3rr13nKgvGuEYSTJly9hw+PtI//+nquOksAQBAVap7tmxxXfawerVrm6enaxRlnz4XHw0BVAEEBUBZysuTjh6VZs06u81mk8aNc03CU1KhodKwYa41c/PFxbnWKC6LTnPPHtfax8uWuX728ZF69JDGjCn9cwEAgMqpqtQ9mZnStGnSjz+6fnZzk2rVkt58s/TPBVRQBAVAWUpKcs3Ym5NzdlvXrlJ09OUP1w8IkEaMKLjt999dywqVJrvdtW7wr7+6EnabTapZ03Xt4ZWYbRgAAFQOVaHukaSJE6XZs11zLUiuyyWee851ySVwlWDVA6AsJSW5Vgk4V8eOrmF0lzNc32ZzdbidOhXcvn372U5Ncq0/vH+/aw3jfNWrS088Ufzhf1OnSkuXSgkJrp9r15buuENq3ZpLDgAAwFnlXfcY41qN4OWXz27z8HB9uVGtWvHOtWWLa7WF/FUOAgOlLl2kW2/lCxJcVQgKgLKUliYdOFBw2zXXlM61bZ6ermv2znXgwNlliCTXEMCTJ12XDeSrX1966KHiBQWHD0vffCPt2+ea1MfX1xUQ3HGH6/IDAACAfOVd9xjjWqng3LrHy0v629+KFxTk5EiffSZt3uz6f5tNatBAuu8+KTj48h8DUIkQFABlKSdHOn264LY6dUqnw3Rzcx3LZjs74+/p065z5l8iYIzr0oFzr98LDnZ96L8QY1z/Jk1yDevLyHBtb9pUuvlmqWHDy28/AACoWsq77pFcX5KcW/d4erqWVCyOuXNdKzzlP4Y6daT+/V0TMZb2KMqsLNccC0lJJbv/vn0lvy9QDAQFQFmy288uqZMvNLR0OhubTfL3dw3ny852bcvJcf1zOFxD7S7H5s3S99+fvfYvIsI1geH111/+sQEAQNVTWeseY6RTp1yrNRw+7AobvL2l9u2lm25yPYbSlp4u/fKLtGRJye6fluYaNQqUEap9oCzl5blmzj2Xv3/pdZj5nWZ+h2mMK6G2210dppubFBQkdehw9n516154QiFjXMcbP9612oHd7vomoGNHacAA14REAAAAf1TedY/kukzy3LrH09N1+cH55I+inDtXio09O4qyQQMpJqbwvAilJSdH2rSpbI4NlAKCAqAsOZ2uzutclzvr7x95exf82W4/O8TOy0tq0kR6++2zt/v6Xvg6O2Nclxt8/fXZjrhGDWnQIKlnz1JtOgAAqELKu+6x2Vw1y7l1j5vbxUcEnDrluk9+SODnJw0c6Kp9LhQyXA6bzVWTlXQkhMPhqtOKe1kFcIkICoCy5OHh6gTOlZVVOsfOT8DzO7V8fn5nOzV3d1fn2KtX8Y+ZmSk9//zZoYPu7tLIkVKfPq6ZfwEAAIpS3nVP/oiDS6l7cnOlDz5wrdaQP4fTdde5goImTUqn7UXx83N9AVOrVsnun5AgrV59dlUqoJQRFABlycPDte7vudLTz07Cc7mMKTjbr7u7q4MuSXqfP3xv8mTXRD75WrSQ/vQn16zFAAAA51PZ6h6HQ9q9W3rrrbMhQWCgNHas1K1b6bT5fEJDpaeecl3eUBJr10qPPEJQgDLDYqBAWfLycs0RcK4TJ0pnmJjT6Roqd+4KBkFBJR8il5cnHTokPffc2W3u7tILL7iWRAQAALiQylT3SK5VAx56yDVfQL7/+z9XSODvX/LjAlUAQQFQloKCCi8luHVr6XSYeXnSli0FtzVq5ErCSzJp0NGj0osvFlxS6P77XZ1lUFDpLwsEAACqlspU9yQmulY5WLXq7La2baVRo1wTP1P34CrHpQdAWQoJkZo3L7jt998LT/RTEnZ74dlyW7RwnfNSJSS4LjeYPbvg9okTpZ9+ck0EdDFnzpxN5J1O1/0iI10/e3pKb74p3XJL8Y4FAAAqn8pS9+TmSrt2Se++W/CyiO3bpd69XSMqLyY9veBSkCtWSP37uy6/cHd31Tz/+U/xjgVUQAQFQFkKCZGaNXN9OM4fKrd6tSvFDgi4vDV/c3KkOXMKbitph5me7rrsIH+Vg3wZGYUnDSqurKyzExh5eBReLgkAAFQtlaXuyc6WTp50/TuX3V5426W079Qp1/+7u0vJySU7DlBB8NUeUJa8vV3D1869xj8pSVqypGAKfamysqRt21zJd77QUKl9+4svAVQUh6Pg9XkAAACXqrLUPfnLOJ473wGAAhhRAJQld3cpIkK64QbX0DvJNcTtu++ka68tuKRPcRnjSqwnTy4482+3blL9+pKPz6W3MyhI6tzZNSdBSc2bJx0/7up4bTbXdYN9+rhuc3eXmjblej8AAKqyylL3eHm56pTLqXs2bXKFF/kjJiMjXW2qVs01oqJzZ+oeVGoEBUBZCw52rcX7xRdnh6StWOGaDyAwUIqKKv71a8a45hOIjZWmTz+73ddXGj7c1TmVpFMKDXUtz9O+/aXfN9/+/a625QcF11wjPfOM6zabTQoLK/mxAQBA5VAZ6h4fH9dcCvl1Skl89JEUF3c2KIiOlsaNc30xIrkutSAoQCVGUACUNR8f1xC8G2+Uvv7aNczfbnd1oIGBru21a188EbfbpdOnXZ3lp59K8fFnb+vYUbr++sJrFzscrqF+u3ef3ebt7eocz030PT1dnW21aiV/nH5+BScqDAhwJf0AAODqUZ51jzGu+Qe2bj27zWZztcfb++w2NzfXff94/0tRrVrBwMPHx/W4qH1QRRAUAGXNZnN1RC+84LpG79AhV0d28KD0xhuuNHrECKlxY1cn5uXl6sBsNlfnmpfnmj8gPl6aMUOaNMk11C1fcLD08suuoX5/TOhzcqSNG6UhQ85ua9hQmjvX1ZkBAACUpvKse4xxHT//0kfJdY7ff3fNnQCg2AgKgCvB3d011O7996W77nJN7GOMdOyY9M470rRpruV4eveW2rRxDdP39nbNvLt/v7R8ubRggbR3b8Elhnx9peeek3r1KnoYnzGuDvfclQsyMgouBQQAAFCayqvukVxhw7l1T24ukxYCJUBQAFwpbm6uyX3mzJFGjXIl3g6Hq+M8fNg1PO/rr4t/vNBQ6d//lu67r+zaDAAAUBLUPUClxvKIwJVks0mdOkmLF0sPPOCaIfdSVa8ujRkjLV0q3XtvqTcRAACgVFD3AJUWIwqAK81mc10n9+ab0iOPSOvWuWYDXrtW2rXLtVbwuTw9XTPpduwode8udekiNWvmGqLndpGsz9dX6tHDldzn8/BwdbqlrWVL1/DAzExXu5o3L/1zAACAyuVK1j02m9SkScG6x2aTatYs9YelunWlrl2lxETXz23aSP7+59+/dm2pQwfXHAuSqxYLCSn5+QMDXatV5U8KWaOG63kDSonNGC5WBsqN3e6anTcz09VRZme7/puc7JrIJyzMNSGQt7frQ7+fn+u/xV2DOP/P+9xr82y2s/9KU0KCq81Op+vYfn4siQgAAM4q67pHctU+f5yTIH+yxNKUmupaWSovz/Wzj4/rg//52pqe7vqXm+v62cPDtXLCxVZ/OJ+cnLPPW/7x/P3PBhHAZSIoACoSY84uI+R0ujobDw/W4QUAAFUPdQ9QYREUAAAAAAAAC5MZAgAAAAAAC5MZAmXN6Tx7TZ6Xl+tatNIaUmeM69q4zEzX5D8+Phef6AcAAKCiMMZ17b7T6ZqfoCzmEwBwyfhEAZQ1u921dvDChdKOHWcnGCwtJ064jr15syuMAAAAqCwcDmn5cmnBAikjo/TrJAAlQlAAlLW0NGnGDOnWW6VPPy08E+/lcDqlZctcx37rrbNL9AAAAFQG2dnSgw9KI0dKR464ggMA5Y6gAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFg8yrsBQFVnithmK81j22yuf6V8bAAAgLJWlnUSgJIjKADKmMNmU7q/v9xatFBeSIhCbaXb/aV6esqtRQvlRkQo2N2dP2oAAFB52GxKadhQ7oGB8vLwkFcp10kASoZLD4AylmKMJmZmqt7OnXopOVkOU1R2XjJOSfPsdtXduVOPJiYq3uEotWMDAACUtQxj1O3gQUXu3KkDDkep1kkASo6gAChjRpLdGKU5ncopg+PbjVGaw6FsY+Qsg+MDAACUpQynU2lOp5zGFHkpAoArj6AAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWDzKuwFAVefu7q6IiAi1atVKderUkc1mK7Vj22w2hYaGqnXr1oqOjpaXl1epHRsAAKCsubm5qXnz5goODpaPj0+p1kkASs5mjDHl3QigKsvJydGhQ4e0ZcsWNWjQQO3atZObW+kM5jHG6PDhw1q7dq0iIyPVtm1b+fn5lcqxAQAAylpeXp7mz5+vjIwMDRw4UAEBAaVWJwEoOYICAAAAAABgIa4DAAAAAAAWggIAAAAAAGBhMkOglBhjlJeXJ6fTKWOM8q/qsdlsstlscnNzk7u7u/VzWXE6nXI4HHI6nXI6ndZ2m80md3f3K9IGAABQtVH3AFUbQQFQQudO72GMUUpKin777Tft379fCQkJSk5Olt1ul5+fnyIiIlSvXj31799fYWFh1v1Kq9PKb4sxRgcPHtTmzZu1a9cuxcfHKzs7W56engoJCVHLli3VpUsX1alTRx4eHqXaBgAAUHVR9wBXF4IC4DLY7XbNmDFD7733ntatWye73a4LzQ/q7u6unj176sEHH9SAAQMUFBR02W3IP9/MmTP13nvvaf369UpLSzvv/m5uburevbv+8pe/qHfv3oqIiLjsNgAAgKqPuge4erDqAVACmZmZ2rhxo55++mlt3LhRdru9wHC3C3Fzc5Onp6f69++vZ555Rp06dZKXl1eJ2pGf6D/88MOaOXOmMjIyitUONzc3eXt7609/+pPuv/9+9ejRo0TnBwAAVR91D3D1YUQBcIkSEhI0b948vfrqq4qLi1Nubm6B20NCQlSrVi2Fh4fL29tbmZmZio+P14EDByS5rqXLycnRokWLlJCQoMcee0zDhg2Tj4/PJbXD6XQqOTlZd955p5YuXaqMjAzrtoiICNWtW1c1a9ZUUFCQ8vLytHfvXm3ZskXGGDmdTmVlZennn3+W3W6Xw+FQ7969L/u5AQAAVQt1D3B1IigALkFSUpJmzpypt956SwcOHJDD4bBu69Wrl2666Sa1bt1aoaGh8vLykpubmxwOh7KyshQXF6c5c+ZoypQpys7OVnZ2tjZv3qy5c+cqOjpaXbt2vaS2ZGRk6M0339SyZcuszrJly5YaPny4unbtqpo1a8rLy0seHh4yxigzM1M7duzQW2+9pW3btslutyszM1OxsbGKiIhQ27ZtFRISUppPFwAAqMSoe4CrF0EBUEwOh0PLli3TlClTtHv3bmuom5eXlx5++GHdeOONat68ucLDw+Xp6VlgshxjjK655ho1atRI11xzjf73v//pxIkTatu2rbp166a6deteUluys7O1Y8cOTZkyRenp6ZKkTp066eGHH1b37t1Vu3btQkm9MUb169dXYGCg/va3v2nPnj2y2+06ffq0tmzZos2bN6tXr16X+SwBAICqgLoHuLoRFADFdODAAS1cuFBr1661OktPT0+NHDlS9913n5o0aSI3N7ci72uz2RQQEKC2bdsqMjJSmZmZ2rNnj/r06aN+/fopMjLyktpis9nk4+Oj5s2b68yZM/Lw8NC9996rwYMHKzQ0tMgZfW02m4KDgzVgwAAtWrRIp06dUkJCgpxOp06dOkWHCQAALNQ9wNWNoAAophUrVmjNmjVKSUmRJHl4eKhp06Z64okn1LRp02Itt2Oz2VS9enXdf//9io+PV1RUlEJDQy+5LV5eXmrSpIn++te/qk6dOnI6nfrTn/6kkJCQi7bDy8tLPXv21Lx585SQkCBJSk9P19GjRy+5HQAAoGqi7gGubgQFQDFkZmbqt99+0759+6xtgYGBGj16tNq2bXvJx6tZs6Zq1qxZ4vbYbDb5+vqqf//+6tmzp06fPq2QkJDzJvt/lD9MMF9eXp4yMzNL3B4AAFB1UPcAICgAimHbtm3auXOnkpOTJbmW2alRo4buueeecm1X/lC8Sx3Cl5KSory8POtnT09PBQQElHbzAABAJUTdA4CgACiGTZs26fTp09bPISEhuu666xQREVGq5zHGSFKhNYHd3NyKNcSvuOfYuHGj0tLSrG0hISFq0qRJqRwfAABUbtQ9AAgKgGLYvn27kpKSrJ+rVaumHj16lPp5nE6nEhISNHXq1ALnGjBggKpXr37ZxzfGKD4+XlOnTtXJkycluTrj6OhoXX/99Zd9fAAAUPlR9wAgKACKYc+ePdZkPpLk7++vZs2alfp58vLydOjQIT3++OPWtiZNmuiaa6657A7TGKO9e/fqrrvuKrAWcteuXTV69OhLHsYHAACqJuoeAAQFQDEkJSUpJyfH+tnX11dRUVFlci5jTIEheE6n0xqaVxyZmZlKTU1VZmamsrKylJSUpP3792vp0qWaPXu2zpw5Y12n17t3bz300EMaMmRIqQ3xAwAAlRt1DwCCAuAicnJylJOTU6AT8/DwUFBQUDm26vzWr1+vCRMmaOXKlTLGyOFwKCcnR+np6UpLS5MxRk2bNtWoUaM0aNAgNW3aVN7e3uXdbAAAUAFQ9wCQCAqAi8rIyLCGq0muGXc9PDzKpJNxc3NTtWrVdOONN1rb6tSpo5CQkGIfIzMzU0ePHtXevXvPu09ERIRq166t2rVrV9iOHwAAXHnUPQAkggLgoowxBYbA5f9/WQxZc3d3V61atfTEE09Y2/z8/FSjRo1iH8Pf31916tRR48aNZYxRbm6uMjIydObMGavtcXFxmjp1qo4cOaKBAweqffv2pOsAAIC6B4AkggLgonx8fOTm5lZgm8PhkN1uL/VOxs3NTUFBQYqJiSnxMerVq6ehQ4eqbdu2cjqdys3NVXp6uuLj47V//35t3LhRR48e1ZEjR7Rz504dP35cOTk56tGjh9zd3Uvx0QAAgMqGugeARFAAXJSPj488PT1ls9msZNrhcCgjI6NCptFRUVGFJhwyxigzM1OrVq3S119/rQULFuj06dM6duyYfv75ZzmdTkVHRys6Orp8Gg0AACoE6h4AkuR28V2Aq5u7u7sCAgLk4XE2V7Pb7Tpz5kw5turS2Gw2+fv7q3///vr44481aNAg6xq906dPa/ny5frpp58uaZZhAABQ9VD3AJAICoBiqVmzpvz8/KyfMzIytG/fvnJsUckFBATob3/7mxo0aGBti4uL0+zZs2W328uxZQAAoCKg7gFAUAAUQ/PmzQvMwJuenq4dO3aUX4MuU9OmTdWsWTMFBwdLknJzc3XixIkLzhgMAACuDtQ9AAgKgGJo2bKlQkNDrZ8TExM1f/78cmzR5atTp06BJYJyc3N15MiRcmwRAACoCKh7ABAUAMXQrVs31a5d21oaKDMzUzt27ND69evLuWUl5+npWWBW4/wlhQAAwNWNugcAQQFQDDVq1FCnTp0KzKp7+vRp/etf/5LD4Si3yXCMMcrOzi7R+U+cOKG0tDTrZw8PjwLfHgAAgKsTdU9Bubm52rdvnxISEpgAEVcNggKgGNzd3dWnTx+1atXK2padna3ly5friy++KJc25eXlaf/+/Ro5cqR+/fXXS0rFDxw4oN27dys1NVXS2dmBGzVqVFbNBQAAlQR1z1nGGCUmJuqFF17QiBEj9OWXXyozM5PAAFWex8V3ASBJrVq1UkxMjPbs2aPdu3fLGKOkpCR98MEH8vb21i233CI/P78Cw9rOJz09XatXr1aDBg1Up04deXl5SXKtU5yQkKBp06ZZ+4aGhur6669X9erVrW12u10HDhzQq6++qkWLFik1NVV//etf1adPHwUEBFhDBYuSkZGhzz//XPv27VNeXp4kKSQkRC1atFBYWFhJnx4AAFCFUPe4OJ1OZWRkKD4+XitXrlRSUpK2bt2qW265RW3btlVAQMBFjwFURgQFQDGFhISoX79+Onz4sE6ePKnk5GQ5nU5t27ZNH3/8sY4dO6YuXbqoYcOGqlGjhnx9fa37Op1Opaen6+jRo9q1a5e2bt2q9evXa8iQIRo0aJDq1q0rydVhnjx5Up988ol13/r166tdu3YFOsycnBzt2bNHM2fOVFZWlpYtWyZPT08dPnxYnTp1UsOGDRUaGip3d3dJrjTcbrdr9+7dmj9/vn744QedPn1akitVj46O1g033GB13AAA4OpG3VPwubj77rsVGRmphQsX6ptvvtGRI0fUu3dv9ejRQ40aNSIwQJVDUABcgmbNmmnYsGE6ceKEFixYoKSkJEnS2rVrtXv3bvXu3Vtt2rRRdHS0QkJC5OHhIZvNppycHJ06dUp79uzR2rVr9fvvvysrK0vu7u5q0qSJ1WEaY6wJg/Ll5eUpKyurQDvyh8y1bt1aq1atktPp1KJFi3Tw4EF1795dbdu2VWRkpHx9feXp6am8vDwlJSVpzZo1mjZtmhISEuR0OiVJkZGR6t27t2JiYq7QswgAACoD6h7XZRgRERG688471aFDB9WoUUMrVqzQggULtHnzZu3cuVN9+vRRmzZtFBUVxZcuqDIICoBL4OnpqS5duljD3JYsWaJTp07J6XQqJSVF06dP1/Tp0+Xm5iY/Pz/5+/vLZrMpLS1NGRkZBY5ls9m0f/9+xcfHyxhzwWFzf+Tv76/rrrtOL730kp588knt3btXmZmZ2r9/v/bv369JkybJ19dXQUFB8vf3V1ZWluLj4+VwOKxjuLm5qUaNGho6dKhGjx6tmjVrltrzBAAAKj/qnrPc3d3VsmVLvfHGG5ozZ47Gjx+vLVu2aOLEiVqyZImGDBmim266SY0bN1Z4ePglPT6gIrIZZuIALpnD4VBiYqK++OILjR8/XomJicrOzlZeXt4FJ7ex2Wzy9PSUt7e3AgIC1KdPH917773q2bOn3NzclJubq82bN2vAgAHWfRo1aqTPPvtMbdu2LXS8/CGAzzzzjNatW6f09HTl5uZesA0eHh7y8fFR9erVNWbMGN16661q0qTJZT0fAACg6qLuKSwlJUXfffedJk+erK1btyo3N1d16tTRn//8Z40ePVqhoaHy9vYmMEClRVAAXAan06lTp07piy++0K+//moNrcsf3nZuYm6z2RQQEKC2bdtq4MCBuuGGG9SgQYNC1/QlJydryZIl1rbAwEB17NjxvEv4GGPkcDg0efJkff/991q/fr2SkpLO24Y2bdpo4MCBuv3229WgQQN5e3uXxVMDAACqGOqewk6cOKGffvpJEydO1KZNmyRJTZs21b/+9S8NHDhQnp6estlsBAaodAgKgFKUmpqq48eP69SpU0pJSVFOTo48PDwUHBysOnXqqGHDhmXeWaSnp+vEiRM6duyYTp8+rdTUVAUGBqpGjRqKjo5WZGRksWYoBgAAuBDqnrOOHTumX375Rf/5z38UFxcnSbrhhhv0+uuvq2nTptbcBQQGqCwICoBSdO6f07n//8dOoSw7ieK0gU4KAABcLuqegu3IHyExfvx4/fOf/1R6erq8vLx066236pFHHlHr1q2Z7BCVBkEB8P+lp6crKSnJmokXLkeOHFG1atXk5+dHwAAAQBVB3VP6jDEyxig7O1uHDh3Sf/7zH02ZMkXGGNWsWVNDhgzRqFGj1KlTJ2spR6CiIijAVS81NVWrV6/W999/r7i4OL322mu69tpry7tZFUJqaqruvfdeZWVlFVr7GAAAVD7UPWXPGKO8vDwlJCRozZo1+uSTT7Ru3ToZY9SgQQP169dPd911l1q0aFHeTQXOi+URcdWy2+1avny5ZsyYod9++02HDh1SYGCgDhw4QIf5/x05ckRxcXHau3ev9u7dqxUrVujGG2/UwIEDFRwcXN7NAwAAxUTdc+Xkr/ZQu3Zt9e3bV/Xq1dOyZcs0Y8YMbd26Vd988402bdqkAQMG6M4771RERASjNlHhEBTgqrRjxw7NmTNHS5Ys0caNG+VwONSmTRsNGjRILVu2LO/mVRhhYWF6+OGHtWjRIq1cuVJz587VwYMHtXHjRvXv31+9evWSh4cHnRsAABUYdU/5CQoKUvv27VWzZk01btxYy5cvV2xsrFauXKmjR49q3759uuGGG9S3b1/5+vpSU6HC4NIDXDWMMUpKStKqVas0b948zZo1S6mpqWrYsKG6du2qvn37qmfPnuddjudqlZubq02bNunXX3/VsmXLtG3bNmVkZKhXr14aOnSounXrpiZNmrCSAgAAFQh1T8WTl5enuLg4LV26VIsWLdK6det0+PBh9erVS8OGDbNqKh8fHwIDlDuCAlR5TqdT2dnZOnjwoH777Td9+umn2rp1qyIjI9WuXTsNHTpUMTExXHt/ESkpKdqwYYOmT5+uJUuW6ODBg/Lx8dHdd9+t4cOHq0mTJgoODmZyHgAAyhF1T8WXm5urw4cPa+rUqZo5c6Z27dolT09PjRgxQjfeeKNatWql6tWry8fHp7ybiqsYQQGqLGOMcnNzlZSUpB07duiLL77Q999/r4CAANWvX19jx47V8OHDVatWLb4NvwTJyclatWqVPvjgA23YsEGnT59W+/bt9eijj6pz586qXbu2fH19eU4BALiCqHsqH4fDoS1btuj111/XypUrlZCQoOjoaN1222264YYb+BIG5YqgAFWSw+FQdna2duzYocmTJ2vChAlKT09XYGCg7rzzTj3zzDOqUaMGb7wlZIxRVlaWpk6dqn/+8586fvy4cnNzFRMTowceeMAaysj8BQAAlD3qnsorf0nFH3/8UW+88Yb27t2rzMxMNWvWTLfffrtGjx6tGjVqyNvbm5oKVxRBAaqU/DfbLVu2aMKECfrxxx917Ngx+fn5KSYmRv/+97/VokUL3mhLiTFGycnJevfdd/XJJ5/o9OnT8vHxUe/evXXPPffoxhtvlJeXF883AABlgLqnasnIyNC3336rzz//XJs3b5YxRlFRUbr//vv18MMPy9fXV5L4feKKIChAlXLgwAH973//048//qijR48qKChInTt31mOPPaYbbrhBEm+upS3/LeTAgQN65513NH36dB07dkyBgYHq3Lmznn32WfXt27ecWwkAQNVD3VO15NdU8fHxmj59ur766iutWbNGHh4eaty4sV555RUNHz6cS0dwRRAUoNIzxigjI0Pvv/++3n//fSUmJsrd3V09e/bUmDFjNHjwYPn5+THcrow5nU7l5eVp69atmjRpkqZNm6b4+Hj5+/urT58+ev3119WsWTMKFgAALgN1T+nIy8tTenq6QkJCyrsphRhj5HA4dOzYMc2ZM0fvvvuu9u3bJ09PT/Xt21f/+Mc/1KpVK3l7e5d3U1GFERSg0jLGKDs7WzNnztRrr72mQ4cOKT09Xb169dK9996rbt26KSIigjVprzC73a60tDRt3bpVkydP1jfffCOHw6Hw8HCNGTNGjzzyiMLCwihgAAC4BNQ9pSclJUWLFy/WCy+8oJEjR+rxxx+Xv79/hXveHA6HMjMzdfz4cU2bNk3vvPOOMjMzFRYWpqFDh+qee+5Ry5YtWR0BZYKgAJVSamqqNm3apC+//FIrV67U4cOH1ahRI91///3q0aOHGjRooMDAQD6MlpP8yQ5PnDihLVu26JtvvtH06dMVERGha665RiNHjtTgwYOZeRkAgGKg7ild+/bt0+eff6533nlH1atXV/v27fXQQw+pT58+FW5upfzRBQkJCdq+fbs+//xzzZgxQ/7+/mrYsKEGDBigm266Se3ataOmQqkiKEClkpGRoS1btmj+/PlasmSJtm3bJkkaPny4hg4dqvbt2ysiIkIeHh7l3FJIrssRMjMztX37di1atEjffPON4uLi1LBhQ3Xr1k2DBg1S165dVaNGjfJuKgAAFQ51T9lISUnRli1bNGfOHE2bNk0nTpxQ27Zt1bNnTw0bNkzXXHONAgICyruZBeQvf7l161YtX75cP/30k7Zv367g4GC1adNGvXv31qBBg9S0adPybiqqCIICVAp2u127du3SypUrtWTJEq1bt06JiYlq3769BgwYoIEDB6p169ay2WwVKgWGizFGJ06c0Lx58xQbG6vY2FjZ7Xa1bNlS3bt3V/fu3dW5c2cFBQWVd1MBACh31D1lLzc3V0eOHNGiRYu0cOFCLVq0SH5+furZs6d69uyp7t27q3HjxtZKAxVJYmKili1bpiVLlmjZsmU6cuSIatSooV69eqlv377q27evQkJCeG3gshAUoMI7ceKENmzYoIULF2rBggWKj49XZGSkOnXqpBtvvFE33ngja8tWIps3b9bUqVO1atUq7dy5UzabTZ06ddKwYcPUsWNHNWvWTB4eHvw+AQBXJeqeK8cYI6fTqXXr1mnixIlavXq19uzZo3r16qlfv37q27ev2rRpo3r16lXIyzr27t2rhQsXavHixdq8ebNOnDihVq1aafTo0erWrZtatGhR4S6lQOVBUIAKyRijzMxMxcXFaeHChfrkk090/PhxBQUFqXXr1rrxxhs1YsQIRURElHdTUQI5OTlatmyZfvjhB61YsULHjx+Xh4eHbr75Zo0dO1aNGzdWeHi4JJZ1AgBUfdQ95St/oshp06bpu+++044dO5ScnKz69etr2LBhGjp0qKKiohQWFlbeTS3E4XBo7969mjFjhqZNm6ZDhw4pNTVVN998sx588EE1adJE4eHh8vT0pKbCJSEoQIVijFFeXp5SU1O1bt06vfTSS9q0aZN8fHzUokUL3XnnnRo8eLDq1atX3k1FKUhNTdXy5cs1ceJELV68WOnp6apWrZqefvppjRw5UtWqVZOXlxeT8wAAqiTqnorn1KlTmjJliqZNm6bNmzfL4XCofv36euihhzRy5EgFBwdXyJGPOTk52r9/vz7//HNNnjxZSUlJCgkJ0X333adbb71VUVFR8vX1ZT4LFBtBASoEY4w1U/7evXv17rvv6ttvv5UkVa9eXQ8++KDuvvtu1alThw+NVVBKSopWrFih1157TWvXrlVeXp7at2+vF198UT169FBwcLDc3NwqXKcMAEBJUPdUfPv379e0adP01Vdfaf/+/ZKk7t27680331Tr1q2tb+grWm2Sl5enTZs26emnn9Zvv/2m3NxcNWnSRPfdd5+GDh2q6Ohoubu7V7h2o+IhKEC5y78+bMeOHfryyy/1xRdfKD09XZ6ennrggQf0/PPPKzw8nI6yijPGKCcnR1OmTNGzzz6rkydPSpIGDhyoRx55RD169LBmIKZzAwBUVtQ9lYcxRsePH9eUKVP0+uuv6/Tp07LZbLrtttv00ksvqX79+vL09JRUsWqT/I93P//8s15++WXt3r1bubm5uuaaa3THHXdo7Nix1opTFandqFgIClDuNmzYoPHjx+vnn3/WyZMn5e3trZiYGL355ptq3ry51VHyRla15b8VGWOUnp6u9957T++//77OnDkjb29v9ejRQ/fee69uuukmeXl5lXNrAQAoGeqeyiO/NsnLy9OpU6f03nvv6X//+5/sdrtCQkJ05513aty4cWrRokWFGtJ/bk2VmZmpSZMm6bPPPtO2bdtks9lUv3593X///XrooYfk4+NTzq1FRUVQgHJz4MABvf/++5o1a5aOHDkiHx8fdezYUY8++qhiYmLk5+fHcPOrUP5bUnZ2tvbs2aPPPvtMM2bMUHx8vMLCwtSpUyc9/PDDGjBgQDm3FACA4qPuqbzOvVRk+/bteuWVVxQbGyun06l69epp2LBhGj16tFq3bl3eTS0kf8TmiRMnNGfOHE2aNEnr1q2Tn5+fGjVqpGeeeUYjR45kBAsKISjAFWWMUWpqqj755BNNnjxZhw8fVl5enrp06aJbb71Vffr0UfXq1eXv71/eTUUFkJubqzNnzmjTpk2aNm2a5s+fr+TkZNWqVUt9+vTRk08+qcaNG1NUAQAqJOqeqsUYo9zcXJ06dUrLli3TZ599pi1btshms6lRo0a6/vrrde+99yo6Orq8m1pIXl6e0tPTtW/fPs2fP18TJkzQ4cOHFRERoX79+umvf/2rmjdvLm9v7/JuKioIggJcEfnLzixevFgTJkzQhg0bdPz4cXXs2FHDhg3Tddddp0aNGik8PJwPfSgkLS1NR48e1fr16zV37lzNnj1bXl5eatGihYYNG6bbb79dYWFhFXKNYwDA1Ye6p+o7c+aM9u7dq8WLF2vmzJnatm2bgoOD1apVKw0bNkyjRo2Sn59fhfv95uTkKCEhQTt27NCcOXP01Vdfyd3dXc2aNVP//v11yy23qHHjxgQGIChA2UtPT9euXbs0e/ZsLV26VKtXr1ZkZKQGDx6snj17qmPHjqpZs6Y1GQxQFGOMkpOTtWPHDi1btkxz5szR6tWr1bRpU/Xq1UsDBw5U165dFRYWVuE6ZQDA1YO65+rhdDp1+PBhbdq0ScuXL9fSpUu1Y8cONW/eXDExMRo8eLC6desmLy+vClWbGGNkt9u1Z88e/frrr5o1a5ZWrFih2rVrq0uXLurdu7d69Oihxo0b8yXMVYygAGUmKytLhw4d0rp167R48WItWLBAWVlZ6ty5swYOHKghQ4aobt26TEyHS+JwOJScnKxFixZp1qxZWr58uRITE9W9e3f169dP3bp1U7NmzVStWrXybioA4CpC3XP1ys3NVVxcnGJjYzVv3jytXr1aqampGjBggG666SZ17dpVdevWrXATB+bPu7Bo0SLNnDlTK1eu1PHjx9WgQQPFxMSod+/eat++vSIjI8u7qSgHBAUodQ6HQydOnNC2bdu0aNEizZs3T/v371d0dLS6du2qMWPGqHv37kyagstijNHRo0f12WefacGCBdqzZ4+8vb3Vu3dv3XDDDercubPq1asnX1/f8m4qAKAKo+5BvuzsbO3evVsTJ07UkiVLtGvXLtWqVUu33XabevXqpRYtWigiIqJChkX5y0AuXLhQ27dvV0ZGhpo0aaJbbrlFMTExatKkifz9/SvUyAiULYIClIr8l1FaWpqOHTumX375RVOnTtXevXsVGBiopk2b6o477tBdd91V4YZfofJbvny5PvzwQ61fv14nTpxQtWrVdP311+vuu+9Ws2bNFBYWxkzSAIBSQ92DC3E4HFq6dKlee+01bd++XadPn1bz5s11++23a9CgQYqKilJQUFCFDI927NihqVOnas6cOdq3b5/sdrs6d+6sxx9/XB07dlR4eDg11VWCoACXzRgjh8OhzMxMff/993rnnXd08OBB2Ww2tWzZUrfffrvuuusuhYWFlXdTUYXZ7XbNnz9f48eP1/Lly5WWlqagoCCNGjVKTz31VIEEn84NAFBS1D0ojvwJLSdOnKiPPvpIBw4cUF5enlq2bKk77rhDt99+u0JDQ625KipSbeJwOLR161Z9++23+vbbb3X69GlJ0pgxY/TMM8+oZs2a8vHxkc1mq1DtRukiKECJ5b908vLytGHDBj3++ONau3atnE6nmjVrpnvuuUfDhw+vkEvEoOpKSUnR0qVL9dlnn2nevHkyxqhmzZp6/vnnNWbMGHl7e9OxAQAuGXUPSurkyZMaP368vvvuO+3cuVOenp5q1qyZnnrqKd16663WhIEVrTbJzc3VgQMH9N///lfjx4+X3W5XtWrV9Nhjj2ns2LGqXbs2owuqMIIClFhOTo727t2r1157TVOmTJExRrVq1dKTTz6p2267TbVq1ZJU8d70ULXlv6WlpqZqxYoVevXVV7VmzRpJUrt27fSPf/xDPXv2VFBQUHk2EwBQyVD3oKTya5PDhw/rhx9+0BdffKFdu3bJ3d1dHTt21LvvvqtOnTpVuBUG8tvtcDj0+++/64knntDy5ctljFGTJk304IMPatiwYYqKiirnlqIsEBTgktntdm3fvl0TJ07UxIkTlZqaKi8vLz3wwAN6/PHHVb169Qo5jApXF2OMnE6nMjIyNH36dD3//PM6evSoPD091adPHz344IPq2bOngoODy7upAIAKjLoHpcXpdMrhcOjIkSOaOnWq3n33XSUmJsrb21u33nqr/va3v6l+/foVbrJDY4y1pOKcOXP08ssva8+ePXI4HGrWrJluu+02jR071grLUDUQFOCSbNmyRd99952mT5+uI0eOyOl0qlevXnrxxRfVvHlzBQYGMgQJFUb+21tWVpbi4+P18ccfa8KECcrMzFRISIh69uyp0aNH6/rrr5e3t3c5txYAUNFQ96C05X+RkZ6ermPHjumzzz7Txx9/LA8PD1WvXl0jRozQHXfcoWuuuUYeHh7l3dxCsrKylJCQoB9//FGTJk3Srl275O3traZNm2r06NEaO3asAgICyruZKAUEBSiW/fv3a9KkSVq8eLF2796tvLw8tW7dWnfddZe6deumqKgo69pvoKLJ75Tzl6+aMmWKFi1apPT0dNWrV0+dOnXSqFGjFBMTw2sYAEDdgzJnjFFeXp5OnDihtWvX6vPPP9fy5csVGBio5s2ba+DAgRo8eLBatmxZ3k0txOl0KjExUXv37tX8+fM1Y8YM7d69WzVq1FCbNm30wAMPaMCAAcwJVckRFOC8jDFKT0/X9OnT9csvv2jDhg1KSkpSkyZNNGDAAPXr108tWrRQeHh4eTcVKLasrCzt3r1bq1at0rx587R27Vo5HA61atVK/fv316hRo1SvXj06NgC4ylD3oDzkv+62bNmixYsX6+eff1ZcXJyqVaumdu3aqW/fvhoyZIhq165d3k0tJCcnR8eOHdPvv/+uJUuWaNasWTp16pRat26tPn366Pbbb1fjxo0r3KUUKB6CAhRijFFubq42btyouXPnatGiRfr999/VsGFDde/eXT169FCnTp3UoEGDCrn+K1AcCQkJ2rx5s5YtW6YVK1Zo/fr1ioiI0PXXX6/+/fsrJiZGwcHBvMYBoIqj7kFFcfjwYS1fvlzLli3TypUrdeLECdWrV099+vRRnz591Lt3b/n5+VW4LzMyMzN14MABxcbGauHChYqNjVVQUJD69eunmJgY9ezZU5GRkVzmWckQFKCA7OxsHT9+XOvXr9eMGTP0008/KSgoSO3bt9eAAQM0aNAg1a9f35q0B6jMjDE6deqUVq5cqV9++UXr1q3T3r171aVLF40aNUqdOnVSkyZNFBAQQHEIAFUQdQ8qGofDoV27dmnhwoVavHixtmzZooSEBHXr1k233367unbtqsaNG8vd3b1CBQbGGOXk5Oi3337T5MmTtXr1au3evVvNmzfX4MGD1aNHD7Vs2VI1a9ascKs7oGgEBZDkmtE3KSlJe/bs0aJFi/TVV1/p5MmTioqKUt++ffXAAw+oSZMmJIGokhwOhxITEzV16lRNmjRJ+/btU15enm666SbdcsstatOmjapXry5/f//ybioAoBRQ96Ciy8nJ0Y4dO/TDDz9o1qxZiouLk9Pp1Lhx43TrrbeqYcOGCg0NrXAhljFG2dnZ+vbbbzVx4kTt2bNH6enp6tChg/70pz+pX79+qlevHstUVwIEBVc5p9OpzMxMxcXFafbs2fruu++0bds2BQcHq1WrVnrhhRcUExPDt6m4KjidTh06dEj//Oc/NW/ePJ05c0bVqlXT4MGDNXLkSHXs2FEBAQEVchZiAMDFUfegssnOztbOnTv11ltvaf78+UpJSVHdunX1l7/8RYMGDVK9evXk4+NTIb+lj4uL04cffqjZs2crLi5Onp6e6tq1q5588kn16dOnvJuHiyAouArl/8odDofOnDmjqVOn6quvvtKWLVvk6empBg0a6MEHH9Sf//xneXp6VqhhTcCVsmzZMr322mtat26dkpOTVbt2bQ0ZMkQPPfSQmjRpwprZAFBJUPegKnA4HJo3b56effZZ7d27V7m5uerQoYPuvfdeDR48WDVq1KhwlyPk27x5s7788kvNmDFDgYGBeuqpp3TnnXeWd7NwEQQFV5n8X7fT6dT333+vf/7zn9q3b5/sdruaN2+uO+64Q/fdd5/CwsLKuaVA+bPb7VqwYIE++eQTxcbGKjMzU2FhYRo9erT+/ve/F5j5uiJ2zABwtaPuQVWSv6Ti+PHj9dZbb1mXI7Rp00Zjx47VmDFjFBAQIKni1SUOh0Pbt2/X8ePH1bdv3wp3yQQKIyi4yjgcDi1btkxPPvmktmzZYq0LfO+99+rmm29WZGQka54C/1/+22N6erqWLl2qzz//XLNmzZIkhYeH6+9//7vuv/9+eXl58TcDABUQdQ+qGmOMNRnzxIkT9c0332j79u3y9vZW06ZN9de//lV33HGH3NzcKtTr+o8fOStS21A0goKrRE5Ojvbv3693331XU6ZMUVZWlqpXr66HHnpII0eOVFRUlDw9PbkmDyhCfoJ/5swZrV27Vv/5z3+0YsUK+fj4qGXLlvrb3/6mPn36KDg4uLybCgAQdQ+qPqfTqdzcXB04cEDTp0/X119/rT179sjf318dO3bUv/71L3Xq1KlCzl2AyuGqCgry18ndu3evjh8/rjNnzig9PV3e3t4KDQ1VRESE6tatq1q1al2xlGv37t3Kzc1V9erVFRoaKi8vr1I9vt1u165duzR16lT98MMPOnHihHJzczV27Fjde++9qlevnoKDgxn+AxSD0+lUVlaW4uPjNXfuXP373/9WYmKiwsLC1LNnT40ZM0ZdunQp9cDA4XAoKSlJBw4c0KlTp5SamqrMzEz5+/srPDxckZGRio6Olp+fX6me94/sdrtOnjypnTt3Kj4+Xjk5OfLy8lJQUJDq16+v6OhowhKgAqHuoe5B1We325Wamqp9+/Zp1qxZ+uKLL5ScnKwaNWpo6NChevjhh1W/fv1S/1srS9Q9FcNVERQkJCRo06ZNWrZsmbZv3674+HhlZ2crNzdXeXl5cnNzk5eXl7y9vVWtWjU1atRInTt31vXXX68aNWqUaluMMXI6nVq9erVmzJihNWvWKCcnRz4+PmrdurWGDh2qvn37lsp5du3apblz52r+/PnasWOHTp8+rS5duuiee+5Rp06dFB0dLQ8PD4b+AJfI6XQqISFB27Zt008//aSpU6fKGKPo6Ghdd911uvHGG9WrVy/5+Phc1nm2b9+u5cuXa9OmTTp06JCSkpKUk5Mju90uh8MhDw8PeXl5ydfXV3Xr1lXr1q3Vo0cPtWvXToGBgaX2WOPi4rR8+XItW7ZMe/fuVWpqqrKzs+V0OuXm5iYPDw8FBASoRo0aatOmjfr166euXbuyOgRQTqh7qHtwdckPBU+cOKGtW7dqypQp+uGHHxQaGqprrrlGgwYN0pAhQ9SwYcMKHZJR91QsVTooSElJ0cqVK7VkyRJt2LBB+/btU3x8vHJzc897Hw8PD4WEhKh+/fpq27atRowYoX79+km6vGtpnE6nMjIytHnzZv36669as2aNtmzZohMnTsjpdEqSOnTooIceekjjxo0r8Xkk6dixY1qwYIGWLl2qtWvXWusCDxo0SL169VLXrl0VEBBARwlchvzid/v27VqyZIkWLFigdevWyd3dXa1bt1avXr3Ur18/dejQ4ZKHth4/flzTp0/X8uXLtWPHDh05ckQpKSlyOBznvY+fn59q1aqlpk2bqlu3bho8eLDatGlzWY8xKSnJemzr16/Xvn37lJycfN79vb29VbNmTbVt21bXX3+9Ro0apZCQkMtqA4Dio+6h7sHVzRijnJwcbd68WbGxsZo+fbq2bNmievXqqVOnToqJiVH37t3VsGHDCvX3QN1TMVWt2OMcp0+f1syZM/XTTz9p9erVSkhIKHC7m5ubfH195evrK7vdrqysLNntduXl5SkxMVGJiYnaunWrTp06JWOM+vTpU6IlR3Jzc3Xy5Ent2rVLW7Zs0bp16xQbG6tTp04V2tfhcCgvL69Ej9cYo4yMDK1fv16//vqr5s6dq4MHDyo8PFyDBg1STEyM+vfvr3r16pXo+AAKstlsVigQHR2txo0b69dff9WKFSu0ceNG7dq1S3v27NGgQYPUo0cP1axZ86LvH8YYHTlyRBMmTNC0adO0b98+5eTkFNjHw8NDvr6+8vHxUVZWlrKysuRwOJSZman9+/fr0KFD2rFjh06ePKm77rpLHTt2LNHjS0pK0vTp0zV58mStXbu2QEfp6empwMBA+fn5yW63Kzk5WTk5OcrJyVFcXJzi4+N18OBBeXh4aMyYMSw3BlwB1D3UPYDNZpOPj486deqkJk2aqG7dupo7d65WrVqlmTNnateuXdq1a5d69eqljh07KiIiolzbS91TsVXJoCAjI0OzZ8/We++9p507d1pJuru7u8LCwlSnTh3VrFlT1apVU1BQkHJycpScnKzExESdPHlS8fHxSk5OVnZ2tqZPn6709HQ1aNDAGrJWXEeOHNHu3bu1fv16LV68WKtWrVJGRoZ1u6enp+x2+2U9VmOMHA6Hjhw5ojVr1uibb75RbGysgoOD1b59e/Xt21cDBw5U69atmbAHKCNBQUEaNGiQ2rZtq/bt22v27NnasGGDpk2bpjVr1ujee+9V79691aRJE/n5+RX5t+h0OnXmzBl9++23evvtt5Wenm7NEBwQEKDq1aurZs2aCg8PV0hIiIKDg5WSkqIzZ87o5MmTOnbsmE6fPq2cnBwdOnRIX3/9tTIzM/Xvf/9bYWFhl9Rh5eXlKTY2Vh9++KE2b94su90ud3d3BQUFKTo6WlFRUapRo4ZCQ0Ot8x04cEBxcXFW57ljxw6988476tChg9q3b19qzzWAwqh7qHuAc9lsNoWGhmr06NFq3769pk2bpmXLlmnnzp36/PPP9dtvv+m2225T9+7d1aRJk3JZvYm6pxIwVYjT6TROp9PExsaaqKgo4+bmZiQZScbX19c0b97cPPTQQ2bRokUmIyPDOJ1O6755eXlm7969ZsKECWbEiBGmRo0axmazGUmmWrVq5ttvvzUZGRmX1J7PP//cdOzYsUA7vLy8TLVq1Ux0dLTp0KGD8fT0tG5r27at+fTTT4t9fLvdbk6fPm02btxonn76aePn52cCAgJMkyZNzOOPP242btxocnJyLqnNAC6P0+k0R48eNV988YXp3r27CQ8PN25ubub66683mzdvLvJv0ul0mszMTPPNN98UeE+w2WwmPDzcDBgwwLz99ttm06ZNhe6fnZ1tVqxYYZ555hnTsWNHExgYaCQZDw8P06hRI7Nq1aoC73XFcfjwYdO1a1fj7e1tHatu3bpm1KhRZt68eSYnJ6fAMTMzM80vv/xibr75ZhMUFGS138PDw9x1110mNzf3ktsA4OKoe6h7gOJwOBxm3bp15qmnnjKtWrUyISEhxt/f3wwbNswsX77cnDx50uTl5V2xvpq6p3KockFBamqq6dq1a4EXnY+Pj+nfv79ZuHBhsY5z8uRJ88knn5iwsDATEBBgnnvuOZOVlXXJv/BvvvnGdOnSxXh5eRl/f38THh5urr32WvPss8+aVatWmT179piaNWuWqMO02+3m4MGD5t133zV16tQx7u7uJiAgwIwYMcJs3rzZZGdnX1JbAZSuvLw8Ex8fb1599VVTvXp14+vra2bOnFlk4W23283u3btN3bp1rfcDSSYiIsK89tpr5uDBgxc9n9PpNMuXLzejR482gYGBpnHjxpdUgJ/rnXfeMdWrV7c67YYNG5oXX3zRpKSkXPB+K1euNHfccYf1YSP/w8qRI0eMw+EoUVsAnB91D3UPcClycnLMhg0bzKOPPmrCwsKMp6en8ff3N0899ZQ5fPiwyczMvCL9NXVP5VBlggKn02kcDod5++23C3SWbm5u5tZbbzXr1q27pOOlpaWZmTNnmokTJ5Y4ETp27Jj5xz/+YcaOHWu++eYbExcXV+D2/fv3l7jDXL9+vbnnnnuMJOPp6Wk6dOhgYmNjq8SLEqhKHA6H2b9/v3n77bdNVlZWodudTqc5fvy4uf/++wt0lr6+vmbChAkmMTHxks63bds2891335lly5aVuM0LFy40N954owkJCTGhoaHmueeeM6mpqRe9n8PhMD/88EOhjn/SpEkmNze3xO0BUBh1D3UPUFKZmZlm/fr15pZbbrFGADVo0MC8//775ujRo8bhcJTZN+LUPZVHlQoK7Ha7adKkSYFUJyYmxkyfPv2SX+z5w/ku54/k3GMUdazL6TAzMjLMzz//bDp27Gg+/vhjY7fbq8QQF6CqudB7gDGukQcbN240wcHBBTqZp59+2hw7dqxc37t+//13s2TJErNv375iHc/pdJr169ebO+64o8Bj+ctf/lJkSAKg5Kh7qHuAksr/+8zLyzMLFiwwrVq1sgKDtm3bmo8++sgcP368TM5N3VN5VJnJDB0Oh2JjY7V3715rIgxPT08NHz5cMTExlzxBR2lM6FGWk4L4+vpq8ODBGjBggDw9PeXu7l5m5wJQchd7Hzh+/LjmzJmjlJQUa1tERIQeeOAB1ahRo1zfu1q1aiVjjGw2W7GOa7PZFBISoujo6ALbExISrPdlAKWDuoe6Byip/L9VNzc39enTR6tWrdLkyZP19ttva+fOnfrrX/+qL7/8UqNHj9bdd9+t0NDQUjs3dU/lUWWCgry8PC1ZsqTAL6Vjx45q2rSp/P39y7FlZSN/abaSLF0EoOI4ffq01qxZU2DbTTfdpJCQkHKfsbsk5/f09JSvr2+BbdnZ2aXVJAD/H3UPgMtls9nk5uYmf39/jRo1Sr1799aUKVM0depUbdmyRQcPHtQPP/yg+++/X6NGjZKbm9tl//1R91QeVSYocDgc2rhxY4FtnTp1UmRkZKm/6NLS0nT8+HFrjWI3NzeFhYWpadOmpXqei6GjBCq/1NRU7dy5s8C2fv36ycfHp9T/xhMTE3Xo0CGrA/P19VWNGjVUp06dUjtHXl5eoQ4yNDSU9yuglFH3ACgN+X9X/v7+ql+/vv785z+rZ8+emjdvnqZPn64NGzbo5Zdf1uzZs/XYY4+pc+fOl/UeQ91TeVSJoMD8/zV1t23bVmB7o0aNSnWoTL6jR4/q+++/1+LFiyVJ3t7e6t69u1588cVSPxeAqis3N1dnzpzRsWPHCmxv3bq1PD09S/18mzdv1meffab4+HhJUmRkpG644QbdcccdpXaO/A8U54qMjKwSHSZQUVD3ACgLHh4eqlmzpoKDgxUZGal27dpp4cKFmj59uubOnauTJ09q4MCBGjlypOrUqXPJtQp1T+VSZYKC7OxsnTx5ssD2unXrKigoqNTPl5aWph07dmjZsmWSJB8fH4WFhZX6eQBUbdnZ2UpKSlJWVpa1LTAwULVq1SqT628TEhK0Zs0axcXFSZIaNmyoli1bluo5EhMTtWvXrgLb2rRpU+7DCYGqhLoHQFny9fVVo0aNVKtWLTVs2FCNGzfW4sWLFRsbq6NHj+rQoUPq2bOnrrvuOtWuXVseHsX7SEndU7lUmaAgIyNDdrvd2ubh4aGwsDD5+PiUY8sA4Pyys7OVnp5eYFtoaKh8fX0rZRKdlpamAwcOaPfu3da2wMBAdejQgYnHgFJE3QPgSvD391f79u11zTXXqFWrVgoPD9e6des0YcIE/f7779q/f7+6d++uZs2aFWsiQuqeyqVKBAVOp1NpaWkFtvn4+MjT07NM0hx3d3f5+/srJCTEOldVnDgIQNnKzc1VZmZmgW2BgYFldj4vLy8FBwdb711BQUGl9qHCGKO9e/dq7dq1SkxMlOS6jrlt27aKioqqlAUAUFFR9wC4Umw2m3x9fTVw4EB17NhREydO1A8//KC9e/fqnXfe0YoVKzRmzBj96U9/kpeX1wWPRd1TuVSZoCAjI6PANl9f3zIb8hEaGqqOHTvK4XBIcr0IO3ToUCbnAlB12e32QhPg+Pn5lVnnEhkZqX79+lkdWo0aNdSkSZNSOXZmZqZiY2O1aNEia5uvr6/uv//+UpklGcBZ1D0ArjSbzaaIiAg9+eSTGjBggD788EPFxsZq06ZNqlevngYMGHDRoIC6p3KpEkGBpELDO5xOZ5mdq0GDBnrooYf00EMPldk5AFR9+csSnass193t3LmzOnfuXOrHdTqdmjVrln7++WcdOnRIkmu5oGuuuUYjRowo9fMBoO4BUH5atWqljz/+WCtXrtTevXsVFRVVrIlUqXsqlyoRFLi7uxeavCc9PV0Oh0PGmCqR6ACoery9veXn51dgW2pqapl2mqXNGKNNmzbp008/1apVqyS5CoGaNWvqgw8+uOi3CwAuHXUPgPJms9nUvXt3de/evdj3oe6pXKpEUODm5laow8zJyVFWVpacTmeVmEwCQNXj5eVV6Drf5OTk8mlMCe3fv19PPPGE1qxZY3X09erV0+OPP65rr722nFsHVE3UPQAqI+qeyqXyr9sgV4rj7+9fqNOMi4tTSkpKObUKAC7M399f1apVK7AtISFBiYmJ1rXAFdm+ffs0evRorVmzxrrmMCoqSnfccYfuu+8+vtUEygh1D4DKiLqncqkyQYG7u7tatWpVYPuBAwcqXUoF4Orh5eWlsLAwRUZGWtuMMdqxY4fy8vLKsWUXZozRvn37dNddd2nLli1WZ9mgQQPdddddeuSRR+Tt7V3OrQSqLuoeAJURdU/lUiWCAsl1vV6bNm0KbNu0aZPi4+PLqUUAcGFubm4KDg5W48aNC2xfvny5cnNzy6lVF5aXl6fNmzfr0Ucf1aZNm6zOsnHjxho3bpzuvvtuRUREVLlUHahoqHsAVDbUPZVLlQoKunTpUmDbpk2bdOjQIeXk5JRTqwDgwkJCQgoV+/Pnz1d6enqFm9wnPT1dCxcu1CuvvKLY2Firs2zZsqX+/Oc/a8SIEapfv36ZLdEG4CzqHgCVEXVP5VFlHlV+hxkeHm5tO3HihNauXav9+/eXY8sA4PzCw8PVtWtX+fj4WNu2bdumdevWKT09vRxbdpYxRqdOndLcuXP10UcfacaMGVZn2aFDB40bN04jRoxQo0aNqmxnCVQ01D0AKiPqnsqjyjwyd3d3NWzYUF26dLFm+3U6nfr1118VGxvLNXsAKqSgoCC1a9dOTZs2tbZlZmbq66+/1oEDB8p9KF5ubq4OHTqkmTNn6r///a9mzZplzarerl07Pfjgg7r77rsVFRVVpTtLoKKh7gFQGVH3VB5VYnnEfDabTXfeeadWrFihlJQUGWO0bds2zZgxQ3Xq1FGfPn0UGBhYrGMZY+R0OpWamqqQkBDr+JJkt9uVkZFhJUs2m03e3t7WfgBQXG5ubgoLC9PIkSO1c+dOq4OcMWOG2rdvLx8fH0VHRxd7khyn0ym73a6cnJxCM6JnZ2crLS3NmlnYw8NDvr6+hZYqypeZman9+/dr2rRpmjRpkg4dOiTJtQ5y48aN9eqrr6pfv35VcgIfoDKg7gFQ2VD3VB5VKiiQpOHDh+vbb7/VwoULlZWVJUlatGiRkpKSlJmZqcGDB8vb21uenp6SVGDiCWOMjDHKy8tTTk6OkpKSNHPmTI0bN67A8JgTJ05o8eLF2rhxoyTJ09NTrVu31t13312oPfnX2hR1zU3++f64zel0Fto3v51VcaIM4GoXGhqqcePGadq0afr9998luQrzV199VWfOnNHo0aPVtGlTeXl5yd3dvdD7gDFGDofDKuYPHjyo/fv3a8SIEQXWU9+xY4dmzpyp06dPS5IiIiLUpUsX9e/fv9DxsrOztWLFCr3//vv69ddfrfdTHx8ftW7dWuPHj1fz5s1ls9mKfM/6I97DgLJB3QOgsqHuqRyqVFCQv1zQe++9p8GDB2vPnj3Ky8uT0+nUunXrdN999+m6667Tgw8+qN69e8vf39/65eW/4E6ePKl169Zpzpw5mjlzphITE1W/fn0NGDDAeuHFx8dr7ty5+v777yW5XkA333xzkR1mXl6esrOzrRT+XElJSQU6zLy8PKWnp1sv5nOFhITIw6NK/boA/H/56foHH3yg/v37W51Tdna23n77bc2aNUuDBw/Wn/70J7Vt21ZeXl4F3rtycnK0Y8cOxcbGas6cOVq3bp1q166tNm3aqHnz5tZ59uzZowkTJiguLk6S1LBhQ7m7uxfqMLOzs/Xmm29q/PjxOnLkiLU9ICBAAwcO1JtvvqmAgACdOXOm2I+R9zCg9FH3AKiMqHsqCVMFOZ1Os2PHDtOvXz/j6+trJBX65+3tbZo3b2569Ohh+vXrZzp16mRq165tPD09C+176623muzsbOv4a9asMSNHjrRu9/HxMbfddluRbdmwYYMZN25ckW24lH/Lli0zWVlZV+opBHCFOZ1O43A4zK+//mrCw8ONzWYr9D7g5uZmgoODTatWrUyfPn1Mv379TOvWrY2/v3+h/cPDw83LL79c4ByTJ082UVFR1j4NGzY0//rXvwq1ZciQIed97yzpv6VLl/IeBpQR6h4AlQ11T8VXyWOOotlsNjVt2lRffvmlxo8fr8mTJ2vPnj0F9snJydGuXbus/U0Rw+Hyb7vciSqKOi4AnCs/Ke/Vq5cWLFigv//971q+fHmBGYCdTqdSUlK0bdu2Asn6H99j8r9lLOl7z4EDB2S320v4SABcadQ9ACob6p6Kr0oGBZJrSEvt2rX1xBNP6KabbtLixYs1Y8YMbdq0SWlpaZLOfx2dj4+P6tevr759+2ro0KHq2LGjvLy8rNtr166tYcOGqVGjRpJcE2O0aNHiCj0yAFVVfkfXsmVLTZo0ScuXL9esWbO0fPlyHTp0yJrw53wFfvXq1dW2bVv17dtXffr0KTCjsCRdc801evjhh63Z0MPCwtShQ4dCx6HIByof6h4AlQ11T8VmM1fBM5OTk6PU1FSdOXNGZ86c0e7du3X8+HElJycrMzNT7u7u8vHxUbVq1RQZGak6deooPDxcISEhqlatmnx9fQtMRJGbm1to9l8fH58iZ/9NTk7WwYMHrWtjSqpHjx4KCQkpMEEHgKotLS1NSUlJSkpK0okTJ3TgwAElJiYqNTVVWVlZ8vLykr+/v2rVqqWoqChVr15doaGhCg4OVnBwcKFZebOysgrM/uvu7i4/Pz8FBAQU2O/XX39Venp6qXac3bt3V2hoKO9hwBVA3QOgMqLuqViuiqAgn/n/M+smJycrIyNDOTk5stvtcnNzszrNgIAA+fv7F0jSL0f+jJx5eXmXdRxfX1+5ublV6pkzAZSM0+m0Cv+srCzl5uYqLy9P7u7u8vDwkL+/vwIDA+Xj41MqHVJmZmaxZvS9FLyHAVcedQ+Ayoi6p2K4qoICAAAAAABwYZc3Ww0AAAAAAKhSCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFgICgAAAAAAgIWgAAAAAAAAWAgKAAAAAACAhaAAAAAAAABYCAoAAAAAAICFoAAAAAAAAFj+H27tYL+M2+B0AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "--- Reactant Structure ---\n" + ] + }, + { + "data": { + "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

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

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "--- Product Structure ---\n" + ] + }, + { + "data": { + "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

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

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "rxn_smiles = \"[C:1]([C:3]([C:2]([H:8])([H:9])[H:10])=[O:4])([H:5])([H:6])[H:7]>>[C:1]([C:3](=[C:2]([H:8])[H:9])[O:4][H:10])([H:5])([H:6])[H:7]\"\n", + "r_smiles, p_smiles = rxn_smiles.split(\">>\")\n", + "params = Chem.SmilesParserParams()\n", + "params.removeHs = False\n", + "r_mol = Chem.MolFromSmiles(r_smiles, params=params)\n", + "p_mol = Chem.MolFromSmiles(p_smiles, params=params)\n", + "\n", + "r_xyz = \"\"\"\\\n", + "10\n", + "\n", + "C -1.460860000000000 -0.144010000000000 -0.053220000000000\n", + "C 1.050380000000000 -0.695950000000000 0.102060000000000\n", + "C -0.033170000000000 0.331050000000000 -0.087970000000000\n", + "O 0.238150000000000 1.516350000000000 -0.275750000000000\n", + "H -2.013860000000000 0.310400000000000 -0.879410000000000\n", + "H -1.910310000000000 0.138980000000000 0.901700000000000\n", + "H -1.508990000000000 -1.230020000000000 -0.166960000000000\n", + "H 1.412030000000000 -1.023800000000000 -0.875360000000000\n", + "H 0.669670000000000 -1.555760000000000 0.659270000000000\n", + "H 1.870170000000000 -0.255950000000000 0.676240000000000\n", + "\"\"\"\n", + "\n", + "p_xyz = \"\"\"\\\n", + "10\n", + "\n", + "C -1.447246410000000 -0.119824230000000 0.000001190000000\n", + "C 0.825515320000000 -1.164474080000000 0.000015810000000\n", + "C 0.048387710000000 -0.090653610000000 -0.000104420000000\n", + "O 0.534209710000000 1.184537660000000 0.000015420000000\n", + "H -1.822704290000000 0.402724090000000 -0.877993030000000\n", + "H -1.822288580000000 0.399353360000000 0.880185840000000\n", + "H -1.811964100000000 -1.141503600000000 -0.001665260000000\n", + "H 1.900576660000000 -1.090080000000000 0.000093490000000\n", + "H 0.410348040000000 -2.153410380000000 0.000030700000000\n", + "H 1.498440620000000 1.164676240000000 0.000145300000000\n", + "\"\"\"\n", + "\n", + "def visualize_molecules_2d_side_by_side(mol1, mol2, title1=\"Molecule 1\", title2=\"Molecule 2\"):\n", + " \"\"\"Display 2D molecular structures side by side.\"\"\"\n", + " fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))\n", + " \n", + " # Remove hydrogens for cleaner visualization\n", + " mol1_no_h = Chem.RemoveHs(mol1)\n", + " mol2_no_h = Chem.RemoveHs(mol2)\n", + " # Create drawing options\n", + " opts = Draw.rdMolDraw2D.MolDrawOptions()\n", + " opts.addAtomIndices = False # Hide atom numbering\n", + " opts.addStereoAnnotation = False # Hide stereo annotations\n", + " \n", + " # Generate 2D images with custom options\n", + " img1 = Draw.MolToImage(mol1_no_h, size=(400, 400), options=opts)\n", + " img2 = Draw.MolToImage(mol2_no_h, size=(400, 400), options=opts)\n", + " \n", + " # Display images\n", + " ax1.imshow(img1)\n", + " ax1.set_title(title1, fontsize=14, pad=10)\n", + " ax1.axis('off')\n", + " \n", + " ax2.imshow(img2)\n", + " ax2.set_title(title2, fontsize=14, pad=10)\n", + " ax2.axis('off')\n", + " \n", + " plt.tight_layout()\n", + " plt.show()\n", + "\n", + "def visualize_molecule(xyz_string: str, name: str=None):\n", + " \"\"\"\n", + " Visualize a molecule from its XYZ string using py3Dmol.\n", + " \n", + " Args:\n", + " xyz_string: XYZ format string of the molecule\n", + " name: Name for the molecule (used as visualization title)\n", + " \"\"\"\n", + " if name:\n", + " print(f\"--- {name.title()} Structure ---\")\n", + " view = py3Dmol.view(width=400, height=400)\n", + " view.addModel(xyz_string, 'xyz')\n", + " view.setStyle({}, {'stick': {'radius': 0.15}, 'sphere': {'scale': 0.25}})\n", + " view.setBackgroundColor('white')\n", + " view.zoomTo()\n", + " view.show()\n", + "\n", + "\n", + "def save_molecule(xyz_string: str, name: str, working_dir: Path):\n", + " \"\"\"\n", + " Save the XYZ string to a file in the specified working directory.\n", + " \n", + " Args:\n", + " xyz_string: XYZ format string of the molecule\n", + " name: Name for the molecule (used for filename)\n", + " working_dir: Directory to save the XYZ file\n", + " \"\"\"\n", + " xyz_file = working_dir / f\"{name}.xyz\"\n", + " with open(xyz_file, \"w\") as f:\n", + " f.write(xyz_string)\n", + "\n", + "visualize_molecules_2d_side_by_side(r_mol, p_mol, \"Reactant Molecule\", \"Product Molecule\")\n", + "visualize_molecule(r_xyz, \"reactant\")\n", + "visualize_molecule(p_xyz, \"product\")\n", + "save_molecule(r_xyz, \"reactant\", base_dir)\n", + "save_molecule(p_xyz, \"product\", base_dir)" + ] + }, + { + "cell_type": "markdown", + "id": "7f287b1f", + "metadata": {}, + "source": [ + "### Step 4: Setup and Run Calculations\n", + "We will run three kinds of calculations: \n", + "1. First, we optimize the geometries of our initial reactant and product structures.\n", + "2. Then we start a nudge elastic band (NEB) TS search starting from the optimized reactant and product geometries.\n", + "3. Finally, we calculate the intrinsic reaction coordinate (IRC) from the found TS. This is not strictly necessary to calculate the barrier height of the reaction but it allows for a nice visualization of the reaction.\n", + "\n", + "For all calculations, we additionally perform vibrational frequency calculations to obtain thermochemistry information including the single point energies needed to determine the barrier height of the reaction.\n", + "\n", + "RGD1 was obtained on a B3LYP-D3/TZVP level of theory and thight SCF. If you would like to reproduce these results, use\n", + "```python\n", + "SIMPLE_KEYWORDS = [\n", + " Dft.B3LYP, # > B3LYP-D3 method\n", + " DispersionCorrection.D3, # > D3\n", + " BasisSet.TZVP, # > TZVP basis set\n", + " Task.FREQ, # > Vibration frequency calculation to get thermochemistry\n", + " Scf.TIGHTSCF, # > tight SCF\n", + "]\n", + "```\n", + "(the commented-out part below). The NEB search will take a few hours then. If you would like to make a dummy calculation with a much lower level of theory, you can use HF/STO-3G as detailed below. The NEB search will then complete in a few minutes." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "ed63e9c6", + "metadata": {}, + "outputs": [], + "source": [ + "# > Define simple keywords for the ORCA calculation, here with a lower level of theory\n", + "SIMPLE_KEYWORDS = [\n", + " Wft.HF, # > Hartree-Fock method\n", + " BasisSet.STO_3G, # > STO3G basis set\n", + " Task.FREQ, # > Vibration frequency calculation to get thermochemistry\n", + "]\n", + "\n", + "'''# > Define simple keywords for the ORCA calculation, consistent with RGD1. Uncomment this part if you like.\n", + "SIMPLE_KEYWORDS = [\n", + " Dft.B3LYP, # > B3LYP-D3 method\n", + " DispersionCorrection.D3, # > D3\n", + " BasisSet.TZVP, # > TZVP basis set\n", + " Task.FREQ, # > Vibration frequency calculation to get thermochemistry\n", + " Scf.TIGHTSCF, # > tight SCF\n", + "]\n", + "'''\n", + "\n", + "# > Set number of cores for the calculation\n", + "NCORES = 2" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "a6453f52", + "metadata": {}, + "outputs": [], + "source": [ + "# > Define setup functions for all three calculator types\n", + "def setup_opt_calc(basename : str, xyz_file: str, base_dir: Path, ncores: int, sk_list: list[SimpleKeyword]) -> Calculator:\n", + " # > Setup subdirectory for geometry optimization calculations\n", + " working_dir = base_dir / \"geometry_optimization\"\n", + " working_dir.mkdir(parents=True, exist_ok=True) # Ensure the subdirectory exists\n", + "\n", + " # > Specify geometry optimization as task\n", + " sk_list = sk_list.copy()\n", + " sk_list.append(Task.OPT)\n", + "\n", + " # > Set up a Calculator object\n", + " calc = Calculator(basename=basename, working_dir=working_dir)\n", + " # > Assign structure to calculator\n", + " calc.structure = Structure.from_xyz(base_dir / xyz_file)\n", + " # > Add simple keywords to calculator\n", + " calc.input.add_simple_keywords(*sk_list)\n", + " # > Set the number of CPUs for the calculation\n", + " calc.input.ncores = ncores\n", + "\n", + " return calc\n", + "\n", + "def setup_neb_calc(basename: str, r_xyz_file: str, p_xyz_file: str, base_dir: Path, ncores: int, sk_list: list[SimpleKeyword]) -> Calculator:\n", + " working_dir = base_dir / \"nudged_elastic_band\"\n", + " working_dir.mkdir(parents=True, exist_ok=True)\n", + " # > Copy the product XYZ file to the calculation directory so ORCA can find it\n", + " p_xyz_path = base_dir / p_xyz_file\n", + " dest_path = working_dir / p_xyz_path.name\n", + " shutil.copy2(p_xyz_path, dest_path)\n", + "\n", + " sk_list = sk_list.copy()\n", + " sk_list.append(Neb.NEB_TS)\n", + "\n", + " calc = Calculator(basename=basename, working_dir=working_dir)\n", + " calc.structure = Structure.from_xyz(base_dir / r_xyz_file)\n", + " # > Define input block for product structure\n", + " options= f\"\"\"\\\n", + " %neb\n", + " neb_end_xyzfile \"{p_xyz_path.name}\"\n", + " end\n", + " \"\"\"\n", + " # > Add options to the calculator\n", + " calc.input.add_arbitrary_string(options, pos=ArbitraryStringPos.BEFORE_COORDS)\n", + " calc.input.add_simple_keywords(*sk_list)\n", + " calc.input.ncores = ncores\n", + "\n", + " return calc\n", + "\n", + "def setup_irc_calc(basename: str, xyz_file: str, base_dir: Path, ncores: int, sk_list: list[SimpleKeyword]) -> Calculator:\n", + " working_dir = base_dir / \"intrinsic_reaction_coordinate\"\n", + " working_dir.mkdir(parents=True, exist_ok=True) # Ensure the subdirectory exists\n", + "\n", + " sk_list = sk_list.copy()\n", + " sk_list.append(Neb.IRC)\n", + "\n", + " calc = Calculator(basename=basename, working_dir=working_dir)\n", + " calc.structure = Structure.from_xyz(base_dir / xyz_file)\n", + " options= \"\"\"\n", + " %scf\n", + " maxiter 200\n", + " end\"\"\"\n", + " calc.input.add_simple_keywords(*sk_list)\n", + " calc.input.add_arbitrary_string(options, pos=ArbitraryStringPos.BEFORE_COORDS)\n", + " calc.input.ncores = ncores\n", + "\n", + " return calc" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "ea1aac5a", + "metadata": {}, + "outputs": [], + "source": [ + "def run_calc(calc: Calculator) -> Output:\n", + " # > Write the ORCA input file\n", + " calc.write_input()\n", + " # > Run the ORCA calculation\n", + " print(\"Running ORCA calculation ...\", end=\"\")\n", + " calc.run()\n", + " print(\" Done\")\n", + " # > Get the output object\n", + " output = calc.get_output()\n", + " \n", + " return output\n", + "\n", + "def check_and_parse_output(output: Output):\n", + " # > Check for proper termination of ORCA\n", + " status = output.terminated_normally()\n", + " if not status:\n", + " # > ORCA did not terminate normally\n", + " raise RuntimeError(f\"ORCA did not terminate normally. Please check {output.basename}.out file\")\n", + " else:\n", + " # > ORCA did terminate normally so we can parse the output\n", + " output.parse()\n", + " # > Now check for convergence of the SCF\n", + " if output.results_properties.geometries[-1].single_point_data.converged:\n", + " print(f\"{output.basename}: SCF CONVERGED\")\n", + " else:\n", + " raise RuntimeError(\"SCF DID NOT CONVERGE\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "0afc66d5", + "metadata": {}, + "source": [ + "First, run the geometry optimization:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20f4eebd", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running ORCA calculation ... Done\n", + "reactant: SCF CONVERGED\n" + ] + } + ], + "source": [ + "r_opt_calc = setup_opt_calc(\"reactant\", \"reactant.xyz\", base_dir, ncores=NCORES, sk_list=SIMPLE_KEYWORDS)\n", + "r_opt_out = run_calc(r_opt_calc)\n", + "check_and_parse_output(r_opt_out)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "438bc289", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running ORCA calculation ... Done\n", + "product: SCF CONVERGED\n" + ] + } + ], + "source": [ + "p_opt_calc = setup_opt_calc(\"product\", \"product.xyz\", base_dir, ncores=NCORES, sk_list=SIMPLE_KEYWORDS)\n", + "p_opt_out = run_calc(p_opt_calc)\n", + "check_and_parse_output(p_opt_out)" + ] + }, + { + "cell_type": "markdown", + "id": "2d436418", + "metadata": {}, + "source": [ + "Then run nudege elastic band search on the optimized reactant and product geometries (this will take a few minutes with HF/STO-3G and several hours for B3LYP-D3/TZVP):" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0af76fd4", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running ORCA calculation ... Done\n", + "ts: SCF CONVERGED\n" + ] + } + ], + "source": [ + "neb_ts_calc = setup_neb_calc(\"ts\", \"geometry_optimization/reactant.xyz\", \"geometry_optimization/product.xyz\", base_dir, ncores=NCORES, sk_list=SIMPLE_KEYWORDS)\n", + "neb_ts_out = run_calc(neb_ts_calc)\n", + "check_and_parse_output(neb_ts_out)" + ] + }, + { + "cell_type": "markdown", + "id": "e0cfe2e4", + "metadata": {}, + "source": [ + "Finally, run the IRC calculation for visualization purposes:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8b01c797", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running ORCA calculation ... Done\n", + "irc: SCF CONVERGED\n" + ] + } + ], + "source": [ + "irc_calc = setup_irc_calc(\"irc\", \"nudged_elastic_band/ts.xyz\", base_dir, ncores=NCORES, sk_list=SIMPLE_KEYWORDS)\n", + "irc_out = run_calc(irc_calc)\n", + "check_and_parse_output(irc_out)" + ] + }, + { + "cell_type": "markdown", + "id": "0039724f", + "metadata": {}, + "source": [ + "### Step 5: Visualize TS and IRC" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "b438883d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "--- Transition State Structure ---\n" + ] + }, + { + "data": { + "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

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

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "ts_xyz_file = base_dir / \"nudged_elastic_band\" / \"ts.xyz\"\n", + "with open(ts_xyz_file, \"r\") as f:\n", + " ts_xyz = f.read()\n", + "visualize_molecule(ts_xyz, \"transition state\")" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "67cad0b4", + "metadata": {}, + "outputs": [ + { + "data": { + "application/3dmoljs_load.v0": "
\n

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

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

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

\n", + "
\n", + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "def animate_trajectory(trajectory_file: Path):\n", + " \"\"\"\n", + " Animate a trajectory from a multi-frame XYZ file using py3Dmol.\n", + " \n", + " Args:\n", + " trajectory_file: Path to the XYZ trajectory file\n", + " name: Name for the animation\n", + " \"\"\"\n", + " with open(trajectory_file, \"r\") as f:\n", + " trajectory_content = f.read()\n", + " view = py3Dmol.view(width=600, height=400)\n", + " view.addModelsAsFrames(trajectory_content, 'xyz')\n", + " view.setStyle({}, {'stick': {'radius': 0.15}, 'sphere': {'scale': 0.25}})\n", + " view.setBackgroundColor('white')\n", + " view.animate({'loop': 'forward'})\n", + " view.zoomTo()\n", + " view.show()\n", + "\n", + "trj_file = base_dir / \"intrinsic_reaction_coordinate\" / \"irc_IRC_Full_trj.xyz\"\n", + "animate_trajectory(trj_file)" + ] + }, + { + "cell_type": "markdown", + "id": "9fa42745", + "metadata": {}, + "source": [ + "### Step 6: Calculate the Barrier Height\n", + "We calculate the barrier height as the energy difference between reactant and transition state.\n", + "We will use electron energy since it will also be used for training in Part II.\n", + "\n", + "The energies can be found in the thermochemistry information of the ORCA output which was obtained via the addtional vibrational frequency calculations." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "cc1b02c1", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "ELECTRON ENERGY (E):\n", + "Activation: 106.655528 kcal/mol, 0.169967 Eh\n", + "Reaction: 19.588445 kcal/mol, 0.031216 Eh\n", + "\n" + ] + } + ], + "source": [ + "E_r = r_opt_out.results_properties.geometries[-1].thermochemistry_energies[0].elenergy\n", + "E_p = p_opt_out.results_properties.geometries[-1].thermochemistry_energies[0].elenergy\n", + "E_ts = neb_ts_out.results_properties.geometries[-1].thermochemistry_energies[0].elenergy\n", + "delta_E_act = E_ts - E_r\n", + "delta_E_rp = E_p - E_r\n", + "\n", + "# > Convert to kcal/mol\n", + "HARTREE_TO_KCAL_MOL = 627.509\n", + "delta_E_act_kcal = delta_E_act * HARTREE_TO_KCAL_MOL\n", + "delta_E_rp_kcal = delta_E_rp * HARTREE_TO_KCAL_MOL\n", + "\n", + "print(f\"\"\"\n", + "ELECTRON ENERGY (E):\n", + "Activation: {delta_E_act_kcal:.6f} kcal/mol, {delta_E_act:.6f} Eh\n", + "Reaction: {delta_E_rp_kcal:.6f} kcal/mol, {delta_E_rp:.6f} Eh\n", + "\"\"\")" + ] + }, + { + "cell_type": "markdown", + "id": "dcdc8ea9", + "metadata": {}, + "source": [ + "Now, the HF/STO-3G values are of course quite off. With B3LYP-D3/TZVP, we obtain `67.935952` kcal/mol. The RGD1 dataset reports an activation energy of `68.35608500000312` kcal/mol (electron energy) for this reaction.\n", + "Given that the RGD1 dataset was constructed using a slightly different approach a small deviation of approximately 0.3 kcal/mol is tolerable.\n", + "\n", + "Now, we will save the reaction SMILES to a CSV file for later use. We will also overwrite `delta_E_act_kcal` to the value at the correct level of theory, since we will compare machine learning model predictions to this value later." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "39e7948a-a8ac-44e1-b767-001d7331e72a", + "metadata": {}, + "outputs": [], + "source": [ + "#Overwrite target to correct level of theory (this is unnecessary if you used to correct level of theory above)\n", + "delta_E_act_kcal = 67.935952" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "ed15a0f8", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
smilesea
0[C:1]([C:3]([C:2]([H:8])([H:9])[H:10])=[O:4])(...67.935952
\n", + "
" + ], + "text/plain": [ + " smiles ea\n", + "0 [C:1]([C:3]([C:2]([H:8])([H:9])[H:10])=[O:4])(... 67.935952" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# > Save the smiles of the reaction to a CSV file \n", + "import pandas as pd\n", + "data = {\n", + " \"smiles\": [rxn_smiles],\n", + " \"ea\": [delta_E_act_kcal]\n", + "}\n", + "\n", + "df = pd.DataFrame(data)\n", + "df.to_csv(\"QM_data.csv\", index=False)\n", + "df" + ] + }, + { + "cell_type": "markdown", + "id": "4af7be35", + "metadata": {}, + "source": [ + "To build a full machine learning dataset we would repeat this for many more reactions.\n", + "Of course, dataset construction also involves other considerations like chemical space coverage and dataset curration.\n", + "However, these topics are beyond the scope of this tutorial." + ] + }, + { + "cell_type": "markdown", + "id": "e586b9c6", + "metadata": {}, + "source": [ + "## Part II: Training a GNN to Predict Activation Free Energies via ChemTorch\n", + "The second part consists of four steps:\n", + "1. Download a small subset of the RGD1 dataset.\n", + "2. Install ChemTorch\n", + "3. Train and Evaluate GNN via ChemTorch CLI.\n", + "4. Use Model to Predict Activation Free Energy\n", + "\n", + "The remainder of this tutorial will run on the command line, please navigate to the folder where you are currently running this tutorial (and thus `QM_data.csv` and `QM_data_precomputed.csv` are already located), and run the commands below on the command line in this folder)." + ] + }, + { + "cell_type": "markdown", + "id": "f3ce895a", + "metadata": {}, + "source": [ + "### Step 1: Download the RGD1 Dataset (a subset)\n", + "\n", + "Since for a tutorial you cannot run hundreds of calculation, we will now donwload a subset of RGD1 to use in our machine learning model instead of the `QM_data.csv` above which we just made. It is exactly the same information, just about 15k reactions instead of the only 1 from above.\n", + "\n", + "```bash\n", + "wget https://raw.githubusercontent.com/heid-lab/reaction_database/main/data/rgd1/barriers/filtered.csv -O QM_data_precomputed.csv\n", + "```\n", + "\n", + "You can inspect the `QM_data_precomputed.csv` to see it is just more reactions than above, but the exact same format." + ] + }, + { + "cell_type": "markdown", + "id": "f853b56a", + "metadata": {}, + "source": [ + "### Step 2: Install ChemTorch\n", + "To begin our ML workflow, clone and install ChemTorch following the [official installation instructions](https://github.com/heid-lab/chemtorch) or the command below (for GPU installation, please refer to the installation instructions). We recommend to set up a separate environment than for the first part of this tutorial. \n", + "\n", + "```bash\n", + "git clone -b tutorial/opi_orca https://github.com/heid-lab/chemtorch.git && \\\n", + "cd chemtorch && \\\n", + "conda deactivate && \\\n", + "conda create -n chemtorch python=3.10 && \\\n", + "conda activate chemtorch && \\\n", + "pip install rdkit numpy==1.26.4 scikit-learn pandas && \\\n", + "pip install torch && \\\n", + "pip install hydra-core && \\\n", + "pip install torch_geometric && \\\n", + "pip install torch_scatter torch_sparse torch_cluster torch_spline_conv -f https://data.pyg.org/whl/torch-2.5.0+cpu.html && \\\n", + "pip install wandb && \\\n", + "pip install ipykernel && \\\n", + "pip install -e .\n", + "```" + ] + }, + { + "cell_type": "markdown", + "id": "e1a57c42", + "metadata": {}, + "source": [ + "### Step 3: Train and Evaluate GNN via ChemTorch CLI\n", + "In this section, we train a Directed Message Passing Neural Network (DMPNN)—a widely used GNN architecture for molecular property prediction.\n", + "ChemTorch provides ready-to-use configuration files, including a complete pipeline for training DMPNNs on the [RGD1 dataset](https://www.nature.com/articles/s41597-023-02043-z).\n", + "> 🔗 For more models, representations, and datasets, visit the [ChemTorch GitHub repository](https://github.com/heid-lab/chemtorch).\n", + "\n", + "> 📚 To define your own GNNs or molecular representations, consult the [ChemTorch documentation](https://github.com/heid-lab/chemtorch)." + ] + }, + { + "cell_type": "markdown", + "id": "fec02883", + "metadata": {}, + "source": [ + "To use our own custom dataset `QM_data_precomputed.csv`, we can either set up our own configuration file, or use an existing one and overwrite the path to the dataset.\n", + "\n", + "To use an existing one, run the following command from the `chemtorch` project root:\n", + "\n", + "```bash\n", + "python chemtorch_cli.py +experiment=graph data_pipeline=rgd1 data_pipeline.data_source.data_path=\"../QM_data_precomputed.csv\"\n", + "```\n", + "\n", + "This tells ChemTorch to use the default graph learning configuration with the RGD1 data pipeline but use our own custom dataset specified via `data_path`.\n", + "Under the hood, this setup will convert each reaction SMILES to a condensed graph of reaction (CGR), train a DMPNN, track metrics of interest and save the best performing model parameters for later.\n", + "to the CLI as well as [Weights & Biases](https://wandb.ai/site/models/) which is a graphical user interface that can be accessed through the browser.\n", + "\n", + "If you would like to make your own configuration file instead, an example is already included in your ChemTorch installation, and can be found in `conf/experiment/opi_tutorial/training.yaml`, so no need to create or change a file. Note that the important lines are setting the `data_pipeline` to `rgd1`, and `data_pipeline/data_source/data_path` to `\"../QM_data_precomputed.csv\"`, just as above. To launch the training process with the config file, run \n", + "```bash \n", + "python chemtorch_cli.py +experiment=opi_tutorial/training\n", + "```\n", + "from the `chemtorch` project root (do not redo this if you already ran ChemTorch above)." + ] + }, + { + "cell_type": "markdown", + "id": "41ec164f", + "metadata": {}, + "source": [ + "Either way, the ChemTorch CLI will output\n", + "1. a link to a detailed real-time view of the ongoing run on [Weights & Biases](https://wandb.ai/site/models/) which can be accessed via the browser (can be disabled by setting `log=false`), \n", + "2. a progress bar, and \n", + "3. a summary of the run including the final training, validation and test metrics." + ] + }, + { + "cell_type": "markdown", + "id": "6b96abe2", + "metadata": {}, + "source": [ + "```txt\n", + " ██████╗██╗ ██╗███████╗███╗ ███╗████████╗ ██████╗ ██████╗ ██████╗██╗ ██╗\n", + "██╔════╝██║ ██║██╔════╝████╗ ████║╚══██╔══╝██╔═══██╗██╔══██╗██╔════╝██║ ██║\n", + "██║ ███████║█████╗ ██╔████╔██║ ██║ ██║ ██║██████╔╝██║ ███████║\n", + "██║ ██╔══██║██╔══╝ ██║╚██╔╝██║ ██║ ██║ ██║██╔══██╗██║ ██╔══██║\n", + "╚██████╗██║ ██║███████╗██║ ╚═╝ ██║ ██║ ╚██████╔╝██║ ██║╚██████╗██║ ██║\n", + " ╚═════╝╚═╝ ╚═╝╚══════╝╚═╝ ╚═╝ ╚═╝ ╚═════╝ ╚═╝ ╚═╝ ╚═════╝╚═╝ ╚═╝\n", + "\n", + "...\n", + "Sanity Checking: | | 0/? [00:00 Load predictions from the CSF file\n", + "df = pd.read_csv(\"../predictions.csv\")\n", + "delta_E_act_kcal_predicted = df[\"prediction\"].values[0]" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "id": "27e875b0", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "PREDICTED ACTIVATION ENERGY: 69.493736 kcal/mol\n", + "ERROR TO ORCA CALCULATION: 1.557785 kcal/mol, 2.29%\n", + "ERROR TO REPORTED VALUE IN RGD1: 1.137651 kcal/mol, 1.66%\n", + "\n" + ] + } + ], + "source": [ + "rgd1_delta_E_act_kcal = 68.35608500000312 # > RGD1 value for the activation energy in kcal/mol\n", + "# > Calculate absolute errors\n", + "error_to_orca_calc = delta_E_act_kcal_predicted - delta_E_act_kcal\n", + "error_to_rgd1 = delta_E_act_kcal_predicted - rgd1_delta_E_act_kcal\n", + "\n", + "# > Calculate relative errors\n", + "relative_error_to_orca_calc = error_to_orca_calc / delta_E_act_kcal\n", + "relative_error_to_rgd1 = error_to_rgd1 / rgd1_delta_E_act_kcal\n", + "\n", + "\n", + "print(f\"\"\"\n", + "PREDICTED ACTIVATION ENERGY: {delta_E_act_kcal_predicted:.6f} kcal/mol\n", + "ERROR TO ORCA CALCULATION: {error_to_orca_calc:.6f} kcal/mol, {relative_error_to_orca_calc:.2%}\n", + "ERROR TO REPORTED VALUE IN RGD1: {error_to_rgd1:.6f} kcal/mol, {relative_error_to_rgd1:.2%}\n", + "\"\"\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "opi", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/index.md b/docs/index.md index 76f75c25..ef0896eb 100644 --- a/docs/index.md +++ b/docs/index.md @@ -72,6 +72,7 @@ contents/notebooks/ir_spectrum.ipynb contents/notebooks/extopt.ipynb contents/notebooks/opencosmors.ipynb contents/notebooks/moplot.ipynb +contents/notebooks/deep_learning_barrier_heights.ipynb ``` ```{toctree} From 75ccdf1edc0069a93f6dc4969e9b8bc00a7e7406 Mon Sep 17 00:00:00 2001 From: Anton Zamyatin Date: Mon, 25 Aug 2025 14:57:53 +0200 Subject: [PATCH 2/6] Update deep learning barrier heights notebook - Fix ckpt_path and prediction_save_path in inference command - Replace arbitrary strings with block inputs instead --- .../deep_learning_barrier_heights.ipynb | 34 +++++++------------ 1 file changed, 13 insertions(+), 21 deletions(-) diff --git a/docs/contents/notebooks/deep_learning_barrier_heights.ipynb b/docs/contents/notebooks/deep_learning_barrier_heights.ipynb index d6ff5626..9c428e4b 100644 --- a/docs/contents/notebooks/deep_learning_barrier_heights.ipynb +++ b/docs/contents/notebooks/deep_learning_barrier_heights.ipynb @@ -5,7 +5,7 @@ "id": "5f7bd60d", "metadata": {}, "source": [ - "# Deep Learning Of Reaction Barrier Heights\n", + "# Deep Learning of Reaction Barrier Heights with ChemTorch\n", "This two-part tutorial showcases how ORCA integrates into downstream deep learning workflows by serving as a data source training and evaluation.\n", "1. First, we will show how to calculate the barrier height of a chemical reaction using ORCA with the ORCA Python inferface (OPI). \n", "2. Second, we use the ChemTorch framework to train and evaluate a graph neural network (GNN) on a curated subset of the popular [RGD1 dataset](https://www.nature.com/articles/s41597-023-02043-z) which contains precomputed barrier heights.\n", @@ -48,7 +48,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "5edf5999", "metadata": {}, "outputs": [], @@ -61,6 +61,7 @@ "from opi.output.core import Output\n", "from opi.input.structures.structure import Structure\n", "from opi.input.arbitrary_string import ArbitraryStringPos\n", + "from opi.input.blocks import BlockNeb, BlockScf\n", "from opi.input.simple_keywords import SimpleKeyword, BasisSet, Scf, Dft, Task, Neb, DispersionCorrection, Wft\n", "\n", "# > for visualization of molecules\n", @@ -382,22 +383,22 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "id": "ed63e9c6", "metadata": {}, "outputs": [], "source": [ "# > Define simple keywords for the ORCA calculation, here with a lower level of theory\n", "SIMPLE_KEYWORDS = [\n", - " Wft.HF, # > Hartree-Fock method\n", - " BasisSet.STO_3G, # > STO3G basis set\n", + " Wft.HF, # > Hartree-Fock method\n", + " BasisSet.STO_3G, # > STO3G basis set\n", " Task.FREQ, # > Vibration frequency calculation to get thermochemistry\n", "]\n", "\n", "'''# > Define simple keywords for the ORCA calculation, consistent with RGD1. Uncomment this part if you like.\n", "SIMPLE_KEYWORDS = [\n", - " Dft.B3LYP, # > B3LYP-D3 method\n", - " DispersionCorrection.D3, # > D3\n", + " Dft.B3LYP, # > B3LYP method\n", + " DispersionCorrection.D3, # > D3 correction\n", " BasisSet.TZVP, # > TZVP basis set\n", " Task.FREQ, # > Vibration frequency calculation to get thermochemistry\n", " Scf.TIGHTSCF, # > tight SCF\n", @@ -410,7 +411,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "id": "a6453f52", "metadata": {}, "outputs": [], @@ -450,13 +451,8 @@ " calc = Calculator(basename=basename, working_dir=working_dir)\n", " calc.structure = Structure.from_xyz(base_dir / r_xyz_file)\n", " # > Define input block for product structure\n", - " options= f\"\"\"\\\n", - " %neb\n", - " neb_end_xyzfile \"{p_xyz_path.name}\"\n", - " end\n", - " \"\"\"\n", + " calc.input.add_blocks(BlockNeb(neb_end_xyzfile=p_xyz_path.name))\n", " # > Add options to the calculator\n", - " calc.input.add_arbitrary_string(options, pos=ArbitraryStringPos.BEFORE_COORDS)\n", " calc.input.add_simple_keywords(*sk_list)\n", " calc.input.ncores = ncores\n", "\n", @@ -471,12 +467,8 @@ "\n", " calc = Calculator(basename=basename, working_dir=working_dir)\n", " calc.structure = Structure.from_xyz(base_dir / xyz_file)\n", - " options= \"\"\"\n", - " %scf\n", - " maxiter 200\n", - " end\"\"\"\n", + " calc.input.add_blocks(BlockScf(maxiter=200))\n", " calc.input.add_simple_keywords(*sk_list)\n", - " calc.input.add_arbitrary_string(options, pos=ArbitraryStringPos.BEFORE_COORDS)\n", " calc.input.ncores = ncores\n", "\n", " return calc" @@ -1125,7 +1117,7 @@ "source": [ "Overwrite the checkpoint path in the command directly, or replace the `ckpt_path: \"lightning_logs/rgd1/dmpnn/seed_0_YYYY-MM-DD_HH-MM-SS/checkpoints/epoch=XX-step=XXXXX.ckpt\"` with your actual log name from the training run in `conf/experiment/opi_tutorial/inference.yaml`. We will overide it so you don't need to open the file:\n", "```bash\n", - "python chemtorch_cli.py +experiment=opi_tutorial/inference ckpt_path=lightning_logs/rgd1/dmpnn/seed_0_2025-08-13_10-28-40/checkpoints/epoch=58-step=15989.ckpt\n", + "python chemtorch_cli.py +experiment=opi_tutorial/inference ckpt_path=\"lightning_logs/rgd1/dmpnn/seed_0_2025-08-13_10-28-40/checkpoints/epoch\\=58-step\\=15989.ckpt prediction_save_path=\"../predictions.csv\"\n", "```\n", "Remember that the file path will be different for you! The predictions will be saved to the specified CSV file, here `../predictions.csv`.\n", "We can then load the model's prediction and compare it to the value from the QM calculations. \n" @@ -1146,7 +1138,7 @@ " ╚═════╝╚═╝ ╚═╝╚══════╝╚═╝ ╚═╝ ╚═╝ ╚═════╝ ╚═╝ ╚═╝ ╚═════╝╚═╝ ╚═╝\n", "\n", "Predicting DataLoader 0: 100%|██████████| 1/1 [00:00<00:00, 8.18it/s]\n", - "Predictions saved to: /your/path/to/chemtorch/predictions/rgd1/dmpnn/seed_0_YYYY-MM-DD_HH-MM-SS/predictions.csv\n", + "Predictions saved to: /your/path/to/opi/docs/contents/notebooks/predictions.csv\n", "```" ] }, From d99407126a3beb7679d40557d8daf49e669359e5 Mon Sep 17 00:00:00 2001 From: Anton Zamyatin Date: Mon, 25 Aug 2025 15:00:09 +0200 Subject: [PATCH 3/6] Move deep learning barrier heights notebook to community example in index --- docs/index.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/index.md b/docs/index.md index ef0896eb..2f2c0505 100644 --- a/docs/index.md +++ b/docs/index.md @@ -72,7 +72,6 @@ contents/notebooks/ir_spectrum.ipynb contents/notebooks/extopt.ipynb contents/notebooks/opencosmors.ipynb contents/notebooks/moplot.ipynb -contents/notebooks/deep_learning_barrier_heights.ipynb ``` ```{toctree} @@ -82,6 +81,7 @@ contents/notebooks/deep_learning_barrier_heights.ipynb contents/notebooks/extrapolate_cps_cbs.ipynb contents/notebooks/extrapolate_ep1.ipynb contents/notebooks/wb97m2.ipynb +contents/notebooks/deep_learning_barrier_heights.ipynb ``` [documentation]: index.html From aba338d94abd39909bd2501281b0921f554914d7 Mon Sep 17 00:00:00 2001 From: Anton Zamyatin Date: Tue, 26 Aug 2025 16:13:31 +0200 Subject: [PATCH 4/6] Update deep_learning_barrier_heights notebook - Fix inference command and comment about escaping special characters - Replace all occurances of activation energy with barrier height for precision - Fix path to predictions.csv --- .../deep_learning_barrier_heights.ipynb | 21 ++++++++++--------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/docs/contents/notebooks/deep_learning_barrier_heights.ipynb b/docs/contents/notebooks/deep_learning_barrier_heights.ipynb index 9c428e4b..900d08f7 100644 --- a/docs/contents/notebooks/deep_learning_barrier_heights.ipynb +++ b/docs/contents/notebooks/deep_learning_barrier_heights.ipynb @@ -10,7 +10,7 @@ "1. First, we will show how to calculate the barrier height of a chemical reaction using ORCA with the ORCA Python inferface (OPI). \n", "2. Second, we use the ChemTorch framework to train and evaluate a graph neural network (GNN) on a curated subset of the popular [RGD1 dataset](https://www.nature.com/articles/s41597-023-02043-z) which contains precomputed barrier heights.\n", "\n", - "The term *barrier height*, also called *activation energy*, refers to the energy difference between reactant and transition state (TS)." + "The term *barrier height* refers to the energy difference between reactant and transition state (TS)." ] }, { @@ -843,7 +843,7 @@ "id": "dcdc8ea9", "metadata": {}, "source": [ - "Now, the HF/STO-3G values are of course quite off. With B3LYP-D3/TZVP, we obtain `67.935952` kcal/mol. The RGD1 dataset reports an activation energy of `68.35608500000312` kcal/mol (electron energy) for this reaction.\n", + "Now, the HF/STO-3G values are of course quite off. With B3LYP-D3/TZVP, we obtain `67.935952` kcal/mol. The RGD1 dataset reports an barrier height of `68.35608500000312` kcal/mol (electron energy) for this reaction.\n", "Given that the RGD1 dataset was constructed using a slightly different approach a small deviation of approximately 0.3 kcal/mol is tolerable.\n", "\n", "Now, we will save the reaction SMILES to a CSV file for later use. We will also overwrite `delta_E_act_kcal` to the value at the correct level of theory, since we will compare machine learning model predictions to this value later." @@ -1105,7 +1105,7 @@ "metadata": {}, "source": [ "### Step 4: Predict Barrier Height for Unseen Reaction\n", - "Now, we can load our trained model and use it to predict the activation energy of new reactions that aren't in the training dataset, for example the datapoint computed in Part I.\n", + "Now, we can load our trained model and use it to predict the barrier height of new reactions that aren't in the training dataset, for example the datapoint computed in Part I.\n", "\n", "We will use the config file in `conf/experiment/opi_tutorial/inference.yaml` (shipped with the ChemTorch repo, no need to make this file again)." ] @@ -1119,8 +1119,9 @@ "```bash\n", "python chemtorch_cli.py +experiment=opi_tutorial/inference ckpt_path=\"lightning_logs/rgd1/dmpnn/seed_0_2025-08-13_10-28-40/checkpoints/epoch\\=58-step\\=15989.ckpt prediction_save_path=\"../predictions.csv\"\n", "```\n", - "Remember that the file path will be different for you! The predictions will be saved to the specified CSV file, here `../predictions.csv`.\n", - "We can then load the model's prediction and compare it to the value from the QM calculations. \n" + "When overriding `ckpt_path` in the CLI make sure to wrap the path in double quotes and escape `=` characters by plycing `\\` in front of them like we did.\n", + "Remember that the file path will be different for you!\n", + "Assuming that the `chemtorch` folder is located in the same folder as this notebook, then the predictions will be saved as `predictions.csv` in the the same folder." ] }, { @@ -1147,7 +1148,7 @@ "id": "2b399a12", "metadata": {}, "source": [ - "Now we can load the prediction and compare it to the activation energy obtained from the ORCA calculation and the value reported by RGD1. Note that model training can differ between architectures (e.g. GPU vs CPU), so you might not get exactly the same numbers." + "Now we can load the prediction and compare it to the barrier height obtained from the ORCA calculation and the value reported by RGD1. Note that model training can differ between architectures (e.g. GPU vs CPU), so you might not get exactly the same numbers." ] }, { @@ -1158,13 +1159,13 @@ "outputs": [], "source": [ "# > Load predictions from the CSF file\n", - "df = pd.read_csv(\"../predictions.csv\")\n", + "df = pd.read_csv(\"predictions.csv\")\n", "delta_E_act_kcal_predicted = df[\"prediction\"].values[0]" ] }, { "cell_type": "code", - "execution_count": 34, + "execution_count": null, "id": "27e875b0", "metadata": {}, "outputs": [ @@ -1181,7 +1182,7 @@ } ], "source": [ - "rgd1_delta_E_act_kcal = 68.35608500000312 # > RGD1 value for the activation energy in kcal/mol\n", + "rgd1_delta_E_act_kcal = 68.35608500000312 # > RGD1 value for the barrier height in kcal/mol\n", "# > Calculate absolute errors\n", "error_to_orca_calc = delta_E_act_kcal_predicted - delta_E_act_kcal\n", "error_to_rgd1 = delta_E_act_kcal_predicted - rgd1_delta_E_act_kcal\n", @@ -1192,7 +1193,7 @@ "\n", "\n", "print(f\"\"\"\n", - "PREDICTED ACTIVATION ENERGY: {delta_E_act_kcal_predicted:.6f} kcal/mol\n", + "PREDICTED BARRIER HEIGHT: {delta_E_act_kcal_predicted:.6f} kcal/mol\n", "ERROR TO ORCA CALCULATION: {error_to_orca_calc:.6f} kcal/mol, {relative_error_to_orca_calc:.2%}\n", "ERROR TO REPORTED VALUE IN RGD1: {error_to_rgd1:.6f} kcal/mol, {relative_error_to_rgd1:.2%}\n", "\"\"\")" From 011f5f5993753ab6df25c5cfc75f08f46aa87f05 Mon Sep 17 00:00:00 2001 From: Anton Zamyatin Date: Wed, 12 Nov 2025 10:28:45 +0100 Subject: [PATCH 5/6] :pencil: Update chemtorch commands in tutorial --- .../deep_learning_barrier_heights.ipynb | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/docs/contents/notebooks/deep_learning_barrier_heights.ipynb b/docs/contents/notebooks/deep_learning_barrier_heights.ipynb index 900d08f7..abf8c9f0 100644 --- a/docs/contents/notebooks/deep_learning_barrier_heights.ipynb +++ b/docs/contents/notebooks/deep_learning_barrier_heights.ipynb @@ -977,16 +977,10 @@ "git clone -b tutorial/opi_orca https://github.com/heid-lab/chemtorch.git && \\\n", "cd chemtorch && \\\n", "conda deactivate && \\\n", - "conda create -n chemtorch python=3.10 && \\\n", + "conda env create -f env/environment.yml && \\\n", "conda activate chemtorch && \\\n", - "pip install rdkit numpy==1.26.4 scikit-learn pandas && \\\n", - "pip install torch && \\\n", - "pip install hydra-core && \\\n", - "pip install torch_geometric && \\\n", - "pip install torch_scatter torch_sparse torch_cluster torch_spline_conv -f https://data.pyg.org/whl/torch-2.5.0+cpu.html && \\\n", - "pip install wandb && \\\n", - "pip install ipykernel && \\\n", - "pip install -e .\n", + "uv sync && \\\n", + "uv pip install torch_scatter torch_sparse torch_cluster torch_spline_conv torch_geometric --no-build-isolation\n", "```" ] }, @@ -1013,14 +1007,16 @@ "To use an existing one, run the following command from the `chemtorch` project root:\n", "\n", "```bash\n", - "python chemtorch_cli.py +experiment=graph data_pipeline=rgd1 data_pipeline.data_source.data_path=\"../QM_data_precomputed.csv\"\n", + "python chemtorch_cli.py +experiment=graph data_module.data_pipeline=rgd1 data_module.data_pipeline.data_source.data_path=\"../QM_data_precomputed.csv\"\n", "```\n", "\n", "This tells ChemTorch to use the default graph learning configuration with the RGD1 data pipeline but use our own custom dataset specified via `data_path`.\n", "Under the hood, this setup will convert each reaction SMILES to a condensed graph of reaction (CGR), train a DMPNN, track metrics of interest and save the best performing model parameters for later.\n", "to the CLI as well as [Weights & Biases](https://wandb.ai/site/models/) which is a graphical user interface that can be accessed through the browser.\n", "\n", - "If you would like to make your own configuration file instead, an example is already included in your ChemTorch installation, and can be found in `conf/experiment/opi_tutorial/training.yaml`, so no need to create or change a file. Note that the important lines are setting the `data_pipeline` to `rgd1`, and `data_pipeline/data_source/data_path` to `\"../QM_data_precomputed.csv\"`, just as above. To launch the training process with the config file, run \n", + "If you would like to make your own configuration file instead, an example is already included in your ChemTorch installation, and can be found in `conf/experiment/opi_tutorial/training.yaml`, so no need to create or change a file.\n", + "Note that the important lines are setting the `data_module.data_pipeline` to `rgd1`, and `data_module.data_pipeline/data_source/data_path` to `\"../QM_data_precomputed.csv\"`, just as above.\n", + "To launch the training process with the config file, run \n", "```bash \n", "python chemtorch_cli.py +experiment=opi_tutorial/training\n", "```\n", From e50b5cc6e2b5ffec211295389bdaff74ef5869b8 Mon Sep 17 00:00:00 2001 From: Anton Zamyatin Date: Sun, 30 Nov 2025 19:09:46 +0100 Subject: [PATCH 6/6] Add version disclaimer and fix typo --- .../notebooks/deep_learning_barrier_heights.ipynb | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/docs/contents/notebooks/deep_learning_barrier_heights.ipynb b/docs/contents/notebooks/deep_learning_barrier_heights.ipynb index abf8c9f0..a95716bc 100644 --- a/docs/contents/notebooks/deep_learning_barrier_heights.ipynb +++ b/docs/contents/notebooks/deep_learning_barrier_heights.ipynb @@ -6,7 +6,15 @@ "metadata": {}, "source": [ "# Deep Learning of Reaction Barrier Heights with ChemTorch\n", - "This two-part tutorial showcases how ORCA integrates into downstream deep learning workflows by serving as a data source training and evaluation.\n", + "\n", + "----\n", + "\n", + "**Disclaimer** \n", + "This community notebook was kindly supplied by [Anton Zamyatin](https://github.com/tonyzamyatin) and [Esther Heid](https://github.com/hesther) for OPI version 1.0.0.\n", + "\n", + "----\n", + "\n", + "This two-part tutorial showcases how ORCA integrates into downstream deep learning workflows by serving as a data source for training and evaluation.\n", "1. First, we will show how to calculate the barrier height of a chemical reaction using ORCA with the ORCA Python inferface (OPI). \n", "2. Second, we use the ChemTorch framework to train and evaluate a graph neural network (GNN) on a curated subset of the popular [RGD1 dataset](https://www.nature.com/articles/s41597-023-02043-z) which contains precomputed barrier heights.\n", "\n",