Using of variable NonbondedSoftcoreForce() during

[Mark Williamson]

Dear OpenMM users,

I'm looking to use NonbondedSoftcoreForce() to "ghost" in an atom (or collection of atoms) during a single MD simulation by changing the value of softcoreLJLambda from 0.0 to 1.0 over a specific amount of simulation time.

I'm using the Python interface and I think I have a clumsy, but correct, first draft. Would it be possible for someone to have a look over it to make sure that I've not missed something fundamental here. I'm worried that I'm misunderstood the concept of context here.

Thanks,

Mark

Code: Select all

# Example illustrating use of Free Energy OpenMM plugin to lambda map in an Atom
# Based upon argon-chemical-potential.py

from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout
import numpy


# =============================================================================
# Specify simulation parameters
# =============================================================================

nparticles = 216 # number of particles

mass = 39.9 * amu # mass
sigma = 3.4 * angstrom # Lennard-Jones sigma
epsilon = 0.238 * kilocalories_per_mole # Lennard-Jones well-depth
charge = 0.0 * elementary_charge # argon model has no charge

nequil_steps = 5000 # number of dynamics steps for equilibration
nprod_steps = 50000 # number of dynamics steps per production iteration

reduced_density = 0.85 # reduced density rho*sigma^3
temperature = 300 * kelvin # temperature
pressure = 1 * atmosphere # pressure
collision_rate = 5.0 / picosecond # collision rate for Langevin thermostat
timestep = 1 * femtosecond # integrator timestep


# =============================================================================
# Compute box size.
# =============================================================================

volume = nparticles*(sigma**3)/reduced_density
box_edge = volume**(1.0/3.0)
cutoff = min(box_edge*0.49, 2.5*sigma) # Compute cutoff
print "sigma = %s" % sigma
print "box_edge = %s" % box_edge
print "cutoff = %s" % cutoff

# =============================================================================
# Build system
# =============================================================================

# Create argon system where first particle is alchemically modified by lambda_value.
# It is currently zero
lambda_value = 0.0


system = System()
system.setDefaultPeriodicBoxVectors(Vec3(box_edge, 0, 0), Vec3(0, box_edge, 0), Vec3(0, 0, box_edge))

topology = Topology()
newChain = topology.addChain()

#nonbondedforce = NonbondedForce()
nonbondedforce = NonbondedSoftcoreForce()

nonbondedforce.setNonbondedMethod(NonbondedForce.CutoffPeriodic)
nonbondedforce.setCutoffDistance(cutoff)


for particle_index in range(nparticles):
  system.addParticle(mass)
  newResidue = topology.addResidue("UNK", newChain)

  if (particle_index == 0 ):
     # Add alchemically-modified particle.
     # I think I'm using addAtom in the wrong way here...
     topology.addAtom("MODD", "Ar", newResidue)
     nonbondedforce.addParticle(charge, sigma, epsilon, lambda_value)
  else:
     # Add normal particle.
     topology.addAtom("Argo", "ar", newResidue)
     nonbondedforce.addParticle(charge, sigma, epsilon)   

system.addForce(nonbondedforce)

print "System.getNumParticles() is %i"  % system.getNumParticles()
print "system.getNumForces() is %i " % system.getNumForces()


# Create random initial positions.
import numpy.random
positions = Quantity(numpy.random.uniform(high=box_edge/angstroms, size=[nparticles,3]), angstrom)


# Create Integrator
integrator = LangevinIntegrator(temperature, collision_rate, timestep)
simulation = Simulation(topology, system, integrator)

# Initiate from Random set of positions
simulation.context.setPositions(positions)

# Minimize energy from coordinates.
simulation.minimizeEnergy()

print "simulation.system.getNumForces() is %i " % simulation.system.getNumForces()

# Equilibrate.
simulation.step(nequil_steps)

# Run dynamics.
simulation.reporters.append( PDBReporter('output.pdb', 100) )
simulation.reporters.append( StateDataReporter(stdout, 100, step=True, potentialEnergy=True, temperature=True) )



# =============================================================================
# Lambda map in the Argon atom
# =============================================================================

for i in range(10):
  print i
  simulation.context.getSystem().getForce(0).setParticleParameters(0, charge, sigma, epsilon, i*0.1)
  print "Force Object status" 
  print simulation.context.getSystem().getForce(0).getParticleParameters(0)
  simulation.step(1000)

[Peter Eastman]

Hi Mark,

First of all, you should know that the free energy plugin (which provides NonbondedSoftcoreForce) will be going away in the next release. It's never been supported with the OpenCL platform, and CustomNonbondedForce provides a much more powerful tool for doing the same thing.

Second, your code won't work correctly. Per-particle parameters are part of the System definition, which means they can't change. When you create a Context, their values are copied over to it, and if you later modify them the Context won't pick up the new values. (OpenMM 5 actually will provide some ability to update them in a running simulation, but for what you want to do there's already a better option.)

So what's the recommended way of doing this? Use a CustomNonbondedForce that gives the exact softcore form you want. It could be the same one currently used by NonbondedSoftcoreForce if you want, or a different one. Then use a global parameter to turn the interaction on and off:

nonbondedforce.addGlobalParameter("lambda", 1.0)

Your energy expression should depend on lambda such that it's fully turned off when lambda=0 and fully turned on when lambda=1. Now you can update its value at any time:

simulation.context.setParameter("lambda", value)

Does that make sense?

Peter

[Aron Broom]

Ahh, I've been wondering about exactly the same thing.

Thanks to both of you for the postings, and great to hear about the new upcoming version

[Mark Williamson]

Dear Peter,

Thank you very much for the pointers and apologies for the delay in replying.

I think I'm very nearly there; I've modified the code as follows:

Code: Select all

# Example illustrating use of CustomNonbondedForce to lambda map in an Atom
# Based upon argon-chemical-potential.py

from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout
import numpy
from array import *


# =============================================================================
# Specify simulation parameters
# =============================================================================

nparticles = 216 # number of particles

mass = 39.9 * amu # mass
sigma = 3.4 * angstrom # Lennard-Jones sigma
epsilon = 0.238 * kilocalories_per_mole # Lennard-Jones well-depth
charge = 0.0 * elementary_charge # argon model has no charge

nequil_steps = 5000 # number of dynamics steps for equilibration
nprod_steps = 50000 # number of dynamics steps per production iteration

reduced_density = 0.85 # reduced density rho*sigma^3
temperature = 300 * kelvin # temperature
pressure = 1 * atmosphere # pressure
collision_rate = 5.0 / picosecond # collision rate for Langevin thermostat
timestep = 1 * femtosecond # integrator timestep


# =============================================================================
# Compute box size.
# =============================================================================

volume = nparticles*(sigma**3)/reduced_density
box_edge = volume**(1.0/3.0)
cutoff = min(box_edge*0.49, 2.5*sigma) # Compute cutoff
print "sigma = %s" % sigma
print "box_edge = %s" % box_edge
print "cutoff = %s" % cutoff

# =============================================================================
# Build system
# =============================================================================

# Create argon system where first particle is alchemically modified by lambda_value.
# It is currently zero
lambda_value = 0.0


system = System()
system.setDefaultPeriodicBoxVectors(Vec3(box_edge, 0, 0), Vec3(0, box_edge, 0), Vec3(0, 0, box_edge))

topology = Topology()
newChain = topology.addChain()


# Normal nonbonded
#nonBondedForce = NonbondedForce()
#nonBondedForce.setNonbondedMethod(NonbondedForce.CutoffPeriodic)
#nonBondedForce.setCutoffDistance(cutoff)


# Modified VdW
customNonBondedForce = CustomNonbondedForce("4*epsilon*lambda*( 1/( (alphaLJ*(1-lambda) + (r/sigma)^6)^2) - 1/( alphaLJ*(1-lambda) + (r/sigma)^6) ) ;sigma=0.5*(sigma1*sigma2); epsilon=sqrt(epsilon1*epsilon2); alphaLJ=0.5");

customNonBondedForce.addPerParticleParameter("sigma")
customNonBondedForce.addPerParticleParameter("epsilon")

customNonBondedForce.addGlobalParameter("lambda", 1.0)

customNonBondedForce.setNonbondedMethod(NonbondedForce.CutoffPeriodic)
customNonBondedForce.setCutoffDistance(cutoff)


# Assign force types
for particle_index in range(nparticles):
  system.addParticle(mass)
  newResidue = topology.addResidue("UNK", newChain)

  if (particle_index == 0 ):
     # Add alchemically-modified particle.
     topology.addAtom("MODD", "Ar", newResidue)
     # Is there a better way to do this?
     customNonBondedForce.addParticle( array('d',[3.4*NmPerAngstrom, 0.238*KJPerKcal ])  )
     
  else:
     # Add normal particle.
     topology.addAtom("Argo", "ar", newResidue)
     # Is there a better way to do this?
     customNonBondedForce.addParticle( array('d',[3.4*NmPerAngstrom, 0.238*KJPerKcal ])  )

system.addForce(customNonBondedForce)

print "system.getNumParticles() is %i"  % system.getNumParticles()
print "system.getNumForces() is %i " % system.getNumForces()


# Create random initial positions.
import numpy.random
positions = Quantity(numpy.random.uniform(high=box_edge/angstroms, size=[nparticles,3]), angstrom)


# Create Integrator
integrator = LangevinIntegrator(temperature, collision_rate, timestep)
simulation = Simulation(topology, system, integrator)

# Initiate from Random set of positions
simulation.context.setPositions(positions)

# Minimize energy from coordinates.
simulation.minimizeEnergy()

print "simulation.system.getNumForces() is %i " % simulation.system.getNumForces()

# Equilibrate.
simulation.step(nequil_steps)

# Run dynamics.
simulation.reporters.append( PDBReporter('output.pdb', 10000) )
simulation.reporters.append( StateDataReporter(stdout, 10000, step=True, potentialEnergy=True, temperature=True) )



# =============================================================================
# Lambda map in the Argon atom
# =============================================================================

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

This seems to work, however, I am applying the lambda value to all Argon atoms in the system. To operate on a specific atom, I did attempt something like this:

Code: Select all

# Normal nonbonded
nonBondedForce = NonbondedForce()
nonBondedForce.setNonbondedMethod(NonbondedForce.CutoffPeriodic)
nonBondedForce.setCutoffDistance(cutoff)


# Modified VdW
customNonBondedForce = CustomNonbondedForce("4*epsilon*lambda*( 1/( (alphaLJ*(1-lambda) + (r/sigma)^6)^2) - 1/( alphaLJ*(1-lambda) + (r/sigma)^6) ) ;sigma=0.5*(sigma1*sigma2); epsilon=sqrt(epsilon1*epsilon2); alphaLJ=0.5");

customNonBondedForce.addPerParticleParameter("sigma")
customNonBondedForce.addPerParticleParameter("epsilon")

customNonBondedForce.addGlobalParameter("lambda", 1.0)

customNonBondedForce.setNonbondedMethod(NonbondedForce.CutoffPeriodic)
customNonBondedForce.setCutoffDistance(cutoff)


# Assign force types
for particle_index in range(nparticles):
  system.addParticle(mass)
  newResidue = topology.addResidue("UNK", newChain)

  if (particle_index == 0 ):
     # Add alchemically-modified particle.
     topology.addAtom("MODD", "Ar", newResidue)
     # Is there a better way to do this?
     customNonBondedForce.addParticle( array('d',[3.4*NmPerAngstrom, 0.238*KJPerKcal ])  )

  else:
     # Add normal particle.
     topology.addAtom("Argo", "ar", newResidue)
     # Is there a better way to do this?
     nonBondedForce.addParticle( charge, sigma, epsilon)

system.addForce(customNonBondedForce)
system.addForce(nonBondedForce)

But I received an error "CustomNonbondedForce must have exactly as many particles as the System it belongs to". Are customNonBondedForce and nonBondedForce mutally exclusive with respect with force objects that can be added to a system or has my poor python coding ability let me down again?

Thanks,

Mark

[Mark Williamson]

Hi,

Is there any update on this?

Thanks,

Mark

[Peter Eastman]

Hi Mark,

Sorry for the delay. I was on vacation, and just got back yesterday.

You can use a per-particle parameter to specify which particles should be modified based on lambda. So you can do something like this:

customNonBondedForce.addGlobalParameter("lambda", 1.0)
customNonBondedForce.addPerParticleParameter("useLambda")

Then in your energy expression, replace every "lambda" by "l12" and specify "l12=1-(1-lambda)*step(useLambda1+useLambda2-0.5)". Finally, set useLambda to 1 for particles that lambda should be applied to and 0 to all others. If useLambda1 and useLambda2 are both 0, then step(useLambda1+useLambda2-0.5) = step(-0.5) = 0, and l12 = 1. If either useLambda1 or useLambda2 is 1, then step(useLambda1+useLambda2-0.5) = 1 and l12 = 1-(1-lambda) = lambda.

Peter

[Mark Williamson]

Hi Peter,

Thank you for your reply... I seem to have it working. For reference, here is my toy example that lambda maps a VdW potential from 1 to 0 for two atoms, thus resulting in their dissociation.

Regards,

Mark

Code: Select all

# Example illustrating use of CustomNonbondedForce to lambda map out an Atom
# Based upon argon-chemical-potential.py

from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout
import numpy
from array import *


# =============================================================================
# Specify simulation parameters
# =============================================================================

nparticles = 2 # number of particles

mass = 39.9 * amu # mass
sigma = 3.4 * angstrom # Lennard-Jones sigma
epsilon = 0.238 * kilocalories_per_mole # Lennard-Jones well-depth

timestep = 1 * femtosecond # integrator timestep


# =============================================================================
# Compute box size.
# =============================================================================

box_edge = 100 * angstrom
cutoff = 8 * angstrom
print "sigma = %s" % sigma
print "box_edge = %s" % box_edge
print "cutoff = %s" % cutoff

# =============================================================================
# Build system
# =============================================================================

# Create argon system where first particle is alchemically modified by lambda_value.
lambda_value = 1.0

system = System()
system.setDefaultPeriodicBoxVectors(Vec3(box_edge, 0, 0), Vec3(0, box_edge, 0), Vec3(0, 0, box_edge))

topology = Topology()
newChain = topology.addChain()


# Modified VdW; see page E in http://dx.doi.org/10.1021/ct300857j
customNonBondedForce = CustomNonbondedForce("4*epsilon*l12*( 1/( (alphaLJ*(1-l12) + (r/sigma)^6)^2) - 1/( alphaLJ*(1-l12) + (r/sigma)^6) ) ;sigma=0.5*(sigma1*sigma2); epsilon=sqrt(epsilon1*epsilon2); alphaLJ=0.5; l12=1-(1-lambda)*step(useLambda1+useLambda2-0.5)");

customNonBondedForce.addPerParticleParameter("sigma")
customNonBondedForce.addPerParticleParameter("epsilon")

customNonBondedForce.addGlobalParameter("lambda", 1.0)
customNonBondedForce.addPerParticleParameter("useLambda")

customNonBondedForce.setNonbondedMethod(NonbondedForce.CutoffPeriodic)
customNonBondedForce.setCutoffDistance(cutoff)


# Assign force types
for particle_index in range(nparticles):
  system.addParticle(mass)
  newResidue = topology.addResidue("UNK", newChain)

  if (particle_index == 0 ):
     # Add alchemically-modified particle.
     topology.addAtom("MODD", "Ar", newResidue)
     # Is there a better way to do this?
     customNonBondedForce.addParticle( array('d',[3.4*NmPerAngstrom, 0.238*KJPerKcal, 1 ])  )
     
  else:
     # Add normal particle.
     topology.addAtom("Argo", "ar", newResidue)
     customNonBondedForce.addParticle( array('d',[3.4*NmPerAngstrom, 0.238*KJPerKcal, 0 ])  )

system.addForce(customNonBondedForce)

print "System.getNumParticles() is %i"  % system.getNumParticles()
print "system.getNumForces() is %i " % system.getNumForces()


# Create initial positions; slightly perturbed from equilibrium.
positions = [  Vec3(0,0,0), Vec3(0,0,0.1) ]


# Create Integrator
integrator = VerletIntegrator(1*femtosecond)
simulation = Simulation(topology, system, integrator)

# Create initial positions.
positions = [  Vec3(0,0,0), Vec3(0,0,0.1) ]
simulation.context.setPositions(positions)

print "simulation.system.getNumForces() is %i " % simulation.system.getNumForces()

# Run dynamics.
simulation.reporters.append( PDBReporter('output.pdb', 100) )
simulation.reporters.append( StateDataReporter(stdout, 100, step=True, potentialEnergy=True) )


# =============================================================================
# Lambda map in the Argon atom
# =============================================================================

# This is NVE; kinetic energy has come from the slight perturbation of the interatomic distance
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]

It's great that's working for you!

Here are a couple of suggestions for improvements to your code. First, instead of multiplying constants by conversion factors, it's better to just specify units for them and let OpenMM take care of the conversion. So

customNonBondedForce.addParticle( array('d',[3.4*NmPerAngstrom, 0.238*KJPerKcal, 1 ]) )

can be changed to

customNonBondedForce.addParticle( array('d',[3.4*angstroms, 0.238*kilocalories, 1 ]) )

That's clearer and easier to read. OpenMM will automatically convert any values with units into its internal unit system. So angstroms will be converted to nm and kcal to kJ.

Second, the second argument to Topology.addAtom() should be an Element object, not a string. So

topology.addAtom("Argo", "ar", newResidue)

should be changed to

topology.addAtom("Argo", element.argon, newResidue)

Peter

[Mark Williamson]

Thanks for the tips. I did try your initial tip, but I received the following error:

Code: Select all

Traceback (most recent call last):
  File "two_atom.py", line 71, in <module>
    customNonBondedForce.addParticle( array('d',[3.4*angstrom, 0.238*kilocalories, 1 ])  )
TypeError: nb_float should return float object

I resolved this by using a plain old list in Python

Code: Select all

 customNonBondedForce.addParticle( [3.4*angstrom, 0.238*kilocalories, 1 ]  )

The other piece of advice worked fine; thanks again.

Mark