19_molecular_dynamics.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyN7soXHJFokcQtIXbExp0Cb",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/francescopatane96/Computer_aided_drug_discovery_kit/blob/main/19_molecular_dynamics.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "AQ0bcJn30a0a",
        "outputId": "e491477d-f182-45d3-8cea-86ebcfea5155"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Collecting condacolab\n",
            "  Downloading condacolab-0.1.7-py3-none-any.whl (7.2 kB)\n",
            "Installing collected packages: condacolab\n",
            "Successfully installed condacolab-0.1.7\n",
            "⏬ Downloading https://github.com/conda-forge/miniforge/releases/download/23.1.0-1/Mambaforge-23.1.0-1-Linux-x86_64.sh...\n",
            "📦 Installing...\n",
            "📌 Adjusting configuration...\n",
            "🩹 Patching environment...\n",
            "⏲ Done in 0:00:17\n",
            "🔁 Restarting kernel...\n"
          ]
        }
      ],
      "source": [
        "try:\n",
        "    import google.colab\n",
        "\n",
        "    !pip install condacolab\n",
        "    import condacolab\n",
        "\n",
        "    condacolab.install()\n",
        "except ModuleNotFoundError:\n",
        "    pass"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "try:\n",
        "    import condacolab\n",
        "    from google.colab import files\n",
        "    from IPython.display import clear_output\n",
        "\n",
        "    condacolab.check()\n",
        "    !conda install -q -y -c conda-forge mdtraj openmm openmmforcefields openff-toolkit pdbfixer pypdb rdkit\n",
        "except ModuleNotFoundError:\n",
        "    on_colab = False\n",
        "else:\n",
        "    # check if installation was succesful\n",
        "    try:\n",
        "        import rdkit\n",
        "\n",
        "        on_colab = True\n",
        "        clear_output()  # clear the excessive installation outputs\n",
        "        print(\"Dependencies successfully installed!\")\n",
        "    except ModuleNotFoundError:\n",
        "        print(\"Error while installing dependencies!\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "49Zn15tn3TW2",
        "outputId": "74750be4-39f0-48f1-94c9-d98b17f1cb16"
      },
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Dependencies successfully installed!\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "import sys\n",
        "\n",
        "if not on_colab and sys.platform.startswith((\"linux\", \"darwin\")):\n",
        "    !mamba install -q -y -c conda-forge openmmforcefields\n",
        "    # Notes:\n",
        "    # - If you do not have mamba installed, install it or use conda instead\n",
        "    # - Under MacOS with an M1 chip you may need to use\n",
        "    #   CONDA_SUBDIR=osx-64 in front of the above command"
      ],
      "metadata": {
        "id": "GMnfN9Lj3sYu"
      },
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import copy\n",
        "from pathlib import Path\n",
        "\n",
        "import requests\n",
        "from IPython.display import display\n",
        "import numpy as np\n",
        "from rdkit import Chem\n",
        "from rdkit.Chem import Draw\n",
        "from rdkit.Chem import AllChem\n",
        "import mdtraj as md\n",
        "import pdbfixer\n",
        "import openmm as mm\n",
        "import openmm.app as app\n",
        "from openmm import unit\n",
        "from openff.toolkit.topology import Molecule, Topology\n",
        "from openmmforcefields.generators import GAFFTemplateGenerator"
      ],
      "metadata": {
        "id": "gDQMA8HG34p0"
      },
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "pdbid = \"3POZ\"\n",
        "ligand_name = \"03P\"\n",
        "pdb_path = f\"{pdbid}.pdb\"\n",
        "pdb_url = f\"https://files.rcsb.org/download/{pdbid}.pdb\""
      ],
      "metadata": {
        "id": "VMUN_T4w4CMA"
      },
      "execution_count": 30,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "r = requests.get(pdb_url)\n",
        "r.raise_for_status()\n",
        "with open(pdb_path, \"wb\") as f:\n",
        "    f.write(r.content)"
      ],
      "metadata": {
        "id": "n262tgb24D0h"
      },
      "execution_count": 31,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "prepare protein"
      ],
      "metadata": {
        "id": "mSd6JJSW6hNK"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def prepare_protein(\n",
        "    pdb_file, ignore_missing_residues=True, ignore_terminal_missing_residues=True, ph=7.0\n",
        "):\n",
        "    \"\"\"\n",
        "    Use pdbfixer to prepare the protein from a PDB file. Hetero atoms such as ligands are\n",
        "    removed and non-standard residues replaced. Missing atoms to existing residues are added.\n",
        "    Missing residues are ignored by default, but can be included.\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    pdb_file: pathlib.Path or str\n",
        "        PDB file containing the system to simulate.\n",
        "    ignore_missing_residues: bool, optional\n",
        "        If missing residues should be ignored or built.\n",
        "    ignore_terminal_missing_residues: bool, optional\n",
        "        If missing residues at the beginning and the end of a chain should be ignored or built.\n",
        "    ph: float, optional\n",
        "        pH value used to determine protonation state of residues\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    fixer: pdbfixer.pdbfixer.PDBFixer\n",
        "        Prepared protein system.\n",
        "    \"\"\"\n",
        "    fixer = pdbfixer.PDBFixer(str(pdb_file))\n",
        "    fixer.removeHeterogens()  # co-crystallized ligands are unknown to PDBFixer\n",
        "    fixer.findMissingResidues()  # identify missing residues, needed for identification of missing atoms\n",
        "\n",
        "    # if missing terminal residues shall be ignored, remove them from the dictionary\n",
        "    if ignore_terminal_missing_residues:\n",
        "        chains = list(fixer.topology.chains())\n",
        "        keys = fixer.missingResidues.keys()\n",
        "        for key in list(keys):\n",
        "            chain = chains[key[0]]\n",
        "            if key[1] == 0 or key[1] == len(list(chain.residues())):\n",
        "                del fixer.missingResidues[key]\n",
        "\n",
        "    # if all missing residues shall be ignored ignored, clear the dictionary\n",
        "    if ignore_missing_residues:\n",
        "        fixer.missingResidues = {}\n",
        "\n",
        "    fixer.findNonstandardResidues()  # find non-standard residue\n",
        "    fixer.replaceNonstandardResidues()  # replace non-standard residues with standard one\n",
        "    fixer.findMissingAtoms()  # find missing heavy atoms\n",
        "    fixer.addMissingAtoms()  # add missing atoms and residues\n",
        "    fixer.addMissingHydrogens(ph)  # add missing hydrogens\n",
        "    return fixer"
      ],
      "metadata": {
        "id": "heG0ehyB4hvm"
      },
      "execution_count": 32,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# prepare protein and build only missing non-terminal residues\n",
        "\n",
        "prepared_protein = prepare_protein(pdb_path, ignore_missing_residues=False)"
      ],
      "metadata": {
        "id": "zt5VDKGl6Y1W"
      },
      "execution_count": 33,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "prepare ligand"
      ],
      "metadata": {
        "id": "K9rPZl1z6eVx"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def prepare_ligand(pdb_file, resname, smiles, depict=True):\n",
        "    \"\"\"\n",
        "    Prepare a ligand from a PDB file via adding hydrogens and assigning bond orders. A depiction\n",
        "    of the ligand before and after preparation is rendered in 2D to allow an inspection of the\n",
        "    results. Huge thanks to @j-wags for the suggestion.\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    pdb_file: pathlib.PosixPath\n",
        "       PDB file containing the ligand of interest.\n",
        "    resname: str\n",
        "        Three character residue name of the ligand.\n",
        "    smiles : str\n",
        "        SMILES string of the ligand informing about correct protonation and bond orders.\n",
        "    depict: bool, optional\n",
        "        show a 2D representation of the ligand\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    prepared_ligand: rdkit.Chem.rdchem.Mol\n",
        "        Prepared ligand.\n",
        "    \"\"\"\n",
        "    # split molecule\n",
        "    rdkit_mol = Chem.MolFromPDBFile(str(pdb_file))\n",
        "    rdkit_mol_split = Chem.rdmolops.SplitMolByPDBResidues(rdkit_mol)\n",
        "\n",
        "    # extract the ligand and remove any already present hydrogens\n",
        "    ligand = rdkit_mol_split[resname]\n",
        "    ligand = Chem.RemoveHs(ligand)\n",
        "\n",
        "    # assign bond orders from template\n",
        "    reference_mol = Chem.MolFromSmiles(smiles)\n",
        "    prepared_ligand = AllChem.AssignBondOrdersFromTemplate(reference_mol, ligand)\n",
        "    prepared_ligand.AddConformer(ligand.GetConformer(0))\n",
        "\n",
        "    # protonate ligand\n",
        "    prepared_ligand = Chem.rdmolops.AddHs(prepared_ligand, addCoords=True)\n",
        "    prepared_ligand = Chem.MolFromMolBlock(Chem.MolToMolBlock(prepared_ligand))\n",
        "\n",
        "    # 2D depiction\n",
        "    if depict:\n",
        "        ligand_2d = copy.deepcopy(ligand)\n",
        "        prepared_ligand_2d = copy.deepcopy(prepared_ligand)\n",
        "        AllChem.Compute2DCoords(ligand_2d)\n",
        "        AllChem.Compute2DCoords(prepared_ligand_2d)\n",
        "        display(\n",
        "            Draw.MolsToGridImage(\n",
        "                [ligand_2d, prepared_ligand_2d], molsPerRow=2, legends=[\"original\", \"prepared\"]\n",
        "            )\n",
        "        )\n",
        "\n",
        "    # return ligand\n",
        "    return prepared_ligand"
      ],
      "metadata": {
        "id": "E_TWK1no4t0u"
      },
      "execution_count": 34,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "smiles = \"CC(C)(O)CC(=O)NCCn1ccc2ncnc(Nc3ccc(Oc4cccc(c4)C(F)(F)F)c(Cl)c3)c12\"\n",
        "rdkit_ligand = prepare_ligand(pdb_path, ligand_name, smiles)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 252
        },
        "id": "pB-3RiPj400o",
        "outputId": "068c34aa-ca4e-42fd-81a5-f507d2b40329"
      },
      "execution_count": 35,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "[17:28:45] WARNING: More than one matching pattern found - picking one\n",
            "\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAADICAIAAABJdyC1AAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO2dd1hTZ9vA75MQwkaEACGAAipuKqDiQFQUpe49aR11AIrW14qftoLVKm3VuqoV25fia2tddVd9XThBBRRRlCEzhCEQdiAheb4/jm9KUSHjhOTA87t69To5nuc+N/rjzhnPIBBCgMFgMHSAoe0EMBgMRlFwwcJgMLQBFywMBkMbcMHCYDC0ARcsDAZDG3DBwmAwtAEXLAwGQxtwwcJgMLQBFywMBkMbcMHCYDC0ARcsDAZDG3DBwmAwtAEXLAwGQxtwwcJgMLQBFywMBkMbcMHCYDC0ARcsDAZDG3DBwmAwtAEXLAwGQxtwwcJgMLQBFywMBkMbcMHCYDC0ARcsDAZDG9pswRKL/96QyaC2FgBAIgGJRItJYdoU2KvWp80WrKAgkMkAAJYtg8xMcHOD2lq4fh0uXNB2Zpi2Avaq9WGGh4drOweNcP48eHpCdTX89ReMHAn5+fD0KTg4gEgEN26AQABVVcBggIkJEIS2c8XQk7Iy7FVro6ftBDTImTNAEMDnAwC4ukJFBaSkAJsNq1f/fYy+PtjZgYdHhYFBMI/Hs7e379SpE7lhY2OjrcwxdKElryrZ7EB7e3vsFVW05YIVFAQMBsTHv/24fj0MGQKBgbByJeTmAp8P+flQWAjZ2cDh1D1+/FuT5p07d167dm1wcHBr542hFaRXy5ZBcDDw+U28Ej1+/HuT4zt37vzFF18EBQVpJVu6QyCEtJ2DRjhzBiZPBoKAb76BmTOhoACGDYPbt8HSEnr3/vuw+nrg86GwsDgz8yqfz+fz+Xl5eXl5ednZ2bW1tdbW1llZWXp6bbmsY1Tj/n2wsQGB4K1X9+5BWBjY2kKnTsDjAZcLHTuKmMzfJZJcsVhcUVFBqpWdnV1TU2NjY4O9UhHUpsnORhYWyN4e5eUp11Amk/Xo0QMAjh8/rpnUMDSG9MrBAfH5b/fs3IkYDATw93+9e1fKf8tYLJajo+PQoUPr6+tJr06cOKHVn4CutNm3hADQ0ADz5oFQCG5uwOMp15YgiBUrVgDAnj17NJIchrbIverbF+zs3u5cswbq6yEnB+7dg2PH4PvvYfp04ZQpU/r378/lchsaGnJzc1++fKmvr4+9UgttV0wNEhqKAJC9PSopUaV5TU2NpaUlAMTFxVGdGobGqOBVfX19ZmZmfHw8auTVw4cPNZhlG6XNdmu4dQuWLwcmEy5cgO7dVYnAYrFKS0vv379fW1s7bdo0qhPE0BK5Vxcvgquroq2YTKaFhYWdnR008kokEk2dOlWDubZF2uYtYWFh4Zo19TIZbN4MQ4eqHic4OFhPT+/+/YcCQT112WHoSmOvhgxRPU4jr8QtH41pRBssWDKZbP78+VlZrosXZ61fr1YoR0fHoKCY/PxXBw6wKcoOQ1eo9uoOn//ywAF9irJrN2j7npR6tm7dCgDW1tYCgUD9aPfvIwDUsSOqqVE/GIbGUOvVvXsIAFlZodpa9YO1I9pawYqLi2OxWAwG4+rVq1TFHDAAAaDDh6mKh6EfmvPq55+pitcuaFO3hEKhcNasWRKJJDQ01M/Pj6qwISEAALt3QxvtY4tpAQ15tXIlAPZKSeja072ysjI/P7+oqCg/P7+4uJjP5xcVFT148CArK2vw4MG3b9+msBuxRALOzsDnw/Xr4OtLVVSMLoK90nF0enCAUCg8fvy4oaFhQSMKCwsFAkEtORHRO+jr669evZraQQ8sFixfDl99BffvY7HaAqp5FRIScurUqdOnT//xxx9MJlP9NFgsWLYMwsIgNhZ7pSi6e4VVW1vbs2fPnJyc9/6psbExj8ezsbEh/29nZ8flcv/73/8ePXp06tSpp0+fpjaZsjI4dgymTwcbG7hyBZycgMUCZ2dIToYOHcDBgdqzYTSIOl69ePEiNTX1zz//nDJlCiXJlJXB8eMwdepbr8aOpSRqW0Z3C9bdu3eHDRump6c3Z84cuUDy/5uYmLzbpLi4uFOnThKJJC0tzdnZmdp8FiwAMzPYuxeWLgU/PzAygo8/hl9/BWdnGDaM2lNhNIg6Xm3cuPHrr7/28fGJiYmhKp/GXkVGUhW1zaK7t4RJSUkAMH/+/KioKAWbWFtbz5w588iRIwcOHNixYwe1+RgaQo8ecP3624//+Q88eADPnsHatdSeB6NZSK/mzZv33XffWVtbK9JE7pVQKDQ3N799+3Z8fLynpycl+TTxCtM8uvuWMC4uDgAGDhyoVKvVq1cDwOHDhysrK1s8WFmWLoXo6LezdwcEwNatgEdW0A7Sq06dOtnY2Pj4+CjYivQqOjp6/vz5APDjjz+qkwNCkJLy98fGXmGaR9cL1qBBg5Rq1a9fv2HDhlVWVkZHR1OeEpMJISGQmkp5YEzrQXpFwuVyFWwl94rD4TCZzGPHjhUWFqqcw5498NFHIC962Csl0HI/sA/w5s0bgiCMjY0lEomybckn7l26dJFKpZQkEx+Pxo5Fr1+//VhdjerqUH09QgjV1iKxmJKTqEVGXYa2U6AHcq8WLVoEAD/88IPibeVeTZo0CQDCwsJUyyE+HrHZiCDQzp1o4kSUnf12f3W1avE0iA56paMF6/z58wAwYsQIFdo2NDSQT9wvXLigfiaVlahLFwSA/u//1A9GJUWSoiU5S0L5oXuL9y7LWabtdOiB3KvevXsDQGxsrOJt5V5t27YNAKytrUUikbIJVFUhV1cEgJYvR127IgC0fr2yMTSLjnulo7eE5HW7l5eXCm2ZTCY5YTYlc6QFBUFGBvTpA199pX4wKtlbvHe19eoIXsQKzgpt50IbSK/69ev38uVLNpvdr18/xdvKvbp586anp2dxcfHx48eVTSAwEFJToU8fqKiA9HTo0wc2bVI2hmbRca/aYMECgCVLlpiZmV2/fv3Zs2fqpPHvf8PRo2BsDCdOgKGhOpGoJ1+S76zvDAAE4PWkFIX0ytzcXCqVuru7s9nKTcIh92rixImg/Dfiv//9a17eD6amaM4cOHYMe6UKuliwpFJpfHw8KP+KUI6Zmdmnn34KAPv27VM5jbS0tws3HTyo4hSAGqWfYb9b1bcAoAE1aDsXeiD3qqamBlT6OpR7lZWVZWtr++TJk9u3byvYNiUlZeXK4Nu312zZcmP7dgCAAwewV8qj7XvS9/D06VMAcHZ2VidIeno6g8Fgs9lFRUUqNK+trR03bl7nzqmffqpOFhqhTlYnEAvqZfXbCrat46/7teTXw2/wVBItI/fq448/BoA//vhDhSByr9auXQsAU6ZMUaSVSCRyc3MDgE8++USXvSoQF+i4V7pYsH766ScAmDdvnppxxo8fDwBbt25Voe2yZcsAwN3ds6pKpmYalLMqbxUniXOj8oa2E6EZcq+srKwAICcnR7U4pFfr16+3srJatWqVTNayIUuWLAGArl27km8nddYr6yTrm5U3tZ1Ic+hiwVqwYAEA7Nu3T804165dAwA7O7t6sg+Cwpw8eRIA2Gx2YmKimjlQzqXyS0QCwUpkxVYr8YYLg/7n1VdffQUAXC5X5Thyr6qqqhQ5PiEhgSAIAwMDcvQF9koddHFojppP3OX4+vra2tqWlZVxOBx7e3s7OztnZ2cul0tuODs7Ozg4sFisJq1yc3PJy6sffvhBqbdIrQBfwv8051ME6Fvet17G6v79tDdIr8iJFpTtkNwYuVfk0vMteuXu7n769OnMzEyyPwT2Sh10rmCVl5enpaUZGBj07dtXnTgNDQ1BQUGFhYUdOnQoLy9PSUlJaTwaAgAAGAyGra2to6Mjj8dzcHBwdHTkcrnfffddWVnZtGnTAgMD1UmAcqRI+kn2JyUNJf5m/qutV2s7HZoh96qgoADU+DpU2avffvsNe6U+OlewTp06JZPJXF1d9fVVn5+/pqZm9uzZFy9eNDAw2L9/f58+fSQSiUAgKCgoyMzMlG/k5uYKBAKBQNC4rbW1tbW19c8//6z2j0Ixu17tuiW6xWPxjnQ+opuvnHUZuVdffPGFl5eXaldYjb06fPiwr69vY52wV62Btu9Jm7Jz5042m62npxcWFqZCT2KEUElJyeDBgwGgY8eO9+7dS01N/frrr997pEQiycvLu3///vHjx3fu3Llq1aphw4YBAI/HU2FIkEa5ceMGU485PHJ4TFWMtnOhJZR71cyR2CvNoXMFq66ujnzrDAAuLi7nzp1TqnlGhrRPHy8AcHZ2TktLk+9X0FGZTNa9e3cAOHnypHJ5a5Li4mJyDc7NmzdrOxe6oiGvFAR7RRU6V7BIbt26JX+G5evr+/z5c0VaPXqErK2Rp2exl9fQwsJC+X4/P7/JkycreOr9+/cDwNChQ1XJuyUyMjICAgIWL1586NChS5cuJScnC4XC5pvIZDLyPbqPj09DQ4Mmsmo/UOuVUmCvKEFHCxZCSCKR7Nmzp0OHDgAwfPh/vvgCVVY2d/zFi8jYGAEgf39UXf2PTi6pqamK/5PU1NR07NgRAB4+fKhy8u9FLBY7Ojq+e1duYGDg7Ow8ZMiQGTNmhIaGHjp06Pz58/Hx8aRzERERAMDhcPLz86nNp31CemVhYQEAw4cfXbtWda+UAntFCbpbsEiKi4uDgtZYWUkBEJeLoqPRe7vpRUUhFgsBoE8/fc98L/Hx8YGBgRMmTFDwpGQn5vnz56uXe1P+9a9/AYCZmVlQUNDixYvHjBnTs2dPU1PTZp4wmpub6+vrEwTx119/UZtMO4f0ytJSUa8WLKBgHiFNexUYGNjmvdLdOd0b8+QJrFwJ9+8DAHh6wt690PglT3g4bN4MALBpE4SHA/HOi46zZ8+mpKTMmDGjS5cuYrG4xSGvOTk5ffu6eXgs+O233QpP8dYCly9fHjduHJPJvH37NvnsVo5IJGrymoncyMjIqKio4HA4JiYmmZmZ1OSBaUQTr/btg8a9HVr0Slk07VVRUdGcOXN4PB7ZKczR0dHIyAgA9PT06uvrS0pK2oJX2q6YiiKToehoxOUiAMRgoE2b0NGjCCF09Spavx4xmejgwQ+2LS0tXbNmDUIoJydn06ZNipxuzpxaAPTVV9QkX1hYaGNjAwDfffedUg3z8vLIm5dHjx5RkwrmnzTxKizsb682bGjBKxUgvfryS2qiNfHq4MGDzfymGxsb9+zZ08/Pb/ny5fT1ijYFi6S6GoWFIQsLFBmJvLxQZib67juUmoqSkpprVVpaumTJkpycnNjYWAUL1t27CABxOKi2Vt2cpVKpr68vAIwZM0aFSVDJ+4iAgAB188B8mA959ewZxSfStFcikej169d37949ceJERERESEjIjBkzPDw8yPJE4uLigmjrFT1uCZtQVQU3boBYDFevQvfuMGkSdOvW3PFlZWWzZ8+ePn16WVmZSCTaTF7ot8SAAfD4MfzyCyxapFa24eHhmzdvtrGxefr0qa2trbLNc3JyunTpQhBEVlYWj8dTKxVMsyjrlWpoy6vKysq8vLzc3FypVDp+/Hi6eqXtiqkiZ86gK1fQkSPIxwelprZwsAq3hAihI0cQAOrd+/2PYxXk9u3bTCaTwWBcu3ZN5SDTpk0DgK+oukHFfBilvFIN7JU60LVgXbqEbtxAMhkaMwZltDRTfm1t7ZUrVxBC5eXlN28qOnuGWIx4PASAbqg6j0tpaSn5vllNJ+7cuQMAHA5HtS7aGMVRyivVwF6pA10LFkIoOxtt3ox+/FGDp/j6a9SrF7p6VZW2MpmMnEjX29tb/QEZ/fv3B4BffvlFzTiYFsnJQZs3owMHNHiKr79GvXuj//5Xlbbt3CsaF6y4OASAPD0pDpuYiMiLsPv30Z07KDgY5eYihNCuXcrFISc/srCwyJYv5KQGR44cAYDevXsrMl0cRh2wV7qMLs7priBkT5aCAorDlpa+jVlcDEIhlJfDN98AALx8CQIBVFcrFCQ+Pn7Dhg0EQURFRXXq1En9rGbPns3j8Z4/fx4TE6N+NEwzYK90GZ2bXkZxbG2BIKCoCKRSYDKpjHzyJKSkQFoazJ8PpqYwYABcuAAAMHMm3L8PHTqAvT04OoK9PdjbQ5cul62t9e3t7R0dHQ0NDQGgqqpq3rx5YrH4888/JxfdVB8Wi7V06dKwsLA9e/aMGDGCkpiY99L6Xs2aBffuNfXKxeWyjQ32qim07NYgh8OBkhIoLAQbG8piXr8OxcUwdy6cPQsMBly+DAcOwMKFwGBATg7ExoJI9I/jO3Z0LStLI7etrKx4PJ5QKMzNzfXw8Hjw4IE6s3o14c2bN46OjmKxOC0tzcXFhaqwmHdpZa9yc+HBg6ZeWVr2KC19RW5jr+TQ+JYQNHP1bmIC5uYAAGZmYGoKVlZAELBmDZiYwI0bUFsLJSXw9ClcuAAHDsCXX6Lx4wcNHz68S5cuBgYGJSUlSUlJb968MTY23rVrF4VWAQCHw5kzZ45MJlNn7TKMIrSyV9evv/UqKemtVxs3onHjBn7Iq4cPHyrYl1ARaOaVth+iqYWfHwJAujN+s7CwMCEhgex//MUXX1AePzk5mSAIU1PT8vJyyoNj5OiyV6tWrSIHNic1P7xDGWjkFb2vsLp3/57D6VNQEKXtRN5iY2Pj7u7+7bffAsDhw4erFXyUqjC9e/cePnx4VVVVVJSu/Mhtku7dd1hbu+mmV9HR0QEBAaD8utPNQCOv6F2wjIzK3rx5XlAgaPnQVsTDw2PIkCHl5eXkO2NqWbJkCUEQ69at69y5s7e399y5c9etW7d37967Fy/C48dQWAhNHkqWlsJ338H330NpKeXJtFWMjEqLi5/prFe2trZMJvO3334rKipSM2ZFRcX48eOTk5M/5NW9S5feetUEbXml7Us8tSC/ZIKDg7WdSFNOnDgBAF27dlVhtHPzzJgx473/jpv790cACADp6yNnZ+TtjebPRzdvotmzUUEBys9Hs2ZRm0kbRve9mjBhAgB8aLECxSF18vb2/pBXWwYMeOsVm42cndGwYWj+fHTrlra8onG3BgDgcrkAUEB5nxm1mTp1qpOTU3p6+pUrV+RTiavPoUOHTp48aWJiEhsba2xszOfzc3Jy8vPz+Xy+J0IglQKfD8XFkJkJmZlw9y4MGQIdOwI5MtbCAqqrwcSEqmTaMLrv1cKFCy9cuLB///5169a1OL/bh/jpp59InUaNGhUWFvZerzwAoKHhH14BwNChWvOqNasj5dy9excABg0apO1E3gP5xMHPz4+qgM+fPyfnYztKztj0IUQilJaGbt1C0dEoPR3J5w+ZPx/RZ+pu7UILr9zc3ADgyJEjqsWR6xQREaGEVzdvatcrehesjIwMAOjcubO2E3kPQqHQxMQEAJ5RMamSSCQiV0/47LPPlGu5axfatAl99RXavVv9NNoJtPCK7Nng7u6uQpDq6uqePXsCwKJFi+jlFY0fuiOEjh49ymQymdT2R6aIDh06fPLJJwDw448/qh8tODj42bNnPXv2VPrd0Oefw9q1sHw5IAS3b6ufSZuHLl7l5eXZ2NgkJibeu3dP2SArV65MSUnp1atXQ0ODul4RRKt61ZrVkUIEAoGfnx8AEAQxb948bafzftLS0hgMhqGhYUlJiTpxjh8/DgAGBgZPnz5VMcT33yMANGaMOmm0B+jlFTlr6PTp05VqLtfp+++/V9erHTsQAKLuuUeL0LJgXblyhZxikcPhnD9/XtvpNIe/vz8AbN++XeUIGRkZZmZmABAZGal6HkIhMjFBANRP+tuGoJ1XGzZs0NfXZzKZmZmZCjZMT08nddq2bRsdvaJZwRKJRCEhIQRBAMCoUaN0f0m1K1euAACPxxOrtEqUWCweOHAgAMyYMUPdVAIDEQBatkzdOG0R+no1d+5cAFi7dq0irerq6tzd3QFg2rRplHkVFNSaXtGpYD1//px8QMhmsyMiIijv4qQh+vTpAwDHjh1Toe2qVasAwNnZmYIxE6mpiMFAhoZIvfvTtgetvfrmm2+YTObChQsVaRISEgIALi4uy5cvp8yrtDTEYCAjo9bxih4FSyaTHTp0iHz52qNHj8TERG1npAQ//fQTAHgqPyPcpUuXCIJgsVixsbHUpOLvjwCQGvenbYy24RWfz1fk+Js3bxIEoa+vv3v3bvp6RYOCVVRUNG7cOPIVQUBAQHV1tbYzUo6amhpTU1OCIHr06DF//vz169fv27fv/PnziYmJRUVFH2qVl5dnZWUFALuUnZKyGa5cKfvoo51jxqh2f9rGaG9eSSSSjRs3hoeHa8Ir4Ucf7Ro7Vv0pm1tE1wvW1atXyW7HVlZW586d03Y6KhIYGPiht7QGBgZdu3YdMWLEJ598snHjxoMHD164cCE2NtbLywsA/P39qZ27Vp3707YE9opCr2QyWY8ePQDgjz/+oCrmh9DdCfzKy8vDwsL27duHEPL19Y2OjqbT6mnvEBcXV11dXVBQkJubm5+fT64Qx+fzy8rK3ns8QRCWlpYvX74kvw+p4tChQ8uXL/fy8oqNjaUwLI3AXmnCq59++ikwMHDQoEEPHjygMOy76G7BcnJyys7OZrPZ33zzzZo1a8g3OG2P2tpa0rD8/PycnBw+n8/n85OSklgsVnh4+IIFCyg/naOjY2lpaVxcHPmeqL2BvaK3V5q+hFMNPp/fsWNHPT29x48fazuXtsb69esBYM6cOdpORAtgrzRH63ilo1dYMpnMwMCgoaFBJBKpPBgd817y8/OdnJxkMtnOnTv79etnb2/P4/HayV8y9kpztI5XOlqwAMDe3p68miUXucVQhVQqdXZ2BoDc3Fz5TgsLCy6Xa2dn5+zsLN9wdnZ2cHBgsVjaS5Z6sFcaonW80t35sLhcbn5+fkFBARaLWsLDw3Nzcy0sLKZPn15YWJiXlycQCIRCoVAoTElJaXIwg8GwtbW1tLTs0KHDyZMnbShcRkZLYK80ROt4pdMFC3RyEjVaExMTs337dgaDcerUqZEjR8r3C4XCzMxMgUBQUFDQeCM3N1cgEAgEAgAYN25cfHy89nKnBuyVJpB7tX379unTp1taWpL7qfdKo0/I1GHp0qUAcODAAW0n0nYoLi4mf103b96sYJOGhoa8vLx9+/YxGAx9ff2CggKNZtgKYK8oR+7VunXruFyug4NDi4OxVfZKd+fDcnHp17//xyIRV9uJtBEQQosWLSooKPDx8dm4caOCrZhMpr29/YoVKyZNmiQWi8nhILQGe0Utjb1KSUkpKChwcnJq8XabIIikpKSEhISJEycq55WaxVVz/PQTAkDKzoOI+RAREREAwOFwVJuKICYmBgCsra1FIhHlubUm2CtqkXtFfgtyOBxFxjY2NDSQT+i3bdumlFe6W7DOnUMAaNw4befRJnj48KG+vj5BEOpM8+Tp6QkAUVFR1OWlBbBXFCL3aufOneSG4uOcdu7cCQC+vr5KeaW7BevRIwSAVJqxGvMPhEKhk5MTqL0Y9a+//goAffr0oXZ4YyuDvaIKuVchISHkhoLTcpFUVlaSMwhu2bJFca90t2Dl5SEAxOVqOw/6M3v2bADo379/fX29OnHq6urICTljYmKoyq31wV5RhdyrmTNnAoCnp6eygq1cuRIAFi5cqLhXuluwxGLEYCAmE69NpRb79+8HAHNzc8Vn0W2GsLAwAJgyZYr6obQF9ooS5F59/fXX5Mbr16+VDZKens5gMNhsNjk5vSJe6W7BQghxOAgA0f9NutZITk42NDQEgN9//52SgEVFRQYGBgwGIyMjg5KAWgF7pSZyr7799lty47ffflMtFLmE9fr16xX0SncLVk0Nio5GFy6gggJUX4/kU5J9eM47zD+orq4mZylavnw5hWHJgf6ff/45hTFbE9KrixexVyoi9+qzzz4jN5apMaH79evXAYDL5ZJrl7Xole4WrL/+QgMGIIkE7dmDnj1DCxa83S/fwDQP+fKlT58+tbW1FIZ99uwZAJiamlZUVFAYttXAXqmJ3Kv58+cDQK9evWpqatQJSE6ov337dkW80t2OowAwciTIl3csL4f4eKD/yJBWQiqVJiYmEgRx4MAB8qKdKvr06ePj41NVVUW+NKQj2CuVkXu1d+/eyspKIyOjEydOkJPiqwz56P3UqVOKeKVGwZIv90puZGTAoUNw44bqARtBTiHh5QW5uZCfDwAgEgGfD3w+6OrsEroFk8ns1KkTQig5OZny4ORaPrt375ZKpZQHbwr2SpeQe/Xy5cuzZ88+evSIXPJeHQICAmxsbBISEsaMGQMteaVGwYqO/nsjLw+2boVx4yAjA37+WfWYAABw4wa4uwM5xetXX8GxYwAAXC5MngyTJ0MbnSGSeshr7N27d8tkMmojT5o0ycXFJSsr69KlS9RGbooGvBIKAbBXaiD3CiHUq1cv9QOy2ewlS5YAQGJiYsteqX7rOXQo+vZb9O23aMECFBWFbtx4u3/+fJVD1tWhNWsQQSAAFByMsrMRQujePVRa+nd4+QameSQSCTmk6/Lly5QHJ8dhkF+2GqSJV0IhGjwYRUUhJfv7NPZqxQrslVpowqv8/HxTU9PFixdv2LChea/UKFgBAaiiAlVUoIUL0R9/ILJLvliMPv0UHTiAwsJQcbFS8V6+RP36IQCkp4fCwnA3GQogx3mNHTuW8siTJk0iCMLNzY3yyP+giVe7dyMABIDs7FBEBBIKFYmBvaIcTXhVVVWFFPBKjYIlX2x24UJUXY3mzEHHj6OVK9H168jODgEgQ0O0bBlKTVUkWGSkyMgIAaCuXdGjR6onhWlMWVmZsbExQRApKSkUhj127BgAGBsba7zLexOvgoPRN98gN7e3ZcvUFK1eLSWvlz7A4cPYK+rRoldqFCyB4B8b9fUoKQmVlSGE0J07aOJExGAgAMRgoMmTxffufSjMmzdvJk2a5OU1DQAFBKCqKtUzwrwLuSh5YGAgVQHT09NNTU2h1UZBy7369de3dWrIEBQRgcaPJ+/xvuzff/z48e8uYiwUCmfNmirHJYMAABYVSURBVIW90hDa8kqT/bAyMlBICDIyQgBnfHzc3d2jo6ObrA177do1Ozs7ALCwsDh3LkuDybRXUlNTGQyGkZFRSUmJ+tHq6ur69esHALNmzVI/mnJkZ6NVq5CJyduy5eaGvvmmMjCQra9PPo0dMWLEpUuXyAG0cq86dOhw9mxWa6faDtCWV5rvOFpUJN20qdf/JvTq0qXLjz/+WFNTIxaLw8LCGAwGAAwaNEiFgUgYBRk7diwAREREqB8qKCiI/EfUWq/RsjK0fTvicsmydWzChC+//HLVqlXm5uakYD179pw6dSrplZeXF/ZKc2jFq1bq6V5TU/Pjjz+6uLiQVllaWpIvGvT09LZu3dqAH4RqksuXLwMAj8cTi8XqxLlw4QJBEGw2OyEhgarcVKSuDv3733XDhpno6wOAmZlZSEhIeHg4OVEv9qp10IpXrTo0RyqVnj9/fsiQIQBgb2/P5XLvffjZFoYqZDIZOebr+PHjKgfJzc3t2LEjAOzdu5fC3NTk7t2748ePJ1dvZjKZtra2LBbLwsICe9UKaMUr7YwlnDFjBgBs375dK2dvhxw4cAAABg8erFpziUQyePBgABg3bpwOzt739OnTgIAA+Tp3GzZs0HZG7YXW90o7BYvsLKvmBJgYxampqSFXXnr48KEKzUNDQ8mLYkqesGqI7OzsoUOHYq9akyZeWVpaOjs7jxo1KiAgIDQ09NChQ9euXXv+/HllZeV7m6vglXYKFjm+cb4afeIxyrJu3ToAmDdvnrINr1y5wmAw9PT0dP8+C3vV+si9Ki0tbWb8jbW1tbu7+8SJE1esWEH23lLNK+0spIoXs2x9goODd+3a9ccffxAE0aVLF0dHR3t7e3t7+06dOjUz2r6oqGjBggUymWzr1q3kw0ddBnvV+jT2KjQ01MTEhMVi1dXVicXikpIScvHUnJyc4uLi4uLixMREAJg7d67KXuGC1V6or6/X09Pr0KHD0aNHm/yRgYGBnZ2ds7Mzl8ttsjFv3rzCwsIRI0aQX6Q6Dvaq9VHEKw8PDzMzM3LJCYlEYmtrq7JXBNLGtBolJSUcDsfCwqKMnJMB0xL/LStLF4kAYL6Njbme0l8z9fX1gwcPTkxMHDFixKxZs/h8fm5uLp/P5/P5eXl5IpHoQw0JgrC0tExOTiaXCdBxsFfKQjuvtHOFZWlpqa+vLxQKRSIRtdPLtVWe1dQE83iGDBWnA1q3bh05d8eZM2fkfSzlCIVCgUBQUFBAXsDLNzIyMths9oYNG2hRrQB7pTy080o7V1gA4OjomJeXl5WV1blzZ60kQC925OXpEYQZk7mIq/Qa65cuXZowYYKent6dO3e8vLyUaku733zslVLQziutTZGMHzcoyzI7u0Vcbo1UGllQIFb4a4bP53/66acIoW+//VZZqwCAXtUKsFfKQ3pVSxOvcMGiB/oMBjkj5nd5eZECwcJXr3Lr6lps1dDQMGfOnNLSUn9//9WrV2s6SV0Ae6UUcq++p4lXahWshoYGldtisZQihMczYDAAYI61tSObnVpbO//ly79aerQcHh5+7949Ho935MgRon3MAYy9Ugq5V7Mae9VsjyrQqlcqFiyxWLxkyZK+ffvW1taqFgGLpRrdjYyO9ujh37FjrUy2KStrU3Z27QembI+JiYmIiGAwGEeOHLGysmrlPFXm7t27774gVxzslWr8w6vsbJ31SsWCNXPmzKioqFevXn3++eeqRcBiqYwRk7nFyemrTp0MGIy/SkvXx8aSawU2pri4eO7cuVKpNDw8fOTIkVrJUwVevnw5YsSIpUuXvnjxQrUI2CuVIb3aJPcqLk4XvVKtP356ejo5V4ytre2zZ89UiHDhwgUA8Pf3Vy0BDELotUi0JDmZ06mTgYHBwYMH5fulUqmfnx8A+Pj40G6KlcWLFwNAr169VFv/FXulPqRX1k5OBgYGBw4ckO/XBa9UH0t4+fJlclLHoUOHqjCCPy4uDgC6dOmicgIYhJBIJAoJCSG/eyZPnlxWVob+N7acw+Hk5+drO0Glqa6udnV1BYDVq1er0Bx7RQlNvCotLUW64ZVag5937NhhZmZmZGS0f/9+pRqWl5f7+/sTBGFra6tOAhiSP//808LCAgAcHR0jIyNZLBZBEOfPn9d2XiqSkJCgr69PEMTFixeVaoi9ohYd9Erd2RoWLFjAYrG6d+8uVGzNJYRQTEyMg4MDABgaGu7cuVPNBDAkr1+/7t+/PwCQb23oPsUKuZCUtbV1QUGBgk2wV5rg9evXAwYM0B2v1C1YUql01KhRoNicHhKJJCwsjMlkAsCAAQPS09PVPDumMRKJRD7cYeTIkQL5skY0RCqV+vr6AsDYsWNbfOCAvdIoEomkR48ebDZbF7yiYD6sysrKvn37cjicuLi4Zg7Lysoi55FgMpmhoaFqzgONeZeGhgYvL69169aRL5u5XO4NOi9nzOfzycnh9u3b18xh2CtNs3//fgAwNjYmbw+16xU1E/hlZGS4uLgMHDhQKpW+94ATJ0506NCBvBm+ffs2JSfFNEYmk02aNOmjjz5CCBUWFpJvcwiCCAkJoe/v8JkzZwCAzWYnJSW99wDslaZ59uwZOYzm999/1wWvKJtx9OrVqxwOZ9u2bU32V1RUzJs3j3zdMH36dPI1FoZyli9fzuFw5CvxNjQ0yO+S5q5cWVRfr930VGbJkiXv7eWAvWoFqquru3fvDo0WTNW6V1ROkbxjx46OHTsWFxfL98TGxjo7OwOAqanpoUOHKDwXpjHrTp82s7dfvnx5k/0xMTHOrq5T4uNHPn16u7xcK7mpSU1NDfk7ExISIt+JvWodFi9bBgDkgJbG+7XoFcVzus+aNcvPzw/98zlo//798XNQzXGquHhYYuKgmJiampp3/7Ssri4kPd0jPt4zPn5nXp5Y99a8aRF5L4fz589jr1qNq6WlExITvaZPl1+2N0ZbXlFcsKRS6aBBg3777Tdy/RIGg0HrZyi6z53ycp/ERK/ExBONLmybIEPo96Iir4QEj/j4eSkpuXV1rZkhJXz//fcA0LFjR/IVO/ZK02SLRN6JiR7x8ec+vJ6NVryifgK/yMjITZs2FRUVOTo6/uc//xk2bBi18TFynlVXL09LYxJED2PjyG7dmj/4RU3Nhqys/Pp6YyYzkMerlUoBYK61tYGqs022JjKZrG/fvi9fvpTJZNgrTSNGaOGrV6m1taMtLLY7Ozd/cCt7ReUUyRUVFYGBgceOHSP7mAUGBmKrNMeLmprgjAwpgBmDsc7BocXjexkb/96jx7bc3Pz6+le1taEODgBAi2pFevXixQvsVeuwOy8vtbbWgc3+slOnFg+We1UsFqfW1q7TtFdUXarFxcWRw6GNjIyWLl1KEASbzX769ClV8TFN+LdAMDopySM+PjwrS6mGNQ0N4VlZ+/j8o4WFmkmNSrBXrcyd8nLP+PhBCQkv3/dItBlqpNJW8IqKgiWRhG3a1KSf8fLlywGgZ8+eqo25x3wIqUwmlEiEEolIKk2qqlr48mWN8uPmN2dnayI3isFetSJyrx5XVIx79uxYUZEKQVrBK7ULVnY28vY+OWwY2Zes/n/9MkQiUZ8+fQAgODhY3VNgGvFaJPo8I+NoYWFCVZXKQWhQsLBXrUtjr6obGlR756fzBevoUWRmhgBkXbo8eKefcXJysoGBAUEQ586dU+ssmEa8Fon28Pmva2tVuLCSo07b1gB71erQxStVC1ZFBQoIQAAIAE2dikpL33vUzp07AYDD4dB6IK5O8Vok2piZeam0tJC2ndebA3ulJejilUoFq6QEOTkhAGRign75pZkDZTLZxx9/DACjR49WYZI/zLu8Folo8bBcFUpLsVfagi5eNfv2MSjo7w2EYO1a+PJLWLoU3ryBIUPA0xMSEmDRomYCEAQRFRVlY2Nz7dq13bt3q/c+EwMA0FFPz83ERNtZqMeHvCouBm9v7JVWoI1XzVWzBQuQTIZkMrRgAbp2DR0+jBBClZVo8WJUXY0U7mf8119/kW+jnzx5onaFxdCfZryqqcFeYZqh2Sus9HTYtAk2bQKJBAoKwNERAMDUFOrrwdgYWCwFa6K/v39gYGB9ff3cuXNVXhYM03ZoxisjI+wVphmaLVhdu8KWLbBlC7BY0L8/XL8OAJCUBAr0f23Czp07yaEVa9euVTVVTFsBe4VRFWZ4ePgH/7CsDPr1e7sxejSIRBAdDXl58MUXin8Nkujp6Q0bNiwqKio2NrZfv37khCGYdgr2CqMq1A9+boY9e/asXr3a39//r7/+arWTYto82Kv2Q6sWLIRQVFRUQEAAS8kvUgymGbBX7YdWLVgYDAajDjSYXQSDwWBIcMHCYDC0ARcsDAZDG3DBwmAwtAEXLAwGQxtwwcJgMLQBFywMBkMbcMHCYDC0ARcsDAZDG3DBwmAwtAEXLAwGQxtwwcJgMLQBFywMBkMbcMHCYDC0ARcsDAZDG3DBwmAwtAEXLAwGQxtwwcJgMLQBFywMBkMbcMHCYDC0ARcsDAZDG3DBwmAwtAEXLAwGQxtwwfoHJSUld+7caX6txtTU1OTkZBWCf/bZZ8OGDVM1NQxGUQ4fPkwQhFQq1XYi1IML1j+4efOmj4+PWCxu5pgtW7b861//arWUMBiMHD1tJ6BbTJ06tbKyks1mN3PMzz//jJfLxmC0QjsqWAih27dvJycnGxoaDh48uGfPngCQlZWVlZXVr1+/Y8eOMRgMf3//O3fuzJs3j8FgAEBNTc3NmzfT09OtrKzMzMyMjIz8/PwSEhLq6up8fX2TkpLEYnGPHj3+/PPP4uJid3f3kSNHkueqrKy8evVqTk6Oubn5qFGjnJyctPmTY6gAIXTq1KmxY8cmJSU9ffrUwcFh1KhRxsbGAHDr1q3OnTtXV1ffuHFjwIABgwcPBoBXr17FxcWJxeKBAwe6ubkBQHJyskgk6tSp07Vr1+rq6gYNGtSrVy8yeENDQ0xMzKtXr/T19eXHC4XCmJiYCRMmnD59uri4+JNPPjE3N5dIJLdu3Xr16hWPx/P19e3QoYM8w0ePHiUlJUml0hcvXmjhL6h1QO0DiUTi7+9vbGw8bdq0kSNHMhiMHTt2IIR++eUXLpfbrVu3Xr16jRkz5vjx4wBQV1eHEEpOTubxeA4ODrNnzx47dqy5ufmUKVMQQvPmzRs9ejRCaM2aNW5ubt26dZs4ceLkyZMZDMbu3bvJ07m6ug4ePDggIGDo0KH6+voxMTEIocWLF3t7e2vtrwCjHg0NDQDg6enp4uIydepUW1tbFxcXgUCAEBoyZIiPj4+JiUnv3r1JryIiIhgMxsiRI/38/AiC2L9/P0IoNDTUxcXF2tr6448/9vb2Jghi165dZPAZM2Z07dp1ypQpo0aNYjKZ+/btQwg9fvwYAHx9fe3s7Lp06VJeXl5aWtq3b19ra+upU6c6OzvzeLzCwkKEkEgk8vf3NzAw8PPzGz9+vL29PQA0NDRo7S9LY7SXghUZGclgMB4/fkx+3LZtG5PJzMrK+uWXXwDg7Nmz5P7GBcvNzc3T07OyspL8o48//vjdgmVsbPz06VPygMWLF/fs2ZPcrq2tlZ/a29t72rRpCBcsmkMWrBEjRojFYoRQUVGRlZXV0qVLEUJDhgxxcnIiawdCKDU1lcFg7Nmzh/wYGhpqZWUlEolCQ0MJgrh9+za5f/369SwWq6CgACFUUVEhP9Hq1avt7e3R/wrW//3f/8n/aMWKFdbW1kVFRQihyspKa2vrjRs3IoS+/PJLNpudkJBAHhYZGdlWC1Z7uSW8f/9+//79PT09yY9BQUEbNmyIjY0lP44fP77J8cXFxUlJScePHzc1NW0mrIODA3n1DgBdu3b9888/yW1DQ0OBQJCSkiIUCtlsdm5uLpU/DEZ7zJ49m8ViAYC1tfWUKVMePnxI7vf29raxsSG379+/L5PJpFIpWTjq6+tLSkpev34NAI6OjvI3xUuXLo2IiEhISBg3bpyZmVltbe2LFy8KCgpEIlF+fr5EIiEPmzZtmvzsd+/edXJyOnv2LPnR1taWTODKlSszZsxwd3dvhb8B7dJeClZhYaG1tbX8o7m5uYGBgUAgsLCweO/x2dnZAKDUsyeCIBBCACCTyYKDg6Oiotzd3W1sbLKysszNzdXKHqMzEAQh37a0tHzz5s27xwgEAhaLJf86BIAZM2aQbpDPRuXNAaCkpAQAIiMjQ0NDLS0tHRwciouLEUIymey9kTkczvXr18mPrq6uDg4OAJCZmTlu3DhKfkAdp70ULHt7+8TERPnH4uLiuro6R0fHqqqq9x5PVrfGOorF4ubfHso5f/78oUOHEhMTP/roIwAICgqSfw9j2hJpaWkuLi7v7ufxeA0NDVFRUeQj+Q+Rnp4OAC4uLrm5uYGBgXv37g0ODgaAqKioRYsWvbeJvb29p6cneeHWGCsrq9LSUvnHNtkDi6S99MMaPXr0kydPrl69Sn6MiIgwNDT09vb+0PEODg7dunU7ePAg+eQiPj7+zp07Cp5LIBCw2exu3boBQGlp6aNHj9ROH6NzxMXFnT9//t2HCQDg7e3NYDB++OEH+Z53L8Tq6+vDw8O5XK6Hh0dubq5MJhs6dCgAyGSyBw8efOikw4cPP3nyZF5eHvmxpqamtrYWAHx8fP7888+ysjIAqKur+/XXX9X98XSV9nKFNXPmzBs3bkyYMMHd3b2ioiInJ+fw4cO2trYfOp7JZEZGRk6fPr1bt248Hi8rK4u89laESZMmbd682cPDw9nZOTk5uXv37o2//TC0Jjw8/MyZM/X19Xfv3p08efLnn3/+7jEuLi67du1as2bNxYsXbW1tX758aWVldf/+fQDIy8vz8vLicrlPnjyprq4+c+aMoaGhh4eHq6vr1KlThw4d+vjxYy6X+6Gzh4WF3bt3r2/fvkOGDBGLxXFxcZGRkbNnzw4LC7t69Wrv3r0HDBjw9OlT8suyTfL2sUs74fXr10+ePDE1NfXw8LCysgIA8tG4r68v+WyiqKgoKSlp9OjR5MeysrJHjx7V1dV5e3sHBAQYGhqePn36+fPnYrHY3d09LS2ttLR00KBBZPCcnJysrKzhw4cDQGlp6cWLF1ks1vDhw1ksVmZm5sCBA1NSUmpqavr376+1nx+jBlKpVE9Pb+XKlX369NHT0+vXrx95yw8At27dMjY2HjBgQOPjMzMz4+LihEKhq6vr8OHD9fT01q9fHx0dvXfv3qKiImdn5yFDhsgfblZWVh47dqyystLb29vNze3ixYvTpk2rqKi4fv36qFGjGj9plUqld+/eTUlJMTU1HTx4sPyetLq6+tatW/n5+Z6enq6urleuXJk+fXrjJ25tg/ZVsNShb9++I0eO3L17t7YTwWgHsmBFRkYuWbJEtQjr168/ceJEZmYmtYm1K9rLLaEKbNmyhcPhkD0/z5079/z58/3792s7KQymXYML1gfh8XhHjhwJCwurrq52dXU9duwYnmsBg9Eu+JYQg8HQhvbSrQGDwbQBcMHCYDC0ARcsDAZDG3DBwmAwtOH/AdILxee0qRQ7AAAEhnpUWHRyZGtpdFBLTCByZGtpdCAyMDIzLjAzLjIAAHiclZZtbFNVGMef23Z9XRlbu27d1q6sQwbK1q5vm1p6yebKxlw7UKajw00H6ECCGl9A1IEQ5Itk2XApopkRJiGKLIo202xnjMSZjA/qhxHEGTRqiC+Q6IwhQb3nOXe66PlyTnJ7/z2999fn9ZxzbeydWVCGFdjQKNcy5VqhXL2SQdNJ56SVOuXmafF7JOVe2WACqKLK49Mqn75AMvKXIvx0anU1/hBgtOux4LP20/ErN/xySKr9Dy3AaHER2rqM+8CPfasoTc9oGkarCzLaGhHaT60KbluIQwsxWqMILXX43FzFMR+HFma0JhHa5FxF6rCmhkOLMNpaYU9vXRA3rZqFuhpGaxbO6SiP5vch7W4R2nIauJNhHi2AtBYR2tG0Mpy8uPlDSEuI0PBOwjxaBGlJ4Zyeq+V4Ws3i1iqeU171VuMLletEaAE63lM9LVfmJN08DR+tXC+c04oFFSJp52lBpN0jQqNN3/pMI6WZNB5KY7CGala894oX753/h3kYbIP4CsKBqUtlmwjMPPPmzmhTlFdstUi7T4SmhCzjvovbCDVIu1+4rdIRXpOGkdYuQqPFkcrBYjNi2FTTEmrYNorAju+Mmme+5a4frHJTIrQtOIKUlrsgo3XNanl0iMCu7KbFFuHlgLXoJuH6MAd4NNZUD4jQFMN223eEOBuzupV2Cm9+JSFeEhitSzxufo6nahYeFN+YMaX6sm5K080bp55BHhI3roq3ijNat3gaeOcj1dXNwrTKAKPR6lWX3RbVtC3i21WQdSk1TVJNS6imbRWhXZhWRjNvdatltIeFD4Ij8r/bleaf7YodGh4RpuU2cDxV49YjQovSRu1Yw9v82BFkmwitlC7kf6JtElTQScmi0YKUBboskPRgMIElGyxWsOboJOtisOaCJQ/MNsizgdkO2XYw5oM9H4wOcBSAoRAchWB0gskJhiKAYigqBiiBYhdo3eByQ0kpuEtBqR2dB1we0C8BfRl4F+mkMi94y8G7FFy36DR26o9eedjtKs7SaN0uj97kNDoKDXqzPduSZyuYVbyQQB3LekZWjP2w72WZfunx9JObbZlxqlu+eJ0Mda3F+cTeAbKjbwD1iwe6SOr69hjVpdODJHKxLUr11PJRctZoJ1ST5v2xAfclfKZ+V698IvQu6rD7s5g8dHIV1R/YfpWjb50eo7rDMSr/tsGH78I3w3L3h2fQhpHB1+Qbg8/h/PSTf5CcTR2odc69cu/2E6g3/vKo/MbRetTBmvPyqaZXUW/99Krsy70Dtblcu3r90JfIfLtzUv7olB//N3HbGXnlRBptu3CoXba/UoU+9k8ekotn9qB+PP57LJm8HbVvLi5PHFuK9h+/eXUsu/8avvt18uD41PkEPhNP1I+/1DuKev8iL9G90I469nQj8UbTqL/ftZgMfzWO+ufLPeSpT4pQX5r9jhS8fxHttD1/ltgHnkD7Hzt4mWzOHGF+TUkTZHgEdX6akMmJj1F/nskhfXtcyHH8DWuyOfpm+XL6AAADgnpUWHRNT0wgcmRraXQgMjAyMy4wMy4yAAB4nH1Wy24kNwy8+yv0Axb4FKVDDuvx7iYIdgZIvPmHADnm/5Gi2u6WASFjH3rY1WyyWCzNU8nPH6+///1vOT/y+vRUCv3P/xij/KVE9PSj5EV5+fr9t3u5vX15+YjcHj/vb38WHcUEz+DvM/bL2+PHR4TLvTxT5dbakGLVXC1GoUrzcz0qCbTqI4KiPHM1I+6+QWq5JXI0lWFFK0d0tQ3QDmAzCiTKCwnXDdATqNVFIzLREHfbAduRMUxbi7zv0kbbACOBXjtHQ/hZqrKS8gbZgeRKqg5+uAYZed/gBnBWia2TAqetd9mxwzQTNtXQvA+S1HYJmQEMEN7YW85IjVlohxQgsxkQNJsJd6XYIXUihZ16zyFqJ0x0h7TZzxALMAnyRxMMfIec04kazRTKAdS897ZjnVv5liSpKVnyz23Y2LIUQGoNzCc0c1Lw0C1NHchWWbk3zuF3l952cuMxCVXMD8yDpiGC/nZazxn1yugZ2waagvHMjnrJIeHt5O48FwlTpZ00RcoD1Du50YAAhDDMXUOSMwKJbA0NSYrZdccRSr/9AzrREEsu2nBCwTukT82hWdSJlIMZWtkB23w30qllu9zQz27oEodvmAdx6r0Ht9hWmQsECskjdaxVrGnfqUNyPsgZreWAIGjpPHasp8flAN0iwNEki2jXkPJhHQxTw1hgDR2V7uaj0+CkigecDW8fw/t2iVQTqRWOGfA1wdu7bXWEu49cDYb/uSULMNmxNST1w5EaeMTrkT43c1/pdLlW+2BHLemg0cN3g9fIAqBkZsNskR8ybbR17X5kham3HtMcR29bleBg+XVSZd6mnqxRyE7LX++vn06b4/x5edxfr/MH8yt2nTKQfGnXUcIw//OrzbtxHQsMe+bF/Tm9kPVy+SNwpu8ZsPSh27fvv/Dl3Uf4zMz0HumLGx+5xuK6nGFZ3ZXTZWR10QNzvoqzaACXSJaNh2TxvnxijfQDcjbGWTGMQi7iaJIhvjjTAbkicmQ56ZRZL60RO7KcRMgsF0u8OATn9i+ROCAnMdLnE9jSZbMZw9Oz5zn4uZ7LomYOrOGykJzhJaJTLHrSoDb1sQT8Pe8VyXoDyGUlZuCkRbPchPlnUWgWjYR6yY8OpJ70GL9HTjZM3isYM9tZR27Cqvv8/vGbDddP/wEAW9ws4CoW4wAAAe56VFh0U01JTEVTIHJka2l0IDIwMjMuMDMuMgAAeJwlUj1rXDEQ/Csp74EstN+7l8agYJLGTpPKuAipDmxsEoMb//iMdPCOd5o3uzM7q8fvT4/z9vaJeJ7mvCf85unhmDgdeB6OewB8P/EB+OS53vN0d+C5OyYfOD0fk758nm64s5p7464+QtpXIKYR3qQrKc7SSaOa9sgMAqDddUS29eawjZQLF0ooImUXuRfaoS2l8kaMJaKNXmymsstCxRc0jL1cN2YVkLmhrjqoAFlPgh8YE5KxCqNTitoiecmmeBDqoZLDarn0nkUmbdmNRP0qC1fJxVLDMJs1KL3aEqn0NUx2IoI/QGbktCNhC9VmvcpyIQOtPMUAJee2OTo5hmbkoiaagAiDwRdyYXW5ktRikDTqCNNBAkfECkAMHb4AF0F0oyPS3Qexsoq2pUAh0FfYxhfq4plIF1aHKQHgQQyCdxpmSAg1IbZSCyTrMANIlAhVmBR7WykGWW6KDEtdURdz7cZJ2NtCwuwqXYz7sUIsZ91e5BoXeWktu4F9bcYIKt1uBDGvXaQx/qAvk41a2pJr6n3b9j3EnbDdFl7gEjI24BK5QHgsRmFB1o72+/315eff17fz6Jd/P17eni9/Lu+d2sevMzEkm6C7fHw7I+PP/6d7qlqArZu9AAACo3pUWHRyZGtpdFBLTDEgcmRraXQgMjAyMy4wMy4yAAB4nHWRa0gUURTH79yZnX372nF1d93ZyVXRKNqgF8W2NxCkB+1CkGAZVygaxIKKCqHHWob5JYjdEnthkBbSw4Ussdi9ywoZ2IfqQ2JmEFT4oQx6EIHQzPGRJl44nN+5c//n/A8zkbw3hrRjR1MHa1GmRYUWUc5IyrXM80ak6hlzIgHg5wAK6DLhXwaJIJL/8wLlgi9I0WC6Ey+YoeTmp8XbmKYG45k6dzHBnGUWWJvveXb96Tw7YnZXMGxCVJ/EzRi3Yl7FnAEJBm28YjSr2GqjVjuyZyN7DrLnImueii0OFec5qEWiNknFpnxFyldMTuosQMZC6izUrlwqNruo0a0gj4rdHgUVUY+X8jL1yiou8qlY9lGsUEFBXoWKS8qxWIyy/KjYj/wlyF+KJEFzI6Iin+z1GDAvexXR7DI5C42iRbJZ8xwFrzjN9fR/R2X1iaXJz2cuEL2oV2JssrovpfP21zdYR91WuA83xdnBi3Hgs+fq2O5vDSGdfUNtbO1wdVDnwYp+1muSmM5sW3MoLo/Am8rGKOlc/QB4jfwyRDrubNT5keM7Cd6+n9S51tlPfuwMgBZ96CJ7H/eAh0TbdfKn7STcDx39zbL31AILriYSbegE3vX1ALl5pRJ41boB0r3lKvD+5+MkkLsB2FLCb9rR8RZ63qUZ8qR7JcwNL+shy9Pt4O1Faw2RLq+AHWOZVuJ5cwL4cNWvUCSyHjjws4qkr5WC/1uT40lbbAK07yMtqcGBMLypClemzkf7gZuz/Ew4XQMcOr6Z+YPtwJ8ac1jXuxTwl9F6duyZG3hk7CMreDgMPh2nepkUPwL+D7WMsn19l6b2GuTSrCsBnN/OWCb9FNj5F5B5vo2k17DHAAADZ3pUWHRNT0wxIHJka2l0IDIwMjMuMDMuMgAAeJx9VltuHDkM/J9T6AIWxLf0GdvZIFhkBki8e4cA+cz9sUW13a0AwraNQQ+7RInFYvXcSl7fX//++bucF7/ebqW0//kfY5R/pbV2+1bypjx//vL1Xl7ePj1/RF4e/9zffhTpRQlr8Pcn9tPb49tHhMq9PLVK7j64aFUTjVFabfO6lnICtdqIaFGeqKo26rZBSnlJ5HDhoUUqRXTRDVAPoGsLJMobDpMN0BIo1VgiMtFgM90B/cgYKu6Rz419+AYYCbTaKRzhJ65C0oQ2yA4k1SZi4IdqNG3WN7gBnNZG2psAJ94779ihNhO6SEg+B0miu4REAAYIdzLPHokScdshGcgsBgTNYsJMWuyQMpFM1nrPJkpv6OgOqbOewRpgEuQPZzR8h5zdiRquAuUAqta771gnL38lSaLSNPknHzq2LAWQUgP9CcmcLWjIlqYOpFcS6k7Z/G7cfSc3GpNQQf/APGgazKhvp/XsUa+EmjFtoCkIa3bUczYJuzczozlI6GrbSZO5PEC9NdM2IABuaOauIM4egURSR0GcYjbZcYSjv/wCnSiIOAdtWMOBd0ibmkOxOCdSDiJoZQf0uTfSiWa55Khn13SOwzfUolHqvQd5bE+ZAwQKm0XqWCqrS9+pg7M/yBnu2SAImjuNHevpcdlA0whwNMlqbVeQ0GEdBFNDW2ANHSfd9UemwXFlCzgbdh/D+naIRBIpFY4Z8DXG7l23OsLTR44Gwf9MkwWY7NgaktjhSA4esT3S52TuTzpdzmsfZDhLOmj0sF3jJfIAUDKRorfID5l627p2P7LC1L3HNMfRfauSz/fXP94hx1vl+XF/vd4q0GPR690BIRe/XhAESz+/6nwal9krTJcWT9d0OJLLu2kGzvQ9AwCdGcd74MxJ7T3SF3c90o7FRTXDvLqlpmvw6ooHhhb30wQuEcdeWMSLl+WKNYI3c0LOkihPjMHni7I2aWBbnEYn5IrwkeUsm+d52xpJbpHlJILtoPvkgfO0mOYlEgfkJIbztPiQdVKTCDlrnl2c47YMniZOzqKF55I1IlMmctIgOpWxBOw97xXx1AqQi8QpAyctMuWALJfekl1ZI5p6QCo5icEvpSNy8pAqXzWd3z9+ZeH+9h8xRM1GQe89mQAAAbx6VFh0U01JTEVTMSByZGtpdCAyMDIzLjAzLjIAAHicJVJLiuUwDLzKLBNwjPWXeMwq0Mt+h8h26BP04afkByFxypKqXOX7Pu7zeJ/3ffx9n9/3/UPP8/BDz/H99PJ4P4L3fXydeL5OQHLi/9/50Pnz/PCf3+OKGa6S45KplhnjdfnMIpNx6fTI8A0tSq9xGfbSDVBOIlJpyIycABnqybVn5bL6QEmYcPEUkiVAYlKK2rhoejWg0yqCGlBdVBsKFY+x5jL2wsgXZhpLNFZsprIh9woZPJVSeTeWC9eQSRH5me66Ikd/OWx3kUYBCBy3NfI0DYiEA9RMBFacA/+sLq1nwZtYJIMmWjzHCzUiVgBi6fIGXARa1gS/KErAwyoKBO6EgEnhInZoimdCLpxbpgSAFzEKfNIygzfoCbH2yyaTrWp3JJuoHYUTbWmYfcYW4wBteznr5pFPMuSl1VICKeyKFVS6mQSJdsZpjMUrpixL7bnFXJCHhOFmMwdZmxlIymEM5IkScQ+GZe5IAIHZpsYU7ECKLTTDF4hbXVG4L7bNJEdujATUtlOAwj3FBg7Hua8AT7ZQBVJlSeP8/Q8i/5obBwHzDAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<IPython.core.display.Image object>"
            ]
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "def rdkit_to_openmm(rdkit_mol, name=\"LIG\"):\n",
        "    \"\"\"\n",
        "    Convert an RDKit molecule to an OpenMM molecule.\n",
        "    Inspired by @hannahbrucemcdonald and @glass-w.\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    rdkit_mol: rdkit.Chem.rdchem.Mol\n",
        "        RDKit molecule to convert.\n",
        "    name: str\n",
        "        Molecule name.\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    omm_molecule: openmm.app.Modeller\n",
        "        OpenMM modeller object holding the molecule of interest.\n",
        "    \"\"\"\n",
        "    # convert RDKit to OpenFF\n",
        "    off_mol = Molecule.from_rdkit(rdkit_mol)\n",
        "\n",
        "    # add name for molecule\n",
        "    off_mol.name = name\n",
        "\n",
        "    # add names for atoms\n",
        "    element_counter_dict = {}\n",
        "    for off_atom, rdkit_atom in zip(off_mol.atoms, rdkit_mol.GetAtoms()):\n",
        "        element = rdkit_atom.GetSymbol()\n",
        "        if element in element_counter_dict.keys():\n",
        "            element_counter_dict[element] += 1\n",
        "        else:\n",
        "            element_counter_dict[element] = 1\n",
        "        off_atom.name = element + str(element_counter_dict[element])\n",
        "\n",
        "    # convert from OpenFF to OpenMM\n",
        "    off_mol_topology = off_mol.to_topology()\n",
        "    mol_topology = off_mol_topology.to_openmm()\n",
        "    mol_positions = off_mol.conformers[0]\n",
        "\n",
        "    # convert units from Ångström to nanometers\n",
        "    # since OpenMM works in nm\n",
        "    mol_positions = mol_positions.to(\"nanometers\")\n",
        "\n",
        "    # combine topology and positions in modeller object\n",
        "    omm_mol = app.Modeller(mol_topology, mol_positions)\n",
        "\n",
        "    return omm_mol"
      ],
      "metadata": {
        "id": "70Ua-pHD5R_T"
      },
      "execution_count": 36,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "omm_ligand = rdkit_to_openmm(rdkit_ligand, ligand_name)"
      ],
      "metadata": {
        "id": "xzSL0hJ05TEn"
      },
      "execution_count": 37,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def merge_protein_and_ligand(protein, ligand):\n",
        "    \"\"\"\n",
        "    Merge two OpenMM objects.\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    protein: pdbfixer.pdbfixer.PDBFixer\n",
        "        Protein to merge.\n",
        "    ligand: openmm.app.Modeller\n",
        "        Ligand to merge.\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    complex_topology: openmm.app.topology.Topology\n",
        "        The merged topology.\n",
        "    complex_positions: openmm.unit.quantity.Quantity\n",
        "        The merged positions.\n",
        "    \"\"\"\n",
        "    # combine topologies\n",
        "    md_protein_topology = md.Topology.from_openmm(protein.topology)  # using mdtraj for protein top\n",
        "    md_ligand_topology = md.Topology.from_openmm(ligand.topology)  # using mdtraj for ligand top\n",
        "    md_complex_topology = md_protein_topology.join(md_ligand_topology)  # add them together\n",
        "    complex_topology = md_complex_topology.to_openmm()\n",
        "\n",
        "    # combine positions\n",
        "    total_atoms = len(protein.positions) + len(ligand.positions)\n",
        "\n",
        "    # create an array for storing all atom positions as tupels containing a value and a unit\n",
        "    # called OpenMM Quantities\n",
        "    complex_positions = unit.Quantity(np.zeros([total_atoms, 3]), unit=unit.nanometers)\n",
        "    complex_positions[: len(protein.positions)] = protein.positions  # add protein positions\n",
        "    complex_positions[len(protein.positions) :] = ligand.positions  # add ligand positions\n",
        "\n",
        "    return complex_topology, complex_positions"
      ],
      "metadata": {
        "id": "yxnPCLmL5uIT"
      },
      "execution_count": 38,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "complex_topology, complex_positions = merge_protein_and_ligand(prepared_protein, omm_ligand)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "GASjYvZZ5wOi",
        "outputId": "0ac6827d-f8a7-4e50-8d99-e83c07c971d8"
      },
      "execution_count": 39,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "/usr/local/lib/python3.10/site-packages/openmm/unit/quantity.py:750: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.\n",
            "  self._value[key] = value / self.unit\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "print(\"Complex topology has\", complex_topology.getNumAtoms(), \"atoms.\")\n",
        "# NBVAL_CHECK_OUTPUT"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "DlNviihM647L",
        "outputId": "2f8001e5-3c54-4f49-b581-68cd03126c65"
      },
      "execution_count": 58,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Complex topology has 5561 atoms.\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "def generate_forcefield(\n",
        "    rdkit_mol=None, protein_ff=\"amber14-all.xml\", solvent_ff=\"amber14/tip3pfb.xml\"\n",
        "):\n",
        "    \"\"\"\n",
        "    Generate an OpenMM Forcefield object and register a small molecule.\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    rdkit_mol: rdkit.Chem.rdchem.Mol\n",
        "        Small molecule to register in the force field.\n",
        "    protein_ff: string\n",
        "        Name of the force field.\n",
        "    solvent_ff: string\n",
        "        Name of the solvent force field.\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    forcefield: openmm.app.Forcefield\n",
        "        Forcefield with registered small molecule.\n",
        "    \"\"\"\n",
        "    forcefield = app.ForceField(protein_ff, solvent_ff)\n",
        "\n",
        "    if rdkit_mol is not None:\n",
        "        gaff = GAFFTemplateGenerator(\n",
        "            molecules=Molecule.from_rdkit(rdkit_mol, allow_undefined_stereo=True)\n",
        "        )\n",
        "        forcefield.registerTemplateGenerator(gaff.generator)\n",
        "\n",
        "    return forcefield"
      ],
      "metadata": {
        "id": "FS0G30OR7Gq9"
      },
      "execution_count": 41,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "forcefield = generate_forcefield(rdkit_ligand)"
      ],
      "metadata": {
        "id": "nGROPvgg7JBF"
      },
      "execution_count": 42,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "modeller = app.Modeller(complex_topology, complex_positions)\n",
        "modeller.addSolvent(forcefield, padding=1.0 * unit.nanometers, ionicStrength=0.15 * unit.molar)"
      ],
      "metadata": {
        "id": "Nq3HtWHb7Sce"
      },
      "execution_count": 43,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "system = forcefield.createSystem(modeller.topology, nonbondedMethod=app.PME)\n",
        "integrator = mm.LangevinIntegrator(\n",
        "    300 * unit.kelvin, 1.0 / unit.picoseconds, 2.0 * unit.femtoseconds\n",
        ")\n",
        "simulation = app.Simulation(modeller.topology, system, integrator)\n",
        "simulation.context.setPositions(modeller.positions)"
      ],
      "metadata": {
        "id": "17hZr0PQ7gky"
      },
      "execution_count": 44,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "simulation.minimizeEnergy()\n",
        "with open(\"topology.pdb\", \"w\") as pdb_file:\n",
        "    app.PDBFile.writeFile(\n",
        "        simulation.topology,\n",
        "        simulation.context.getState(getPositions=True, enforcePeriodicBox=True).getPositions(),\n",
        "        file=pdb_file,\n",
        "        keepIds=True,\n",
        "    )"
      ],
      "metadata": {
        "id": "8BblL0fX73kH"
      },
      "execution_count": 45,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# output settings\n",
        "if on_colab:\n",
        "    steps = 500  # corresponds to 100 ps\n",
        "    write_interval = 50  # write every 10 ps\n",
        "    log_interval = 25  # log progress to stdout every 5 ps\n",
        "else:\n",
        "    steps = 100  # corresponds to 20 fs\n",
        "    write_interval = 10  # write every 2 fs\n",
        "    log_interval = 10  # log progress to stdout every 2 fs\n",
        "simulation.reporters.append(\n",
        "    md.reporters.XTCReporter(file=str(\"trajectory.xtc\"), reportInterval=write_interval)\n",
        ")\n",
        "simulation.reporters.append(\n",
        "    app.StateDataReporter(\n",
        "        sys.stdout,\n",
        "        log_interval,\n",
        "        step=True,\n",
        "        potentialEnergy=True,\n",
        "        temperature=True,\n",
        "        progress=True,\n",
        "        remainingTime=True,\n",
        "        speed=True,\n",
        "        totalSteps=steps,\n",
        "        separator=\"\\t\",\n",
        "    )\n",
        ")"
      ],
      "metadata": {
        "id": "F6Lcxnf4A604"
      },
      "execution_count": 55,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "simulation.context.setVelocitiesToTemperature(300 * unit.kelvin)\n",
        "simulation.step(steps)  # perform the simulation"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "iJyno0HEE8TS",
        "outputId": "239eb3af-c820-4b95-9178-0d7f4eb16239"
      },
      "execution_count": 56,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "#\"Progress (%)\"\t\"Step\"\t\"Potential Energy (kJ/mole)\"\t\"Temperature (K)\"\t\"Speed (ns/day)\"\t\"Time Remaining\"\n",
            "275.0%\t1375\t-820069.5917694856\t300.70801789383444\t0\t--\n",
            "280.0%\t1400\t-819983.1464936219\t300.0918382591453\t0.81\t23:56:49\n",
            "285.0%\t1425\t-819263.0368271468\t298.8138801747044\t0.866\t23:56:56\n",
            "290.0%\t1450\t-818919.0744361738\t298.8305427291525\t0.725\t23:56:14\n",
            "295.0%\t1475\t-819748.4363111855\t300.729636904383\t0.647\t23:55:40\n",
            "300.0%\t1500\t-819660.1624241127\t299.66761221053065\t0.687\t23:55:49\n",
            "305.0%\t1525\t-819323.6770758838\t299.71689710004654\t0.683\t23:55:41\n",
            "310.0%\t1550\t-819539.2441378905\t300.3312551362554\t0.657\t23:55:24\n",
            "315.0%\t1575\t-819116.2467783571\t299.9009869911623\t0.683\t23:55:29\n",
            "320.0%\t1600\t-818403.6001041392\t298.94589120878936\t0.706\t23:55:31\n",
            "325.0%\t1625\t-818228.2611690544\t299.15140107215495\t0.68\t23:55:15\n",
            "330.0%\t1650\t-818716.549979551\t300.31481806182217\t0.68\t23:55:08\n",
            "335.0%\t1675\t-818214.4521257979\t298.8141444495686\t0.697\t23:55:09\n",
            "340.0%\t1700\t-818554.5553893867\t299.4424987187032\t0.702\t23:55:05\n",
            "345.0%\t1725\t-819182.7033567722\t300.57025441743275\t0.683\t23:54:51\n",
            "350.0%\t1750\t-818565.5740845099\t299.88551212408834\t0.694\t23:54:49\n",
            "355.0%\t1775\t-819000.0796485256\t301.0857574413104\t0.704\t23:54:48\n",
            "360.0%\t1800\t-818265.8023260854\t299.6084377600189\t0.696\t23:54:38\n",
            "365.0%\t1825\t-818490.5172593745\t300.6368918015818\t0.69\t23:54:29\n",
            "370.0%\t1850\t-818362.2791066584\t300.80767221479937\t0.7\t23:54:27\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "if on_colab:\n",
        "    files.download(\"topology.pdb\")\n",
        "    files.download(\"trajectory.xtc\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 335
        },
        "id": "JOR1aLbSFCXh",
        "outputId": "3d8b5072-1105-4c95-9c13-8eaa65483588"
      },
      "execution_count": 57,
      "outputs": [
        {
          "output_type": "error",
          "ename": "FileNotFoundError",
          "evalue": "ignored",
          "traceback": [
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
            "\u001b[0;32m<ipython-input-57-12ff510b6faf>\u001b[0m in \u001b[0;36m<cell line: 1>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mon_colab\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m     \u001b[0mfiles\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdownload\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"topology.pdb\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m     \u001b[0mfiles\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdownload\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"trajectory.xtc\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m/usr/local/lib/python3.10/dist-packages/google/colab/files.py\u001b[0m in \u001b[0;36mdownload\u001b[0;34m(filename)\u001b[0m\n\u001b[1;32m    223\u001b[0m   \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0m_os\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexists\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    224\u001b[0m     \u001b[0mmsg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'Cannot find file: {}'\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 225\u001b[0;31m     \u001b[0;32mraise\u001b[0m \u001b[0mFileNotFoundError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# pylint: disable=undefined-variable\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    226\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    227\u001b[0m   \u001b[0mcomm_manager\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_IPython\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_ipython\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkernel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcomm_manager\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;31mFileNotFoundError\u001b[0m: Cannot find file: topology.pdb"
          ]
        }
      ]
    }
  ]
}