minimization-with-nagl.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Minimization of multiple conformers against single target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "315aa6402c3b4482a5895f07cb1a7891",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": []
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from openff.toolkit import Molecule, Topology, ForceField\n",
    "from openff.toolkit.utils.nagl_wrapper import NAGLToolkitWrapper\n",
    "from openff.interchange import Interchange\n",
    "from openmm import LangevinMiddleIntegrator, Platform\n",
    "from openmm.unit import *\n",
    "from openmm.app import PDBFile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in input files:\n",
    "- conformers: SDF with multiple conformers of the same ligann\n",
    "- ligand: reference pose which is used to relax the protein initially\n",
    "- receptor: protein without bound ligand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "conformers = Molecule.from_file(\"conformers.sdf\")\n",
    "ligand_ref = Molecule.from_file(\"ligand.sdf\")\n",
    "receptor_path = 'receptor.pdb'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build protein topology"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "top = Topology.from_pdb(receptor_path)\n",
    "sage_ff14sb = ForceField(\"openff-2.2.0.offxml\", \"ff14sb_off_impropers_0.0.4.offxml\")\n",
    "protein_intrcg = Interchange.from_smirnoff(force_field=sage_ff14sb,topology=top)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build ligand topology using NAGL charges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "ligand_topology = ligand_ref.to_topology()\n",
    "NAGLToolkitWrapper().assign_partial_charges(ligand_ref, 'openff-gnn-am1bcc-0.1.0-rc.3.pt')\n",
    "ligand_intrcg = sage_ff14sb.create_interchange(ligand_topology, charge_from_molecules=[ligand_ref])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combine protein+ligand topologies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "env: INTERCHANGE_EXPERIMENTAL=1\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/homezvol3/ekubinco/miniforge3/envs/ssbind/lib/python3.11/site-packages/openff/interchange/components/interchange.py:950: UserWarning: Interchange object combination is experimental and likely to produce strange results. Any workflow using this method is not guaranteed to be suitable for production. Use with extreme caution and thoroughly validate results!\n",
      "  return _combine(self, other)\n",
      "/data/homezvol3/ekubinco/miniforge3/envs/ssbind/lib/python3.11/site-packages/openff/interchange/operations/_combine.py:25: InterchangeDeprecationWarning: The `handlers` attribute is deprecated. Use `collections` instead.\n",
      "  \"vdW\" in interchange1.handlers\n",
      "/data/homezvol3/ekubinco/miniforge3/envs/ssbind/lib/python3.11/site-packages/openff/interchange/operations/_combine.py:26: InterchangeDeprecationWarning: The `handlers` attribute is deprecated. Use `collections` instead.\n",
      "  and \"vdW\" in interchange2.handlers\n",
      "/data/homezvol3/ekubinco/miniforge3/envs/ssbind/lib/python3.11/site-packages/openff/interchange/operations/_combine.py:27: InterchangeDeprecationWarning: The `handlers` attribute is deprecated. Use `collections` instead.\n",
      "  and \"Electrostatics\" in interchange1.handlers\n",
      "/data/homezvol3/ekubinco/miniforge3/envs/ssbind/lib/python3.11/site-packages/openff/interchange/operations/_combine.py:28: InterchangeDeprecationWarning: The `handlers` attribute is deprecated. Use `collections` instead.\n",
      "  and \"Electrostatics\" in interchange2.handlers\n"
     ]
    }
   ],
   "source": [
    "%env INTERCHANGE_EXPERIMENTAL=1\n",
    "docked_intrcg = protein_intrcg.combine(ligand_intrcg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize the complex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "05456d9131024724a70e3649bb98a6ff",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "NGLWidget()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "w = docked_intrcg.visualize()\n",
    "w.clear_representations()\n",
    "w.add_representation(\n",
    "    \"licorice\",\n",
    "    radius=0.1,\n",
    "    selection=[*range(protein_intrcg.topology.n_atoms)],\n",
    ")\n",
    "w.add_representation(\n",
    "    \"spacefill\",\n",
    "    selection=[*range(protein_intrcg.topology.n_atoms, docked_intrcg.topology.n_atoms)],\n",
    ")\n",
    "w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up OpenMM simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)\n",
    "platform = Platform.getPlatformByName(\"CPU\")\n",
    "openmm_simulation = docked_intrcg.to_openmm_simulation(integrator, platform=platform)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the PDB does not have a CRYST1 entry so default periodic box vectors are used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Quantity(value=[Vec3(x=2.0, y=0.0, z=0.0), Vec3(x=0.0, y=2.0, z=0.0), Vec3(x=0.0, y=0.0, z=2.0)], unit=nanometer)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "openmm_simulation.context.getState().getPeriodicBoxVectors()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check platform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'CPU'"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "openmm_simulation.context.getPlatform().getName()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Relax protein against reference ligand\n",
    "This speeds up subsequent minimization of conformers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "openmm_simulation.minimizeEnergy()\n",
    "\n",
    "# save protein positions for reference\n",
    "new_positions = openmm_simulation.context.getState(getPositions=True).getPositions()\n",
    "protein_pos = new_positions[0:protein_intrcg.topology.n_atoms]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Minimize conformers and structure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i,conformer in enumerate(conformers):\n",
    "    new_positions = openmm_simulation.context.getState(getPositions=True).getPositions()\n",
    "    new_positions[0:protein_intrcg.topology.n_atoms] = protein_pos\n",
    "    new_positions[protein_intrcg.topology.n_atoms:docked_intrcg.topology.n_atoms] = conformer.to_topology().get_positions().to_openmm()\n",
    "    openmm_simulation.context.setPositions(new_positions)\n",
    "    openmm_simulation.minimizeEnergy()\n",
    "    state = openmm_simulation.context.getState(getPositions=True, getEnergy=True)\n",
    "    with open('result.pdb', 'a+') as output:\n",
    "        PDBFile.writeModel(openmm_simulation.topology, state.getPositions(), output, modelIndex=i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract minimized conformers to SDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/data/homezvol3/ekubinco/miniforge3/envs/ssbind/lib/python3.11/site-packages/MDAnalysis/converters/RDKit.py:473: UserWarning: No `bonds` attribute in this AtomGroup. Guessing bonds based on atoms coordinates\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "import MDAnalysis as mda\n",
    "from rdkit import Chem\n",
    "\n",
    "u = mda.Universe(\"result.pdb\")\n",
    "elements = mda.topology.guessers.guess_types(u.atoms.names)\n",
    "u.add_TopologyAttr(\"elements\", elements)\n",
    "atoms = u.select_atoms(\"resname UNK\")\n",
    "\n",
    "sdwriter = Chem.SDWriter(\"minimized_conformers.sdf\")\n",
    "for _ in u.trajectory:\n",
    "    sdwriter.write(atoms.convert_to(\"RDKIT\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}