{ "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) solvated in water is assigned AMBER and TIP3P-FB parameters through the ParmEd library.\n", "\n", "We'll need two PDB files. One for the ligand in vacuum, and one for the solvated protein without ligand. The coordinates of the protein-ligand complex will be determined by these PDB files, so their positions needs to be consistent with the ligand being positioned in the binding pocket if this is the desired initial configuration of your simulation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Parametrize a molecule with the Parsley force field\n", "\n", "First, we parametrize the ligand (toluene) with the SMIRNOFF-format Parsley 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 openff.toolkit.utils import get_data_file_path\n", "from openff.toolkit.topology import Molecule, Topology\n", "from openff.toolkit.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_file_path('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 Parsley force field from disk.\n", "force_field = ForceField('openff_unconstrained-1.0.0.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", "
\n", " Warning: 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", "
" ] }, { "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 in TIP3P-FB water\n", "\n", "We have to create a ParmEd `Structure` of the receptor (T4 Lysozyme) to combine to the toluene `Structure`. Here we assign `amber99sbildn` to a T4 Lysozyme receptor solvated in TIP3P-FB water using OpenMM.\n", "\n", "First, we load the parameters and the PDB file including positions for the protein, water, and ion atoms.\n", "\n", "
\n", " Note: If you already have AMBER (prmtop/inpcrd), GROMACS (top/gro), or any other file supported by ParmEd specifying the parameters for the solvated protein, you can simply load the files directly into a Structure using ParmEd's functionalities. See https://parmed.github.io/ParmEd/html/readwrite.html .\n", "
" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Load the AMBER protein force field parameters through OpenMM.\n", "from simtk.openmm import app\n", "omm_forcefield = app.ForceField('amber99sbildn.xml', 'tip3pfb.xml')\n", "\n", "# Load the solvated receptor from a PDB file.\n", "t4_pdb_file_path = get_data_file_path('systems/test_systems/T4_lysozyme_water_ions.pdb')\n", "t4_pdb_file = PDBFile(t4_pdb_file_path)\n", "\n", "# Obtain the updated OpenMM Topology and positions.\n", "omm_topology = t4_pdb_file.getTopology()\n", "positions = t4_pdb_file.getPositions()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We then create a parametrized OpenMM `System` and convert it to a `Structure`.\n", "\n", "Note the `rigidWater=False` argument in `ForceField.createSystem()`. This is necessary to work around a problem araising with ParmEd in reading the parameters of constrained bonds from an OpenMM `System` (see https://github.com/openforcefield/openff-toolkit/issues/259 for more details). We'll re-add the hydrogen bond constraints when we'll create the `System` for the complex.\n", "\n", "
\n", " Note: 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", "
" ] }, { "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. Note that the coordinates of protein and ligand in this example are determined by the PDB file, and they are already consistent with the ligand being positioned in the binding pocket.\n", "\n", "
\n", " Note: If you want to include water molecules in the binding site, you will have to be careful to place them so that they won't create steric clashes once the ligand is inserted.\n", "
" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "complex_structure = toluene_structure + t4_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", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }