Host-guest systems

This is a refresh of a notebook written by David Mobley with the following modifications:

  • A pre-docked guest is used docked_guest.mol2 was saved and converted to guest.sdf)

  • Mol2 files converted to SDF (OA.mol2 was converted to host.sdf)

In this notebook, a prepared host-guest complex is loaded

from openff.toolkit import ForceField, Molecule, Topology

Just like any other molecule, we begin by loading representations (in this case, files on disk) into Molecule objects.

guest = Molecule.from_file("guest.sdf")
host = Molecule.from_file("host.sdf")

Again, like many other workflows, we “combine” multiple Molecules into a single Topology object. We can also visualize the result to ensure the guest looks reasonably docked into the host. (Use your cursor to move the complex around.)

docked_topology = Topology.from_molecules([guest, host])
docked_topology.visualize()

We can safely use Sage or a similar small molecule force field to parameterize the guest. But the host is large (128 heavy atoms, 184 atoms total) and charge assignment using AM1-BCC may take tens of minutes to hours. We can instead use OpenFF NAGL, which uses a graph-convolutional neural network (GCNN or GNN) to mimic AM1-BCC partial charges (specifically the ELF10 variant). Not counting the time it takes to import the Python module and load the model, charge assignment for the host should take on the order of hundreds of milliseconds. (Larger polymers may take slightly longer but the scaling is sub-linear with number of atoms!)

For more on GCNNs, see Espaloma and its associated paper.

from openff.toolkit.utils.nagl_wrapper import NAGLToolkitWrapper

NAGLToolkitWrapper().assign_partial_charges(
    molecule=host,
    partial_charge_method="openff-gnn-am1bcc-0.1.0-rc.1.pt",
)

host.partial_charges.round(3)

DGL backend not selected or invalid.  Assuming PyTorch for now.



Setting the default backend to "pytorch". You can change it in the ~/.dgl/config.json file or export the DGLBACKEND environment variable.  Valid options are: pytorch, mxnet, tensorflow (all lowercase)
Magnitude
[-0.093 0.013 -0.113 0.013 -0.283 -0.093 -0.441 0.047 0.225 -0.441 0.013
 -0.093 -0.113 -0.283 -0.093 0.013 -0.441 0.047 -0.093 0.013 0.225 -0.441
 -0.441 0.047 -0.093 -0.441 0.013 -0.283 0.013 -0.283 0.013 -0.441 -0.441
 0.225 -0.113 -0.093 -0.113 0.047 -0.093 0.225 -0.179 -0.12 0.098 -0.474
 -0.179 -0.12 0.098 -0.474 -0.179 -0.12 0.098 -0.474 -0.179 -0.12 0.098
 -0.474 -0.171 -0.316 -0.278 -0.278 0.095 0.095 -0.171 -0.316 -0.278
 -0.278 0.095 0.095 -0.171 -0.316 -0.278 -0.278 0.095 0.095 -0.171 -0.316
 -0.278 -0.278 0.095 0.095 -0.259 -0.259 -0.259 -0.259 -0.259 -0.259
 -0.259 -0.073 -0.183 -0.168 -0.183 -0.207 -0.073 -0.259 -0.073 -0.183
 -0.183 -0.168 -0.073 -0.207 -0.073 -0.183 -0.183 -0.168 -0.073 -0.207
 0.285 -0.635 0.285 -0.635 0.285 -0.635 -0.474 -0.474 -0.474 0.285 -0.635
 -0.635 -0.635 -0.635 -0.635 -0.474 -0.073 -0.183 -0.168 -0.183 -0.207
 -0.073 0.023 0.057 0.047 -0.047 0.023 0.057 0.047 -0.047 0.047 0.057
 0.057 -0.047 0.023 0.023 0.047 -0.047 -0.213 -0.213 -0.071 -0.071 -0.213
 -0.213 -0.071 -0.071 -0.213 -0.213 -0.071 -0.071 -0.213 -0.213 -0.071
 -0.071 0.058 0.046 0.046 0.058 0.046 0.046 0.058 0.046 0.046 0.058 0.046
 0.046 -0.049 -0.008 -0.049 -0.049 -0.049 -0.008 -0.049 -0.049 -0.008
 -0.049 -0.008 -0.049]
Unitselementary_charge

Now that we have partial charges for the host, we can wire it through parameterization using the charge_from_molecules argument. Interchange will recognize these charges and not attempt running AM1-BCC for the host, though it will use to assign charges to the guest molecule.

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

out = sage.create_interchange(topology=docked_topology, charge_from_molecules=[host])

Now that the Interchange object is created, you can run simulations in a number of engines. Here we’ll run a quick energy minimization and then a thirty-second OpenMM simulation. The result is a trajectory, viewable with NGLview, that shows a few tens or hundreds of frames of this host-guest complex dancing around in vacuo.

None of this workflow required OpenMM until now - you can swap these steps out for analogous operations in GROMACS, Amber, or LAMMPS!

import openmm
import openmm.app
import openmm.unit

simulation = out.to_openmm_simulation(
    openmm.LangevinMiddleIntegrator(
        300.0 * openmm.unit.kelvin,
        1.0 / openmm.unit.picosecond,
        2.0 * openmm.unit.femtosecond,
    )
)

simulation.minimizeEnergy()

dcd_reporter = openmm.app.DCDReporter("trajectory.dcd", 1000)
simulation.reporters.append(dcd_reporter)

simulation.context.setVelocitiesToTemperature(300.0 * openmm.unit.kelvin)

simulation.runForClockTime(0.5 * openmm.unit.minute)

import mdtraj
import nglview

nglview.show_mdtraj(
    mdtraj.load(
        "trajectory.dcd",
        top=mdtraj.Topology.from_openmm(
            out.to_openmm_topology(),
        ),
    )
)