CustomNonbondedForce + NBFix

The functionality of OpenMM will (eventually) include everything that one would need to run modern molecular simulation.
POST REPLY
User avatar
Daniel Roston
Posts: 5
Joined: Mon Sep 12, 2016 5:03 pm

CustomNonbondedForce + NBFix

Post by Daniel Roston » Thu Feb 14, 2019 3:16 pm

I'm running into some problems using CustomNonbondedForce for alchemical FEP in a system that uses charmm with NBFix. From what I can tell, openmm's solution to NBFix automatically uses CustomNonbondedForce. But then the FEP using CustomNonbondedForce doesn't do what it's supposed to do. I've been using this example (http://openmm.org/tutorials/alchemical-free-energy/) as the basis for the FEP in similar systems that don't involve NBFix and it works great. And that script works fine if i comment out the NBFix lines in the charmm parameter files. but I kind of want to use NBFix. Is there some workaround?

User avatar
Peter Eastman
Posts: 1851
Joined: Thu Aug 09, 2007 1:25 pm

Re: CustomNonbondedForce + NBFix

Post by Peter Eastman » Thu Feb 14, 2019 4:08 pm

You can do it in a similar way, but the details are a little different. Quoting from the page you linked, it defines the energy function as

Code: Select all

energy_function = 'lambda*4*epsilon*x*(x-1.0); x = (sigma/reff_sterics)^6;'
energy_function += 'reff_sterics = sigma*(0.5*(1.0-lambda) + (r/sigma)^6)^(1/6);'
energy_function += 'sigma = 0.5*(sigma1+sigma2); epsilon = sqrt(epsilon1*epsilon2);'
Note the third line: it computes sigma and epsilon by applying combination rules to per-particle parameters.

Now take a look at how CharmmPsfFile defines its CustomNonbondedForce:

Code: Select all

cforce = mm.CustomNonbondedForce('(a/r6)^2-b/r6; r6=r^6;'
                                 'a=acoef(type1, type2);'
                                 'b=bcoef(type1, type2)')
It writes the Lennard-Jones potential in a slightly different form, using the coefficients A and B instead of sigma and epsilon. It's trivial to convert between the two though. See https://en.wikipedia.org/wiki/Lennard-J ... xpressions. More important, it doesn't use combination rules. Instead, it looks up the coefficients in a pair of tables based on the types of the two atoms. That's necessary when using NBFix.

So you can mostly use the same energy function as in the example, but you'll need to modify the last line. Instead of computing sigma and epsilon with combination rules, you'll instead look up A and B from the tables and use those to compute sigma and epsilon.

User avatar
Daniel Roston
Posts: 5
Joined: Mon Sep 12, 2016 5:03 pm

Re: CustomNonbondedForce + NBFix

Post by Daniel Roston » Fri Feb 15, 2019 1:21 pm

Thanks, this helps. But how do I look up the A and B parameters? when i use:

cforce.getFunctionParameters(0)

i get

Exception: CustomNonbondedForce: function is not a Continuous1DFunction

what's the proper way to do it?

User avatar
Peter Eastman
Posts: 1851
Joined: Thu Aug 09, 2007 1:25 pm

Re: CustomNonbondedForce + NBFix

Post by Peter Eastman » Fri Feb 15, 2019 1:49 pm

getFunctionParameters() is deprecated. Use getTabulatedFunction() instead.

User avatar
Daniel Roston
Posts: 5
Joined: Mon Sep 12, 2016 5:03 pm

Re: CustomNonbondedForce + NBFix

Post by Daniel Roston » Tue Feb 19, 2019 12:41 pm

Hi,

Still having some trouble here...I'm pretty new to openmm (and python!) so thanks for bearing with me. When I try:

custom_force.addTabulatedFunction('acoef',openmm.Discrete2DFunction(cforce.getTabulatedFunction(0)))

I get:

TypeError: __init__() missing 2 required positional arguments: 'ysize' and 'values'

how should i be using getTabulatedFunction?

User avatar
Peter Eastman
Posts: 1851
Joined: Thu Aug 09, 2007 1:25 pm

Re: CustomNonbondedForce + NBFix

Post by Peter Eastman » Tue Feb 19, 2019 1:03 pm

Try this:

Code: Select all

from copy import deepcopy
custom_force.addTabulatedFunction('acoef',deepcopy(cforce.getTabulatedFunction(0)))

User avatar
Daniel Roston
Posts: 5
Joined: Mon Sep 12, 2016 5:03 pm

Re: CustomNonbondedForce + NBFix

Post by Daniel Roston » Thu Feb 21, 2019 2:14 pm

Hi,

I think deepcopy is working, but there still seems to be some numerical issue. When I turn off VDW for a *single* tip3p water, the PE of the system changes by ca. 5000 kj/mol. Should there be some sort of unit conversion factor somewhere? here's the relevant section of the script that I'm using:

# Retrieve the NonbondedForce
forces = { force.__class__.__name__ : force for force in system_nbfix.getForces() }
nbforce = forces['NonbondedForce']
cforce = forces['CustomNonbondedForce']

# Add a CustomNonbondedForce to handle only alchemically-modified interactions
# set charge of alcemical particles to zero, so don't need that component of energy_function
alchemical_particles = set([0,1,2])
chemical_particles = set(range(system_nbfix.getNumParticles())) - alchemical_particles
energy_function = 'lambda*4*epsilon*x*(x-1.0); x = (sigma/reff_sterics)^6;'
energy_function += 'reff_sterics = sigma*(0.5*(1.0-lambda) + (r/sigma)^6)^(1/6);'
energy_function += 'sigma = (a/b)^(1/6) ; epsilon = (b^2)/(4*a);'
energy_function += 'a=acoef(type1, type2); b=bcoef(type1, type2)'

custom_force = openmm.CustomNonbondedForce(energy_function)

lambda_sim=1.00

from copy import deepcopy
custom_force.addTabulatedFunction('acoef',deepcopy(cforce.getTabulatedFunction(0)))
custom_force.addTabulatedFunction('bcoef',deepcopy(cforce.getTabulatedFunction(1)))
custom_force.addGlobalParameter('lambda',lambda_sim)
custom_force.addPerParticleParameter('type')

custom_force.setNonbondedMethod(custom_force.CutoffNonPeriodic)
custom_force.setCutoffDistance(1.2*nanometer)

for atom in psfsim.topology.atoms():
#get atom types from cforce
par1 = cforce.getParticleParameters(atom.index)
custom_force.addParticle(par1)
#turn off nbforce for alchemical particles
if atom.index in alchemical_particles:
[charge, sigma, epsilon] = nbforce.getParticleParameters(atom.index)
nbforce.setParticleParameters(atom.index,0*charge,sigma,0*epsilon)

custom_force.addInteractionGroup(alchemical_particles, chemical_particles)

system_nbfix.addForce(custom_force)

User avatar
Daniel Roston
Posts: 5
Joined: Mon Sep 12, 2016 5:03 pm

Re: CustomNonbondedForce + NBFix

Post by Daniel Roston » Tue Feb 26, 2019 3:00 pm

Ok, I think I got the numerical issue with the alchemical force sorted, but now I need to turn off the non-alchemical LJ force between the alchemical particle and chemical particles. but when i do

cforce.addExclusion(atom1.index,atom2.index)

I end up getting:

Exception: All Forces must have identical exceptions

is there another way to do it?

POST REPLY