Merge pull request #1605 from peastman/hbondbug

Bug fixes to building CustomHbondForce

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

@@ -2760,7 +2760,7 @@ class CustomHbondGenerator(object):
         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)
+            types = ff._findAtomTypes(donor.attrib, 3)[:generator.particlesPerDonor]
             if None not in types:
                 generator.donorTypes1.append(types[0])
                 if len(types) > 1:
@@ -2769,7 +2769,7 @@ class CustomHbondGenerator(object):
                     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)
+            types = ff._findAtomTypes(acceptor.attrib, 3)[:generator.particlesPerAcceptor]
             if None not in types:
                 generator.acceptorTypes1.append(types[0])
                 if len(types) > 1:
@@ -2813,8 +2813,10 @@ class CustomHbondGenerator(object):
                 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])
+                    if type1 in types1 and type2 in types2:
+                        force.addDonor(bond.atom1, bond.atom2, -1, self.donorParamValues[i])
+                    elif type1 in types2 and type2 in types1:
+                        force.addDonor(bond.atom2, bond.atom1, -1, self.donorParamValues[i])
         else:
             for angle in data.angles:
                 type1 = data.atomType[data.atoms[angle[0]]]
@@ -2843,8 +2845,10 @@ class CustomHbondGenerator(object):
                 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):
+                    if type1 in types1 and type2 in types2:
                         force.addAcceptor(bond.atom1, bond.atom2, -1, self.acceptorParamValues[i])
+                    elif type1 in types2 and type2 in types1:
+                        force.addAcceptor(bond.atom2, bond.atom1, -1, self.acceptorParamValues[i])
         else:
             for angle in data.angles:
                 type1 = data.atomType[data.atoms[angle[0]]]

wrappers/python/tests/TestGenerators.py

@@ -67,5 +67,38 @@ class TestGenerators(unittest.TestCase):
                 self.assertTrue(tuple(hbond.getExclusionParticles(i)) in expectedExclusions)
 
 
+    def test_CustomHbondGenerator2(self):
+        """Test the generator for CustomHbondForce with different parameters."""
+
+        xml = """
+<ForceField>
+ <CustomHbondForce energy="a*b*distance(a1,d1)" particlesPerDonor="2" particlesPerAcceptor="1" bondCutoff="0">
+  <PerDonorParameter name="a"/>
+  <PerAcceptorParameter name="b"/>
+  <Donor class1="N" class2="H" a="3"/>
+  <Acceptor class1="O" b="2"/>
+ </CustomHbondForce>
+</ForceField>"""
+        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))
+        expectedDonors = [(6,7,-1), (16,17,-1)]
+        expectedAcceptors = [(5,-1,-1), (15,-1,-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)
+
+
 if __name__ == '__main__':
     unittest.main()