<commit_before>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a system mixing SMIRNOFF and non-SMIRNOFF-formatted force fields\n",
    "\n",
    "This example shows how to create a receptor-ligand `System` where the ligand (toluene) is parametrized with a SMIRNOFF force field and the protein (T4 Lysozyme) is assigned AMBER parameters through the ParmEd library."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parametrize a molecule with smirnoff99Frosst\n",
    "\n",
    "First, we parametrize the ligand (toluene) with the smirnoff99Frosst force field through the usual route"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simtk.openmm.app import PDBFile\n",
    "\n",
    "from openforcefield.utils import get_data_filename\n",
    "from openforcefield.topology import Molecule, Topology\n",
    "from openforcefield.typing.engines.smirnoff import ForceField"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an OpenFF Topology of toluene from a pdb file.\n",
    "toluene_pdb_file_path = get_data_filename('molecules/toluene.pdb')\n",
    "toluene_pdbfile = PDBFile(toluene_pdb_file_path)\n",
    "toluene = Molecule.from_smiles('Cc1ccccc1')\n",
    "off_topology = Topology.from_openmm(openmm_topology=toluene_pdbfile.topology,\n",
    "                                    unique_molecules=[toluene])\n",
    "\n",
    "# Load the smirnoff99Frosst system from disk.\n",
    "force_field = ForceField('smirnoff99Frosst.offxml')\n",
    "\n",
    "# Parametrize the toluene molecule.\n",
    "toluene_system = force_field.create_openmm_system(off_topology)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and we convert the OpenMM `System` to a ParmEd `Structure` that we'll be able to mix with the protein\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "  <b>Warning:</b> ParmEd's Structure model is inspired by AMBER. Some information in an OpenMM System are not directly translatable into a Structure. In particular, long-range interaction treatment method (e.g., PME, CutoffPeriodic) and parameters (e.g., cutoff and cutoff switching distance, PME error tolerance) are known to be lost during the conversion.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import parmed\n",
    "\n",
    "# Convert OpenMM System into a ParmEd Structure.\n",
    "toluene_structure = parmed.openmm.load_topology(toluene_pdbfile.topology,\n",
    "                                                toluene_system,\n",
    "                                                xyz=toluene_pdbfile.positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a ParmEd `Structure` of an AMBER-parametrized receptor\n",
    "\n",
    "We have to create a ParmEd `Structure` of the receptor (T4 Lysozyme) to combine to the toluene `Structure`. Here we assign `amber99sbildn` parameters using OpenMM.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "    <b>Note:</b> If you already have AMBER (prmtop/inpcrd), GROMACS (top/gro), or any other file specifying the protein parameters supported by ParmEd, you can simply load the files directly into a Structure using ParmEd's functionalities. See https://parmed.github.io/ParmEd/html/readwrite.html .\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "t4_pdb_file_path = get_data_filename('proteins/T4-protein.pdb')\n",
    "t4_pdbfile = PDBFile(t4_pdb_file_path)\n",
    "\n",
    "# Load the AMBER protein force field through OpenMM.\n",
    "from simtk.openmm import app\n",
    "omm_forcefield = app.ForceField('amber99sbildn.xml')\n",
    "\n",
    "# Parameterize the protein.\n",
    "t4_system = omm_forcefield.createSystem(t4_pdbfile.topology)\n",
    "\n",
    "# Convert the protein System into a ParmEd Structure.\n",
    "t4_structure = parmed.openmm.load_topology(t4_pdbfile.topology,\n",
    "                                           t4_system,\n",
    "                                           xyz=t4_pdbfile.positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Combine receptor and ligand structures\n",
    "\n",
    "We can then merge the receptor and ligand `Structure` objects into a single one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "complex_structure = t4_structure + toluene_structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert back the structure into an OpenMM System or export to AMBER/GROMACS/CHARMM/OpenMM files\n",
    "\n",
    "Once we have the `Structure` of the complex, we can chose to create a `System` object that we can simulate with OpenMM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simtk.openmm.app import NoCutoff, HBonds\n",
    "from simtk import unit\n",
    "\n",
    "# Convert the Structure to an OpenMM System in vacuum.\n",
    "complex_system = complex_structure.createSystem(nonbondedMethod=NoCutoff,\n",
    "                                                nonbondedCutoff=9.0*unit.angstrom,\n",
    "                                                constraints=HBonds,\n",
    "                                                removeCMMotion=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or export it in some file format to be run later or by a different simulation engine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export the System to an OpenMM XML file.\n",
    "from simtk.openmm import XmlSerializer\n",
    "with open('complex.xml', 'w') as f:\n",
    "    f.write(XmlSerializer.serialize(complex_system))\n",
    "\n",
    "# Export the Structure to AMBER files.\n",
    "complex_structure.save('complex.prmtop', overwrite=True)\n",
    "complex_structure.save('complex.inpcrd', overwrite=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",<commit_msg>Fix #259: Add example showing how to workaround ParmEd problems with TIP3P water<commit_after>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a system mixing SMIRNOFF and non-SMIRNOFF-formatted force fields\n",
    "\n",
    "This example shows how to create a receptor-ligand `System` where the ligand (toluene) is parametrized with a SMIRNOFF force field and the protein (T4 Lysozyme) is assigned AMBER parameters through the ParmEd library."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parametrize a molecule with smirnoff99Frosst\n",
    "\n",
    "First, we parametrize the ligand (toluene) with the smirnoff99Frosst force field through the usual route to create an OpenMM `System`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simtk.openmm.app import PDBFile\n",
    "\n",
    "from openforcefield.utils import get_data_filename\n",
    "from openforcefield.topology import Molecule, Topology\n",
    "from openforcefield.typing.engines.smirnoff import ForceField"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an OpenFF Topology of toluene from a pdb file.\n",
    "toluene_pdb_file_path = get_data_filename('molecules/toluene.pdb')\n",
    "toluene_pdbfile = PDBFile(toluene_pdb_file_path)\n",
    "toluene = Molecule.from_smiles('Cc1ccccc1')\n",
    "off_topology = Topology.from_openmm(openmm_topology=toluene_pdbfile.topology,\n",
    "                                    unique_molecules=[toluene])\n",
    "\n",
    "# Load the smirnoff99Frosst system from disk.\n",
    "force_field = ForceField('smirnoff99Frosst.offxml')\n",
    "\n",
    "# Parametrize the toluene molecule.\n",
    "toluene_system = force_field.create_openmm_system(off_topology)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and we convert the OpenMM `System` to a ParmEd `Structure` that we'll be able to mix with the protein.\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "  <b>Warning:</b> ParmEd's Structure model is inspired by AMBER. Some information in an OpenMM System are not directly translatable into a Structure. In particular, long-range interaction treatment method (e.g., PME, CutoffPeriodic) and parameters (e.g., cutoff and cutoff switching distance, PME error tolerance) are known to be lost during the conversion.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import parmed\n",
    "\n",
    "# Convert OpenMM System into a ParmEd Structure.\n",
    "toluene_structure = parmed.openmm.load_topology(toluene_pdbfile.topology,\n",
    "                                                toluene_system,\n",
    "                                                xyz=toluene_pdbfile.positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a ParmEd `Structure` of an AMBER-parametrized receptor\n",
    "\n",
    "We have to create a ParmEd `Structure` of the receptor (T4 Lysozyme) to combine to the toluene `Structure`. Here we assign AMBER ff14SB parameters using OpenMM.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "    <b>Note:</b> If you already have AMBER (prmtop/inpcrd), GROMACS (top/gro), or any other file specifying the protein parameters supported by ParmEd, you can simply load the files directly into a Structure using ParmEd's functionalities. See https://parmed.github.io/ParmEd/html/readwrite.html .\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the (unsolvated) receptor from a PDB file.\n",
    "t4_pdb_file_path = get_data_filename('proteins/T4-protein.pdb')\n",
    "t4_pdbfile = PDBFile(t4_pdb_file_path)\n",
    "\n",
    "# Load the AMBER protein force field through OpenMM.\n",
    "from simtk.openmm import app\n",
    "omm_forcefield = app.ForceField('amber14/protein.ff14SB.xml', 'amber14/tip3p.xml')\n",
    "\n",
    "# Solvate the protein in TIP3P water.\n",
    "from simtk import unit\n",
    "from simtk.openmm.app.modeller import Modeller\n",
    "modeller = Modeller(t4_pdbfile.topology, t4_pdbfile.positions)\n",
    "modeller.addSolvent(omm_forcefield, model='tip3p', padding=11*unit.angstrom, neutralize=True)\n",
    "\n",
    "# Obtain the updated OpenMM Topology and positions including water molecules.\n",
    "omm_topology = modeller.getTopology()\n",
    "positions = modeller.getPositions()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When creating an OpenMM `System` with water molecules that has to be converted to a ParmEd `Structure`, you should always specify `rigidWater=False` in `ForceField.createSystem()`. This is to work around a problem with ParmEd in reading parameters off a `System` for constrained bonds (see https://github.com/openforcefield/openforcefield/issues/259 for more details). We'll re-add the hydrogen bonds constraints when we'll create the `System` for the complex.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "    <b>Note:</b> If you don't solvate the system or if you load it directly from AMBER, GROMACS, or other files directly to ParmEd, you won't need extra precautions.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameterize the protein.\n",
    "t4_system = omm_forcefield.createSystem(omm_topology, rigidWater=False)\n",
    "\n",
    "# Convert the protein System into a ParmEd Structure.\n",
    "t4_structure = parmed.openmm.load_topology(omm_topology,\n",
    "                                           t4_system,\n",
    "                                           xyz=positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Combine receptor and ligand structures\n",
    "\n",
    "We can then merge the receptor and ligand `Structure` objects to form the complex."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "complex_structure = t4_structure + toluene_structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert back the structure into an OpenMM System\n",
    "\n",
    "Once we have the `Structure` of the complex, we can chose to create a `System` object that we can simulate with OpenMM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from simtk.openmm.app import NoCutoff, HBonds\n",
    "from simtk import unit\n",
    "\n",
    "# Convert the Structure to an OpenMM System in vacuum.\n",
    "complex_system = complex_structure.createSystem(nonbondedMethod=NoCutoff,\n",
    "                                                nonbondedCutoff=9.0*unit.angstrom,\n",
    "                                                constraints=HBonds,\n",
    "                                                removeCMMotion=False)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
{'fraction_non_alphanumeric': 0.14709560867148416, 'fraction_numerical': 0.006461923290717065, 'mean_word_length': 2.7609093284557096, 'pattern_counts': {'":': 141, '<': 23, '<?xml version=': 0, '>': 23, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 6}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 32, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7593477', 'n_tokens_mistral': 4900, 'n_tokens_neox': 4708, 'n_words': 1205}