Implementing a TI protocol using CustomNonbondedForce

[Mark Williamson]

Dear Forum,

I am attempting to develop a TI protocol using OpenMM 5.0’s Python wrappers and the CustomNonbondedForce object, for lambda mapping in a (non-bonded) ligand into an existing simulation.

I am looking to use/test a range of bespoke EE and VDW softcore potentials that I have successfully created using CustomNonbondedForce object. In addition, this can also be scaled as a function of lamda, (i.e. simulation.context.setParameter("lambda", value ) during a simulation (please see this post for more information https://simtk.org/forums/viewtopic.php?f=161&t=3781 ). The flexibility of the CustomNonbondedForce object allows rapid prototyping of different softcore potential forms.

The argon-chemical-potential.py example provides a great template for this, but my issue at the moment is developing the best strategy for tagging/marking a subset of ligand atoms in an existing system, for these bespoke potentials to be applied to. Ironically, I assumed this part would be the easiest, however, I am not clear on the best route. I am looking for some general guidance on which approach to use. Let me enumerate the various approaches I see to this and perhaps solicit advice on the best one:

1) Build my final system object normally and then substitute all ligan NonbondedForce objects for a pair of CustomNonbondedForce objects for a softcore VDW and EE term. The general approach would be to copy the existing system object to a newSystem object and then replace certain NonbondedForce with CustomNonbondedForce terms.

However, how I do select the specific NonbondedForce objects that pertain to the ligand since this information has already been “lost” once Modeller has generated a system object. Could I perhaps use forcegroups in Modeller to tag NonbondedForces that need to be replaced? What about exclusions (i.e. bonded+angle) and 1-4 interactions ? Unlike NonbondedForce, CustomNonbondedForce does not have an addException() method. I assume I am going to have to code methods to work these out again?

2) Add another method to Modeller, similar to add(), say ghost(), that assigns softcore CustomNonbondedForces for the NB of the method’s Topology parameter. (Maybe not, this could be tricky since point 3 might need to be implemented as well.)

3) Create a CustomNonbondedForceGenerator in ForceField, giving a set of CustomNonbondedForces instead of NonbondedForces for the ligand?

Also, are there any other things that I may have overlooked here?

Thanks,

Mark

[Peter Eastman]

Hi Mark,

However, how I do select the specific NonbondedForce objects that pertain to the ligand

Under what circumstances would you have more than one NonbondedForce in a System? The ForceField would never do that, and I'm not sure what it would be for. I could see having multiple CustomNonbondedForces (if they had different functions), but not NonbondedForces.

Another important point is that every NonbondedForce or CustomNonbondedForce is applied to every atom in the system, not just a subset of them. You can, of course, make it so the force on a particular atom is 0: just multiply by a per-particle parameter that is 1 for the atoms you want to affect, 0 for others. This will give the result you want, but won't help the performance at all: it will still calculate the interaction for every particle, even if that often comes out to 0.

since this information has already been “lost” once Modeller has generated a system object.

I'm not sure what you mean by this. What information has been lost?

Peter

[Justin MacCallum]

Hi Mark,

Why not make lambda a per-particle parameter? Then you replace the existing non-bonded function with this one and set lambda to zero for the particles that have the normal interactions.

Justin

[Justin MacCallum]

Just to elaborate, this would be kind of like your method 1.

• Build your list normally
• Create your custom non-bonded force object
• Loop over all particles
  • Copy the parameters from the old non-bonded force object to your custom one
  • Set lambda appropriately for each particle
• Remove the old non-bonded force object from the system
• Add your new one

I haven't tried anything like this yet (still learning openMM myself), but I believe John Chodera outlines something like this in one of the talks available on the downloads page.

[Mark Williamson]

Dear all,

First, thanks for your responses and suggestions, it is very much appreciated.
Responding to Peter first:

Under what circumstances would you have more than one NonbondedForce in a System?

Apologies, I have been lapse in my terminology. What I was attempting to convey, is, how do I select a specific subset of pairs within a NonBondedForce object?

I'm not sure what you mean by this. What information has been lost?

This relates to both points above. Given an existing System object, I know that I can iterate through a NonBondedForce object using its methods, getNumParticles() and getParticleParameters(), and from this I can obtain a particle’s, index, mass and VDW+EE parameters, however I no longer know the particle’s name, type or even which molecule/residue it is in. Hence, for instance, I cannot change NonBondedForce parameters for pairs in a specific ligand. Now, as I write this, I realise that maybe, I could query the Modeller object to return a list of particle indices of a subset of atoms, say within a ligand or residue. This would work if the indices in the Modeller and System objects are the same... would this be a solution?

(ninja edit; I've attempted this in the code below)

Justin, thanks for the pointers. Let me show you what I have at the moment.... I think I’m very nearly there. This is applying the lambda mapping to all CustomNonbondedForces, however, given what I have said above, I may be able to apply it to a subset.

(2nd ninja edit; I'm also attempting this in the code below)

Code: Select all

# WIP Example use of OpenMM 5.0 to carry out a TI (like) calculation
# with bespoke softcore potentials.
#
# Mark J. Williamson

from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *

from sys import stdout
from copy import deepcopy

#platform = openmm.Platform_getPlatformByName("Reference")


# Set up platform
platform = openmm.Platform_getPlatformByName("CUDA")
platformProperties = {"CudaPrecision":"mixed"}
platformProperties = {"CudaDeviceIndex":"0"}




print "Building system"
# Add ketoconazole which has been parameterised using GAFF
# and converted using JAMBER: https://bitbucket.org/mjw99/jamber
ketoconazolePdb = PDBFile('./leap/ketoconazole.pdb')
forceField = ForceField('tip3p.xml', 'ketoconazole_ff.xml')
modeller = Modeller(ketoconazolePdb.topology, ketoconazolePdb.positions)

# Add missing hydrogens
modeller.addHydrogens(forceField)
# Add TIP3P solvent
#modeller.addSolvent(forceField, model='tip3p', padding=1*angstrom)
modeller.addSolvent(forceField, model='tip3p', padding=8*angstrom)
# Dump modeller structure
PDBFile.writeFile(modeller.getTopology(), modeller.getPositions(), open('modeller.pdb', 'w'))



# Now, define two softcore potentials:


######################################################################################
# caseSoftCoreVDW; see see Eq. 3 in http://dx.doi.org/10.1002/jcc.21909
caseSoftCoreVDW = CustomNonbondedForce("4*epsilon*(l12)*( (sigma / (alphaLJ*(1-l12)*sigma + r^6)^(1/6) )^12 - (sigma / ( alphaLJ*(1-l12)*sigma + r^6)^(1/6) )^6) ;"
"sigma=0.5*(sigma1+sigma2);"
"epsilon=sqrt(epsilon1*epsilon2);"
"alphaLJ=0.5;"
"l12=1-(1-lambda)*step(useLambda1+useLambda2-0.5)");

# Note, the Lorentz-Bertelot rules are being invoked here....

caseSoftCoreVDW.addPerParticleParameter("sigma")
caseSoftCoreVDW.addPerParticleParameter("epsilon")
caseSoftCoreVDW.addPerParticleParameter("useLambda")

# "useLamba" is a per particle parameter that is a function of the global parameter, lambda
# This enables a subset of particles to be affected by the lambda value.
caseSoftCoreVDW.addGlobalParameter("lambda", 1.0)

caseSoftCoreVDW.setNonbondedMethod(CustomNonbondedForce.CutoffPeriodic)
caseSoftCoreVDW.setCutoffDistance(8*angstrom)

####################################################################################
# caseSoftCoreEE; see Eq. 5 in http://dx.doi.org/10.1002/jcc.21909
# Note the reversal of (1-lambda) here... I've got to check this
caseSoftCoreEE = CustomNonbondedForce("l12 * ((q1 * q2) / ( 4 * PI * EPSILON0 * (((1-l12) * BETAC) + r^M)^(1/M) ))  ;"
"PI=3.141592653589;"
"EPSILON0=5.72765E-4;"
"M=6;"
"BETAC=2.5;"
"l12=1-(1-lambda)*step(useLambda1+useLambda2-0.5)");

# EPSILON0 is from ./SimTKUtilities/SimTKOpenMMRealType.h , hence in the correct internal units
# For values of M and B, see page 3260 of 10.1002/jcc.21909


caseSoftCoreEE.addPerParticleParameter("q")
caseSoftCoreEE.addPerParticleParameter("useLambda")

caseSoftCoreEE.addGlobalParameter("lambda", 1.0)

caseSoftCoreEE.setNonbondedMethod(CustomNonbondedForce.CutoffPeriodic)
caseSoftCoreEE.setCutoffDistance(8*angstrom)






# Get ligand index
LigandAtomIndex = []

for residue in modeller.topology.residues():
  if (residue.name == "KTN"):
     for atom in residue.atoms():
        #print atom.name, atom.index
        LigandAtomIndex.append(atom.index)





# Create the original system
system = forceField.createSystem(modeller.topology, nonbondedMethod=PME, nonbondedCutoff=8*angstrom, constraints=HBonds)




# Create new System which will have NonbondedForce terms replaced with CustomNonbondedForce
newSystem = System()

# Copy across particles
for j in range (system.getNumParticles()):
  newSystem.addParticle( system.getParticleMass(j) )

# Copy across constraints
for j in range (system.getNumConstraints()):
  particle1, particle2, distance = system.getConstraintParameters(j)
  newSystem.addConstraint(particle1, particle2, distance)

# Copy across Box vectors
a, b, c = system.getDefaultPeriodicBoxVectors()
newSystem.setDefaultPeriodicBoxVectors(a,b,c)

# Decompose and copy across respective force terms
for i in range(system.getNumForces()):

   force = system.getForce(i)
   print force

   if isinstance( force, openmm.HarmonicBondForce ):
     print "Found " + str(force.getNumBonds()) + " HarmonicBondForce terms"
     # Deep copy these forces into out new state
     copyOfForce = deepcopy(force)
     # Add this to our HarmonicBondSystem
     newSystem.addForce(copyOfForce)

     # Exceptions
     for bondIndex in range(force.getNumBonds()):
       particle1, particle2, length, k = force.getBondParameters(bondIndex)
       caseSoftCoreVDW.addExclusion(particle1, particle2)
       caseSoftCoreEE.addExclusion(particle1, particle2)


   if isinstance( force, openmm.HarmonicAngleForce ):
     print "Found " + str(force.getNumAngles()) + " HarmonicAngleForce terms"
     copyOfForce = deepcopy(force)
     newSystem.addForce(copyOfForce)

     # Exceptions
     for angleIndex in range(force.getNumAngles()):
       particle1, particle2, particle3, length, k = force.getAngleParameters(angleIndex)
       caseSoftCoreVDW.addExclusion(particle1, particle3)
       caseSoftCoreEE.addExclusion(particle1, particle3)


   if isinstance( force, openmm.PeriodicTorsionForce ):
     print "Found " + str(force.getNumTorsions()) + " PeriodicTorsionForce terms"
     copyOfForce = deepcopy(force)
     newSystem.addForce(copyOfForce)


   # We want to substitute NonbondedForce for caseSoftCoreVDW & caseSoftCoreEE
   # TODO, how does one select a subset?
   if isinstance( force, openmm.NonbondedForce ):
     print "Found " + str(force.getNumParticles()) + " NonbondedForce terms"
     for j in range(force.getNumParticles()):
       charge, sigma, epsilon =  force.getParticleParameters(j)
       #print charge
       #print sigma
       #print epsilon

       if j in LigandAtomIndex:
         caseSoftCoreVDW.addParticle( [sigma, epsilon, 1 ]  )
         caseSoftCoreEE.addParticle( [charge, 1 ]  )

       else:

       # Else, a "normal" particle
         caseSoftCoreVDW.addParticle( [sigma, epsilon, 0 ]  )
         caseSoftCoreEE.addParticle( [charge, 0 ]  )


     # Poor initial attempt Exception terms
     print "Found " + str(force.getNumExceptions()) + " NonbondedForce exception terms"
     for k in range(force.getNumExceptions()):

       print force.getExceptionParameters(k)
       particle1, particle2, chargeProd, sigma, epsilon  = force.getExceptionParameters(k)


       # hmmm has this already been done?
       if ( (chargeProd == 0*elementary_charge) & (epsilon == 0*kilojoule/mole) ):
         caseSoftCoreVDW.addExclusion(particle1, particle2)
         caseSoftCoreEE.addExclusion(particle1, particle2)


     # TODO
     # What about 1-4 scaling?

     newSystem.addForce(caseSoftCoreVDW)
     newSystem.addForce(caseSoftCoreEE)




# Print out info about the new system
print "For newSystem"
for i in range(newSystem.getNumForces()):

   force = system.getForce(i)
   print force

   if isinstance( force, openmm.NonbondedForce ):
     print "Found " + str(force.getNumParticles()) + " NonbondedForce terms"
     print "Found " + str(force.getNumExceptions()) + " NonbondedForce exception terms"

   if isinstance( force, openmm.CustomNonbondedForce ):
     print "Found " + str(force.getNumParticles()) + " CustomNonbondedForce terms"
     print "Found " + str(force.getNumExceptions()) + " CustomNonbondedForce exception terms"






######################
# 1) Minimisation    #
######################
print "Minimising system"

integrator = VerletIntegrator(1*femtosecond)

simulation = Simulation(modeller.topology, newSystem, integrator, platform, platformProperties)
#simulation.context.setParameter("lambda",  1  )

simulation.context.setPositions(modeller.positions)
simulation.minimizeEnergy(maxIterations=1000)

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



################################
# 2) Thermalisation under NVT  #
################################
print "Heating system under NVT"
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 1*femtoseconds)


# Set the COM Removal to something sensible
for i in range(newSystem.getNumForces()):
   if (type(newSystem.getForce(i)) == openmm.CMMotionRemover):
      newSystem.getForce(i).setFrequency(1000)

simulation = Simulation(modeller.topology, newSystem, integrator, platform, platformProperties)
simulation.context.setPositions(positions)

simulation.reporters.append(StateDataReporter("heating.txt", 1000, step=True, potentialEnergy=True, temperature=True, density=True))
simulation.reporters.append(PDBReporter('heating.pdb', 1000))

# Debug
#simulation.reporters.append(StateDataReporter("heating.txt", 1, step=True, potentialEnergy=True, temperature=True, density=True))
#simulation.reporters.append(PDBReporter('heating.pdb', 1))


simulation.step(35000) 

# Save the positions and velocities
positions = simulation.context.getState(getPositions=True).getPositions()
velocities = simulation.context.getState(getVelocities=True).getVelocities()


#clear reporters
simulation.reporters = []


####################################
# 4) Density correction under NPT  #
####################################
print "Density correction under NPT"
#TODO, Why does the density decrease here? 
#      Have I missed something?

newSystem.addForce(MonteCarloBarostat(1*bar, 300*kelvin))
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 2*femtoseconds)

simulation = Simulation(modeller.topology, newSystem, integrator, platform, platformProperties)

simulation.context.setPositions(positions)
simulation.context.setVelocities(velocities)

simulation.reporters.append(StateDataReporter("density.txt", 1000, step=True, potentialEnergy=True, temperature=True, density=True))
simulation.reporters.append(PDBReporter('density.pdb', 1000))

simulation.step(35000) # i.e. 20,000 fs == 20 ps == 0.02 ns

# Save the positions and velocities
positions = simulation.context.getState(getPositions=True).getPositions()
velocities = simulation.context.getState(getVelocities=True).getVelocities()

#clear reporters
simulation.reporters = []



####################################
# 5) lambda scaling under NPT      #
####################################

print "lambda scaling under NPT"

simulation.context.setPositions(positions)
simulation.context.setVelocities(velocities)


# Report every 0.1 ns / 100 ps
simulation.reporters.append(StateDataReporter("lambda.txt", 100, step=True, potentialEnergy=True, temperature=True, density=True))
simulation.reporters.append(PDBReporter('lambda.pdb',100))

for i in range(10):
  #1, 0.9, 0.8
  simulation.context.setParameter("lambda",  1 - (i*0.1)  )
  #simulation.context.setParameter("lambda",  1 )
  print "Current lambda value is " + str(simulation.context.getParameter("lambda"))
  simulation.step(10000)

[Peter Eastman]

Ok, got it. When you create a System from a Topology, the particles remain in the same order. So you can just match up per-particle information between them by index.

Peter

[Mark Williamson]

Just an update on this. As part of my ongoing testing, I ran a pretty basic MD protocol on a box of TIP3P water. This involved an initial minimisation, thermalisation under NVT and then density correction under NPT. Using a periodic system of just NonbondedForce objects with NonbondedForce.CutoffPeriodic, this resulted in an acceptable density ~ 1 g/cm^3 @300K.

I then attempted to reproduce this by using a CustomNonbondedForce() representation of the NonbondedForce() object, i.e.:

Code: Select all

  combinedVDWEE = CustomNonbondedForce("4*eps*((sigma/r)^12-(sigma/r)^6)+138.935456*q/r;"
  "q=q1*q2;"
  "sigma=0.5*(sig1+sig2);"
  "eps=sqrt(eps1*eps2)");
  combinedVDWEE.addPerParticleParameter("eps")
  combinedVDWEE.addPerParticleParameter("q")
  combinedVDWEE.addPerParticleParameter("sig")

This was achieved by copying the System object and then replacing all instances of NonbondedForce with CustomNonbondedForce, pretty much as outlined in this post’s code. In theory this should have given the same result, however, this resulted in a final density that was around 0.2 g/cm^3. Investigations into this revealed that the lack of the long range dispersion correction in CustomNonbondedForce, intrinsic in the NonbondedForce, was the cause of this incorrect density.

Now, my question is, how do I work around this? I still want to essentially carry out a TI style calculation using a periodic system.

I’ve not looked too closely at this yet, but how difficult would it be to add the long range dispersion correction and/or PME to CustomNonbondedForce?