Create a system mixing SMIRNOFF and non-SMIRNOFF-formatted force fields

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.

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.

Parametrize a molecule with the Parsley force field

First, we parametrize the ligand (toluene) with the SMIRNOFF-format Parsley force field through the usual route to create an OpenMM System.



In [1]:

    
from simtk.openmm.app import PDBFile

from openforcefield.utils import get_data_file_path
from openforcefield.topology import Molecule, Topology
from openforcefield.typing.engines.smirnoff import ForceField



In [2]:

    
# 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 Parsley force field from disk.
force_field = ForceField('openff_unconstrained-1.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.

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.



In [3]:

    
import parmed

# Convert OpenMM System into a ParmEd Structure.
toluene_structure = parmed.openmm.load_topology(toluene_pdbfile.topology,
                                                toluene_system,
                                                xyz=toluene_pdbfile.positions)

Create a ParmEd `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.

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 .



In [4]:

    
# Load the AMBER protein force field parameters through OpenMM.
from simtk.openmm import app
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/openforcefield/issues/259 for more details). We'll re-add the hydrogen bond constraints when we'll create the System for the complex.

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.



In [5]:

    
# 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)

Combine receptor and ligand structures

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.

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.



In [6]:

    
complex_structure = toluene_structure + t4_structure

Convert back the structure into an OpenMM System

Once we have the Structure of the complex, we can chose to create a System object that we can simulate with OpenMM.



In [7]:

    
from simtk.openmm.app import NoCutoff, HBonds
from simtk import unit

# 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)

Create a system mixing SMIRNOFF and non-SMIRNOFF-formatted force fields

Parametrize a molecule with the Parsley force field

Create a ParmEd Structure of an AMBER-parametrized receptor in TIP3P-FB water

Combine receptor and ligand structures

Convert back the structure into an OpenMM System

Create a ParmEd `Structure` of an AMBER-parametrized receptor in TIP3P-FB water