Skip to content

GitLab

Commit b2212517 authored by peastman's avatar peastman
Browse files

Merge pull request #506 from peastman/forcefield 

Began creating API for programmatically building and editing ForceFields
parents e62739d2 9bce462d
Hide whitespace changes
Inline Side-by-side
Showing with 220 additions and 135 deletions
wrappers/python/simtk/openmm/app/forcefield.py
View file @ b2212517
...@@ -40,6 +40,11 @@ import simtk.unit as unit ...@@ -40,6 +40,11 @@ import simtk.unit as unit
import element as elem import element as elem
from simtk.openmm.app import Topology from simtk.openmm.app import Topology
def _convertParameterToNumber(param):
if unit.is_quantity(param):
return mm.stripUnits((param,))[0]
return float(param)
# Enumerated values for nonbonded method # Enumerated values for nonbonded method
class NoCutoff(object): class NoCutoff(object):
...@@ -115,13 +120,7 @@ class ForceField(object): ...@@ -115,13 +120,7 @@ class ForceField(object):
if tree.getroot().find('AtomTypes') is not None: if tree.getroot().find('AtomTypes') is not None:
for type in tree.getroot().find('AtomTypes').findall('Type'): for type in tree.getroot().find('AtomTypes').findall('Type'):
element = None self.registerAtomType(type.attrib)
if 'element' in type.attrib:
element = elem.get_by_symbol(type.attrib['element'])
name = type.attrib['name']
if name in self._atomTypes:
raise ValueError('Found multiple definitions for atom type: '+name)
self._atomTypes[name] = (type.attrib['class'], float(type.attrib['mass']), element)
# Load the residue templates. # Load the residue templates.
...@@ -129,38 +128,17 @@ class ForceField(object): ...@@ -129,38 +128,17 @@ class ForceField(object):
for residue in root.find('Residues').findall('Residue'): for residue in root.find('Residues').findall('Residue'):
resName = residue.attrib['name'] resName = residue.attrib['name']
template = ForceField._TemplateData(resName) template = ForceField._TemplateData(resName)
self._templates[resName] = template
for atom in residue.findall('Atom'): for atom in residue.findall('Atom'):
template.atoms.append(ForceField._TemplateAtomData(atom.attrib['name'], atom.attrib['type'], self._atomTypes[atom.attrib['type']][2])) template.atoms.append(ForceField._TemplateAtomData(atom.attrib['name'], atom.attrib['type'], self._atomTypes[atom.attrib['type']][2]))
for site in residue.findall('VirtualSite'): for site in residue.findall('VirtualSite'):
template.virtualSites.append(ForceField._VirtualSiteData(site)) template.virtualSites.append(ForceField._VirtualSiteData(site))
for bond in residue.findall('Bond'): for bond in residue.findall('Bond'):
b = (int(bond.attrib['from']), int(bond.attrib['to'])) template.addBond(int(bond.attrib['from']), int(bond.attrib['to']))
template.bonds.append(b)
template.atoms[b[0]].bondedTo.append(b[1])
template.atoms[b[1]].bondedTo.append(b[0])
for bond in residue.findall('ExternalBond'): for bond in residue.findall('ExternalBond'):
b = int(bond.attrib['from']) b = int(bond.attrib['from'])
template.externalBonds.append(b) template.externalBonds.append(b)
template.atoms[b].externalBonds += 1 template.atoms[b].externalBonds += 1
for template in self._templates.values(): self.registerResidueTemplate(template)
signature = _createResidueSignature([atom.element for atom in template.atoms])
if signature in self._templateSignatures:
self._templateSignatures[signature].append(template)
else:
self._templateSignatures[signature] = [template]
# Build sets of every atom type belonging to each class
for type in self._atomTypes:
atomClass = self._atomTypes[type][0]
if atomClass in self._atomClasses:
typeSet = self._atomClasses[atomClass]
else:
typeSet = set()
self._atomClasses[atomClass] = typeSet
typeSet.add(type)
self._atomClasses[''].add(type)
# Load force definitions # Load force definitions
...@@ -171,12 +149,53 @@ class ForceField(object): ...@@ -171,12 +149,53 @@ class ForceField(object):
# Load scripts # Load scripts
for node in tree.getroot().findall('Script'): for node in tree.getroot().findall('Script'):
self._scripts.append(node.text) self.registerScript(node.text)
def getGenerators(self):
"""Get the list of all registered generators."""
return self._forces
def registerGenerator(self, generator):
"""Register a new generator."""
self._forces.append(generator)
def registerAtomType(self, parameters):
"""Register a new atom type."""
name = parameters['name']
if name in self._atomTypes:
raise ValueError('Found multiple definitions for atom type: '+name)
atomClass = parameters['class']
mass = _convertParameterToNumber(parameters['mass'])
element = None
if 'element' in parameters:
element = parameters['element']
if not isinstance(element, elem.Element):
element = elem.get_by_symbol(element)
self._atomTypes[name] = (atomClass, mass, element)
if atomClass in self._atomClasses:
typeSet = self._atomClasses[atomClass]
else:
typeSet = set()
self._atomClasses[atomClass] = typeSet
typeSet.add(name)
self._atomClasses[''].add(name)
def registerResidueTemplate(self, template):
"""Register a new residue template."""
self._templates[template.name] = template
signature = _createResidueSignature([atom.element for atom in template.atoms])
if signature in self._templateSignatures:
self._templateSignatures[signature].append(template)
else:
self._templateSignatures[signature] = [template]
def registerScript(self, script):
"""Register a new script to be executed after building the System."""
self._scripts.append(script)
def _findAtomTypes(self, node, num): def _findAtomTypes(self, attrib, num):
"""Parse the attributes on an XML tag to find the set of atom types for each atom it involves.""" """Parse the attributes on an XML tag to find the set of atom types for each atom it involves."""
types = [] types = []
attrib = node.attrib
for i in range(num): for i in range(num):
if num == 1: if num == 1:
suffix = '' suffix = ''
...@@ -186,7 +205,7 @@ class ForceField(object): ...@@ -186,7 +205,7 @@ class ForceField(object):
typeAttrib = 'type'+suffix typeAttrib = 'type'+suffix
if classAttrib in attrib: if classAttrib in attrib:
if typeAttrib in attrib: if typeAttrib in attrib:
raise ValueError('Tag specifies both a type and a class for the same atom: '+etree.tostring(node)) raise ValueError('Specified both a type and a class for the same atom: '+str(attrib))
if attrib[classAttrib] not in self._atomClasses: if attrib[classAttrib] not in self._atomClasses:
types.append(None) # Unknown atom class types.append(None) # Unknown atom class
else: else:
...@@ -198,18 +217,17 @@ class ForceField(object): ...@@ -198,18 +217,17 @@ class ForceField(object):
types.append([attrib[typeAttrib]]) types.append([attrib[typeAttrib]])
return types return types
def _parseTorsion(self, node): def _parseTorsion(self, attrib):
"""Parse the node defining a torsion.""" """Parse the node defining a torsion."""
types = self._findAtomTypes(node, 4) types = self._findAtomTypes(attrib, 4)
if None in types: if None in types:
return None return None
torsion = PeriodicTorsion(types) torsion = PeriodicTorsion(types)
attrib = node.attrib
index = 1 index = 1
while 'phase%d'%index in attrib: while 'phase%d'%index in attrib:
torsion.periodicity.append(int(attrib['periodicity%d'%index])) torsion.periodicity.append(int(attrib['periodicity%d'%index]))
torsion.phase.append(float(attrib['phase%d'%index])) torsion.phase.append(_convertParameterToNumber(attrib['phase%d'%index]))
torsion.k.append(float(attrib['k%d'%index])) torsion.k.append(_convertParameterToNumber(attrib['k%d'%index]))
index += 1 index += 1
return torsion return torsion
...@@ -234,6 +252,11 @@ class ForceField(object): ...@@ -234,6 +252,11 @@ class ForceField(object):
self.virtualSites = [] self.virtualSites = []
self.bonds = [] self.bonds = []
self.externalBonds = [] self.externalBonds = []
def addBond(self, atom1, atom2):
self.bonds.append((atom1, atom2))
self.atoms[atom1].bondedTo.append(atom2)
self.atoms[atom2].bondedTo.append(atom1)
class _TemplateAtomData: class _TemplateAtomData:
"""Inner class used to encapsulate data about an atom in a residue template definition.""" """Inner class used to encapsulate data about an atom in a residue template definition."""
...@@ -660,23 +683,27 @@ def _findMatchErrors(forcefield, res): ...@@ -660,23 +683,27 @@ def _findMatchErrors(forcefield, res):
class HarmonicBondGenerator: class HarmonicBondGenerator:
"""A HarmonicBondGenerator constructs a HarmonicBondForce.""" """A HarmonicBondGenerator constructs a HarmonicBondForce."""
def __init__(self): def __init__(self, forcefield):
self.ff = forcefield
self.types1 = [] self.types1 = []
self.types2 = [] self.types2 = []
self.length = [] self.length = []
self.k = [] self.k = []
def registerBond(self, parameters):
types = self.ff._findAtomTypes(parameters, 2)
if None not in types:
self.types1.append(types[0])
self.types2.append(types[1])
self.length.append(_convertParameterToNumber(parameters['length']))
self.k.append(_convertParameterToNumber(parameters['k']))
@staticmethod @staticmethod
def parseElement(element, ff): def parseElement(element, ff):
generator = HarmonicBondGenerator() generator = HarmonicBondGenerator(ff)
ff._forces.append(generator) ff.registerGenerator(generator)
for bond in element.findall('Bond'): for bond in element.findall('Bond'):
types = ff._findAtomTypes(bond, 2) generator.registerBond(bond.attrib)
if None not in types:
generator.types1.append(types[0])
generator.types2.append(types[1])
generator.length.append(float(bond.attrib['length']))
generator.k.append(float(bond.attrib['k']))
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args): def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
existing = [sys.getForce(i) for i in range(sys.getNumForces())] existing = [sys.getForce(i) for i in range(sys.getNumForces())]
...@@ -707,25 +734,29 @@ parsers["HarmonicBondForce"] = HarmonicBondGenerator.parseElement ...@@ -707,25 +734,29 @@ parsers["HarmonicBondForce"] = HarmonicBondGenerator.parseElement
class HarmonicAngleGenerator: class HarmonicAngleGenerator:
"""A HarmonicAngleGenerator constructs a HarmonicAngleForce.""" """A HarmonicAngleGenerator constructs a HarmonicAngleForce."""
def __init__(self): def __init__(self, forcefield):
self.ff = forcefield
self.types1 = [] self.types1 = []
self.types2 = [] self.types2 = []
self.types3 = [] self.types3 = []
self.angle = [] self.angle = []
self.k = [] self.k = []
def registerAngle(self, parameters):
types = self.ff._findAtomTypes(parameters, 3)
if None not in types:
self.types1.append(types[0])
self.types2.append(types[1])
self.types3.append(types[2])
self.angle.append(_convertParameterToNumber(parameters['angle']))
self.k.append(_convertParameterToNumber(parameters['k']))
@staticmethod @staticmethod
def parseElement(element, ff): def parseElement(element, ff):
generator = HarmonicAngleGenerator() generator = HarmonicAngleGenerator(ff)
ff._forces.append(generator) ff.registerGenerator(generator)
for angle in element.findall('Angle'): for angle in element.findall('Angle'):
types = ff._findAtomTypes(angle, 3) generator.registerAngle(angle.attrib)
if None not in types:
generator.types1.append(types[0])
generator.types2.append(types[1])
generator.types3.append(types[2])
generator.angle.append(float(angle.attrib['angle']))
generator.k.append(float(angle.attrib['k']))
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args): def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
existing = [sys.getForce(i) for i in range(sys.getNumForces())] existing = [sys.getForce(i) for i in range(sys.getNumForces())]
...@@ -789,23 +820,29 @@ class PeriodicTorsion: ...@@ -789,23 +820,29 @@ class PeriodicTorsion:
class PeriodicTorsionGenerator: class PeriodicTorsionGenerator:
"""A PeriodicTorsionGenerator constructs a PeriodicTorsionForce.""" """A PeriodicTorsionGenerator constructs a PeriodicTorsionForce."""
def __init__(self): def __init__(self, forcefield):
self.ff = forcefield
self.proper = [] self.proper = []
self.improper = [] self.improper = []
def registerProperTorsion(self, parameters):
torsion = self.ff._parseTorsion(parameters)
if torsion is not None:
self.proper.append(torsion)
def registerImproperTorsion(self, parameters):
torsion = self.ff._parseTorsion(parameters)
if torsion is not None:
self.improper.append(torsion)
@staticmethod @staticmethod
def parseElement(element, ff): def parseElement(element, ff):
generator = PeriodicTorsionGenerator() generator = PeriodicTorsionGenerator(ff)
generator.ff = ff ff.registerGenerator(generator)
ff._forces.append(generator)
for torsion in element.findall('Proper'): for torsion in element.findall('Proper'):
torsion = ff._parseTorsion(torsion) generator.registerProperTorsion(torsion.attrib)
if torsion is not None:
generator.proper.append(torsion)
for torsion in element.findall('Improper'): for torsion in element.findall('Improper'):
torsion = ff._parseTorsion(torsion) generator.registerImproperTorsion(torsion.attrib)
if torsion is not None:
generator.improper.append(torsion)
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args): def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
existing = [sys.getForce(i) for i in range(sys.getNumForces())] existing = [sys.getForce(i) for i in range(sys.getNumForces())]
...@@ -888,21 +925,21 @@ class RBTorsion: ...@@ -888,21 +925,21 @@ class RBTorsion:
class RBTorsionGenerator: class RBTorsionGenerator:
"""An RBTorsionGenerator constructs an RBTorsionForce.""" """An RBTorsionGenerator constructs an RBTorsionForce."""
def __init__(self): def __init__(self, forcefield):
self.ff = forcefield
self.proper = [] self.proper = []
self.improper = [] self.improper = []
@staticmethod @staticmethod
def parseElement(element, ff): def parseElement(element, ff):
generator = RBTorsionGenerator() generator = RBTorsionGenerator(ff)
generator.ff = ff ff.registerGenerator(generator)
ff._forces.append(generator)
for torsion in element.findall('Proper'): for torsion in element.findall('Proper'):
types = ff._findAtomTypes(torsion, 4) types = ff._findAtomTypes(torsion.attrib, 4)
if None not in types: if None not in types:
generator.proper.append(RBTorsion(types, [float(torsion.attrib['c'+str(i)]) for i in range(6)])) generator.proper.append(RBTorsion(types, [float(torsion.attrib['c'+str(i)]) for i in range(6)]))
for torsion in element.findall('Improper'): for torsion in element.findall('Improper'):
types = ff._findAtomTypes(torsion, 4) types = ff._findAtomTypes(torsion.attrib, 4)
if None not in types: if None not in types:
generator.improper.append(RBTorsion(types, [float(torsion.attrib['c'+str(i)]) for i in range(6)])) generator.improper.append(RBTorsion(types, [float(torsion.attrib['c'+str(i)]) for i in range(6)]))
...@@ -984,15 +1021,15 @@ class CMAPTorsion: ...@@ -984,15 +1021,15 @@ class CMAPTorsion:
class CMAPTorsionGenerator: class CMAPTorsionGenerator:
"""A CMAPTorsionGenerator constructs a CMAPTorsionForce.""" """A CMAPTorsionGenerator constructs a CMAPTorsionForce."""
def __init__(self): def __init__(self, forcefield):
self.ff = forcefield
self.torsions = [] self.torsions = []
self.maps = [] self.maps = []
@staticmethod @staticmethod
def parseElement(element, ff): def parseElement(element, ff):
generator = CMAPTorsionGenerator() generator = CMAPTorsionGenerator(ff)
generator.ff = ff ff.registerGenerator(generator)
ff._forces.append(generator)
for map in element.findall('Map'): for map in element.findall('Map'):
values = [float(x) for x in map.text.split()] values = [float(x) for x in map.text.split()]
size = sqrt(len(values)) size = sqrt(len(values))
...@@ -1000,7 +1037,7 @@ class CMAPTorsionGenerator: ...@@ -1000,7 +1037,7 @@ class CMAPTorsionGenerator:
raise ValueError('CMAP must have the same number of elements along each dimension') raise ValueError('CMAP must have the same number of elements along each dimension')
generator.maps.append(values) generator.maps.append(values)
for torsion in element.findall('Torsion'): for torsion in element.findall('Torsion'):
types = ff._findAtomTypes(torsion, 5) types = ff._findAtomTypes(torsion.attrib, 5)
if None not in types: if None not in types:
generator.torsions.append(CMAPTorsion(types, int(torsion.attrib['map']))) generator.torsions.append(CMAPTorsion(types, int(torsion.attrib['map'])))
...@@ -1064,28 +1101,32 @@ parsers["CMAPTorsionForce"] = CMAPTorsionGenerator.parseElement ...@@ -1064,28 +1101,32 @@ parsers["CMAPTorsionForce"] = CMAPTorsionGenerator.parseElement
class NonbondedGenerator: class NonbondedGenerator:
"""A NonbondedGenerator constructs a NonbondedForce.""" """A NonbondedGenerator constructs a NonbondedForce."""
def __init__(self, coulomb14scale, lj14scale): def __init__(self, forcefield, coulomb14scale, lj14scale):
self.ff = forcefield
self.coulomb14scale = coulomb14scale self.coulomb14scale = coulomb14scale
self.lj14scale = lj14scale self.lj14scale = lj14scale
self.typeMap = {} self.typeMap = {}
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
@staticmethod @staticmethod
def parseElement(element, ff): def parseElement(element, ff):
existing = [f for f in ff._forces if isinstance(f, NonbondedGenerator)] existing = [f for f in ff._forces if isinstance(f, NonbondedGenerator)]
if len(existing) == 0: if len(existing) == 0:
generator = NonbondedGenerator(float(element.attrib['coulomb14scale']), float(element.attrib['lj14scale'])) generator = NonbondedGenerator(ff, float(element.attrib['coulomb14scale']), float(element.attrib['lj14scale']))
ff._forces.append(generator) ff.registerGenerator(generator)
else: else:
# Multiple <NonbondedForce> tags were found, probably in different files. Simply add more types to the existing one. # Multiple <NonbondedForce> tags were found, probably in different files. Simply add more types to the existing one.
generator = existing[0] generator = existing[0]
if generator.coulomb14scale != float(element.attrib['coulomb14scale']) or generator.lj14scale != float(element.attrib['lj14scale']): 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') raise ValueError('Found multiple NonbondedForce tags with different 1-4 scales')
for atom in element.findall('Atom'): for atom in element.findall('Atom'):
types = ff._findAtomTypes(atom, 1) generator.registerAtom(atom.attrib)
if None not in types:
values = (float(atom.attrib['charge']), float(atom.attrib['sigma']), float(atom.attrib['epsilon']))
for t in types[0]:
generator.typeMap[t] = values
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args): def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
methodMap = {NoCutoff:mm.NonbondedForce.NoCutoff, methodMap = {NoCutoff:mm.NonbondedForce.NoCutoff,
...@@ -1149,24 +1190,28 @@ parsers["NonbondedForce"] = NonbondedGenerator.parseElement ...@@ -1149,24 +1190,28 @@ parsers["NonbondedForce"] = NonbondedGenerator.parseElement
class GBSAOBCGenerator: class GBSAOBCGenerator:
"""A GBSAOBCGenerator constructs a GBSAOBCForce.""" """A GBSAOBCGenerator constructs a GBSAOBCForce."""
def __init__(self): def __init__(self, forcefield):
self.ff = forcefield
self.typeMap = {} self.typeMap = {}
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
@staticmethod @staticmethod
def parseElement(element, ff): def parseElement(element, ff):
existing = [f for f in ff._forces if isinstance(f, GBSAOBCGenerator)] existing = [f for f in ff._forces if isinstance(f, GBSAOBCGenerator)]
if len(existing) == 0: if len(existing) == 0:
generator = GBSAOBCGenerator() generator = GBSAOBCGenerator(ff)
ff._forces.append(generator) ff.registerGenerator(generator)
else: else:
# Multiple <GBSAOBCForce> tags were found, probably in different files. Simply add more types to the existing one. # Multiple <GBSAOBCForce> tags were found, probably in different files. Simply add more types to the existing one.
generator = existing[0] generator = existing[0]
for atom in element.findall('Atom'): for atom in element.findall('Atom'):
types = ff._findAtomTypes(atom, 1) generator.registerAtom(atom.attrib)
if None not in types:
values = (float(atom.attrib['charge']), float(atom.attrib['radius']), float(atom.attrib['scale']))
for t in types[0]:
generator.typeMap[t] = values
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args): def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
methodMap = {NoCutoff:mm.NonbondedForce.NoCutoff, methodMap = {NoCutoff:mm.NonbondedForce.NoCutoff,
...@@ -1206,7 +1251,6 @@ class GBVIGenerator: ...@@ -1206,7 +1251,6 @@ class GBVIGenerator:
"""A GBVIGenerator constructs a GBVIForce.""" """A GBVIGenerator constructs a GBVIForce."""
def __init__(self,ff): def __init__(self,ff):
self.ff = ff self.ff = ff
self.fixedParameters = {} self.fixedParameters = {}
self.fixedParameters['soluteDielectric'] = 1.0 self.fixedParameters['soluteDielectric'] = 1.0
...@@ -1224,9 +1268,9 @@ class GBVIGenerator: ...@@ -1224,9 +1268,9 @@ class GBVIGenerator:
if (key in element.attrib): if (key in element.attrib):
generator.fixedParameters[key] = float(element.attrib[key]) generator.fixedParameters[key] = float(element.attrib[key])
ff._forces.append(generator) ff.registerGenerator(generator)
for atom in element.findall('Atom'): for atom in element.findall('Atom'):
types = ff._findAtomTypes(atom, 1) types = ff._findAtomTypes(atom.attrib, 1)
if None not in types: if None not in types:
values = (float(atom.attrib['charge']), float(atom.attrib['radius']), float(atom.attrib['gamma'])) values = (float(atom.attrib['charge']), float(atom.attrib['radius']), float(atom.attrib['gamma']))
for t in types[0]: for t in types[0]:
...@@ -1292,7 +1336,8 @@ parsers["GBVIForce"] = GBVIGenerator.parseElement ...@@ -1292,7 +1336,8 @@ parsers["GBVIForce"] = GBVIGenerator.parseElement
class CustomBondGenerator: class CustomBondGenerator:
"""A CustomBondGenerator constructs a CustomBondForce.""" """A CustomBondGenerator constructs a CustomBondForce."""
def __init__(self): def __init__(self, forcefield):
self.ff = forcefield
self.types1 = [] self.types1 = []
self.types2 = [] self.types2 = []
self.globalParams = {} self.globalParams = {}
...@@ -1301,15 +1346,15 @@ class CustomBondGenerator: ...@@ -1301,15 +1346,15 @@ class CustomBondGenerator:
@staticmethod @staticmethod
def parseElement(element, ff): def parseElement(element, ff):
generator = CustomBondGenerator() generator = CustomBondGenerator(ff)
ff._forces.append(generator) ff.registerGenerator(generator)
generator.energy = element.attrib['energy'] generator.energy = element.attrib['energy']
for param in element.findall('GlobalParameter'): for param in element.findall('GlobalParameter'):
generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue']) generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
for param in element.findall('PerBondParameter'): for param in element.findall('PerBondParameter'):
generator.perBondParams.append(param.attrib['name']) generator.perBondParams.append(param.attrib['name'])
for bond in element.findall('Bond'): for bond in element.findall('Bond'):
types = ff._findAtomTypes(bond, 2) types = ff._findAtomTypes(bond.attrib, 2)
if None not in types: if None not in types:
generator.types1.append(types[0]) generator.types1.append(types[0])
generator.types2.append(types[1]) generator.types2.append(types[1])
...@@ -1339,7 +1384,8 @@ parsers["CustomBondForce"] = CustomBondGenerator.parseElement ...@@ -1339,7 +1384,8 @@ parsers["CustomBondForce"] = CustomBondGenerator.parseElement
class CustomAngleGenerator: class CustomAngleGenerator:
"""A CustomAngleGenerator constructs a CustomAngleForce.""" """A CustomAngleGenerator constructs a CustomAngleForce."""
def __init__(self): def __init__(self, forcefield):
self.ff = forcefield
self.types1 = [] self.types1 = []
self.types2 = [] self.types2 = []
self.types3 = [] self.types3 = []
...@@ -1349,15 +1395,15 @@ class CustomAngleGenerator: ...@@ -1349,15 +1395,15 @@ class CustomAngleGenerator:
@staticmethod @staticmethod
def parseElement(element, ff): def parseElement(element, ff):
generator = CustomAngleGenerator() generator = CustomAngleGenerator(ff)
ff._forces.append(generator) ff.registerGenerator(generator)
generator.energy = element.attrib['energy'] generator.energy = element.attrib['energy']
for param in element.findall('GlobalParameter'): for param in element.findall('GlobalParameter'):
generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue']) generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
for param in element.findall('PerAngleParameter'): for param in element.findall('PerAngleParameter'):
generator.perAngleParams.append(param.attrib['name']) generator.perAngleParams.append(param.attrib['name'])
for angle in element.findall('Angle'): for angle in element.findall('Angle'):
types = ff._findAtomTypes(angle, 3) types = ff._findAtomTypes(angle.attrib, 3)
if None not in types: if None not in types:
generator.types1.append(types[0]) generator.types1.append(types[0])
generator.types2.append(types[1]) generator.types2.append(types[1])
...@@ -1401,7 +1447,8 @@ class CustomTorsion: ...@@ -1401,7 +1447,8 @@ class CustomTorsion:
class CustomTorsionGenerator: class CustomTorsionGenerator:
"""A CustomTorsionGenerator constructs a CustomTorsionForce.""" """A CustomTorsionGenerator constructs a CustomTorsionForce."""
def __init__(self): def __init__(self, forcefield):
self.ff = forcefield
self.proper = [] self.proper = []
self.improper = [] self.improper = []
self.globalParams = {} self.globalParams = {}
...@@ -1409,20 +1456,19 @@ class CustomTorsionGenerator: ...@@ -1409,20 +1456,19 @@ class CustomTorsionGenerator:
@staticmethod @staticmethod
def parseElement(element, ff): def parseElement(element, ff):
generator = CustomTorsionGenerator() generator = CustomTorsionGenerator(ff)
generator.ff = ff ff.registerGenerator(generator)
ff._forces.append(generator)
generator.energy = element.attrib['energy'] generator.energy = element.attrib['energy']
for param in element.findall('GlobalParameter'): for param in element.findall('GlobalParameter'):
generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue']) generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
for param in element.findall('PerTorsionParameter'): for param in element.findall('PerTorsionParameter'):
generator.perTorsionParams.append(param.attrib['name']) generator.perTorsionParams.append(param.attrib['name'])
for torsion in element.findall('Proper'): for torsion in element.findall('Proper'):
types = ff._findAtomTypes(torsion, 4) types = ff._findAtomTypes(torsion.attrib, 4)
if None not in types: if None not in types:
generator.proper.append(CustomTorsion(types, [float(torsion.attrib[param]) for param in generator.perTorsionParams])) generator.proper.append(CustomTorsion(types, [float(torsion.attrib[param]) for param in generator.perTorsionParams]))
for torsion in element.findall('Improper'): for torsion in element.findall('Improper'):
types = ff._findAtomTypes(torsion, 4) types = ff._findAtomTypes(torsion.attrib, 4)
if None not in types: if None not in types:
generator.improper.append(CustomTorsion(types, [float(torsion.attrib[param]) for param in generator.perTorsionParams])) generator.improper.append(CustomTorsion(types, [float(torsion.attrib[param]) for param in generator.perTorsionParams]))
...@@ -1491,7 +1537,8 @@ parsers["CustomTorsionForce"] = CustomTorsionGenerator.parseElement ...@@ -1491,7 +1537,8 @@ parsers["CustomTorsionForce"] = CustomTorsionGenerator.parseElement
class CustomNonbondedGenerator: class CustomNonbondedGenerator:
"""A CustomNonbondedGenerator constructs a CustomNonbondedForce.""" """A CustomNonbondedGenerator constructs a CustomNonbondedForce."""
def __init__(self, energy, bondCutoff): def __init__(self, forcefield, energy, bondCutoff):
self.ff = forcefield
self.energy = energy self.energy = energy
self.bondCutoff = bondCutoff self.bondCutoff = bondCutoff
self.typeMap = {} self.typeMap = {}
...@@ -1501,14 +1548,14 @@ class CustomNonbondedGenerator: ...@@ -1501,14 +1548,14 @@ class CustomNonbondedGenerator:
@staticmethod @staticmethod
def parseElement(element, ff): def parseElement(element, ff):
generator = CustomNonbondedGenerator(element.attrib['energy'], int(element.attrib['bondCutoff'])) generator = CustomNonbondedGenerator(ff, element.attrib['energy'], int(element.attrib['bondCutoff']))
ff._forces.append(generator) ff.registerGenerator(generator)
for param in element.findall('GlobalParameter'): for param in element.findall('GlobalParameter'):
generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue']) generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
for param in element.findall('PerParticleParameter'): for param in element.findall('PerParticleParameter'):
generator.perParticleParams.append(param.attrib['name']) generator.perParticleParams.append(param.attrib['name'])
for atom in element.findall('Atom'): for atom in element.findall('Atom'):
types = ff._findAtomTypes(atom, 1) types = ff._findAtomTypes(atom.attrib, 1)
if None not in types: if None not in types:
values = [float(atom.attrib[param]) for param in generator.perParticleParams] values = [float(atom.attrib[param]) for param in generator.perParticleParams]
for t in types[0]: for t in types[0]:
...@@ -1589,7 +1636,8 @@ parsers["CustomNonbondedForce"] = CustomNonbondedGenerator.parseElement ...@@ -1589,7 +1636,8 @@ parsers["CustomNonbondedForce"] = CustomNonbondedGenerator.parseElement
class CustomGBGenerator: class CustomGBGenerator:
"""A CustomGBGenerator constructs a CustomGBForce.""" """A CustomGBGenerator constructs a CustomGBForce."""
def __init__(self): def __init__(self, forcefield):
self.ff = forcefield
self.typeMap = {} self.typeMap = {}
self.globalParams = {} self.globalParams = {}
self.perParticleParams = [] self.perParticleParams = []
...@@ -1599,14 +1647,14 @@ class CustomGBGenerator: ...@@ -1599,14 +1647,14 @@ class CustomGBGenerator:
@staticmethod @staticmethod
def parseElement(element, ff): def parseElement(element, ff):
generator = CustomGBGenerator() generator = CustomGBGenerator(ff)
ff._forces.append(generator) ff.registerGenerator(generator)
for param in element.findall('GlobalParameter'): for param in element.findall('GlobalParameter'):
generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue']) generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
for param in element.findall('PerParticleParameter'): for param in element.findall('PerParticleParameter'):
generator.perParticleParams.append(param.attrib['name']) generator.perParticleParams.append(param.attrib['name'])
for atom in element.findall('Atom'): for atom in element.findall('Atom'):
types = ff._findAtomTypes(atom, 1) types = ff._findAtomTypes(atom.attrib, 1)
if None not in types: if None not in types:
values = [float(atom.attrib[param]) for param in generator.perParticleParams] values = [float(atom.attrib[param]) for param in generator.perParticleParams]
for t in types[0]: for t in types[0]:
...@@ -1729,7 +1777,7 @@ class AmoebaBondGenerator: ...@@ -1729,7 +1777,7 @@ class AmoebaBondGenerator:
generator = AmoebaBondGenerator(float(element.attrib['bond-cubic']), float(element.attrib['bond-quartic'])) generator = AmoebaBondGenerator(float(element.attrib['bond-cubic']), float(element.attrib['bond-quartic']))
forceField._forces.append(generator) forceField._forces.append(generator)
for bond in element.findall('Bond'): for bond in element.findall('Bond'):
types = forceField._findAtomTypes(bond, 2) types = forceField._findAtomTypes(bond.attrib, 2)
if None not in types: if None not in types:
generator.types1.append(types[0]) generator.types1.append(types[0])
generator.types2.append(types[1]) generator.types2.append(types[1])
...@@ -1842,7 +1890,7 @@ class AmoebaAngleGenerator: ...@@ -1842,7 +1890,7 @@ class AmoebaAngleGenerator:
generator = AmoebaAngleGenerator(forceField, float(element.attrib['angle-cubic']), float(element.attrib['angle-quartic']), float(element.attrib['angle-pentic']), float(element.attrib['angle-sextic'])) generator = AmoebaAngleGenerator(forceField, float(element.attrib['angle-cubic']), float(element.attrib['angle-quartic']), float(element.attrib['angle-pentic']), float(element.attrib['angle-sextic']))
forceField._forces.append(generator) forceField._forces.append(generator)
for angle in element.findall('Angle'): for angle in element.findall('Angle'):
types = forceField._findAtomTypes(angle, 3) types = forceField._findAtomTypes(angle.attrib, 3)
if None not in types: if None not in types:
generator.types1.append(types[0]) generator.types1.append(types[0])
...@@ -2327,7 +2375,7 @@ class AmoebaTorsionGenerator: ...@@ -2327,7 +2375,7 @@ class AmoebaTorsionGenerator:
# where ti=[amplitude_i,angle_i] # where ti=[amplitude_i,angle_i]
for torsion in element.findall('Torsion'): for torsion in element.findall('Torsion'):
types = forceField._findAtomTypes(torsion, 4) types = forceField._findAtomTypes(torsion.attrib, 4)
if None not in types: if None not in types:
generator.types1.append(types[0]) generator.types1.append(types[0])
...@@ -2438,7 +2486,7 @@ class AmoebaPiTorsionGenerator: ...@@ -2438,7 +2486,7 @@ class AmoebaPiTorsionGenerator:
forceField._forces.append(generator) forceField._forces.append(generator)
for piTorsion in element.findall('PiTorsion'): for piTorsion in element.findall('PiTorsion'):
types = forceField._findAtomTypes(piTorsion, 2) types = forceField._findAtomTypes(piTorsion.attrib, 2)
if None not in types: if None not in types:
generator.types1.append(types[0]) generator.types1.append(types[0])
generator.types2.append(types[1]) generator.types2.append(types[1])
...@@ -2557,7 +2605,7 @@ class AmoebaTorsionTorsionGenerator: ...@@ -2557,7 +2605,7 @@ class AmoebaTorsionTorsionGenerator:
# <TorsionTorsion class1="3" class2="1" class3="2" class4="3" class5="1" grid="0" nx="25" ny="25" /> # <TorsionTorsion class1="3" class2="1" class3="2" class4="3" class5="1" grid="0" nx="25" ny="25" />
for torsionTorsion in element.findall('TorsionTorsion'): for torsionTorsion in element.findall('TorsionTorsion'):
types = forceField._findAtomTypes(torsionTorsion, 5) types = forceField._findAtomTypes(torsionTorsion.attrib, 5)
if None not in types: if None not in types:
generator.types1.append(types[0]) generator.types1.append(types[0])
...@@ -2795,7 +2843,7 @@ class AmoebaStretchBendGenerator: ...@@ -2795,7 +2843,7 @@ class AmoebaStretchBendGenerator:
# <StretchBend class1="2" class2="1" class3="4" k1="3.14005676385" k2="3.14005676385" /> # <StretchBend class1="2" class2="1" class3="4" k1="3.14005676385" k2="3.14005676385" />
for stretchBend in element.findall('StretchBend'): for stretchBend in element.findall('StretchBend'):
types = forceField._findAtomTypes(stretchBend, 3) types = forceField._findAtomTypes(stretchBend.attrib, 3)
if None not in types: if None not in types:
generator.types1.append(types[0]) generator.types1.append(types[0])
...@@ -2952,7 +3000,7 @@ class AmoebaVdwGenerator: ...@@ -2952,7 +3000,7 @@ class AmoebaVdwGenerator:
# sigma is modified based on radiustype and radiussize # sigma is modified based on radiustype and radiussize
for atom in element.findall('Vdw'): for atom in element.findall('Vdw'):
types = forceField._findAtomTypes(atom, 1) types = forceField._findAtomTypes(atom.attrib, 1)
if None not in types: if None not in types:
values = [float(atom.attrib['sigma']), float(atom.attrib['epsilon']), float(atom.attrib['reduction'])] values = [float(atom.attrib['sigma']), float(atom.attrib['epsilon']), float(atom.attrib['reduction'])]
...@@ -3225,7 +3273,7 @@ class AmoebaMultipoleGenerator: ...@@ -3225,7 +3273,7 @@ class AmoebaMultipoleGenerator:
# set type map: [ kIndices, multipoles, AMOEBA/OpenMM axis type] # set type map: [ kIndices, multipoles, AMOEBA/OpenMM axis type]
for atom in element.findall('Multipole'): for atom in element.findall('Multipole'):
types = forceField._findAtomTypes(atom, 1) types = forceField._findAtomTypes(atom.attrib, 1)
if None not in types: if None not in types:
# k-indices not provided default to 0 # k-indices not provided default to 0
...@@ -3282,7 +3330,7 @@ class AmoebaMultipoleGenerator: ...@@ -3282,7 +3330,7 @@ class AmoebaMultipoleGenerator:
# polarization parameters # polarization parameters
for atom in element.findall('Polarize'): for atom in element.findall('Polarize'):
types = forceField._findAtomTypes(atom, 1) types = forceField._findAtomTypes(atom.attrib, 1)
if None not in types: if None not in types:
classIndex = atom.attrib['type'] classIndex = atom.attrib['type']
...@@ -3810,7 +3858,7 @@ class AmoebaWcaDispersionGenerator: ...@@ -3810,7 +3858,7 @@ class AmoebaWcaDispersionGenerator:
# typeMap[] = [ radius, epsilon ] # typeMap[] = [ radius, epsilon ]
for atom in element.findall('WcaDispersion'): for atom in element.findall('WcaDispersion'):
types = forceField._findAtomTypes(atom, 1) types = forceField._findAtomTypes(atom.attrib, 1)
if None not in types: if None not in types:
values = [float(atom.attrib['radius']), float(atom.attrib['epsilon'])] values = [float(atom.attrib['radius']), float(atom.attrib['epsilon'])]
...@@ -4154,7 +4202,7 @@ class AmoebaUreyBradleyGenerator: ...@@ -4154,7 +4202,7 @@ class AmoebaUreyBradleyGenerator:
generator = AmoebaUreyBradleyGenerator() generator = AmoebaUreyBradleyGenerator()
forceField._forces.append(generator) forceField._forces.append(generator)
for bond in element.findall('UreyBradley'): for bond in element.findall('UreyBradley'):
types = forceField._findAtomTypes(bond, 3) types = forceField._findAtomTypes(bond.attrib, 3)
if None not in types: if None not in types:
generator.types1.append(types[0]) generator.types1.append(types[0])
...@@ -4205,20 +4253,21 @@ parsers["AmoebaUreyBradleyForce"] = AmoebaUreyBradleyGenerator.parseElement ...@@ -4205,20 +4253,21 @@ parsers["AmoebaUreyBradleyForce"] = AmoebaUreyBradleyGenerator.parseElement
class DrudeGenerator: class DrudeGenerator:
"""A DrudeGenerator constructs a DrudeForce.""" """A DrudeGenerator constructs a DrudeForce."""
def __init__(self): def __init__(self, forcefield):
self.ff = forcefield
self.typeMap = {} self.typeMap = {}
@staticmethod @staticmethod
def parseElement(element, ff): def parseElement(element, ff):
existing = [f for f in ff._forces if isinstance(f, DrudeGenerator)] existing = [f for f in ff._forces if isinstance(f, DrudeGenerator)]
if len(existing) == 0: if len(existing) == 0:
generator = DrudeGenerator() generator = DrudeGenerator(ff)
ff._forces.append(generator) ff.registerGenerator(generator)
else: else:
# Multiple <DrudeForce> tags were found, probably in different files. Simply add more types to the existing one. # Multiple <DrudeForce> tags were found, probably in different files. Simply add more types to the existing one.
generator = existing[0] generator = existing[0]
for particle in element.findall('Particle'): for particle in element.findall('Particle'):
types = ff._findAtomTypes(particle, 5) types = ff._findAtomTypes(particle.attrib, 5)
if None not in types[:2]: if None not in types[:2]:
aniso12 = 0.0 aniso12 = 0.0
aniso34 = 0.0 aniso34 = 0.0
......
wrappers/python/tests/TestForceField.py
View file @ b2212517
...@@ -4,6 +4,7 @@ from simtk.openmm.app import * ...@@ -4,6 +4,7 @@ from simtk.openmm.app import *
from simtk.openmm import * from simtk.openmm import *
from simtk.unit import * from simtk.unit import *
import simtk.openmm.app.element as elem import simtk.openmm.app.element as elem
import simtk.openmm.app.forcefield as forcefield
class TestForceField(unittest.TestCase): class TestForceField(unittest.TestCase):
"""Test the ForceField.createSystem() method.""" """Test the ForceField.createSystem() method."""
...@@ -149,6 +150,41 @@ class TestForceField(unittest.TestCase): ...@@ -149,6 +150,41 @@ class TestForceField(unittest.TestCase):
if diff > 0.1 and diff/norm(f1) > 1e-3: if diff > 0.1 and diff/norm(f1) > 1e-3:
numDifferences += 1 numDifferences += 1
self.assertTrue(numDifferences < system.getNumParticles()/20) # Tolerate occasional differences from numerical error self.assertTrue(numDifferences < system.getNumParticles()/20) # Tolerate occasional differences from numerical error
def test_ProgrammaticForceField(self):
"""Test building a ForceField programmatically."""
# Build the ForceField for TIP3P programmatically.
ff = ForceField()
ff.registerAtomType({'name':'tip3p-O', 'class':'OW', 'mass':15.99943*daltons, 'element':elem.oxygen})
ff.registerAtomType({'name':'tip3p-H', 'class':'HW', 'mass':1.007947*daltons, 'element':elem.hydrogen})
residue = ForceField._TemplateData('HOH')
residue.atoms.append(ForceField._TemplateAtomData('O', 'tip3p-O', elem.oxygen))
residue.atoms.append(ForceField._TemplateAtomData('H1', 'tip3p-H', elem.hydrogen))
residue.atoms.append(ForceField._TemplateAtomData('H2', 'tip3p-H', elem.hydrogen))
residue.addBond(0, 1)
residue.addBond(0, 2)
ff.registerResidueTemplate(residue)
bonds = forcefield.HarmonicBondGenerator(ff)
bonds.registerBond({'class1':'OW', 'class2':'HW', 'length':0.09572*nanometers, 'k':462750.4*kilojoules_per_mole/nanometer})
ff.registerGenerator(bonds)
angles = forcefield.HarmonicAngleGenerator(ff)
angles.registerAngle({'class1':'HW', 'class2':'OW', 'class3':'HW', 'angle':1.82421813418*radians, 'k':836.8*kilojoules_per_mole/radian})
ff.registerGenerator(angles)
nonbonded = forcefield.NonbondedGenerator(ff, 0.833333, 0.5)
nonbonded.registerAtom({'type':'tip3p-O', 'charge':-0.834, 'sigma':0.31507524065751241*nanometers, 'epsilon':0.635968*kilojoules_per_mole})
nonbonded.registerAtom({'type':'tip3p-H', 'charge':0.417, 'sigma':1*nanometers, 'epsilon':0*kilojoules_per_mole})
ff.registerGenerator(nonbonded)
# Build a water box.
modeller = Modeller(Topology(), [])
modeller.addSolvent(ff, boxSize=Vec3(3, 3, 3)*nanometers)
# Create a system using the programmatic force field as well as one from an XML file.
system1 = ff.createSystem(modeller.topology)
ff2 = ForceField('tip3p.xml')
system2 = ff2.createSystem(modeller.topology)
self.assertEqual(XmlSerializer.serialize(system1), XmlSerializer.serialize(system2))
class AmoebaTestForceField(unittest.TestCase): class AmoebaTestForceField(unittest.TestCase):
"""Test the ForceField.createSystem() method with the AMOEBA forcefield.""" """Test the ForceField.createSystem() method with the AMOEBA forcefield."""
......
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