Merge pull request #1593 from peastman/hbond

ForceField supports CustomHbondForce

docs-source/usersguide/application.rst

@@ -2696,6 +2696,62 @@ must include an attribute called :code:`radius`\ .
 CustomGBForce also allows you to define tabulated functions.  See section
 :ref:`tabulated-functions` for details.
 
+<CustomHbondForce>
+=========================
+
+To add a CustomHbondForce to the System, include a tag that looks like this:
+
+.. code-block:: xml
+
+    <CustomHbondForce particlesPerDonor="3" particlesPerAcceptor="2" bondCutoff="2"
+        energy="scale*k*(distance(a1,d1)-r0)^2*(angle(a1,d1,d2)-theta0)^2">
+     <GlobalParameter name="scale" defaultValue="1"/>
+     <PerDonorParameter name="theta0"/>
+     <PerAcceptorParameter name="k"/>
+     <PerAcceptorParameter name="r0"/>
+     <Donor class1="H" class2="N" class3="C" theta0="2.1"/>
+     <Acceptor class1="O" class2="C" k="115.0" r0="0.2"/>
+     ...
+    </CustomHbondForce>
+
+The energy expression for the CustomHbondForce is specified by the
+:code:`energy` attribute.  This is a mathematical expression that gives the
+energy of each donor-acceptor interaction as a function of various particle coordinates,
+distances, and angles.  See the API documentation for details.  :code:`particlesPerDonor`
+specifies the number of particles that make up a donor group, and :code:`particlesPerAcceptor`
+specifies the number of particles that make up an acceptor group.
+
+The expression may depend on an arbitrary list of global, per-donor, or
+per-acceptor parameters.  Use a :code:`<GlobalParameter>` tag to define a global
+parameter, a :code:`<PerDonorParameter>` tag to define a per-donor parameter,
+and a :code:`<PerAcceptorParameter>` tag to define a per-acceptor parameter.
+
+Exclusions are created automatically based on the :code:`bondCutoff` attribute.
+If any atom of a donor is within the specified distance (measured in bonds) of
+any atom of an acceptor, an exclusion is added to prevent them from interacting
+with each other.  If a donor and an acceptor share any atom in common, that is a
+bond distance of 0, so they are always excluded.
+
+Every :code:`<Donor>` or :code:`<Acceptor>` tag defines a rule for creating donor
+or acceptor groups.  The number of atoms specified in each one must match the
+value of :code:`particlesPerDonor` or :code:`particlesPerAcceptor` specified in the
+parent tag. Each tag may identify the atoms either by type (using the attributes
+:code:`type1`\ , :code:`type2`\ , ...) or by class (using the attributes
+:code:`class1`\ , :code:`class2`\ , ...).  The force field considers every atom
+in the system (if the number of atoms is 1), every pair of bonded atoms (if the number
+of atoms is 2), or every set of three atoms where the first is bonded to the second
+and the second to the third (if the number of atoms is 3).  For each one, it searches
+for a rule whose atom types or atom classes match the atoms.  If it finds one,
+it calls :code:`addDonor()` or :code:`addAcceptor()` on the CustomHbondForce.
+Otherwise, it ignores that set and continues. The remaining attributes are the
+values to use for the per-donor and per-acceptor parameters. All parameters must
+be specified for every tag, and the attribute name must match the name of the
+parameter.  For instance, if there is a per-donor parameter with the name “k”,
+then every :code:`<Donor>` tag must include an attribute called :code:`k`\ .
+
+CustomHbondForce also allows you to define tabulated functions.  See section
+:ref:`tabulated-functions` for details.
+
 <CustomManyParticleForce>
 =========================
 
@@ -2716,7 +2772,7 @@ To add a CustomManyParticleForce to the System, include a tag that looks like th
 
 The energy expression for the CustomManyParticleForce is specified by the
 :code:`energy` attribute.  This is a mathematical expression that gives the
-energy of each pairwise interaction as a function of various particle coordinates,
+energy of each interaction as a function of various particle coordinates,
 distances, and angles.  See the API documentation for details.  :code:`particlesPerSet`
 specifies the number of particles involved in the interaction and
 :code:`permutationMode` specifies the permutation mode.

wrappers/python/simtk/openmm/app/forcefield.py

@@ -53,6 +53,40 @@ def _convertParameterToNumber(param):
         return param.value_in_unit_system(unit.md_unit_system)
     return float(param)
 
+def _parseFunctions(element):
+    """Parse the attributes on an XML tag to find any tabulated functions it defines."""
+    functions = []
+    for function in element.findall('Function'):
+        values = [float(x) for x in function.text.split()]
+        if 'type' in function.attrib:
+            functionType = function.attrib['type']
+        else:
+            functionType = 'Continuous1D'
+        params = {}
+        for key in function.attrib:
+            if key.endswith('size'):
+                params[key] = int(function.attrib[key])
+            elif key.endswith('min') or key.endswith('max'):
+                params[key] = float(function.attrib[key])
+        functions.append((function.attrib['name'], functionType, values, params))
+    return functions
+
+def _createFunctions(force, functions):
+    """Add TabulatedFunctions to a Force based on the information that was recorded by _parseFunctions()."""
+    for (name, type, values, params) in functions:
+        if type == 'Continuous1D':
+            force.addTabulatedFunction(name, mm.Continuous1DFunction(values, params['min'], params['max']))
+        elif type == 'Continuous2D':
+            force.addTabulatedFunction(name, mm.Continuous2DFunction(params['xsize'], params['ysize'], values, params['xmin'], params['xmax'], params['ymin'], params['ymax']))
+        elif type == 'Continuous3D':
+            force.addTabulatedFunction(name, mm.Continuous2DFunction(params['xsize'], params['ysize'], params['zsize'], values, params['xmin'], params['xmax'], params['ymin'], params['ymax'], params['zmin'], params['zmax']))
+        elif type == 'Discrete1D':
+            force.addTabulatedFunction(name, mm.Discrete1DFunction(values))
+        elif type == 'Discrete2D':
+            force.addTabulatedFunction(name, mm.Discrete2DFunction(params['xsize'], params['ysize'], values))
+        elif type == 'Discrete3D':
+            force.addTabulatedFunction(name, mm.Discrete2DFunction(params['xsize'], params['ysize'], params['zsize'], values))
+
 # Enumerated values for nonbonded method
 
 class NoCutoff(object):
@@ -462,7 +496,7 @@ class ForceField(object):
             torsion.k.append(_convertParameterToNumber(attrib['k%d'%index]))
             index += 1
         return torsion
-
+        
     class _SystemData(object):
         """Inner class used to encapsulate data about the system being created."""
         def __init__(self):
@@ -2600,6 +2634,7 @@ class CustomNonbondedGenerator(object):
             generator.perParticleParams.append(param.attrib['name'])
         generator.params = ForceField._AtomTypeParameters(ff, 'CustomNonbondedForce', 'Atom', generator.perParticleParams)
         generator.params.parseDefinitions(element)
+        generator.functions += _parseFunctions(element)
 
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
         methodMap = {NoCutoff:mm.CustomNonbondedForce.NoCutoff,
@@ -2612,19 +2647,7 @@ class CustomNonbondedGenerator(object):
             force.addGlobalParameter(param, self.globalParams[param])
         for param in self.perParticleParams:
             force.addPerParticleParameter(param)
-        for (name, type, values, params) in self.functions:
-            if type == 'Continuous1D':
-                force.addTabulatedFunction(name, mm.Continuous1DFunction(values, params['min'], params['max']))
-            elif type == 'Continuous2D':
-                force.addTabulatedFunction(name, mm.Continuous2DFunction(params['xsize'], params['ysize'], values, params['xmin'], params['xmax'], params['ymin'], params['ymax']))
-            elif type == 'Continuous3D':
-                force.addTabulatedFunction(name, mm.Continuous2DFunction(params['xsize'], params['ysize'], params['zsize'], values, params['xmin'], params['xmax'], params['ymin'], params['ymax'], params['zmin'], params['zmax']))
-            elif type == 'Discrete1D':
-                force.addTabulatedFunction(name, mm.Discrete1DFunction(values))
-            elif type == 'Discrete2D':
-                force.addTabulatedFunction(name, mm.Discrete2DFunction(params['xsize'], params['ysize'], values))
-            elif type == 'Discrete3D':
-                force.addTabulatedFunction(name, mm.Discrete2DFunction(params['xsize'], params['ysize'], params['zsize'], values))
+        _createFunctions(force, self.functions)
         for atom in data.atoms:
             values = self.params.getAtomParameters(atom, data)
             force.addParticle(values)
@@ -2671,19 +2694,7 @@ class CustomGBGenerator(object):
             generator.computedValues.append((value.attrib['name'], value.text, computationMap[value.attrib['type']]))
         for term in element.findall('EnergyTerm'):
             generator.energyTerms.append((term.text, computationMap[term.attrib['type']]))
-        for function in element.findall("Function"):
-            values = [float(x) for x in function.text.split()]
-            if 'type' in function.attrib:
-                type = function.attrib['type']
-            else:
-                type = 'Continuous1D'
-            params = {}
-            for key in function.attrib:
-                if key.endswith('size'):
-                    params[key] = int(function.attrib[key])
-                elif key.endswith('min') or key.endswith('max'):
-                    params[key] = float(function.attrib[key])
-            generator.functions.append((function.attrib['name'], type, values, params))
+        generator.functions += _parseFunctions(element)
 
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
         methodMap = {NoCutoff:mm.CustomGBForce.NoCutoff,
@@ -2700,19 +2711,7 @@ class CustomGBGenerator(object):
             force.addComputedValue(value[0], value[1], value[2])
         for term in self.energyTerms:
             force.addEnergyTerm(term[0], term[1])
-        for (name, type, values, params) in self.functions:
-            if type == 'Continuous1D':
-                force.addTabulatedFunction(name, mm.Continuous1DFunction(values, params['min'], params['max']))
-            elif type == 'Continuous2D':
-                force.addTabulatedFunction(name, mm.Continuous2DFunction(params['xsize'], params['ysize'], values, params['xmin'], params['xmax'], params['ymin'], params['ymax']))
-            elif type == 'Continuous3D':
-                force.addTabulatedFunction(name, mm.Continuous2DFunction(params['xsize'], params['ysize'], params['zsize'], values, params['xmin'], params['xmax'], params['ymin'], params['ymax'], params['zmin'], params['zmax']))
-            elif type == 'Discrete1D':
-                force.addTabulatedFunction(name, mm.Discrete1DFunction(values))
-            elif type == 'Discrete2D':
-                force.addTabulatedFunction(name, mm.Discrete2DFunction(params['xsize'], params['ysize'], values))
-            elif type == 'Discrete3D':
-                force.addTabulatedFunction(name, mm.Discrete2DFunction(params['xsize'], params['ysize'], params['zsize'], values))
+        _createFunctions(force, self.functions)
         for atom in data.atoms:
             values = self.params.getAtomParameters(atom, data)
             force.addParticle(values)
@@ -2723,6 +2722,167 @@ class CustomGBGenerator(object):
 parsers["CustomGBForce"] = CustomGBGenerator.parseElement
 
 
+## @private
+class CustomHbondGenerator(object):
+    """A CustomHbondGenerator constructs a CustomHbondForce."""
+
+    def __init__(self, forcefield):
+        self.ff = forcefield
+        self.donorTypes1 = []
+        self.donorTypes2 = []
+        self.donorTypes3 = []
+        self.acceptorTypes1 = []
+        self.acceptorTypes2 = []
+        self.acceptorTypes3 = []
+        self.globalParams = {}
+        self.perDonorParams = []
+        self.perAcceptorParams = []
+        self.donorParamValues = []
+        self.acceptorParamValues = []
+        self.functions = []
+
+    @staticmethod
+    def parseElement(element, ff):
+        generator = CustomHbondGenerator(ff)
+        ff.registerGenerator(generator)
+        generator.energy = element.attrib['energy']
+        generator.bondCutoff = int(element.attrib['bondCutoff'])
+        generator.particlesPerDonor = int(element.attrib['particlesPerDonor'])
+        generator.particlesPerAcceptor = int(element.attrib['particlesPerAcceptor'])
+        if generator.particlesPerDonor < 1 or generator.particlesPerDonor > 3:
+            raise ValueError('Illegal value for particlesPerDonor for CustomHbondForce')
+        if generator.particlesPerAcceptor < 1 or generator.particlesPerAcceptor > 3:
+            raise ValueError('Illegal value for particlesPerAcceptor for CustomHbondForce')
+        for param in element.findall('GlobalParameter'):
+            generator.globalParams[param.attrib['name']] = float(param.attrib['defaultValue'])
+        for param in element.findall('PerDonorParameter'):
+            generator.perDonorParams.append(param.attrib['name'])
+        for param in element.findall('PerAcceptorParameter'):
+            generator.perAcceptorParams.append(param.attrib['name'])
+        for donor in element.findall('Donor'):
+            types = ff._findAtomTypes(donor.attrib, generator.particlesPerDonor)
+            if None not in types:
+                generator.donorTypes1.append(types[0])
+                if len(types) > 1:
+                    generator.donorTypes2.append(types[1])
+                if len(types) > 2:
+                    generator.donorTypes3.append(types[2])
+                generator.donorParamValues.append([float(donor.attrib[param]) for param in generator.perDonorParams])
+        for acceptor in element.findall('Acceptor'):
+            types = ff._findAtomTypes(acceptor.attrib, generator.particlesPerAcceptor)
+            if None not in types:
+                generator.acceptorTypes1.append(types[0])
+                if len(types) > 1:
+                    generator.acceptorTypes2.append(types[1])
+                if len(types) > 2:
+                    generator.acceptorTypes3.append(types[2])
+                generator.acceptorParamValues.append([float(acceptor.attrib[param]) for param in generator.perAcceptorParams])
+        generator.functions += _parseFunctions(element)
+
+    def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
+        methodMap = {NoCutoff:mm.CustomHbondForce.NoCutoff,
+                     CutoffNonPeriodic:mm.CustomHbondForce.CutoffNonPeriodic,
+                     CutoffPeriodic:mm.CustomHbondForce.CutoffPeriodic}
+        if nonbondedMethod not in methodMap:
+            raise ValueError('Illegal nonbonded method for CustomNonbondedForce')
+        force = mm.CustomHbondForce(self.energy)
+        sys.addForce(force)
+        for param in self.globalParams:
+            force.addGlobalParameter(param, self.globalParams[param])
+        for param in self.perDonorParams:
+            force.addPerDonorParameter(param)
+        for param in self.perAcceptorParams:
+            force.addPerAcceptorParameter(param)
+        _createFunctions(force, self.functions)
+        force.setNonbondedMethod(methodMap[nonbondedMethod])
+        force.setCutoffDistance(nonbondedCutoff)
+
+        # Add donors.
+
+        if self.particlesPerDonor == 1:
+            for atom in data.atoms:
+                type1 = data.atomType[atom]
+                for i in range(len(self.donorTypes1)):
+                    types1 = self.donorTypes1[i]
+                    if type1 in self.donorTypes1[i]:
+                        force.addDonor(atom.index, -1, -1, self.donorParamValues[i])
+        elif self.particlesPerDonor == 2:
+            for bond in data.bonds:
+                type1 = data.atomType[data.atoms[bond.atom1]]
+                type2 = data.atomType[data.atoms[bond.atom2]]
+                for i in range(len(self.donorTypes1)):
+                    types1 = self.donorTypes1[i]
+                    types2 = self.donorTypes2[i]
+                    if (type1 in types1 and type2 in types2) or (type1 in types2 and type2 in types1):
+                        force.addDonor(bond[0], bond[1], -1, self.donorParamValues[i])
+        else:
+            for angle in data.angles:
+                type1 = data.atomType[data.atoms[angle[0]]]
+                type2 = data.atomType[data.atoms[angle[1]]]
+                type3 = data.atomType[data.atoms[angle[2]]]
+                for i in range(len(self.donorTypes1)):
+                    types1 = self.donorTypes1[i]
+                    types2 = self.donorTypes2[i]
+                    types3 = self.donorTypes3[i]
+                    if (type1 in types1 and type2 in types2 and type3 in types3) or (type1 in types3 and type2 in types2 and type3 in types1):
+                        force.addDonor(angle[0], angle[1], angle[2], self.donorParamValues[i])
+
+        # Add acceptors.
+
+        if self.particlesPerAcceptor == 1:
+            for atom in data.atoms:
+                type1 = data.atomType[atom]
+                for i in range(len(self.acceptorTypes1)):
+                    types1 = self.acceptorTypes1[i]
+                    if type1 in self.acceptorTypes1[i]:
+                        force.addAcceptor(atom.index, -1, -1, self.acceptorParamValues[i])
+        elif self.particlesPerAcceptor == 2:
+            for bond in data.bonds:
+                type1 = data.atomType[data.atoms[bond.atom1]]
+                type2 = data.atomType[data.atoms[bond.atom2]]
+                for i in range(len(self.acceptorTypes1)):
+                    types1 = self.acceptorTypes1[i]
+                    types2 = self.acceptorTypes2[i]
+                    if (type1 in types1 and type2 in types2) or (type1 in types2 and type2 in types1):
+                        force.addAcceptor(bond.atom1, bond.atom2, -1, self.acceptorParamValues[i])
+        else:
+            for angle in data.angles:
+                type1 = data.atomType[data.atoms[angle[0]]]
+                type2 = data.atomType[data.atoms[angle[1]]]
+                type3 = data.atomType[data.atoms[angle[2]]]
+                for i in range(len(self.acceptorTypes1)):
+                    types1 = self.acceptorTypes1[i]
+                    types2 = self.acceptorTypes2[i]
+                    types3 = self.acceptorTypes3[i]
+                    if (type1 in types1 and type2 in types2 and type3 in types3) or (type1 in types3 and type2 in types2 and type3 in types1):
+                        force.addAcceptor(angle[0], angle[1], angle[2], self.acceptorParamValues[i])
+
+        # Add exclusions.
+
+        for donor in range(force.getNumDonors()):
+            (d1, d2, d3, params) = force.getDonorParameters(donor)
+            outerAtoms = set((d1, d2, d3))
+            if -1 in outerAtoms:
+                outerAtoms.remove(-1)
+            excludedAtoms = set(outerAtoms)
+            for i in range(self.bondCutoff):
+                newOuterAtoms = set()
+                for atom in outerAtoms:
+                    for bond in data.atomBonds[atom]:
+                        b = data.bonds[bond]
+                        bondedAtom = (b.atom2 if b.atom1 == atom else b.atom1)
+                        if bondedAtom not in excludedAtoms:
+                            newOuterAtoms.add(bondedAtom)
+                            excludedAtoms.add(bondedAtom)
+                outerAtoms = newOuterAtoms
+            for acceptor in range(force.getNumAcceptors()):
+                (a1, a2, a3, params) = force.getAcceptorParameters(acceptor)
+                if a1 in excludedAtoms or a2 in excludedAtoms or a3 in excludedAtoms:
+                    force.addExclusion(donor, acceptor)
+
+parsers["CustomHbondForce"] = CustomHbondGenerator.parseElement
+
+
 ## @private
 class CustomManyParticleGenerator(object):
     """A CustomManyParticleGenerator constructs a CustomManyParticleForce."""

wrappers/python/tests/TestGenerators.py

+import unittest
+from validateConstraints import *
+from simtk.openmm.app import *
+from simtk.openmm import *
+from simtk.unit import *
+import simtk.openmm.app.element as elem
+import simtk.openmm.app.forcefield as forcefield
+import math
+try:
+    from cStringIO import StringIO
+except ImportError:
+    from io import StringIO
+import os
+import warnings
+
+class TestGenerators(unittest.TestCase):
+    """Test the various generators found in forcefield.py."""
+
+    def setUp(self):
+        # alanine dipeptide with implicit water
+        self.pdb1 = PDBFile('systems/alanine-dipeptide-implicit.pdb')
+
+
+    def test_CustomHbondGenerator(self):
+        """Test the generator for CustomHbondForce."""
+
+        for bondCutoff in range(4):
+            xml = """
+<ForceField>
+ <CustomHbondForce energy="a*b*distance(a1,d1)" particlesPerDonor="3" particlesPerAcceptor="2" bondCutoff="%d">
+  <PerDonorParameter name="a"/>
+  <PerAcceptorParameter name="b"/>
+  <Donor class1="C" class2="N" class3="H" a="3"/>
+  <Acceptor class1="C" class2="O" b="2"/>
+  <Function name="test" min="1" max="2" type="Continuous1D">
+   0 1 2 3 4 5
+  </Function>
+ </CustomHbondForce>
+</ForceField>""" % bondCutoff
+            ff = ForceField('amber99sb.xml', StringIO(xml))
+            system = ff.createSystem(self.pdb1.topology)
+            hbond = [f for f in system.getForces() if isinstance(f, CustomHbondForce)][0]
+            self.assertEqual(1, hbond.getNumPerDonorParameters())
+            self.assertEqual(1, hbond.getNumPerAcceptorParameters())
+            self.assertEqual('a', hbond.getPerDonorParameterName(0))
+            self.assertEqual('b', hbond.getPerAcceptorParameterName(0))
+            self.assertEqual(1, hbond.getNumTabulatedFunctions())
+            expectedDonors = [(4,6,7), (14,16,17)]
+            expectedAcceptors = [(4,5,-1), (14,15,-1)]
+            self.assertEqual(len(expectedDonors), hbond.getNumDonors())
+            self.assertEqual(len(expectedAcceptors), hbond.getNumAcceptors())
+            for i in range(hbond.getNumDonors()):
+                atom1, atom2, atom3, params = hbond.getDonorParameters(i)
+                self.assertTrue((atom1, atom2, atom3) in expectedDonors)
+                self.assertEqual((3.0,), params)
+            for i in range(hbond.getNumAcceptors()):
+                atom1, atom2, atom3, params = hbond.getAcceptorParameters(i)
+                self.assertTrue((atom1, atom2, atom3) in expectedAcceptors)
+                self.assertEqual((2.0,), params)
+            expectedExclusions = [(0,0), (1,1)]
+            if bondCutoff >= 2:
+                expectedExclusions.append((0,1))
+            if bondCutoff >= 3:
+                expectedExclusions.append((1,0))
+            self.assertEqual(len(expectedExclusions), hbond.getNumExclusions())
+            for i in range(hbond.getNumExclusions()):
+                self.assertTrue(tuple(hbond.getExclusionParticles(i)) in expectedExclusions)
+
+
+if __name__ == '__main__':
+    unittest.main()