Merge pull request #1359 from jchodera/master

Add a method for retrieving list of matching templates to ForceField

docs-source/usersguide/application.rst

@@ -1976,7 +1976,7 @@ residue template defined in the :class:`ForceField`.
     forcefield = ForceField('amber99sb.xml', 'tip3p.xml')
     unmatched_residues = forcefield.getUnmatchedResidues(topology)
 
-This is useful for idenfitying issues with prepared systems, debugging issues
+This is useful for identifying issues with prepared systems, debugging issues
 with residue template definitions, or identifying which additional residues need
 to be parameterized.
 
@@ -1994,6 +1994,17 @@ residues and empty residue templates using :method:`generateTemplatesForUnmatche
         # Register the template with the forcefield.
         forcefield.registerResidueTemplate(template)
 
+If you find that templates seem to be incorrectly matched, another useful
+function :method:`getMatchingTemplates()` can help you identify which templates
+are being matched:
+::
+
+    pdb = PDBFile('input.pdb')
+    forcefield = ForceField('amber99sb.xml', 'tip3p.xml')
+    templates = forcefield.getMatchingTemplates(topology)
+    for (residue, template) in zip(pdb.topology.residues(), templates):
+        print "Residue %d %s matched template %s" % (residue.id, residue.name, template.name)
+
 <HarmonicBondForce>
 ===================
 

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

@@ -616,16 +616,47 @@ class ForceField(object):
         # Find the template matching each residue, compiling a list of residues for which no templates are available.
         bondedToAtom = self._buildBondedToAtomList(topology)
         unmatched_residues = list() # list of unmatched residues
-        for chain in topology.chains():
-            for res in chain.residues():
-                # Attempt to match one of the existing templates.
-                [template, matches] = self._getResidueTemplateMatches(res, bondedToAtom)
-                if matches is None:
-                    # No existing templates match.
-                    unmatched_residues.append(res)
+        for res in topology.residues():
+            # Attempt to match one of the existing templates.
+            [template, matches] = self._getResidueTemplateMatches(res, bondedToAtom)
+            if matches is None:
+                # No existing templates match.
+                unmatched_residues.append(res)
 
         return unmatched_residues
 
+    def getMatchingTemplates(self, topology):
+        """Return a list of forcefield residue templates matching residues in the specified topology.
+
+        .. CAUTION:: This method is experimental, and its API is subject to change.
+
+        Parameters
+        ----------
+        topology : Topology
+            The Topology whose residues are to be checked against the forcefield residue templates.
+
+        Returns
+        -------
+        templates : list of _TemplateData
+            List of forcefield residue templates corresponding to residues in the topology.
+            templates[index] is template corresponding to residue `index` in topology.residues()
+
+        This method may be of use in debugging issues related to parameter assignment.
+        """
+        # Find the template matching each residue, compiling a list of residues for which no templates are available.
+        bondedToAtom = self._buildBondedToAtomList(topology)
+        templates = list() # list of templates matching the corresponding residues
+        for res in topology.residues():
+            # Attempt to match one of the existing templates.
+            [template, matches] = self._getResidueTemplateMatches(res, bondedToAtom)
+            # Raise an exception if we have found no templates that match.
+            if matches is None:
+                raise ValueError('No template found for residue %d (%s).  %s' % (res.index+1, res.name, _findMatchErrors(self, res)))
+            else:
+                templates.append(template)
+
+        return templates
+
     def generateTemplatesForUnmatchedResidues(self, topology):
         """Generate forcefield residue templates for residues in specified topology for which no forcefield templates are available.
 

wrappers/python/tests/TestForceField.py

@@ -431,6 +431,22 @@ class TestForceField(unittest.TestCase):
         system = forcefield.createSystem(pdb.topology, nonbondedMethod=NoCutoff)
         # TODO: Test energies are finite?
 
+    def test_getMatchingTemplates(self):
+        """Test retrieval of list of templates that match residues in a topology."""
+
+        # Load the PDB file.
+        pdb = PDBFile(os.path.join('systems', 'ala_ala_ala.pdb'))
+        # Create a ForceField object.
+        forcefield = ForceField('amber99sb.xml')
+        # Get list of matching residue templates.
+        templates = forcefield.getMatchingTemplates(pdb.topology)
+        # Check results.
+        residues = [ residue for residue in pdb.topology.residues() ]
+        self.assertEqual(len(templates), len(residues))
+        self.assertEqual(templates[0].name, 'NALA')
+        self.assertEqual(templates[1].name, 'ALA')
+        self.assertEqual(templates[2].name, 'CALA')
+
 class AmoebaTestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method with the AMOEBA forcefield."""