Mixing Sage and Amber force fields: BRD4 benchmark

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.

# Retrieve protein and ligand files for BRD4 and a docked inhibitor from the benchmark systems GitHub repository
# https://github.com/MobleyLab/benchmarksets
import requests

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)

from openff.toolkit import ForceField, Molecule

Later we will use Interchange.__add__, which is an experimental feature. It needs to be turned on by setting the environment variable below, and by doing so we accept some stability and accuracy risks.

%env INTERCHANGE_EXPERIMENTAL=1

env: INTERCHANGE_EXPERIMENTAL=1

ligand_molecule = Molecule("ligand.sdf")
sage = ForceField("openff-2.1.0.offxml")

ligand = sage.create_interchange(topology=ligand_molecule.to_topology())

receptor_molecule = Molecule.from_polymer_pdb("receptor.pdb")

ff14sb = ForceField("ff14sb_off_impropers_0.0.3.offxml")

receptor = ff14sb.create_interchange(topology=receptor_molecule.to_topology())

complex_system = receptor + ligand

# TODO
# complex.box = pdbfile box vectors ...
# complex.positions = np.vstack([receptor.positions, ligand.positions])

/home/runner/micromamba/envs/openff-docs-examples/lib/python3.10/site-packages/openff/utilities/utilities.py:80: MoleculeDeprecationWarning: `Molecule.from_polymer_pdb` is deprecated in favor of `Topology.from_pdb`, the recommended method for loading PDB files. This method will be removed in a future release of the OpenFF Toolkit.
  return function(*args, **kwargs)

/home/runner/micromamba/envs/openff-docs-examples/lib/python3.10/site-packages/openff/interchange/_experimental.py:35: InterchangeDeprecationWarning: The `+` operator is deprecated. Use `Interchange.combine` instead.
  return func(*args, **kwargs)
/home/runner/micromamba/envs/openff-docs-examples/lib/python3.10/site-packages/openff/interchange/components/interchange.py:954: UserWarning: Interchange object combination is experimental and likely to produce strange results. Any workflow using this method is not guaranteed to be suitable for production. Use with extreme caution and thoroughly validate results!
  return _combine(self, other)

Export to OpenMM

complex_system.to_openmm()

<openmm.openmm.System; proxy of <Swig Object of type 'OpenMM::System *' at 0x7fdf5566fd20> >

Export to Amber

# TODO: Fix inferring residue information with mixed topology
if False:
    complex_system.to_inpcrd("complex.inpcrd")
    complex_system.to_prmtop("complex.prmtop")

Export to GROMACS

complex_system.to_gro("complex.gro")
complex_system.to_top("complex.top")

/home/runner/micromamba/envs/openff-docs-examples/lib/python3.10/site-packages/openff/interchange/interop/gromacs/export/_export.py:48: UserWarning: WARNING: System defined with no box vectors, which GROMACS does not offically support in versions 2020 or newer (see https://gitlab.com/gromacs/gromacs/-/issues/3526). Setting box vectors to a 5  nm cube.
  self._write_gro(gro, decimal)