Merge pull request #856 from swails/fix_hangles

Fix HAngles with Amber prmtop files

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

@@ -73,9 +73,8 @@ class AmberPrmtopFile(object):
 
     def __init__(self, file):
         """Load a prmtop file."""
-        top = Topology()
         ## The Topology read from the prmtop file
-        self.topology = top
+        self.topology = top = Topology()
         self.elements = []
 
         # Load the prmtop file
@@ -229,7 +228,7 @@ class AmberPrmtopFile(object):
         elif implicitSolvent is None:
             implicitSolventKappa = 0.0
 
-        sys = amber_file_parser.readAmberSystem(prmtop_loader=self._prmtop, shake=constraintString,
+        sys = amber_file_parser.readAmberSystem(self.topology, prmtop_loader=self._prmtop, shake=constraintString,
                         nonbondedCutoff=nonbondedCutoff, nonbondedMethod=methodMap[nonbondedMethod],
                         flexibleConstraints=False, gbmodel=implicitString, soluteDielectric=soluteDielectric,
                         solventDielectric=solventDielectric, implicitSolventKappa=implicitSolventKappa,

wrappers/python/simtk/openmm/app/internal/amber_file_parser.py

@@ -122,6 +122,7 @@ class PrmtopLoader(object):
 
         with open(inFilename, 'r') as fIn:
             for line in fIn:
+                if line[0] == '%':
                     if line.startswith('%VERSION'):
                         tag, self._prmtopVersion = line.rstrip().split(None, 1)
                     elif line.startswith('%FLAG'):
@@ -134,7 +135,7 @@ class PrmtopLoader(object):
                         index1=format.index(')')
                         format = format[index0+1:index1]
                         m = FORMAT_RE_PATTERN.search(format)
-                    self._raw_format[self._flags[-1]] = (format, m.group(1), m.group(2), m.group(3), m.group(4))
+                        self._raw_format[self._flags[-1]] = (format, m.group(1), m.group(2), int(m.group(3)), m.group(4))
                     elif line.startswith('%COMMENT'):
                         continue
                 elif self._flags \
@@ -144,8 +145,7 @@ class PrmtopLoader(object):
                 else:
                     flag=self._flags[-1]
                     (format, numItems, itemType,
-                     itemLength, itemPrecision) = self._getFormat(flag)
-                    iLength=int(itemLength)
+                     iLength, itemPrecision) = self._getFormat(flag)
                     line = line.rstrip()
                     for index in range(0, len(line), iLength):
                         item = line[index:index+iLength]
@@ -651,7 +651,7 @@ class PrmtopLoader(object):
 # AMBER System builder (based on, but not identical to, systemManager from 'zander')
 #=============================================================================================
 
-def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmodel=None,
+def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=None, gbmodel=None,
           soluteDielectric=1.0, solventDielectric=78.5,
           implicitSolventKappa=0.0*(1/units.nanometer), nonbondedCutoff=None,
           nonbondedMethod='NoCutoff', scee=None, scnb=None, mm=None, verbose=False,
@@ -659,6 +659,9 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
     """
     Create an OpenMM System from an Amber prmtop file.
 
+    REQUIRED ARGUMENT
+      topology (forcefield.Topology) The topology for the system that is about
+      to be created
     ARGUMENTS (specify  one or the other, but not both)
       prmtop_filename (String) - name of Amber prmtop file (new-style only)
       prmtop_loader (PrmtopLoader) - the loaded prmtop file
@@ -739,7 +742,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
         system.addParticle(mass)
 
     # Add constraints.
-    isWater = [prmtop.getResidueLabel(i) in ('WAT', 'TP4', 'TP5', 'T4E') for i in range(prmtop.getNumAtoms())]
+    isWater = [prmtop.getResidueLabel(i) in ('WAT', 'HOH', 'TP4', 'TP5', 'T4E') for i in range(prmtop.getNumAtoms())]
     if shake in ('h-bonds', 'all-bonds', 'h-angles'):
         for (iAtom, jAtom, k, rMin) in prmtop.getBondsWithH():
             system.addConstraint(iAtom, jAtom, rMin)
@@ -774,13 +777,14 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
             distance = c[2].value_in_unit(units.nanometer)
             atomConstraints[c[0]].append((c[1], distance))
             atomConstraints[c[1]].append((c[0], distance))
+    topatoms = list(topology.atoms())
     for (iAtom, jAtom, kAtom, k, aMin) in prmtop.getAngles():
         if shake == 'h-angles':
-            type1 = prmtop.getAtomType(iAtom)
-            type2 = prmtop.getAtomType(jAtom)
-            type3 = prmtop.getAtomType(kAtom)
-            numH = len([type for type in (type1, type3) if type.startswith('H')])
-            constrained = (numH == 2 or (numH == 1 and type2.startswith('O')))
+            atomI = topatoms[iAtom]
+            atomJ = topatoms[jAtom]
+            atomK = topatoms[kAtom]
+            numH = ((atomI.element.atomic_number == 1) + (atomK.element.atomic_number == 1))
+            constrained = (numH == 2 or (numH == 1 and atomJ.element is elem.oxygen))
         else:
             constrained = False
         if constrained:

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -12,6 +12,7 @@ prmtop2 = AmberPrmtopFile('systems/alanine-dipeptide-implicit.prmtop')
 prmtop3 = AmberPrmtopFile('systems/ff14ipq.parm7')
 prmtop4 = AmberPrmtopFile('systems/Mg_water.prmtop')
 prmtop5 = AmberPrmtopFile('systems/tz2.truncoct.parm7')
+prmtop6 = AmberPrmtopFile('systems/gaffwat.parm7')
 inpcrd3 = AmberInpcrdFile('systems/ff14ipq.rst7')
 inpcrd4 = AmberInpcrdFile('systems/Mg_water.inpcrd')
 
@@ -200,6 +201,24 @@ class TestAmberPrmtopFile(unittest.TestCase):
         # Amber using this force field.
         self.assertAlmostEqual(-7042.3903307/ene, 1, places=3)
 
+    def test_HAngle(self):
+        """ Test that HAngle constraints are properly handled for all hydrogens """
+        system = prmtop6.createSystem(nonbondedMethod=PME,
+                                      nonbondedCutoff=1*nanometers,
+                                      constraints=HBonds)
+        self.assertEqual(system.getForce(0).getNumBonds(), 0)
+        self.assertEqual(system.getNumParticles(), 3000)
+        self.assertEqual(system.getNumConstraints(), 2000)
+        self.assertEqual(system.getForce(1).getNumAngles(), 1000)
+
+        system = prmtop6.createSystem(nonbondedMethod=PME,
+                                      nonbondedCutoff=1*nanometers,
+                                      constraints=HAngles)
+        self.assertEqual(system.getForce(0).getNumBonds(), 0)
+        self.assertEqual(system.getNumParticles(), 3000)
+        self.assertEqual(system.getNumConstraints(), 3000)
+        self.assertEqual(system.getForce(1).getNumAngles(), 0)
+
     def test_LJ1264(self):
         """Test prmtop with 12-6-4 vdW potential implemented"""
         system = prmtop4.createSystem(nonbondedMethod=PME,