Specific atom type pair nonboned forces in XML? (Interaction groups?)

[George Pantelopulos]

Hi all,

Is it possible to write NonbondedForces which are specific to a pair of atom types or more in force field XML files? In particular, I am trying to set different 12-6 LJ parameters for self-self type interactions and self-other type interactions.

From reading up on github it seems that InteractionGroups with a CustomNonbondedForce is how achieving this kind of goal has been implemented, but I haven't identified how this is done with XML force field files. Is this possible?

Thank you for the help!

[Peter Eastman]

Generally the best way to do this is to have a per-particle atom type parameter, then use those to look up force parameters in a table. To do that in an XML file, you currently need to include a <Script> tag. For an example of doing this, take a look at https://raw.githubusercontent.com/pande ... r_2013.xml and scroll down to the very end.

OpenMM 7.1 will include a new <LennardJonesForce> tag to simplify this case.

Peter

[George Pantelopulos]

Dear Peter,

I realize this might be a bit of a silly question to ask to revive this thread, but I'd been using the custom nonbonded forces with 2D tabulated functions successfully for this purpose, however I was using this for LJ particles I was initializing within a python script, and not as a force field xml file.

I want to use a force field file to assign the parameters and to be able to build a Simulation object, instead of running the system through a context and integrator, but strangely it seems that, upon trying to build the system, the residues in an input PDB file are not recognized (error: No template found for residue 1 (AYY). This might mean your input topology is missing some atoms or bonds, or possibly that you are using the wrong force field.), even though these residues are very clearly defined in the xml file. I am clueless as to the cause of this.

Trying to load a pdb file which looks something like so:

ATOM 1 A AYY A 1 7.671 23.131 4.217
ATOM 2 B BEE A 2 17.532 26.381 2.292

And trying to build a system from a force field which looks like so:

Code: Select all

<ForceField>
 <AtomTypes>
  <Type name="0" class="LJA" element="Ar" mass="39.948"/>
  <Type name="1" class="LJB" element="Ar" mass="39.948"/>
 </AtomTypes>
 <Residues>
  <Residue name="AYY">
   <Atom name="A" type="0"/>
  </Residue>
  <Residue name="BEE">
   <Atom name="B" type="1"/>
  </Residue>
 </Residues>
 
 <Script>
import simtk.openmm as mm
nonbondedForce = [f for f in [sys.getForce(i) for i in range(sys.getNumForces())] if type(f) == mm.NonbondedForce][0]

# making self-interaction well 4x deeper
epsilon = [
3.9832, 0.9958,   # A-A, A-B
0.9958, 3.9832]   # B-A, B-B

sigma = [
0.3405, 0.3405,   # A-A, A-B
0.3405, 0.3405]   # B-A, B-B

# Create a CustomNonbondedForce to compute Lennard-Jones interactions.

customNonbondedForce = mm.CustomNonbondedForce('4*eps*((sig/r)^12-(sig/r)^6); eps=epsilon(type1, type2); sig=sigma(type1, type2)')
customNonbondedForce.setNonbondedMethod(min(nonbondedForce.getNonbondedMethod(), 2))
customNonbondedForce.setCutoffDistance(nonbondedForce.getCutoffDistance())
customNonbondedForce.addTabulatedFunction('epsilon', mm.Discrete2DFunction(2, 2, epsilon))
customNonbondedForce.addTabulatedFunction('sigma', mm.Discrete2DFunction(2, 2, sigma))
customNonbondedForce.addPerParticleParameter('type')
for atom in data.atoms:
    customNonbondedForce.addParticle([data.atomType[atom]])
sys.addForce(customNonbondedForce)

 </Script>
</ForceField>

[Peter Eastman]

The problem is that your PDB file is missing the element column. So it doesn't know those atoms are supposed to be argon, and therefore can't match them against the template. See http://docs.openmm.org/7.0.0/userguide/ ... -templates for information on how template matching is done, and http://www.wwpdb.org/documentation/file ... .html#ATOM for the format of ATOM records in PDB files.

Peter

[George Pantelopulos]

Hi Peter,

Oh, I see thank you very much for pointing that out for me! Fixing that bit of business

Another question on using scripts in force field xml files, as must be done for pair-specific LJ interactions:
*How do I add particles to the force based on their type as interpreted by templating?
*How do I ensure that the force's parameters will be set to a particular nonbonded method and nonbonded cutoff specified when using the loaded force field to create a simulation (to set the method to CutoffPeriodic and the cutoff to 1.0 nm when doing forcefield.createSystem(pdb.topology, nonbondedMethod=app.CutoffPeriodic, nonbondedCutoff=1.0*unit.nanometers)?

As of right now, trying something like the following:

Code: Select all

<ForceField>
 <AtomTypes>
  <Type name="0" class="LJA" element="Ar" mass="39.948"/>
  <Type name="1" class="LJB" element="Ar" mass="39.948"/>
 </AtomTypes>
 <Residues>
  <Residue name="AYY">
   <Atom name="A" type="0"/>
  </Residue>
  <Residue name="BEE">
   <Atom name="B" type="1"/>
  </Residue>
 </Residues>
 
 <Script>
import simtk.openmm as mm
import simtk.openmm.app as app
import simtk.unit as u
cutoff = 2.5 * u.nanometers

# making self-interaction well 4x deeper
epsilon = [
3.9832, 0.9958,   # A-A, A-B
0.9958, 3.9832]   # B-A, B-B

sigma = [
0.3405, 0.3405,   # A-A, A-B
0.3405, 0.3405]   # B-A, B-B

# Create a CustomNonbondedForce to compute Lennard-Jones interactions.

customNonbondedForce = mm.CustomNonbondedForce('4*eps*((sig/r)^12-(sig/r)^6); eps=epsilon(type1, type2); sig=sigma(type1, type2)')
customNonbondedForce.setNonbondedMethod(mm.NonbondedForce.CutoffNonPeriodic)
customNonbondedForce.setCutoffDistance(cutoff)
customNonbondedForce.addTabulatedFunction('epsilon', mm.Discrete2DFunction(2, 2, epsilon))
customNonbondedForce.addTabulatedFunction('sigma', mm.Discrete2DFunction(2, 2, sigma))
customNonbondedForce.addPerParticleParameter('type')
for atom in data.atoms:
    customNonbondedForce.addParticle([data.atomType[atom]])
sys.addForce(customNonbondedForce)

 </Script>
</ForceField>

Gives the following error:

Traceback (most recent call last):
File "teststruct.py", line 7, in <module>
forcefield.createSystem(pdb.topology, nonbondedMethod=app.CutoffNonPeriodic, nonbondedCutoff=1.0*u.nanometers)
File "/Users/gpantel/miniconda2/lib/python2.7/site-packages/simtk/openmm/app/forcefield.py", line 956, in createSystem
exec(script, locals())
File "<string>", line 25, in <module>
File "/Users/gpantel/miniconda2/lib/python2.7/site-packages/simtk/openmm/openmm.py", line 11554, in addParticle
return _openmm.CustomNonbondedForce_addParticle(self, *args)
NotImplementedError: Wrong number or type of arguments for overloaded function 'CustomNonbondedForce_addParticle'.
Possible C/C++ prototypes are:
OpenMM::CustomNonbondedForce::addParticle(std::vector< double,std::allocator< double > > const &)
OpenMM::CustomNonbondedForce::addParticle(

[Peter Eastman]

Per-particle parameters for a CustomNonbondedForce are always numeric. So that's what addParticle() is expecting for its argument: a list of numbers. But data.atomType[atom] is an object, not a number. That's why it's giving you the error:

NotImplementedError: Wrong number or type of arguments for overloaded function 'CustomNonbondedForce_addParticle'.

So you need to map from atom type objects to the indices used to look up parameters in your tables.

How do I ensure that the force's parameters will be set to a particular nonbonded method and nonbonded cutoff specified when using the loaded force field

The nonbondedMethod and nonbondedCutoff variables will hold the corresponding arguments that were passed to createSystem(). You can write something along these lines:

Code: Select all

if nonbondedMethod in [CutoffPeriodic, Ewald, PME]:
    force.setNonbondedMethod(CustomNonbondedForce.CutoffPeriodic)
elif nonbondedMethod is NoCutoff:
    force.setNonbondedMethod(CustomNonbondedForce.NoCutoff)
elif nonbondedMethod is CutoffNonPeriodic:
    force.setNonbondedMethod(CustomNonbondedForce.CutoffNonPeriodic)
else:
    raise ValueError('Unrecognized nonbonded method [%s]' % nonbondedMethod)
force.setCutoffDistance(nonbondedCutoff)

[George Pantelopulos]

Thank you once again, Peter!