simulating pdbmols.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3ba2789b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import prody\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f94747cc",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'temp.pdb'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#load pdb using prody:\n",
    "pdbid = '3ET7'\n",
    "chA = prody.parsePDB(pdbid, chain='A',verbose=False)\n",
    "ligand =chA.select('resname 349')\n",
    "prody.writePDB('temp.pdb', ligand)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b53a87c",
   "metadata": {},
   "source": [
    "# Parameterize for simulation in TIP3P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a27e8e57",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning: Unable to load toolkit 'OpenEye Toolkit'. The Open Force Field Toolkit does not require the OpenEye Toolkits, and can use RDKit/AmberTools instead. However, if you have a valid license for the OpenEye Toolkits, consider installing them for faster performance and additional file format support: https://docs.eyesopen.com/toolkits/python/quickstart-python/linuxosx.html OpenEye offers free Toolkit licenses for academics: https://www.eyesopen.com/academic-licensing\n"
     ]
    }
   ],
   "source": [
    "from simtk.openmm import app\n",
    "from simtk import openmm\n",
    "from simtk import unit\n",
    "\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import AllChem\n",
    "from rdkit.Chem.rdPartialCharges import ComputeGasteigerCharges\n",
    "\n",
    "from openff.toolkit.topology import Molecule\n",
    "from openmmforcefields.generators import GAFFTemplateGenerator\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "TEMPERATURE = 298*unit.kelvin\n",
    "FRICTION = 1/unit.picosecond\n",
    "TIMESTEP = 2*unit.femtosecond\n",
    "\n",
    "pdb = Chem.MolFromPDBFile('./temp.pdb')\n",
    "mol = Chem.MolFromSmiles('c1ccc(c(c1)CNc2c(cnc(n2)NC3=CC4=CC(=O)N=C4C=C3)C(F)(F)F)S(=O)(=O)N5CCCC5')\n",
    "pdb = AllChem.AssignBondOrdersFromTemplate(mol, pdb)\n",
    "pdb.SetProp('_Name', 'UNL')\n",
    "\n",
    "pdb = Chem.AddHs(pdb, addCoords=True)\n",
    "ComputeGasteigerCharges(pdb)\n",
    "charges = [a.GetProp('_GasteigerCharge') for a in pdb.GetAtoms()]\n",
    "\n",
    "with open('mol.sdf', 'w') as f:\n",
    "    f.write(Chem.MolToMolBlock(pdb))\n",
    "\n",
    "with open('mol.pdb', 'w') as f:\n",
    "    f.write(Chem.MolToPDBBlock(pdb))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "5e2586ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "# smiles = 'Cc1cccc(c1NC(=O)c2cnc(s2)Nc3cc(nc(n3)C)N4CCN(CC4)CCO)Cl'\n",
    "# dasatinib = Chem.MolFromSmiles(smiles)\n",
    "# dasatinib = Chem.AddHs(dasatinib)\n",
    "# AllChem.EmbedMolecule(dasatinib)\n",
    "\n",
    "# with open('dasatinib.sdf', 'w') as f:\n",
    "#     f.write(Chem.MolToMolBlock(pdb))\n",
    "\n",
    "# Create an OpenFF Molecule object for benzene from SMILES\n",
    "molecule = Molecule.from_file('mol.sdf')\n",
    "\n",
    "#warning - using gasteiger!\n",
    "molecule.assign_partial_charges('gasteiger')\n",
    "\n",
    "\n",
    "# Create the GAFF template generator\n",
    "gaff = GAFFTemplateGenerator(molecules=molecule)\n",
    "# Create an OpenMM ForceField object with AMBER ff14SB and TIP3P with compatible ions\n",
    "from simtk.openmm.app import ForceField\n",
    "forcefield = ForceField( 'amber/tip3p_standard.xml', )\n",
    "# Register the GAFF template generator\n",
    "forcefield.registerTemplateGenerator(gaff.generator)\n",
    "# You can now parameterize an OpenMM Topology object that contains the specified molecule.\n",
    "# forcefield will load the appropriate GAFF parameters when needed, and antechamber\n",
    "# will be used to generate small molecule parameters on the fly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "001e9c35",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAACWCAIAAADCEh9HAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3deVyN6f8/8PdZWp32jRZbQovKUOKkULb5FIUwKMaSwTh8zODDIM2MrJ/RDGOmYfhkGt+RfqKG0IKKjHaSVFJatahTp86ps9y/Py7ONBVaTkW9n3/d5z7XfZ3rfkzzct33fd3XRaMoChBCCHUWvbcbgBBCHzaMUYQQ6hKMUYQQ6hKMUYQQ6hKMUYQQ6hKMUYS6S3FxcX19PdkuKChobGzs3fagboIxilB32b9/f3JyMtneunVrXl5e77YHdROMUYQQ6hJmbzcAob7s0KFDZ8+eBYCkpKTebgvqLhijCL3DxYqKeC4XaDQAODh8uByN1v5jP/vss4kTJwLAqlWrWnz1xdOn5A3CEYqK6w0MZNVa1PMwRhF6h2cCwVI9vXEqKu0/JC4u7tSpUywWS1VVVUtLCwDk5eVblCltbDxnZibLhqJegjGK0LsVNjaqMBgDGAwDBYW3l5RIJH5+fr6+viKRaM6cOYqKimS/rq5uiySVANzhcgFgmJKSfquQRR8QjFGE3i2FxyttajJSUHh7jJaXly9fvvzatWs0Go3D4Rw5ckROTo58dfz4cWmxapFIIJFQACVNTQCggxn6gcMYRejd5mppkYv6mzU192trNxsaKtBbjnKJiYlZtmxZaWmpjo7O2bNnZ82a1boeIUWFVFQElJRYsVgMAA8dnZ5oPepmGKMItVcTRR16/rxCKEzm8fYNG2aipET2i0Sib7/99ptvvpFIJFOnTg0KCtLX129xrATgWlXVjyUlL5qayB6cobLPoOF8owi9XTafryMnp8FkAkAen//Vs2c5fL48jbZWX99z4MDiwsIlS5bEx8czGIxdu3bt3r2bwWC0qCGxru77oqKshgYAGK6ktMnAgK2mlt3QMFJZuRfOB8kaxihCHdMokXxXVPT/KioAwFggiHZ3r6qoMDQ0PHfu3OTJk1sUzszMDCotvaGuDgB68vLr9fVna2nhSy99DP4HRagDxGLx+jVrZr98uZlOp9XV3QsPr6qosLGxSUtLa5GhlZWVmzZtsrKyurJzpzKN5j1o0EVz839hhvZJFEKo3Xbu3AkAysrKNBpNTlt7gLo6ANBotI0bN/L5fFKGx+N9/fXXLBYLAOTk5DZs2FDF4/Vus1G3wot6hNpQWAgDB4KcHEgk8OIFDBoEPB7s2BH9448zKUoMADQajaKonTt36urqbt++vbGx0czMLCgoKDc3d+vWrQUFBQDg7Ox89OhRCwuL3j4b1L0wRhFqw6hR4OkJu3ZBTQ3MmQMTJ0JAQBGXOxag0tR0EZ//V35+vouLS1hYGI1GS0lJ+eSTT7Kzs5lMpkgkAgBbW9vDhw87ODj09nmgnoA3ahBqg7o63LkDubkAABIJ+PuLuNzFAJWWlrOnTNHOz88fPHjw//73PxqNlpmZGRwcXFFRAQAMBkNOTu7nn39OSEjADO0/MEYRatv+/bBpEwCAggJs3pysoJBiZGS0adOCn376UU5O7syZM6Ghoba2tubm5gcPHqyurra3t1dRUREKhcbGxvRWg/NRH4b/sRFqm7U1jBgBISEAAAcPTkhM/Mvf33/Lli0AMGXKFA8PjzVr1iQmJqqpqXl7e6enp8fFxXE4HAD45ZdferflqIfhvVGE/iExERQUYM0a+Osv4HJhyhTQ1IToaKirq7O0tMzPz5eWtLe39/b29vDwkM4/UlpaOmTIEAAoKCgYNGhQr7Qf9TzsjSL0t99+AwcHmDsXSDCqqcG2bcBiwaNHjyZPntzQ0AAAqqqq3t7eaWlpcXFxnp6e0gwFgEGDBv3rX/8SCoWBgYG9dQqo52GMIgQAIBbDf/4DXl4gEMCMGRAZ+Wr/okUSNvvQuHHj0tPTtbW1f/jhh7KysoCAACsrqzbr8fb2BoBffvlFIpH0WONR78KpSRCC2lpYtgzCw4HJhG+/he3bX+3Pz89fsWLF7du3AcDT0/PEiRNkUP1bzJw5cwKbrWNnd7+qyg4ncOof8N4o6u+ys2HuXMjKAm1tCA6GqVNf7T979uznn39eV1enp6d38uRJV1fXdlZ4urT0REnJNA2NQ8OHd1ej0fsEL+pRvxYeHu7tHZ6VBWPHQnIyTJ0KFAW7dvHY7K+WL19eV1fn4eHx6NGj9mcoAMzV1mbSaLdraiqEwu5rOXp/YIyifoqiqH379rm5uSUnL/3ss+fx8TB4MDQ0gL194b59rLt3N6mpGZw9ezY4OJgsptR+WnJyk9XUxBT1Z1VVNzUevVcwRlGfsn79erIhEonIKM42CQQCLy+vXbt2URT1xRdbTpwwUlaGx4/rhg4tuHvXCKDO2vqnjIx7np6enWvGPB0dAAitrMTHTP0BxijqU+Lj48mGRCK5e/dum2WKiors7e2DgoJUVFRCQ0P37t1Lo9GOHg0YM4ZXUTGERsvZseNSSsoeQ0PDTjfDTlXVUEGhpLHxr9raTleCPhT4pB71KSKR6MmTJ3Q6XfiG+5KxsbELFy588eKFiYnJpUuXzMzM+Hy+u7v79evXAU6qq/8YGak9fnwnO6FSNIA52to/Fxdn1tdPVFXtYm3oPYcxivoUHo83e/ZseXl5R0fH5vslEklMTMyePXuys7OrqqocHBwuX76srq4eHx/v6upaU1MDAGy2YkzMOHl52fxPsUBbewCdXtjYeKSwkA6wxchIJtWi9xBe1KM+RVFRXUkpWF/f/ttvvyV7KiuFhw4dGj58xPTp0xMSEqqqqgBg4sSJcnJyHA7HwcGhpqaGTqf7+fnFx8fLKkMBQJXJzObzbVVVl+jpfaKnJ6tq0XsIe6OoL+DxIDoa5s4FJhNmzRrf1HRKTa2poQHWroWgIDld3aCCgmcARgDFuro6fn77ampqdHR09PX1KYpisVgRERH29vbd0TA1BkOLyaTTaN1ROXpP4PB71Bfk5ICFBQQGgqbmjYqKGXFxkJQkSU6OAXCm08HV9WpYGJ2iZsydm2dtfe/Ikc/q6+sBQEFBYd++fcuWLdPT0wOA58+fDx48WIat+qagoLSpSZXBMB8wwBM7pH0XXtSjPsLFBb7/HiZMmAEAFAXJyfQBA5w5HMjOBgbjY2XlWSNG0K9cGXHunFV9fQOLZauiohUdHR0YGEgylKIoFxcXmbfqcwODA8OHY4b2bXhRj/oIJSX48kv46iuYOBFMTGDECNDTg337QEkJsrOhvh5yc0FODtTVxwA8FYuHXbvWwGbjMvFIBjBGUd8xfz6cOQMDBoCuLmhqwsaNsHcvHDwI+fkwahS4usKmTeDqCrNmDVu3TqSk9CgsrLS2tvYto/S7aIq6uq6cXDdVjt4fGKPow8bnw44dsHDhq4/+/mBjA7t2AQAsWgSnT8PDhwAAxsawZg0YGoJEAmIxbN1alJ1tq62traent3XrVnLsrVu3ZNu2yWpqsq0QvZ8wRtEHrKgI3N0hKQnS08HdHQBgxAg4eBD09V8V8PeHzZsBXi+sFBEBTCYkJ8PLl0MnTZpiamp87949IyMjAMBnrajT8BETeu9ERESkpqaS7eDg4FyyPmcrCQlgawtJSWBsDMeOgfTS3NsbpM+KNDWBTFY3ZgyYmsL58wAAYWFQVgZ37tw8deqUzuspQWk0WkenIEGIwBhFPerEiRMvXrwg2wcPHiTLcrSQmJgojc64uLji4uLWZX7/HZycoLQUHBwgIQEsLFoW2LED9u6FkSMhLw/IC0S+vvDDD0BRwGaD9Mn5zZs3pYc030ao/TBGUY+KiYnhcrlk+8qVK42NjW0Wk0gkYrFYLBa3vtYmq30sWwZ8Pnh7Q1QUtJhjvrYWfH1hxQrw9YW6OlBWBhUVAAAVFdi4EfDaHckc3htFPY3P55PR729Zrcjf3//ChQsAkJaW5uHhId1fVVW1Y0fpyZMWCgrw00/w6af/OKqhAU6ehP37gfR3nZ3Bzw9sbP4usHgxLF4s27NBCGMU9bg9e/aoqKgAQE5OzpvKbNmyhaTnxo0bJRJJeXm5rq5uRkaGm5vb8+el06blf/utzsSJf5dvbGz87bfK3bsNysoAABwcYN8+6J7XOxFqCWMU9bTDhw+PHDkSABwcHFp/u379+oEDB0q3GQzGmTNnrl27tnnz5kOHDnG53LFjx545w5e+tCkUCv/v//7P19dXTW1YWVmUrS3s2gUdWfIDoa7Ce6Po/RIfH+/s7GxpaUm2XV1dnz17VlFR8dVXX3G53MWLF8fHx5M33yUSSVBQkJmZ2fLly/Py8sTiiqtXa/76CzMU9TTsjaIe9emnn+bm5t6/f9/e3p7D4Sgrt3wdUyKR6OrqAgB5WK+mpqatrS39ls/nC4VCiqL+/PPPPXv2pKWlAcDo0aN37NixdOlSBoPRg6eC0Cs4wxPqaUuXLj137lxQUNDSpUtbf2tiYuLm5iYWi8vLy6Ojo6uqqsg89jo6OjU1NUKhUEtLa+DAgY8ePQKAYcOG+fj4LFu2DAMU9SLsjaKeVlFRAQA6/xymJBaLs7KykpOTq6urExISkpOTBQIBACgoKDAYDLFYTI4CgOrq6qqqKkNDwy+//PKzzz5TUFDo+VNAqDnsjaKeNnbs2LS0tJSUFHV19fv37ycmJiYmJqakpPB4PGkZOp0+evTovLw8gUDg5+fH4XByc3Nzc3MzMjL279/f2NiYm5trbGzci2eBkBT2RlFPKykpAQBnZ+eXL1823z906FAbGxsVFRU2my0QCIRCYXh4eHR0tL6+/oABA6ysrKysrObPn5+SkhIWFhYfH48xit4TGKOo54hEok2bNtXV1SkqKr58+VJdXd3c3Nze3t7U1JROpz99+jQ5OTkyMvL06dMAoKamNmHCBADQ++ecx05OTmFhYdHR0cuXL++d00DonzBGUQ95+fLlggULbt68qaysvHv37iVLlhQWFu7bt+/UqVNkmTmp4cOH29jY2NjYnD17FlrFqLOzMwBERkZSFEXDNY7QewDvjaKe8PTpUxcXl6ysrEGDBl2+fNnGxgYAbt26NXXqVABQV1cfP348m80eN27chAkTyIAnABg0aFBZWVlxcbG+dOY7AAAwMjIqKip6+PChhYVFQ0PD6tWrz5071/MnhRCBvVHU7e7cuePu7l5RUWFpaRkeHi5dNm78+PEhISE2NjZtLiQnkUgqKytpNFqLZ/oA4OTkFBgYGBUVZWFhIRaLCwoKuv0cEHozjFHUvU6fPr1u3bqmpqbZs2f/8ccfqqqq0q9YLNb8+fOlHymKEovFTOarv8mqqiqRSKSlpSXXbB0OiUSSl5dH3iX9448/zM3N+Xx+T50KQm3DGEXdhaIoX19fX19fAOBwOEePHqXT3/jysUAgWLNmjby8/K+//kr2kGlJtbW14+PjMzMzHz16lJycnJ6ezuPxyMv46enpWVlZIpGoR84GoTejEOoGPB7Pzc0NAJhM5o8//vj2ws+fPx87diwAqKmp5efnk52HDx9u8y928ODBnp6eioqKABAXF1dbWztp0qTuPyGE3gh7o0j2SkpK5syZk5ycrKGhERISMm3atLcUTkhImDdvXllZmbGx8eXLl4cMGQIAx48f37Fjh5qamkgkGj16tJmZ2bhx48zNzS0tLckDqIiICIFAEBUVZWVl1UNnhdCb9HaOo74mNTXVyMhoGECIikqdkxP1+edURcWbCv/++++kX+ng4FBeXk5RlFAo/Pzzz8kf58aNG8kE+K0dP34cANhstlgszszM7K6TQagdMEaRLIWEhCgrKysCPBkwoOrOHYqiqOhoaurU1iVFItH27dtJXHp7ezc1NVEUVVlZSYZAKSoqBgYGvuWHamtr5eTkmEwml8vtnlNBqL0wRpHM+Pv7k4dIBz7+WLxq1d9fODpSJSXNS3K5XFdXVwBgMpnHjx8nOx88eDBs2DAA0NfX/+uvv975c/b29gAQFhYm05NAqMNw2mYEZI5kAKAoqnO3GpuamlasWLF582aKonx8fLZ7e9ObL1asoQHV1dJPubm5dnZ24eHh2trakZGRGzZsAIArV67Y29s/e/bso48+unfvnq2t7Tt/1MnJCQCio6M70WCEZKm3cxz1vjFjxpANiURiaWnZoWMzMjJ8fHyGDBmioKDAYrEuX75MURSVlUXNnPmqhFhMWVhQ9fXk07Vr19TV1QHA0tLy2bNn5EcPHDhAurGLFy9uaGho50/Hx8cDgJmZWYca3E63b98mG0Kh8A65O4HQG2CMImr48OHfv9bOGE1LS9u5c2fzOZZoNFpERMTfJby9qS++oEJDqWXLqMOHyb6oqCgyv/L8+fN5PB5FUXw+f9myZeRwHx+fDjVbKBSSwfyFhYUdOrA9pP+0cLlcNpst8/pRX4IDnhDIycmNGTOmPSUzMzODg4ODg4MfP35M9mhra7u5uZ05cwZeX2UDAHA4IBAAmw1ZWfDFF2BtTXY7OjpOmTJl/Pjxfn5+dDq9uLjYzc0tKSlJRUXlt99+mzt3boeazWQyHR0dw8PDY2JivLy8OnQsQjKEMYpAXl6ePB+n3jBPTX5+/uXLly9cuHDnzh2yR1NTc8GMGV5jx9rNnFmkpnbq1ClDQ8NXb21SFJw/D+XlMHUq7NgBsbFw9So5islkXrt2jbzu2Xy4aFhYmJmZWSda7uTkROYklXmMlpSUuLu7AwC+JYXeCWMUtZSQkGBra8tgMJ4/fx4aGnrhwoW7d++ShNXQ0HBxcfHw8JitpcXcvh3Gj4dff2WlpDAApNOLiIuyGZQY9PWBLEP/+vkVQTL03Llzq1atEggEDg4OISEhrScfaScyaV5UVBQl60nz9PX1Q0NDAaC2tvbjjz+WYc2o78EYRXDv3j2yQaPRvLy87O3tHR0deTxeYmIi2a+uru7u7r5w4UInJ6dXXc6pUwUBAZFPn15IT7dLS/MAgNcxWq0UV3D15SBYpZRVqryILT/DsXm8icXir7766uDBgwDg7e19/Pjx5jOPdJS5ubmenl5JSclPP/20Zs2arlQlFRgY2GKGU4TeoZfvzaL3zLVr1+h0uqamJgAoKyu7uLgEBwcLBAJpgYaGhpCQkGIVFenayPNotP0A2zZsoB48oCjq2TOvpCQoKzuakqKYlEQTCqvIgTU1NREREZMnTwYAeXn5kydPdr21sbGxioqKAwcOBAANDQ1vb++4uLhO1yYSibZs2QIA6urqixcvJjvr6+s3bNjQ9aaiPgxjFP3D6tWrAcDExOTixYt8Pl+6XyAQhIWFeXp6kofjmQAMOp3NZvv7+/9ib78T4NiECRJF+ZIn28rKjubkfFxZGZiUBNevjwoMDORwOOPGjSNDmoYOHaqmpnbz5s2uNzUqKopEuZ2dnYWFhbRnYGFhcejQoeLi4g7V1vyNgGPHjnW9eaj/wBhFf7tw4QIAKCoqpqWlkT1CofDKlSteXl5qamokpGg0mrq6eoypaZW/P0VRlECQpqVlBXDZ1bXi+ymxsXD9+owrV0zPnh1obq7S/LpHQUGB3D/dsmVL15t6+/btAQMGAMCqVavIe/cZGRnbt2+XXo/T6XQ2mx0QEFBXV/fO2rKzs01NTQFAW1tbJhGP+hWMUfTK8+fPybX8iRMnpDvr6+tZLBYJJjMzMy8vL/I4aPyoUZK1a6k5c6h//esLQ0MAWGCgZ2ExxMiIduoUyMmBmhooKICBgcGCBQu+++67u3fvCgQC8hL9/v37u9jUGzduKCkpAcDq1atbzF0iEokiIyM9PT2l9xyUlJQ8PDzCwsJEIlGbtUVERJA3AqysrMgbAQh1CMYooiiKEgqFkyZNAgB3d/cWX3311Vf79u3Lysry8fEhg+cdHR2LiorItyEhITqamtKHSNrajIgIOR8f18DA/2VkZLSoirz32cVL5uvXr5MMXbNmjTRDuVzumTNnamtrpcVqamoCAwOdnZ2lT/ANDQ05HE56enrz2vz9/clJeXh41L9+1QqhDsEYRRRFUTt37iRBU1lZ2frbggJqyZKtAMBgMHx9fUm3jsfjrVy5kiTUjGnTdu10Cg52bmho4PHuP3u2PDNzXGNjwYsXRysqfhWJXk3CRJZEPn36dKfb2WaGUhR16tQpcjuCdDyFQmGzxhccOHBgxIgR0tsL48aN8/f3LyoqIu2h0Wjbt2+XSCSdbhXq5zBGEXXr1i0Gg0Gn02NiYlp/e+kSpalJaWu/sLa2lb5p/ujRI/Lik6Kioj+5SfpP5eU/lZTsS03VSEqClBTFhw9np6X9GBh4xtvbOyYmpnOvb167do1MTurt7d0i9SIiIhwdHaUdTwMDg23btjXvDkskktu3b69atUp6k5eMYFVRUXk1DwBCnYUx2t+Vl5eT5Yv37t3b4is+n9qwgQKgAKg5cyhpPzUwMJDceTQ1NW1xjcznPy4r+29NzZ/Z2dOrqy9XVv7vyZNpSUn0c+fg448hPX1gYeG/w8JOdOIpU0REBMnQtWvXvqnnWFhYeODAAbLgnfR+7oEDB8rKypq1kB8WFjZv3jxVVVUlJaWHDx92tCUItYAx2q9JJJI5c+YAwOTJk1s8gcnNpaysKABKQYH6/nuKBNfLl7Xz5s0jCeXt7d1iNqaqqt+bmoq43BsVFb/U1ydJ9zc1FUVHb50zRzMpCZKS4PhxWL5cq7jYp7Exv53tbE+GNpeUlMThcKQvRzEYDGdn58DAwOZ3P8loU+nSTwh1GsZov3b06FEycL11mrx4QQ0cSI0cSaWkvNpz/z5lYiIcM4atqqp67ty5FuXr61OSkxlpaToiUS3VSmpq6vDhw5cs+Xj+/NF2dvJLl0JSEvz888SZM2cGBQW9/dnO1atXSYZu3ry5Q3cwBQJBSEiIq6ur9O0mGxsb6bfkRdLw8PD2V4hQmzBG+68HDx6QeAoNDW2zQGoqRcZcisWUnx/FZFIAlItLQV5eXquy4sePJyYlQWHhF2+oKnXx4sVcLpfL5V68GLJhw9ynT+e7uk4i6aaiorJixYqYmJjWKy9JM/Tf//53p8/05cuXAQEBbDb766+/lu7897//DQB+fn6drhYhAmO0n+LxeKNHjwaAd77p+OIFNWsWBUDRaBSHQzU2tlGmvPynpCRITx8kfSjfQmpqqqenJ9m+fv06uTfaekySgYEBh8ORDv6/cuWKgoJCFzO0ueY3Lk6fPg0AS5YskUnNqD/DGO13yHzJK1asAABzc/O3zzYfFUUNGkQBUDo61NWrbZcRCivT0rSSkuDly5A31fP48eNt27aR7fj4+BYj8HNycvbs2TN8+HDpoyFbW9vPPvuMZKhM3npq7f79+wAgnZ4ZoU7DGO1HmpqaAgICdHR0fH19yVilFs/Z09OpW7debUdGUtu2UTQaBUBNn041e9bdEpmLJDt75htLtBt5NKStrU3ePqLRaFu3bu16tW2qr6+n0+lycnKNbXawEWo3jNF+QSKRnD9/XtrdI5M0KysrT5s2bc+ek4GBVFoa1dREHTlC6ehQ5AWlhQup3bspJpPy8aHesFY8RVFUWlpNaurglBQlgSBXVq3l8/lBQUFkKCvpO3cTMiYfxzyhLsKVQfu+xMRER0fHRYsW5eXlmZqanjlzRk1NbcCAAQ0NDTExMTEx9OXLwdoaWCxITYXly2HLllcHenpCejrs3Qv0N/yZCIWwZIna9OmZjY1hCgrGbRfqOEVFxaVLl1pYWEgkkvT0dLLz+fPnaWlpsvoJgrxBkJGRIdtqUX+DMdqnFRaCl9fvW7fGxcVpa2sfPHjQ3d19/fr1ly5dGjx4cGFh4Z9//rloEdvDA0xMQCQCdXWwsQEW69WqHzQaSJf2qKiokEgkZLusrIxs/Pe/kJkJGhoDrK2dZd728ePHA0BSUhIAxMbGDhkyxNvbW7Y/QWL04cOHsq0W9Tu93R1G3aO6mvryS0pBgQKo0tffu2dPQEAAGXBOo9E8PDyePn3a4oi6OurQIer8eaqigpowgXJzo3Jy/v52+vTppaWlZNva2lokEhUUUCwWBUBFRnbLGZw4cQIAvLy8KIqqq6uj0+ny8vLNJ5DuuuDgYABwdXWVYZ2oH8LeaJ8jFMIvv8CoUXDkCDQ1gYeHpp+fQXb22rVry8rK2Gx2QkJCcHBw88fiBIv16uJdWxtWroSIiHf8zsaNwOPBJ5+As+x7ogD/7I2yWKxRo0Y1NTXJ9gIce6NIJnAtpr5CLIaqKtDVhfBwWLsWAGDKFFi3Dk6fhhUrvOTlL9jbr9m0af78+W9Z+m3UKNDSAgBYvRrS0uD1RKOvfPPNN+RV+pKSkmvXICwMVFXhyJHuOiErKysFBYWsrKy6ujoVFZXx48c/fvw4MTFx3LhxsvoJExMTJSWlgoKC2tpaMqs/Qp2AvdE+4eefwdkZdu0Ce3vQ1YWlS+G338DSEpYsgevXQUND/uuvr0dGLliw4O3LZ7q4wMSJAAB0Opw4AQMH/uNbT0/PtWvXrl27Vltb298fAODbb0Ffv7vOSV5enjxlSk1NBQCSnsnJyTL8CQaDMXr0aIqiMjMzZVgt6m+wN/rhy8uDc+fg1i2g0+HFC5g5EyIiYNQoqKsDOTn4/HPYvftVJ7NruNyhWloDTUxAXl5+0yZwd3/V6+0+NjY2ycnJSUlJDg4Oza/xZWjMmDGpqakPHz60s7OTbc2o/8AY/fAlJsLUqa/ua+rpgZoa0GgwcybU1IC/P5iby+RHWCzW2bP0rCy4dw9UVFR8fGgy7Re2jfRASXSOHTuWyWRmZGTw+Xwyc7NMkLXw8PYo6gq8qP/wMRjweigSAIBYDEwm/P47REbKKkMB4OLFixIJq7z8kpXV/2JjY3vmL6d5D1RZWdnU1FQkEklHksoEDh1FXYcx+uGzs4OoKBAKAQAKC6GxEbS1QV4eAFJSUvLz80mpu3fvSsd7dg6dTi8q8snJMS0o6GqT28nCwkJJSSk3N7empgYAbGxsQBbX9aQ2Ah/Wo67DGP3wGRrChg0wfTosXQrLl8Mvv0i/CQ8Plz6T+f3339XcODkAAAyuSURBVHNycrryO/Ly8jSaRCL5YuvWLrW3/ZhMppWVFUVR5CyaX+N32vXr14cPHx4ZGUk+GhgYaGhoVFZWXrp0SdK8U49Qu2GM9gnLlsGtW3DyJMTEwNixMq/++XMICgI6na6kpCSR3FFWFublyfxH2tb8ur7rT5lu3Ljh5uZWXV198+ZNsuf69ev19fU6Ojru7u6DBw/etGlTfHw8RVGyaDvqLzBG+5DXK7M3d+jQIXd3d3d394h3jqd/g6QksLOD5cuBwQAVFQYArFv3wt6+Sy1tPxKdpDdqZWUlLy//+PHjurq6TlQVHx8/b948gUCwceNGPz8/ALh+/frcuXObmpqMjY2HDRtWXFz8ww8/TJ482djYeMeOHTJ/hR/1Wb38FhXqTnv37g0JeTUH6Pr162NjY5uamjpUQ2gopaxMAVBOTlR1NaWnZwEAa9akv/tIGSEPf4YNG0Y+fvTRRwAQGxvb0Xru3LnDYrEAYOXKlWQlkri4uAEDBgDAxo0bpb/l4+PT/P2uoUOHcjic5ORkiqJiYmKkC6Nevnw5Rbq4Cur3sDfaj6Smpo4aNeru3bvtLP/99zB/PjQ0wKefQkQEqKuDoiILAOTkyruzmf9gamqqoqLy7NmziooKAFi2bNmXX36pq6vboUru3bs3a9YsHo+3YsWKkydP0mi0u3fvzp49u76+ftWqVd9//z0pZm5uvnfv3pycnNjY2A0bNujp6eXn5//www/jxo2zsrK6dOlSVlYWKZmcnCx9docQ9kb7spCQkL/++otsnz592tHREQDk5eX/+9//vn1tOJHo1dLKNBrl4/P3/qFDhwLAmjVrurPVLU2ePBkAIiIiOnd4SkqKhoYGACxcuJAsIpKQkKCiogIAK1asaL36k5RYLI6Li+NwOCS1//Of/xw/fryhoaGhoWHnzp0XL17s5PmgPgdjtB8RCoU+Pj50Oh0AnJ2dy94woz2Xy507d4mpab2SEhUc/Pf+58+fq6ury8nJjRw5sifXJSZrz33zzTedODYtLU1LSwsA5s+fLxQKqbZS9Z2ampqioqKOHTtma2v7ySeffPLJJxYWFhijSApjtN8JCwsjyaKnp3fjxo0W3xYUFJChlGZmTgkJf/dYL126pKmpCQBkfSQWi3X48GESTN2NrD1nZWVVWFjYoQMfP36sp6cHAO7u7uSmsDRVFyxY0NHGHzt2LCAggGzv2bMHYxRJYYz2R4WFheRKmcFg+Pj4SDtlaWlphoaGAGBiYpKdnU12CoXC7du3kzlNXF1dHz9+7OnpSe4IjRkz5s6dO93a1JKSEmtrazL9Ep1OZ7PZAQEBdWTd57d68uTJoEGDAGDWrFlkltIHDx6QVZ6kqdohGKPoTTBG+5To6Gjp4kVXrlx5y0WrSCTy8fFhMBgAMGXKlOLi4tDQUDIPnpOTU3V1NSmWn59P5uxgMpk+Pj7SO4kxMTFkfWYajebp6VleXt4dp5OammpkZAQAhoaGs2fPJgvWA4CqquqqVatiY2Pfcod37dq1zTP0yZMnZNZq6Z6OSktLe/DgAdlOSEjIzZXZ2lPoQ4cx2qe4uLhI71pOmDChvr7+7eUjIyNJuKioqJB7pitXrpT21C5evEhuIw4ZMuTu3bstjuXz+T4+PiTaNDU1/f393/K4phOuXr1KOqGTJk168eIF1da69kZGRtu3b89pPk3/a01NTX5+fmT56OzsbH19fQCYMWMGn8+XYSMRojBG+5iOxihFUeXl5Ww2m6xm7OzsTHYKBAIOh0Oiys3N7eXLl286PCcnZ8aMGaTk5MmTMzIymodamwHXHv7+/iTWFy9e3Dr4njx54uPjM2zYMOmAk3HjxgUEBNTW1rauKj8/f8iQIaSXTVIVIdnCGO1TXFxcZs6c6ebm5ubmpqGh0Z4YpShq6dKl0jzy8vJ68OCBtbU1eZrk7+/fnhqCgoLIwxwVFRUzMzPpfgsLi46eglAoXLduHbld4OPj85bLdjIgydvbm/V6mn5FRUUPD4+wsDDp46Pnz5+TtGWz2e25o4pQJ2CM9ikteqMPHz5sfTHeQlFRkby8PIPBOHXqFLmIJjdMR44cmZqa2v6frq6uXrdu3ddffz1mzBjpzo7GaFVV1dSpU0kgnjt3rp1HNTQ0BAcHN7/YNzAw4HA4N27cMDY2BoCJEye22VFFSCYwRvuUFjE6c+ZM8giI3Fts07Zt2wBg0aJFfD5/yZIlJIaUlZVPnz5NCiQnJ3do6PvQoUO9X3tnjGZlZUnvqMbGxo4YMQIA9PX1k5KS2v+LUnl5eb6+viQ6SX8WAGxtbblcbidqQ6idMEb7lH379kkfmm/btm379u1kmKeGhsaPP/7Y+sE9j8cjo0HPnz9vaWlJLuRNTEy2bdvG4XBImfPnz3do6LuFhcWL194Zo3Z2dtKhBWPGjHFzc7OysiooKGj/z7UmkUji4uJWrlyppqZGp9O7WBtC74Qx2sfl5OR8/PHHpHdmbW1NZoGT8vf3B4BRo0aRSTpGjx6dnp4uFotjYmI4HI5YLBaLxX/88UeHYrRDF/UtYpTH40k/dh2bzQaA8PBwWVWIUJtwLaY+bsSIEVeuXAkPD+dwOGlpaZMnT162bNmRI0d0dXXFYjGZlePJkycA4Onp+dNPP5E8BYArV64UFRUBQHFxsYuLS/t/kQx6b739JmvXrmUymQBQXl4u/XWZmDZt2p07d27evNmh9iPUYb2d46iHNDQ0SId5qqur+/v7H3m9xrySktLJkyebFya9UbLd0Yv6DrGzsyspKampqampqTE3N5dt5TExMQBgbW0t22oRagF7o/2FkpLS3r17ly5dyuFwrl27tnnzZvIExsDAIDo6etSoUb3VMFVVVdIJJQNFZWjixIlKSkoPHjyorKwkr4Ei1B1wvtH+xcTEJCIiIiwsjMFgUBTFZDLJJKQtihkbG7u6upJtS0vLadOmdVN7jIyMpOlJBsnLkKKiop2dnUQiiY2NlW3NCDWHMdofkRujAGBqaqqjo9O6wODBg52dncn26NGjJ02a1E0tCQ4Oli46Hx4eLvP6yShU6cpLCHUHjNH+6MSJE2Sjc4safUBIP5rcJEWom9AoXASxnxGJRGpqag0NDUwmUyQSlZaWktlJ+iSRSKSpqcnj8UpKSvrwaaLehb3RfufChQsNDQ3wetn3hISE3m5RN2IymWw2m6IovK5H3QdjtN/Zv38/AIwcOXLWrFkAcOfOnd5uUffC26Oou+GAp/6lqqqKLFm8Z88eslJbn49RZ2dne3t7Kyur3m4I6rPw3mh/kZeXFxkZeezYsUePHikpKdXX1/N4PA0NDTqdXlNTQ+a975PWr1/v7e1Npv7z8PA4evQoWSgFIVnB3mhfxuPx7t27FxUVFR4enpmZSXYqKyuzWCw+n0+j0by8vAYPHiwUCrlcLoPBkE7c2Zfw+XwyugsAGhoaJBJJ77YH9T0Yox+2uXPnBgYGqqurA4Cjo+Pt27ebmpoSEhIiIyOjoqKSkpKkCaKlpTVt2jQnJ6fjx49nZGTs3bvXwsLCzMzsyy+/BADylr2Xl1dvnky3+fPPP8mtjOLi4t5uC+qDMEY/bDU1NdLuFZk4+bfffquvryd7FBQUHBwcpk+fPn369I8++oi8L2RnZ2djY/Pdd9/t3r1btlOBvLdYLJaamhoAkDlQEJIt/Kv64F28eJGkYW1tLQDU19cPHz7c2dnZ2dl55syZZEL75qysrDZv3nz48OFff/11yJAhpIOWnJy8evXqnm98z5gyZQoZ3XXy5MnebgvqgzBGP3gMBoMs+wEAO3bs8PHxeec4c19f39DQ0Nzc3BEjRmzatAkwXxDqAnxS/2FzdHQMDQ0lM9hbWlo+ePCgnQfeunVr2rRpDAbjwYMHpqamffjeaHV1NYvFkpOTA4Cqqip1dXXpvzoIyQQOv++npkyZMmXKFJFItG7dOoqilJWVyVSkfY+GhgbJUADQ0tLCDEUyh73RD1tpaamenh55dlRcXGxgYND+Y7lcrrm5eXFxcUBAgLe3d7e1EaE+DmO0X7tw4cLChQtVVVUfPXqEg9IR6hy8qO/XPDw83Nzchg0bRp7yI4Q6AWO0vzt9+vSSJUtKS0vJxx07dqSlpfVukxD6sGCM9ncaGhpcLpdMnQcAlZWVAoGgd5uE0IcFx40iAIDExESywl1hYWFvtwWhDwzGKAIAqKmpefHiBQBIu6UIoXbCGEUAANOnTydLgd67d6+324LQBwbvjSKEUJfguFEEZWVlysrKZBKToqIiTU3NPjyLM0IyhzGKEEJdghf1CCHUJRijCCHUJRijCCHUJRijCCHUJRijCCHUJf8foB2uiC1aVC0AAAQAelRYdHJka2l0UEtMIHJka2l0IDIwMjEuMDMuMwAAeJzNkn1M1HUcx7/3u0eOA47jOPCOg4MTPSJEzYfUVbw/OXKWZlarqNhVDm9uiYsGCF5wHtzxoAYYytMdGEfyYApkiMhQBkED21AzZg4f0tyaNlxQc01X38H4Q5etP/tun732+e79efx+J/u+vML4UXETsdmzmtuz3PJFCmbhFAtqZuKUzbnyGVcsm4Hkn5E0I/3vnA2Tz1YQz13Lk/7dfyTcbwaih/FozsfxsTr2MOe2IJqhSPR/p7+IKQRBpBDEIiYRM6mEyaRMLmNyMeMP7KewCUo/q7/SJqj8rQEqmxCgsAYGsKBApg4yBattgibYGqKxCdoQa6jWJoSqrbpQpgtj4To2T8f0GpNBb4owWI0RNiHSaI2KtAimSGu0yRptsAkxJotgjjKZY0zzzRZBEsskC5h0IZNamCyOyZ5g8ngmf5L5JTDlIuafyFSLWeASFriUBT3FtMuYfjmLWMGMK1n00yxmFdNKpRKxXCaV+Sn9VQEKWbAmRBuqlkUYI03RBmlUpCnGHCaa+8ar77ccQe+yAoqVHUO6pg2xrt1kvtyEhM4v8MsCJyV01mLj+z44KvOoansrXlrpw/4qOzUPfc31Nfjt/E5ynunCxewqdHqyaPjFNtxoqkPd2l3UoPfhG1M37r3qoJQtjVie2oM/LHbSHm2CzDGI5uRcmifxYfrBKDqNeRRV2IAG/QiKsZvuvOJFctdpJK9wUlZfPc/dyzUu3psH+8OHsTnbTd7JWrR/NIKGlCJac7UcyvFBHuOmYp7BmPItEj1FdKMpDzXVY3hmfQnPuQt3D/+Agd9LqNPjgm/HBaTdKqa4d8p43CXEu0uo/NdqpPacR/xwKe/tEO95AvbniigtrR5hl35C/0gpn+sAr9WP6LtuyrClw607iQ8Xuqj34wxeqxtrRwtoicLJZ2nHepWTzsUUoGvqMPYe3M17y0HihBeld/Lp6su5XN+MzfeclK7ZCu9kO673FFLh7Q/4buuxM85J0w/exIUN1TAkOvjsW/k+y9GvzOea93B6Xxscr+dSvPsEtn3XiLey8ulg3imEV9Zi095MWrOjF+o/K1A3lUv6MydhX+TB2FQGXVa2Ytxchml7Jo1XtOCvzga8UZtDNQoPKn+uxM2BT+hi9iEUBBzHrbwcOtXdjMWtQ8gMyqJ1m5qgT/oet7+y01K+k8aKMaz6zEmjZbXY8/YJDG1wkXL8AFrmd+PafRc5Kn04EnEWjkE3Jb77ObYfn0DR83vIe6UQhsgBdLQX0dEXrLhW1wdLUCF5U4vx6c1j6LjuoPiqUsg3dkD9o4u2nH0NvdN1cJ5z0rbeddD9DcRhbEqaMZmvAAAEK3pUWHRNT0wgcmRraXQgMjAyMS4wMy4zAAB4nJVXTY41NQzczyn6BFH8F8dr5kNIwCCB4A7sub+wE6e7FyB9Hj2Nuvqp65Xtsp3+8+uXjyv+fv/8+e9/rvuPPj8+LpnXwOvq//kxs+sv7L3H8whNJvnta4wGfV35vamzX79d/0fx/iQLTV4s0nDgYulNya/++H6W3iZpsoiOuAJrHHwFLb0Jz8VCbdrIiKi7vq8SC2qykNliwaboWn74bhZXj7KeHdhYRrJM5iJLTxZo1jkjYqhqkQkZUR+aNRKgCgv6s7DiGNzG4GRhLEXkLLgrM+KqH79AL7Fwk0GZ3Ql2IuqlvKA03vV1LaqHpQ8psmBqkQa4/OcZV6tFRM3G6aPRTwfYmEUW0e1ddefcWqiWF25zHhYlSRY2LPWReOftSs8GWxXMhlTzi3utS2qZu49AG4gWtXTm1EJ9a5FmWsvuaF23FmvAkiweWknL8FmHyeLjarH4vZg5JS1MPSMy6ZndzjXvakPg1ILbfzG5GIpaiCFZ5jx9BBHRjxUtQ+H4BY8WQymxeH23S0LLnjSeF40OLbCw+7Sn62ysiID9ikqVpgY0kqXrcR3GFK3N3XSd59k0/UIxrwos6NP7TCnJ7A7fATXX9WZ2b/vMrs+rml+8qnPPJt8kvM8M0Ue91tNw58U7eUrWiLk267ABaGZ3jsOCtT6KiO5TUO5a8Fk8qnsacpP4fjM8EVktu9492Y3cOkKy+JGheK5jZs0NSyQ7VxbMP1X8wiq0tSARbD8zKRVYYvewzmudggYt58Sm88NehcUngyLsiGyq5WkEZonF3O268uK9MCfsXKFIKSKf1KCSEeE5e4MNrWW3T+Ud0fBWTn1oWNOipCsv1IYB5fQOPxe0oJ++xkgtyJYdar2iJU4e2yWeF/WQthY/bZZYfK7RisMjYhQ9m8RVlVi6wtg9rdPkOFFKNXKvgeneASBgO+NTpJTdOHsvxzoLG9nOiyiVWHw25bb3s9SkuVnQRCosfiaEVSN1l9COaLiHqJQX9vfGnRfvS16V9im1uqKUl4nbdV5zX0jXPh0qlPzi71Zzs/iO5z43C3uia5X26Q1ZI3/52zsA46W7Nl9geS22WvdO9nveWzwKlfaHYP0PQOsqAZ8vA8jCCcbCCXThBHPhBzzPOJUdEHvG1pcBYOEEGPgACnwAB37Ai0AWTuoR+AANfMAMfJ6xwAl8ecLN5lMfbzmIgR/wECAFTuoAtwKUwAeMN9D1Owlm4AMscFJ7n+Itx01CtxzCwFkfogXyGw78gBebrG/ymREc+aMB7lR5Qzo+pfdn5qv0ZK/Sx1B4ALx8wPgG9PIB8xvIywf81AcuvhPidedbgdedH+/IJf0N4FV6eUoCl7wTL/zKtTxh90tuBZ43uRV4auQo+PW6vn19fvwLqZILDNSVuv0AAASeelRYdFNNSUxFUyByZGtpdCAyMDIxLjAzLjMAAHicjVZNaxxHFPwrOUqwGfp9dvcsOgmMT7YhR+ODMMEJJHYwPvrHp6pHlmd2OpBF2pVeV1e/V6/6zb5//eFRHj7qx7v3rz/cf7x7s33a5x+Bx7tX9/h59bL2ON63Vf+BOr37b3cPb+/5+8Z3O/7jTx8f/P1sRx59eIP83v7y/U76ktLjkrH00vUiuoTnFeHWnNFochFbIvWqspS6QVtcxBd3BHURqZfMpeUIqscIpskl6+IanQTae7+qLWLJcKn1IrGo5FWxxxjr2UnQ0xiz1kHbFvMQhotUEMdS3Ee4DIJe21Vzy6sv4sFgrYw1U8a8XQT/KGJ1qYKsakEW0hmuZgPrRk7oIA0HBKEqzu2afoEa6kKcuTDYWr8oWDbSrMLdogPZFduRX+RA9iRnZU2BpHMgWXxD1c5CW6usvlpwu3clMnh6xZo1nqQ9gupFHdgSjdhqQ/8ufUQlWFRb2oaFYH10oNvINqox3HP0NZEKs83G4zSNQdGNt2+aSjJYW98y86jEegetQeTBWjIIDhNjFO2umwreWVyQWJcmWwpFR8VGbwhUQWsRzqSR6JEkhw5sFdYWJZMECb4LDdiLbsUh7NxJrKvTnbZZVlMHgQ3LRjO6U4pxf6NfCjKow9x1iOA0N7KDVXA8W0ay4iOolUFD0Uin6uB09zrKMjSNHaA7kEYNY1bK+pCkGySneXAFgHZYAuDmfh1y1kYRkh1GZ3H6Fc2zqkJwb9QXYkhrhBdyCC6pMzkXZTCdfigA45YirBFGEt6Vwa2jaOlZCeddRi0l62iGGJMurfpQOGR0LhUOpGWtDnR2dJZmb8OtWWIYE7UILVdHN3gcbgd8iEtlddz0pmwdrMG2Eltl2LXI8FrwumBVKhXtLUdDu5dGsMO1zK73raWkgGmhM0cT7A8s2jZuexnFSQ7tg0MMa11YhPC+YnboGCGi95enb1/+fvf1yz9rWd49ff3259Nfj388ff30+0IrN8hbfr7qKmeQoG991fNCqZl1t7voajeoXwEzwexbfXK698Yp+fJa4wwyM+D2oDyfoc0CF+XlJWs9MUlEL7bHtBNGYRzc3N1hfXKYYkQeiGQiraf7ngnSTrRtaFvsFPRVzkKLeLFc5VZcjoFSuXIrLrKUjC5+yPJW3TlqIu8ENdHXwmKVW1Gx3fAsKflTirbKrazY7pj8hz7rVNUeuz4jE52oGqqGlYlln1fOSpZuqnI8f6IqHpXmWJoYNQ1fJ3b7c9WJlCWl6P7K6apnLUuqMoOJmAWPFNUdga16VhMEGGcHoWx2/7VXXNqJhDVE8yCHneU0N0yzw/2fDADBzax5oJooiyHjuTNAW23i1wlqJnIVVH/ohU0MW7Qx5fMYKPjmoHFUb6LxGeW3GjMZjRZ2kNIngmvDc3D1iWefV2ajdUI9mbPPBBPPPq/cSviijU+sOan6fynok1GKC4fZsMZMNkNp+/vU1pjIhkl8GLO5xq2Eg8z4DIvJvU/H94A1JrJVfLvZl6BrTBypmBw19oXGWU6txa0dQLfKkiqr+w6l3/8FH8lHhQPbcoUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<rdkit.Chem.rdchem.Mol at 0x121b8d220>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "molecule.to_rdkit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3a70e242",
   "metadata": {},
   "outputs": [],
   "source": [
    "top = molecule.to_topology().to_openmm()\n",
    "pos = pdb.GetConformer(0).GetPositions()/10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5a22e7c0",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "#and add salt water:\n",
    "modeller = app.Modeller(top, pos)\n",
    "modeller.addSolvent(forcefield, \n",
    "                    ionicStrength = 0.15*unit.molar,\n",
    "                   padding=1*unit.nanometer)\n",
    "app.PDBFile.writeFile(modeller.topology, modeller.positions, open('system.pdb', 'w'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "472dc2c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "system =forcefield.createSystem(modeller.topology,\n",
    "                                nonbondedCutoff=0.9*unit.nanometer,\n",
    "                               constraints=app.HBonds,\n",
    "                               hydrogenMass=4*unit.amu,\n",
    "                               rigidWater=True,\n",
    "                               nonbondedMethod=app.PME)\n",
    "#'constraints' : app.HBonds, 'rigidWater' : True, 'removeCMMotion' : False, 'hydrogenMass' : 4*unit.amu \n",
    "integrator = openmm.LangevinIntegrator(TEMPERATURE, FRICTION, TIMESTEP)\n",
    "platform = openmm.Platform.getPlatformByName('OpenCL')\n",
    "simulation = app.Simulation(modeller.topology, system, integrator, platform)\n",
    "simulation.context.setPositions(modeller.positions)\n",
    "simulation.minimizeEnergy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ecd024a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "\n",
    "simulation.reporters.append(app.DCDReporter('output.dcd', 500))\n",
    "simulation.reporters.append(app.StateDataReporter(sys.stdout, 5000, \n",
    "                                              step=True, \n",
    "                                              potentialEnergy=True, \n",
    "                                              temperature=True,\n",
    "                                             speed=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "0b9db12c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#\"Step\",\"Potential Energy (kJ/mole)\",\"Temperature (K)\",\"Speed (ns/day)\"\n",
      "5000,-34554.1703483302,295.72005133727146,0\n",
      "10000,-34134.8734733302,295.1417582956899,69.7\n",
      "15000,-34450.9984733302,293.1647432752746,70.1\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m/var/folders/jh/02165y2n7kq2y5ychxtzcjm40000gn/T/ipykernel_9068/2982974720.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0msimulation\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstep\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1e6\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/miniconda3/envs/sci/lib/python3.9/site-packages/simtk/openmm/app/simulation.py\u001b[0m in \u001b[0;36mstep\u001b[0;34m(self, steps)\u001b[0m\n\u001b[1;32m    130\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mstep\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msteps\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    131\u001b[0m         \u001b[0;34m\"\"\"Advance the simulation by integrating a specified number of time steps.\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 132\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_simulate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mendStep\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcurrentStep\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0msteps\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    133\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    134\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mrunForClockTime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcheckpointFile\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstateFile\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcheckpointInterval\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\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[0;32m~/miniconda3/envs/sci/lib/python3.9/site-packages/simtk/openmm/app/simulation.py\u001b[0m in \u001b[0;36m_simulate\u001b[0;34m(self, endStep, endTime)\u001b[0m\n\u001b[1;32m    195\u001b[0m             \u001b[0mstepsToGo\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnextSteps\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    196\u001b[0m             \u001b[0;32mwhile\u001b[0m \u001b[0mstepsToGo\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m10\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 197\u001b[0;31m                 \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mintegrator\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstep\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# Only take 10 steps at a time, to give Python more chances to respond to a control-c.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    198\u001b[0m                 \u001b[0mstepsToGo\u001b[0m \u001b[0;34m-=\u001b[0m \u001b[0;36m10\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    199\u001b[0m                 \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcurrentStep\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m10\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/miniconda3/envs/sci/lib/python3.9/site-packages/simtk/openmm/openmm.py\u001b[0m in \u001b[0;36mstep\u001b[0;34m(self, steps)\u001b[0m\n\u001b[1;32m  12922\u001b[0m             \u001b[0mthe\u001b[0m \u001b[0mnumber\u001b[0m \u001b[0mof\u001b[0m \u001b[0mtime\u001b[0m \u001b[0msteps\u001b[0m \u001b[0mto\u001b[0m \u001b[0mtake\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m  12923\u001b[0m         \"\"\"\n\u001b[0;32m> 12924\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0m_openmm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mLangevinIntegrator_step\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msteps\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  12925\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m  12926\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__init__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\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[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "simulation.step(1e6)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}