This example shows how to create a receptor-ligand OpenMM 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.
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.
First, we parametrize the ligand (toluene) with the SMIRNOFF-format Sage force field through the usual route to create an OpenMM System
.
try:
from openmm import app, unit
from openmm.app import HBonds, NoCutoff, PDBFile
except ImportError:
from simtk import unit
from simtk.openmm import app
from simtk.openmm.app import PDBFile, HBonds, NoCutoff
from openff.toolkit.topology import Molecule, Topology
from openff.toolkit.typing.engines.smirnoff import ForceField
from openff.toolkit.utils import get_data_file_path
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
# Create an OpenFF Topology of toluene from a pdb file.
toluene_pdb_file_path = get_data_file_path("molecules/toluene.pdb")
toluene_pdbfile = PDBFile(toluene_pdb_file_path)
toluene = Molecule.from_smiles("Cc1ccccc1")
off_topology = Topology.from_openmm(
openmm_topology=toluene_pdbfile.topology, unique_molecules=[toluene]
)
# Load the Sage force field from disk.
force_field = ForceField("openff_unconstrained-2.0.0.offxml")
# Parametrize the toluene molecule.
toluene_system = force_field.create_openmm_system(off_topology)
and we convert the OpenMM System
to a ParmEd Structure
that we'll be able to mix with the protein.
import parmed
# Convert OpenMM System into a ParmEd Structure.
toluene_structure = parmed.openmm.load_topology(
toluene_pdbfile.topology, toluene_system, xyz=toluene_pdbfile.positions
)
Structure
of an AMBER-parametrized receptor in TIP3P-FB water¶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.
First, we load the parameters and the PDB file including positions for the protein, water, and ion atoms.
# Load the AMBER protein force field parameters through OpenMM.
omm_forcefield = app.ForceField("amber99sbildn.xml", "tip3pfb.xml")
# Load the solvated receptor from a PDB file.
t4_pdb_file_path = get_data_file_path("systems/test_systems/T4_lysozyme_water_ions.pdb")
t4_pdb_file = PDBFile(t4_pdb_file_path)
# Obtain the updated OpenMM Topology and positions.
omm_topology = t4_pdb_file.getTopology()
positions = t4_pdb_file.getPositions()
We then create a parametrized OpenMM System
and convert it to a Structure
.
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.
# Parameterize the protein.
t4_system = omm_forcefield.createSystem(omm_topology, rigidWater=False)
# Convert the protein System into a ParmEd Structure.
t4_structure = parmed.openmm.load_topology(omm_topology, t4_system, xyz=positions)
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.
complex_structure = toluene_structure + t4_structure
Once we have the Structure
of the complex, we can chose to create a System
object that we can simulate with OpenMM.
# Convert the Structure to an OpenMM System in vacuum.
complex_system = complex_structure.createSystem(
nonbondedMethod=NoCutoff,
nonbondedCutoff=9.0 * unit.angstrom,
constraints=HBonds,
removeCMMotion=False,
)