Obtaining CustomNonbondedForce's shortest pairwise interaction

[Mark Williamson]

Dear Forum,

This question is cast within the context of using the rich OpenMM framework as a basic docking engine, in which there is only one class of force, CustomNonbondedForce, and making use of State.getPotentialEnergy() as a scoring metric.

If one should set the following cutoff scheme with a CustomNonbondedForce object :

Code: Select all

  exampleCustomNonbondedForce.setNonbondedMethod(mm.CustomNonbondedForce.CutoffNonPeriodic)
  exampleCustomNonbondedForce.setCutoffDistance(cutoff)

Then the sum of pairwise interactions will be calculated according the object's given energy expression, that satisfy the the distance cutoff.

My question is, is it possible to select only one of these interactions pairs, based the minimum distance of the set of all pairs? Current ideas include perhaps using setExclusionParticles() in an iterative manner to whittle the pair list down to one, or perhaps making use of the select() and/or step() functions within the energy expression. I sense I am only scratching at the surface of Lepton's functionality.

Regards,

Mark

[Peter Eastman]

You might be able to hack something together like that, but it would be a really inefficient way of implementing it. Let me make sure I understand exactly what you want to do. You want to calculate an interaction between two molecules by identifying the very closest distance between them, right? Given all atoms in molecule 1 and all atoms in molecule 2, you want to find the very closest pair. And then you want to compute a function of that distance.

I also gather that you aren't trying to do dynamics? So you only need to compute energy, not force. Is all of that correct?

[Mark Williamson]

Thank you the reply on this.

You want to calculate an interaction between two molecules by identifying the very closest distance between them, right? Given all atoms in molecule 1 and all atoms in molecule 2, you want to find the very closest pair. And then you want to compute a function of that distance.

Correct, I also attach a diagram of a contrived example of what we hope to achieve with this.

The following code attempts to match what is going on in the diagram, but currently includes all pairwise interactions, hence the energy is 49.0 kJ/mol and not 16.0 kJ/mol.

Code: Select all

import simtk.openmm.app as app
import simtk.openmm as mm
import simtk.unit as unit
import simtk.openmm.vec3 as v
import sys



# Energy rule between two points
exampleCustomNonbondedForce = mm.CustomNonbondedForce("ssip1 * ssip2")
exampleCustomNonbondedForce.addPerParticleParameter("ssip")

# Set the nature of the Cutoff and its distance
exampleCustomNonbondedForce.setNonbondedMethod(mm.CustomNonbondedForce.CutoffNonPeriodic)
exampleCustomNonbondedForce.setCutoffDistance(20 * unit.angstrom)

system = mm.System()
topology = app.Topology()
newChain = topology.addChain()

for i in range(2):
    system.addParticle(15.99943*unit.amu)
    system.addParticle(1.007947*unit.amu)
    system.addParticle(1.007947*unit.amu)
    h2o = topology.addResidue("WAT", newChain)

    topology.addAtom("O" , app.element.oxygen, h2o)
    topology.addAtom("H1", app.element.hydrogen, h2o)
    topology.addAtom("H2", app.element.hydrogen, h2o)

    # Contrived SSIP terms for the example
    exampleCustomNonbondedForce.addParticle([3])
    exampleCustomNonbondedForce.addParticle([2])
    exampleCustomNonbondedForce.addParticle([2])


# Exclude intra molecular terms
exampleCustomNonbondedForce.addExclusion(0,1)
exampleCustomNonbondedForce.addExclusion(0,2)
exampleCustomNonbondedForce.addExclusion(1,2)

exampleCustomNonbondedForce.addExclusion(3,4)
exampleCustomNonbondedForce.addExclusion(3,5)
exampleCustomNonbondedForce.addExclusion(4,5)

system.addForce(exampleCustomNonbondedForce)

# Checking
print("system.getNumParticles() is ", system.getNumParticles())
print("system.getNumForces() is ", system.getNumForces())

# Create Integrator and simulations
integrator = mm.VerletIntegrator(1 * unit.femtosecond)
simulation = app.Simulation(topology, system, integrator)

positions = [v.Vec3(x=-0.6197, y=-0.4387, z=0.0044), v.Vec3(x=-0.5227, y=-0.45439999999999997, z=0.0061), v.Vec3(x=-0.6508, y=-0.5183, z=-0.0414), 
             v.Vec3(x=-0.3885, y=-0.5031, z=0.0077), v.Vec3(x=-0.2502, y=-0.5304000000000001, z=-0.05310000000000001), v.Vec3(x=-0.49829999999999997, y=-0.5980000000000001, z=0.0054)]

simulation.context.setPositions(positions)


# Obtain the "energy"
print(simulation.context.getState(getEnergy=True).getPotentialEnergy())


# Saving minimised positions
positions = simulation.context.getState(getPositions=True).getPositions()
app.PDBFile.writeFile(simulation.topology, positions, open('wat.pdb', 'w'))

I also gather that you aren't trying to do dynamics? So you only need to compute energy, not force. Is all of that correct?

This is correct as well, only energy, which is being used as a score within a docking context. No forces, no dynamics.

As a complete aside, we also discovered this paper in which they may be doing something similar via their reference of a "custom knowledge-based force object in OpenMM".

Thanks,

Mark

[Peter Eastman]

Since you described this as docking, I assume you aren't really planning to do it with water, right? That is, you only have two molecules and your potential only depends on a single distance between them? This isn't anything like a water box where there's a separate energy term based on the distance between every pair of molecules?

Assuming that's correct, doing this with OpenMM and a custom force seems a bit of an overkill. All you really need to do is identify the closest distance, and then it's trivial for you to compute the corresponding energy yourself. Something like the compute_contacts() function from MDTraj with scheme='closest' seems like it might be much easier (and likely faster)? See http://mdtraj.org/development/api/gener ... e_contacts.

[Mark Williamson]

Thank you for the reply on this.

Since you described this as docking, I assume you aren't really planning to do it with water, right? That is, you only have two molecules and your potential only depends on a single distance between them?

Correct; the above was just a simple example. We plan to do this with a set of protein/ligand systems.

We have a coarse graining approach using the QM molecular electrostatic potential (MEP) on the Vdw surface which manifests as a set of VirtualSites around a molecule. These VirtualSites (termed Surface Site Interaction Points (SSIPs) within that paper) also have a scalar value associated with them, which is related to the magnitude of the MEP at that point.

We specifically make use of OutOfPlaneSite in order that the correct orientation of this VirtualSites was obtained as a function of the relative heavy atom positions; i.e. a back projection of the critical points of the MEP to obtain the correct weight12, weight13 and weightCross values. The idea is to develop a VirtualSite description for both the protein residues and the ligand and use this to compute the docking score. Please see the image attached.

Something like the compute_contacts() function from MDTraj with scheme='closest' seems like it might be much easier (and likely faster)? See http://mdtraj.org/development/api/gener ... e_contacts.

OK, I did not know about this functionality; thanks for the useful pointer. Do you know if it supports virtual sites and any values associated with them?

[Peter Eastman]

Do you know if it supports virtual sites and any values associated with them?

MDTraj is designed for trajectory analysis, so it doesn't really "know" about virtual sites. That function just takes a set of positions which might have been loaded from a file, or might be something you computed yourself, and calculates the distances between residues.