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, Quantity, Topology

Combining parametrization results from different force fields requires the use of Interchange.combine, which is a newer feature and has not been tested for every conceivable combination of inputs. It may not work as expected for all systems. If extending this workflow, please validate results to avoid surprises in simulation results!

ligand_molecule = Molecule.from_file("ligand.sdf")

sage = ForceField("openff-2.2.0.offxml")

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

receptor_topology = Topology.from_pdb("receptor.pdb")

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

receptor_system = ff14sb.create_interchange(topology=receptor_topology)

complex_system = receptor_system.combine(ligand_system)

# Make pseudo-vacuum because GROMACS does not support proper vacuum simulations
complex_system.box = Quantity([5, 5, 5], "nanometer")

/home/runner/micromamba/envs/openff-docs-examples/lib/python3.11/site-packages/openff/interchange/operations/_combine.py:77: InterchangeCombinationWarning: Interchange object combination is complex and likely to produce strange results outside of a limited set of use cases it has been tested in. Any workflow using this method is not guaranteed to be suitable for production or stable between versions. Use with extreme caution and thoroughly validate results!
  warnings.warn(

/home/runner/micromamba/envs/openff-docs-examples/lib/python3.11/site-packages/openff/interchange/operations/_combine.py:57: InterchangeCombinationWarning: Found electrostatics 1-4 scaling factors of 5/6 with slightly different rounding (0.833333 and 0.8333333333). This likely stems from OpenFF using more digits in rounding 1/1.2. The value of 0.8333333333 will be used, which may or may not introduce small errors. 
  warnings.warn(

Export to OpenMM

By default, Interchange.to_openmm() returns an openmm.System object. There is also Interchange.to_openmm_simulation() which requires an openmm.Integrator as an argument and returns a openmm.app.Simulation object containing information about the topology, box vectors, and particle positions.

See the Interchange user docs for more optional features.

complex_system.to_openmm()

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

Export to Amber

Amber parameter/topology (.prmtop), coordinate (.inpcrd) and run control (.in) files can each be written out separately with individual Interchange.to_x methods. These are wrapped up into a convenient Interchange.to_amber method which handles them all at once.

See the Interchange user docs for more.

complex_system.to_amber("complex")

/home/runner/micromamba/envs/openff-docs-examples/lib/python3.11/site-packages/openff/interchange/components/mdconfig.py:473: UserWarning: Ambiguous failure while processing constraints. Constraining h-bonds as a stopgap.
  warnings.warn(
/home/runner/micromamba/envs/openff-docs-examples/lib/python3.11/site-packages/openff/interchange/components/mdconfig.py:401: SwitchingFunctionNotImplementedWarning: A switching distance 8.0 angstrom was specified by the force field, but Amber does not implement a switching function. Using a hard cut-off instead. Non-bonded interactions will be affected.
  warnings.warn(

Export to GROMACS

Gromacs topology (.top), coordinate (.gro) and run control (.mdp) files can each be written out separately with individual Interchange.to_x methods. These are wrapped up into a convenient Interchange.to_gromacs method which handles them all at once.

See the Interchange user docs for more.

complex_system.to_gromacs(prefix="complex")

/home/runner/micromamba/envs/openff-docs-examples/lib/python3.11/site-packages/openff/interchange/components/mdconfig.py:473: UserWarning: Ambiguous failure while processing constraints. Constraining h-bonds as a stopgap.
  warnings.warn(