Mixing Sage and Amber: BRD4 benchmark with ParmEd

This example applies SMIRNOFF-format parameters to BRD4 inhibitors from the living review on binding free energy benchmark systems by Mobley and Gilson. The BRD4 system comes from the accompanying GitHub repository.

This example uses ParmEd to combine a protein parameterized with Amber with a ligand parameterized with SMIRNOFF. This example is meant to illustrate how to apply parameters to a single ligand, but it’s also easy to process many ligands.

# Retrieve protein and ligand files for BRD4 and a docked inhibitor from the benchmark systems GitHub repository
# https://github.com/MobleyLab/benchmarksets
import parmed
import requests
from openmm import XmlSerializer, app, unit
from openmm.app import HBonds, NoCutoff, PDBFile

repo_url = (
    "https://raw.githubusercontent.com/MobleyLab/benchmarksets/master/input_files/"
)
sources = {
    "receptor.pdb": repo_url + "BRD4/pdb/BRD4.pdb",
    "ligand.pdb": repo_url + "BRD4/pdb/ligand-1.pdb",
    "ligand.sdf": repo_url + "BRD4/sdf/ligand-1.sdf",
}
for filename, url in sources.items():
    r = requests.get(url)
    open(filename, "w").write(r.text)

Parametrize a molecule with the SMIRNOFF-format “Sage” force field

First, we parametrize the ligand with the Sage force field.

We use both a PDB file and an SDF file for the ligand.

# Create an OpenFF Molecule object from the ligand SDF fiel
from openff.toolkit import ForceField, Molecule

ligand_off_molecule = Molecule("ligand.sdf")

# Load the SMIRNOFF-format Sage force field
force_field = ForceField("openff_unconstrained-2.2.0.offxml")

# Parametrize the ligand molecule by creating a Topology object from it
ligand_system = force_field.create_openmm_system(ligand_off_molecule.to_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. 
# Read in the coordinates of the ligand from the PDB file
ligand_pdbfile = PDBFile("ligand.pdb")

# Convert OpenMM System object containing ligand parameters into a ParmEd Structure.
ligand_structure = parmed.openmm.load_topology(
    ligand_pdbfile.topology, ligand_system, xyz=ligand_pdbfile.positions
)

Create a ParmEd Structure of an AMBER-parametrized receptor

We have to create a ParmEd Structure containing positions and parameters of the receptor (BRD4) that we will then combine with the positions and parameters in the ligand Structure we created above.

First, we use OpenMM to assign Amber14 parameters using OpenMM, but you can also create a Structure from Amber prmtop and inpcrd files.

Note: 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 . 
receptor_pdbfile = PDBFile("receptor.pdb")

# Load the AMBER protein force field through OpenMM.
omm_forcefield = app.ForceField("amber14-all.xml")

# Parameterize the protein.
receptor_system = omm_forcefield.createSystem(receptor_pdbfile.topology)

# Convert the protein System into a ParmEd Structure.
receptor_structure = parmed.openmm.load_topology(
    receptor_pdbfile.topology, receptor_system, xyz=receptor_pdbfile.positions
)

Combine receptor and ligand structures

We can then merge the receptor and ligand Structure objects into a single Structure containing both positions and parameters.

complex_structure = receptor_structure + ligand_structure

Export to OpenMM

Once we have the Structure for the complex, containing both positions and parameters, we can chose to create an OpenMM System object that we can simulate using 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,
)

# Export the System to an OpenMM System XML and PDB file
complex_structure.save("complex.pdb", overwrite=True)

with open("complex.xml", "w") as f:
    f.write(XmlSerializer.serialize(complex_system))

Export to Amber

We can also export the Structure to Amber prmtop and inpcrd files using ParmEd:

# Export the Structure to AMBER files
complex_structure.save("complex.prmtop", overwrite=True)
complex_structure.save("complex.inpcrd", overwrite=True)

Export to gromacs

We can export the Structure to gromacs gro and top files using ParmEd:

# Export the Structure to Gromacs files
complex_structure.save("complex.gro", overwrite=True)
complex_structure.save("complex.top", overwrite=True)