Skip to content

GitLab

Commit 65a581bf authored by peastman's avatar peastman
Browse files

Per-atom parameters can be taken from residue templates

parent d1afc976
Show whitespace changes
Inline Side-by-side
Showing with 220 additions and 119 deletions
wrappers/python/simtk/openmm/app/forcefield.py
View file @ 65a581bf
......@@ -150,7 +150,11 @@ class ForceField(object):
resName = residue.attrib['name']
template = ForceField._TemplateData(resName)
for atom in residue.findall('Atom'):
template.atoms.append(ForceField._TemplateAtomData(atom.attrib['name'], atom.attrib['type'], self._atomTypes[atom.attrib['type']][2]))
params = {}
for key in atom.attrib:
if key not in ('name', 'type'):
params[key] = _convertParameterToNumber(atom.attrib[key])
template.atoms.append(ForceField._TemplateAtomData(atom.attrib['name'], atom.attrib['type'], self._atomTypes[atom.attrib['type']][2], params))
for site in residue.findall('VirtualSite'):
template.virtualSites.append(ForceField._VirtualSiteData(site))
for bond in residue.findall('Bond'):
......@@ -256,6 +260,7 @@ class ForceField(object):
"""Inner class used to encapsulate data about the system being created."""
def __init__(self):
self.atomType = {}
self.atomParameters = {}
self.atoms = []
self.excludeAtomWith = []
self.virtualSites = {}
......@@ -282,10 +287,11 @@ class ForceField(object):
class _TemplateAtomData:
"""Inner class used to encapsulate data about an atom in a residue template definition."""
def __init__(self, name, type, element):
def __init__(self, name, type, element, parameters={}):
self.name = name
self.type = type
self.element = element
self.parameters = parameters
self.bondedTo = []
self.externalBonds = 0
......@@ -325,6 +331,77 @@ class ForceField(object):
else:
self.excludeWith = self.atoms[0]
class _AtomTypeParameters:
"""Inner class used to record parameter values for atom types."""
def __init__(self, forcefield, forceName, atomTag, paramNames):
self.ff = forcefield
self.forceName = forceName
self.atomTag = atomTag
self.paramNames = paramNames
self.paramsForType = {}
self.extraParamsForType = {}
def registerAtom(self, parameters, expectedParams=None):
if expectedParams is None:
expectedParams = self.paramNames
types = self.ff._findAtomTypes(parameters, 1)
if None not in types:
values = {}
extraValues = {}
for key in parameters:
if key in expectedParams:
values[key] = _convertParameterToNumber(parameters[key])
else:
extraValues[key] = parameters[key]
if len(values) < len(expectedParams):
for key in expectedParams:
if key not in values:
raise ValueError('%s: No value specified for "%s"' % (self.forceName, key))
for t in types[0]:
self.paramsForType[t] = values
self.extraParamsForType[t] = extraValues
def parseDefinitions(self, element):
""""Load the definitions from an XML element."""
expectedParams = list(self.paramNames)
excludedParams = [node.attrib['name'] for node in element.findall('UseAttributeFromResidue')]
for param in excludedParams:
if param not in expectedParams:
raise ValueError('%s: <UseAttributeFromResidue> specified an invalid attribute: %s' % (self.forceName, param))
expectedParams.remove(param)
for atom in element.findall(self.atomTag):
for param in excludedParams:
if param in atom.attrib:
raise ValueError('%s: The attribute "%s" appeared in both <%s> and <UseAttributeFromResidue> tags' % (self.forceName, param, self.atomTag))
self.registerAtom(atom.attrib, expectedParams)
def getAtomParameters(self, atom, data):
"""Get the parameter values for a particular atom."""
t = data.atomType[atom]
p = data.atomParameters[atom]
if t in self.paramsForType:
values = self.paramsForType[t]
result = [None]*len(self.paramNames)
for i, name in enumerate(self.paramNames):
if name in values:
result[i] = values[name]
elif name in p:
result[i] = p[name]
else:
raise ValueError('%s: No value specified for "%s"' % (self.forceName, name))
return result
else:
raise ValueError('%s: No parameters defined for atom type %s' % (self.forceName, t))
def getExtraParameters(self, atom, data):
"""Get extra parameter values for an atom that appeared in the <Atom> tag but were not included in paramNames."""
t = data.atomType[atom]
if t in self.paramsForType:
return self.extraParamsForType[t]
else:
raise ValueError('%s: No parameters defined for atom type %s' % (self.forceName, t))
def createSystem(self, topology, nonbondedMethod=NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
constraints=None, rigidWater=True, removeCMMotion=True, hydrogenMass=None, **args):
"""Construct an OpenMM System representing a Topology with this force field.
......@@ -385,6 +462,7 @@ class ForceField(object):
matchAtoms = dict(zip(matches, res.atoms()))
for atom, match in zip(res.atoms(), matches):
data.atomType[atom] = template.atoms[match].type
data.atomParameters[atom] = template.atoms[match].parameters
for site in template.virtualSites:
if match == site.index:
data.virtualSites[atom] = (site, [matchAtoms[i].index for i in site.atoms], matchAtoms[site.excludeWith].index)
......@@ -1139,13 +1217,10 @@ class NonbondedGenerator:
self.coulomb14scale = coulomb14scale
self.lj14scale = lj14scale
self.typeMap = {}
self.params = ForceField._AtomTypeParameters(forcefield, 'NonbondedForce', 'Atom', ('charge', 'sigma', 'epsilon'))
def registerAtom(self, parameters):
types = self.ff._findAtomTypes(parameters, 1)
if None not in types:
values = (_convertParameterToNumber(parameters['charge']), _convertParameterToNumber(parameters['sigma']), _convertParameterToNumber(parameters['epsilon']))
for t in types[0]:
self.typeMap[t] = values
self.params.registerAtom(parameters)
@staticmethod
def parseElement(element, ff):
......@@ -1158,8 +1233,7 @@ class NonbondedGenerator:
generator = existing[0]
if generator.coulomb14scale != float(element.attrib['coulomb14scale']) or generator.lj14scale != float(element.attrib['lj14scale']):
raise ValueError('Found multiple NonbondedForce tags with different 1-4 scales')
for atom in element.findall('Atom'):
generator.registerAtom(atom.attrib)
generator.params.parseDefinitions(element)
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
methodMap = {NoCutoff:mm.NonbondedForce.NoCutoff,
......@@ -1171,12 +1245,8 @@ class NonbondedGenerator:
raise ValueError('Illegal nonbonded method for NonbondedForce')
force = mm.NonbondedForce()
for atom in data.atoms:
t = data.atomType[atom]
if t in self.typeMap:
values = self.typeMap[t]
values = self.params.getAtomParameters(atom, data)
force.addParticle(values[0], values[1], values[2])
else:
raise ValueError('No nonbonded parameters defined for atom type '+t)
force.setNonbondedMethod(methodMap[nonbondedMethod])
force.setCutoffDistance(nonbondedCutoff)
if 'ewaldErrorTolerance' in args:
......@@ -1226,13 +1296,10 @@ class GBSAOBCGenerator:
def __init__(self, forcefield):
self.ff = forcefield
self.typeMap = {}
self.params = ForceField._AtomTypeParameters(forcefield, 'GBSAOBCForce', 'Atom', ('charge', 'radius', 'scale'))
def registerAtom(self, parameters):
types = self.ff._findAtomTypes(parameters, 1)
if None not in types:
values = (_convertParameterToNumber(parameters['charge']), _convertParameterToNumber(parameters['radius']), _convertParameterToNumber(parameters['scale']))
for t in types[0]:
self.typeMap[t] = values
self.params.registerAtom(parameters)
@staticmethod
def parseElement(element, ff):
......@@ -1243,8 +1310,7 @@ class GBSAOBCGenerator:
else:
# Multiple <GBSAOBCForce> tags were found, probably in different files. Simply add more types to the existing one.
generator = existing[0]
for atom in element.findall('Atom'):
generator.registerAtom(atom.attrib)
generator.params.parseDefinitions(element)
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
methodMap = {NoCutoff:mm.NonbondedForce.NoCutoff,
......@@ -1254,12 +1320,8 @@ class GBSAOBCGenerator:
raise ValueError('Illegal nonbonded method for GBSAOBCForce')
force = mm.GBSAOBCForce()
for atom in data.atoms:
t = data.atomType[atom]
if t in self.typeMap:
values = self.typeMap[t]
values = self.params.getAtomParameters(atom, data)
force.addParticle(values[0], values[1], values[2])
else:
raise ValueError('No GBSAOBC parameters defined for atom type '+t)
force.setNonbondedMethod(methodMap[nonbondedMethod])
force.setCutoffDistance(nonbondedCutoff)
if 'soluteDielectric' in args:
......@@ -1590,12 +1652,8 @@ class CustomNonbondedGenerator:
generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
for param in element.findall('PerParticleParameter'):
generator.perParticleParams.append(param.attrib['name'])
for atom in element.findall('Atom'):
types = ff._findAtomTypes(atom.attrib, 1)
if None not in types:
values = [float(atom.attrib[param]) for param in generator.perParticleParams]
for t in types[0]:
generator.typeMap[t] = values
generator.params = ForceField._AtomTypeParameters(ff, 'CustomNonbondedForce', 'Atom', generator.perParticleParams)
generator.params.parseDefinitions(element)
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
methodMap = {NoCutoff:mm.CustomNonbondedForce.NoCutoff,
......@@ -1622,11 +1680,8 @@ class CustomNonbondedGenerator:
elif type == 'Discrete3D':
force.addTabulatedFunction(name, mm.Discrete2DFunction(params['xsize'], params['ysize'], params['zsize'], values))
for atom in data.atoms:
t = data.atomType[atom]
if t in self.typeMap:
force.addParticle(self.typeMap[t])
else:
raise ValueError('No CustomNonbonded parameters defined for atom type '+t)
values = self.params.getAtomParameters(atom, data)
force.addParticle(values)
force.setNonbondedMethod(methodMap[nonbondedMethod])
force.setCutoffDistance(nonbondedCutoff)
sys.addForce(force)
......@@ -1686,12 +1741,8 @@ class CustomGBGenerator:
generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
for param in element.findall('PerParticleParameter'):
generator.perParticleParams.append(param.attrib['name'])
for atom in element.findall('Atom'):
types = ff._findAtomTypes(atom.attrib, 1)
if None not in types:
values = [float(atom.attrib[param]) for param in generator.perParticleParams]
for t in types[0]:
generator.typeMap[t] = values
generator.params = ForceField._AtomTypeParameters(ff, 'CustomGBForce', 'Atom', generator.perParticleParams)
generator.params.parseDefinitions(element)
computationMap = {"SingleParticle" : mm.CustomGBForce.SingleParticle,
"ParticlePair" : mm.CustomGBForce.ParticlePair,
"ParticlePairNoExclusions" : mm.CustomGBForce.ParticlePairNoExclusions}
......@@ -1742,11 +1793,8 @@ class CustomGBGenerator:
elif type == 'Discrete3D':
force.addTabulatedFunction(name, mm.Discrete2DFunction(params['xsize'], params['ysize'], params['zsize'], values))
for atom in data.atoms:
t = data.atomType[atom]
if t in self.typeMap:
force.addParticle(self.typeMap[t])
else:
raise ValueError('No CustomGB parameters defined for atom type '+t)
values = self.params.getAtomParameters(atom, data)
force.addParticle(values)
force.setNonbondedMethod(methodMap[nonbondedMethod])
force.setCutoffDistance(nonbondedCutoff)
sys.addForce(force)
......@@ -3083,8 +3131,6 @@ class AmoebaVdwGenerator:
self.vdw14Scale = vdw14Scale
self.vdw15Scale = vdw15Scale
self.typeMap = {}
#=============================================================================================
@staticmethod
......@@ -3097,31 +3143,10 @@ class AmoebaVdwGenerator:
generator = AmoebaVdwGenerator(element.attrib['type'], element.attrib['radiusrule'], element.attrib['radiustype'], element.attrib['radiussize'], element.attrib['epsilonrule'],
float(element.attrib['vdw-13-scale']), float(element.attrib['vdw-14-scale']), float(element.attrib['vdw-15-scale']))
forceField._forces.append(generator)
generator.params = ForceField._AtomTypeParameters(forceField, 'AmoebaVdwForce', 'Vdw', ('sigma', 'epsilon', 'reduction'))
generator.params.parseDefinitions(element)
two_six = 1.122462048309372
# types[] = [ sigma, epsilon, reductionFactor, class ]
# sigma is modified based on radiustype and radiussize
for atom in element.findall('Vdw'):
types = forceField._findAtomTypes(atom.attrib, 1)
if None not in types:
values = [float(atom.attrib['sigma']), float(atom.attrib['epsilon']), float(atom.attrib['reduction'])]
if (generator.radiustype == 'SIGMA'):
values[0] *= two_six
if (generator.radiussize == 'DIAMETER'):
values[0] *= 0.5
for t in types[0]:
generator.typeMap[t] = values
else:
outputString = "AmoebaVdwGenerator: error getting type: %s" % (atom.attrib['class'])
raise ValueError(outputString)
#=============================================================================================
# Return a set containing the indices of particles bonded to particle with index=particleIndex
......@@ -3199,28 +3224,25 @@ class AmoebaVdwGenerator:
force = existing[0]
# add particles to force
# throw error if particle type not available
sigmaScale = 1
if self.radiustype == 'SIGMA':
sigmaScale = 1.122462048309372
if self.radiussize == 'DIAMETER':
sigmaScale = 0.5
for (i, atom) in enumerate(data.atoms):
t = data.atomType[atom]
if t in self.typeMap:
values = self.typeMap[t]
values = self.params.getAtomParameters(atom, data)
# ivIndex = index of bonded partner for hydrogens; otherwise ivIndex = particle index
ivIndex = i
mass = sys.getParticleMass(i)/unit.dalton
if (mass < 1.9 and len(data.atomBonds[i]) == 1):
if atom.element == elem.hydrogen and len(data.atomBonds[i]) == 1:
bondIndex = data.atomBonds[i][0]
if (data.bonds[bondIndex].atom1 == i):
ivIndex = data.bonds[bondIndex].atom2
else:
ivIndex = data.bonds[bondIndex].atom1
force.addParticle(ivIndex, values[0], values[1], values[2])
else:
raise ValueError('No vdw type for atom %s' % (atom.name))
force.addParticle(ivIndex, values[0]*sigmaScale, values[1], values[2])
# set combining rules
......@@ -3937,8 +3959,6 @@ class AmoebaWcaDispersionGenerator:
self.dispoff = dispoff
self.shctd = shctd
self.typeMap = {}
#=========================================================================================
@staticmethod
......@@ -3957,19 +3977,8 @@ class AmoebaWcaDispersionGenerator:
element.attrib['dispoff'],
element.attrib['shctd'])
forceField._forces.append(generator)
# typeMap[] = [ radius, epsilon ]
for atom in element.findall('WcaDispersion'):
types = forceField._findAtomTypes(atom.attrib, 1)
if None not in types:
values = [float(atom.attrib['radius']), float(atom.attrib['epsilon'])]
for t in types[0]:
generator.typeMap[t] = values
else:
outputString = "AmoebaWcaDispersionGenerator: error getting type: %s" % (atom.attrib['class'])
raise ValueError(outputString)
generator.params = ForceField._AtomTypeParameters(forceField, 'AmoebaWcaDispersionForce', 'WcaDispersion', ('radius', 'epsilon'))
generator.params.parseDefinitions(element)
#=========================================================================================
......@@ -3997,14 +4006,9 @@ class AmoebaWcaDispersionGenerator:
force.setAwater( float(self.awater ))
force.setShctd( float(self.shctd ))
for (i, atom) in enumerate(data.atoms):
t = data.atomType[atom]
if t in self.typeMap:
values = self.typeMap[t]
for atom in data.atoms:
values = self.params.getAtomParameters(atom, data)
force.addParticle(values[0], values[1])
else:
raise ValueError('No WcaDispersion type for atom %s of %s %d' % (atom.name, atom.residue.name, atom.residue.index))
parsers["AmoebaWcaDispersionForce"] = AmoebaWcaDispersionGenerator.parseElement
......
wrappers/python/tests/TestForceField.py
View file @ 65a581bf
......@@ -6,6 +6,10 @@ from simtk.unit import *
import simtk.openmm.app.element as elem
import simtk.openmm.app.forcefield as forcefield
import math
if sys.version_info >= (3, 0):
from io import StringIO
else:
from cStringIO import StringIO
class TestForceField(unittest.TestCase):
"""Test the ForceField.createSystem() method."""
......@@ -144,7 +148,8 @@ class TestForceField(unittest.TestCase):
context = Context(system, integrator)
context.setPositions(pdb.positions)
state1 = context.getState(getForces=True)
state2 = XmlSerializer.deserialize(open('systems/lysozyme-implicit-forces.xml').read())
with open('systems/lysozyme-implicit-forces.xml') as input:
state2 = XmlSerializer.deserialize(input.read())
numDifferences = 0
for f1, f2, in zip(state1.getForces().value_in_unit(kilojoules_per_mole/nanometer), state2.getForces().value_in_unit(kilojoules_per_mole/nanometer)):
diff = norm(f1-f2)
......@@ -198,6 +203,50 @@ class TestForceField(unittest.TestCase):
self.assertEqual(vectors[i], self.pdb1.topology.getPeriodicBoxVectors()[i])
self.assertEqual(vectors[i], system.getDefaultPeriodicBoxVectors()[i])
def test_ResidueAttributes(self):
"""Test a ForceField that gets per-particle parameters from residue attributes."""
xml = """
<ForceField>
<AtomTypes>
<Type name="tip3p-O" class="OW" element="O" mass="15.99943"/>
<Type name="tip3p-H" class="HW" element="H" mass="1.007947"/>
</AtomTypes>
<Residues>
<Residue name="HOH">
<Atom name="O" type="tip3p-O" charge="-0.834"/>
<Atom name="H1" type="tip3p-H" charge="0.417"/>
<Atom name="H2" type="tip3p-H" charge="0.417"/>
<Bond from="0" to="1"/>
<Bond from="0" to="2"/>
</Residue>
</Residues>
<NonbondedForce coulomb14scale="0.833333" lj14scale="0.5">
<UseAttributeFromResidue name="charge"/>
<Atom type="tip3p-O" sigma="0.315" epsilon="0.635"/>
<Atom type="tip3p-H" sigma="1" epsilon="0"/>
</NonbondedForce>
</ForceField>"""
ff = ForceField(StringIO(xml))
# Build a water box.
modeller = Modeller(Topology(), [])
modeller.addSolvent(ff, boxSize=Vec3(3, 3, 3)*nanometers)
# Create a system and make sure all nonbonded parameters are correct.
system = ff.createSystem(modeller.topology)
nonbonded = [f for f in system.getForces() if isinstance(f, NonbondedForce)][0]
atoms = list(modeller.topology.atoms())
for i in range(len(atoms)):
params = nonbonded.getParticleParameters(i)
if atoms[i].element == elem.oxygen:
self.assertEqual(params[0], -0.834*elementary_charge)
self.assertEqual(params[1], 0.315*nanometers)
self.assertEqual(params[2], 0.635*kilojoule_per_mole)
else:
self.assertEqual(params[0], 0.417*elementary_charge)
self.assertEqual(params[1], 1.0*nanometers)
self.assertEqual(params[2], 0.0*kilojoule_per_mole)
class AmoebaTestForceField(unittest.TestCase):
"""Test the ForceField.createSystem() method with the AMOEBA forcefield."""
......@@ -270,6 +319,22 @@ class AmoebaTestForceField(unittest.TestCase):
self.assertAlmostEqual(constraints[(0,2)], hoDist)
self.assertAlmostEqual(constraints[(1,2)], hohDist)
def test_Forces(self):
"""Compute forces and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""
pdb = PDBFile('systems/alanine-dipeptide-implicit.pdb')
forcefield = ForceField('amoeba2013.xml', 'amoeba2013_gk.xml')
system = forcefield.createSystem(pdb.topology, polarization='direct')
integrator = VerletIntegrator(0.001)
context = Context(system, integrator)
context.setPositions(pdb.positions)
state1 = context.getState(getForces=True)
with open('systems/alanine-dipeptide-amoeba-forces.xml') as input:
state2 = XmlSerializer.deserialize(input.read())
for f1, f2, in zip(state1.getForces().value_in_unit(kilojoules_per_mole/nanometer), state2.getForces().value_in_unit(kilojoules_per_mole/nanometer)):
diff = norm(f1-f2)
self.assertTrue(diff < 0.1 or diff/norm(f1) < 1e-3)
if __name__ == '__main__':
unittest.main()
wrappers/python/tests/systems/alanine-dipeptide-amoeba-forces.xml 0 → 100644
View file @ 65a581bf
<?xml version="1.0" ?>
<State time="0" type="State" version="1">
<PeriodicBoxVectors>
<A x="2" y="0" z="0"/>
<B x="0" y="2" z="0"/>
<C x="0" y="0" z="2"/>
</PeriodicBoxVectors>
<Forces>
<Force x="-60.899006905648534" y="-704.8991281062977" z=".0364180169671772"/>
<Force x="1179.7263311404506" y="629.2168542831286" z="1.563821009034657"/>
<Force x="-484.4401097634541" y="135.94210610096957" z="558.0155391926069"/>
<Force x="-484.1625627429897" y="136.02915381525668" z="-557.7203977673528"/>
<Force x="-804.8403408476064" y="-857.579131840898" z="118.43755474210333"/>
<Force x="84.24441920455223" y="-421.86322192608037" z="120.734742628977"/>
<Force x="666.4942549112259" y="-479.03339298154526" z="19.66166555862896"/>
<Force x="-518.9062723002917" y="686.3850843706331" z=".7311137123231983"/>
<Force x="-105.62094199225953" y="710.9845370176521" z="-508.0677463920413"/>
<Force x="285.6690654282444" y="-111.35548608152155" z="847.8096733172159"/>
<Force x="-396.4822698774795" y="24.462223992557014" z="297.9581483103431"/>
<Force x="-168.22340166944585" y="110.99779797739818" z="-824.0916382722021"/>
<Force x="805.7023041833917" y="15.738288953519556" z="-340.2263402925698"/>
<Force x="366.3025998383222" y="-647.7535895818237" z="-273.83047817845636"/>
<Force x="-649.7310636705703" y="-495.52797252280067" z="317.36013856851315"/>
<Force x="-480.120788018345" y="172.44477988154952" z="-16.235455767965004"/>
<Force x="-486.20230228601736" y="982.6376366390087" z="127.20818643942287"/>
<Force x="820.5730252078641" y="-78.71234329038225" z="109.31395566365498"/>
<Force x="355.03388135028814" y="-1391.2776654673507" z=".9942543092954484"/>
<Force x="-511.3467906691061" y="677.7717086552518" z=".008253291392112616"/>
<Force x="293.66157422204157" y="452.68343054060256" z="502.83463917630775"/>
<Force x="293.5683952568295" y="452.70832957117113" z="-502.49604726620026"/>
</Forces>
</State>
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment