Add methods to get unmatched residues without forcefield templates.

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

@@ -569,6 +569,89 @@ class ForceField(object):
                     break
         return [template, matches]
 
+    def getUnmatchedResidues(self, topology):
+        """Return a list of Residue objects from specified topology for which no forcefield templates are available.
+
+        Parameters
+        ----------
+        topology : Topology
+            The Topology for which to create a System
+
+        Returns
+        -------
+        unmatched_residues : list of Residue
+            List of residue templates from `topology` for which no forcefield residue templates are available.
+            Note that multiple instances of the same residue appearing at different points in the topology may be returned.
+
+        This method may be of use in generating missing residue templates or diagnosing parameterization failures.
+        """
+        # Find the template matching each residue, compiling a list of residues for which no templates are available.
+        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)
+
+        return unmatched_residues
+
+    def getUniqueUnmatchedResidues(self, topology):
+        """Returns a unique list of Residue objects from specified topology for which no forcefield templates are available.
+
+        Parameters
+        ----------
+        topology : Topology
+            The Topology for which to create a System
+
+        Returns
+        -------
+        unmatched_residues : list of Residue
+            List of residue templates from `topology` for which no forcefield residue templates are available.
+            Note that only a single instances each missing residue type will be returned.
+
+        This method may be of use in generating missing residue templates.
+        """
+        # Get a non-unique list of unmatched residues.
+        unmatched_residues = self.getUnmatchedResidues(topology)
+        # Record which atoms are bonded to each other atom
+        bondedToAtom = []
+        for atom in topology.atoms():
+            bondedToAtom.append(set())
+        for bond in topology.bonds():
+            bondedToAtom[bond.atom1.index].add(bond.atom2.index)
+            bondedToAtom[bond.atom2.index].add(bond.atom1.index)
+        # Generate a unique list of unmatched residues by comparing fingerprints.
+        unique_unmatched_residues = list()
+        signatures = set()
+        for residue in unmatched_residues:
+            signature = _createResidueSignature([ atom.element for atom in residue.atoms() ])
+            is_unique = True
+            if signature in signatures:
+                # Signature is the same as an existing residue; check connectivity.
+                template = ForceField._TemplateData(residue.name)
+                for atom in residue.atoms():
+                    template.atoms.append(ForceField._TemplateAtomData(atom.name, None, atom.element))
+                for (atom1,atom2) in residue.internal_bonds():
+                    template.addBondByName(atom1.name, atom2.name)
+                residue_atoms = [ atom for atom in residue.atoms() ]
+                for (atom1,atom2) in residue.external_bonds():
+                    if atom1 in residue_atoms:
+                        template.addExternalBondByName(atom1.name)
+                    elif atom2 in residue_atoms:
+                        template.addExternalBondByName(atom2.name)
+                for check_residue in unique_unmatched_residues:
+                    matches = _matchResidue(check_residue, template, bondedToAtom)
+                    if matches is not None:
+                        is_unique = False
+            if is_unique:
+                # Residue is unique.
+                unique_unmatched_residues.append(residue)
+                signatures.add(signature)
+
+        return unique_unmatched_residues
+
     def createSystem(self, topology, nonbondedMethod=NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
                      constraints=None, rigidWater=True, removeCMMotion=True, hydrogenMass=None, **args):
         """Construct an OpenMM System representing a Topology with this force field.

wrappers/python/tests/TestForceField.py

@@ -13,9 +13,9 @@ else:
 
 class TestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method."""
- 
+
     def setUp(self):
-        """Set up the tests by loading the input pdb files and force field 
+        """Set up the tests by loading the input pdb files and force field
         xml files.
 
         """
@@ -33,16 +33,16 @@ class TestForceField(unittest.TestCase):
     def test_NonbondedMethod(self):
         """Test all five options for the nonbondedMethod parameter."""
 
-        methodMap = {NoCutoff:NonbondedForce.NoCutoff, 
-                     CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic, 
-                     CutoffPeriodic:NonbondedForce.CutoffPeriodic, 
+        methodMap = {NoCutoff:NonbondedForce.NoCutoff,
+                     CutoffNonPeriodic:NonbondedForce.CutoffNonPeriodic,
+                     CutoffPeriodic:NonbondedForce.CutoffPeriodic,
                      Ewald:NonbondedForce.Ewald, PME: NonbondedForce.PME}
         for method in methodMap:
             system = self.forcefield1.createSystem(self.pdb1.topology,
                                                   nonbondedMethod=method)
             forces = system.getForces()
-            self.assertTrue(any(isinstance(f, NonbondedForce) and 
-                                f.getNonbondedMethod()==methodMap[method] 
+            self.assertTrue(any(isinstance(f, NonbondedForce) and
+                                f.getNonbondedMethod()==methodMap[method]
                                 for f in forces))
 
     def test_DispersionCorrection(self):
@@ -50,7 +50,7 @@ class TestForceField(unittest.TestCase):
 
         for useDispersionCorrection in [True, False]:
             system = self.forcefield1.createSystem(self.pdb1.topology,
-                                                   nonbondedCutoff=2*nanometer, 
+                                                   nonbondedCutoff=2*nanometer,
                                                    useDispersionCorrection=useDispersionCorrection)
 
             for force in system.getForces():
@@ -62,8 +62,8 @@ class TestForceField(unittest.TestCase):
 
         for method in [CutoffNonPeriodic, CutoffPeriodic, Ewald, PME]:
             system = self.forcefield1.createSystem(self.pdb1.topology,
-                                                   nonbondedMethod=method, 
-                                                   nonbondedCutoff=2*nanometer, 
+                                                   nonbondedMethod=method,
+                                                   nonbondedCutoff=2*nanometer,
                                                    constraints=HBonds)
             cutoff_distance = 0.0*nanometer
             cutoff_check = 2.0*nanometer
@@ -81,29 +81,29 @@ class TestForceField(unittest.TestCase):
 
     def test_RigidWaterAndConstraints(self):
         """Test all eight options for the constraints and rigidWater parameters."""
-        
+
         topology = self.pdb1.topology
         for constraints_value in [None, HBonds, AllBonds, HAngles]:
             for rigidWater_value in [True, False]:
                 system = self.forcefield1.createSystem(topology,
-                                                       constraints=constraints_value, 
+                                                       constraints=constraints_value,
                                                        rigidWater=rigidWater_value)
-                validateConstraints(self, topology, system, 
+                validateConstraints(self, topology, system,
                                     constraints_value, rigidWater_value)
 
     def test_ImplicitSolvent(self):
-        """Test the four types of implicit solvents using the implicitSolvent 
+        """Test the four types of implicit solvents using the implicitSolvent
         parameter.
 
         """
-        
-        topology = self.pdb2.topology 
+
+        topology = self.pdb2.topology
         system = self.forcefield2.createSystem(topology)
         forces = system.getForces()
         self.assertTrue(any(isinstance(f, GBSAOBCForce) for f in forces))
 
     def test_ImplicitSolventParameters(self):
-        """Test that solventDielectric and soluteDielectric are passed correctly 
+        """Test that solventDielectric and soluteDielectric are passed correctly
         for the different types of implicit solvent.
 
         """
@@ -121,12 +121,12 @@ class TestForceField(unittest.TestCase):
                     found_matching_solute_dielectric = True
             if isinstance(force, NonbondedForce):
                 self.assertEqual(force.getReactionFieldDielectric(), 1.0)
-        self.assertTrue(found_matching_solvent_dielectric and 
+        self.assertTrue(found_matching_solvent_dielectric and
                         found_matching_solute_dielectric)
 
     def test_HydrogenMass(self):
         """Test that altering the mass of hydrogens works correctly."""
-        
+
         topology = self.pdb1.topology
         hydrogenMass = 4*amu
         system1 = self.forcefield1.createSystem(topology)
@@ -138,10 +138,10 @@ class TestForceField(unittest.TestCase):
         totalMass1 = sum([system1.getParticleMass(i) for i in range(system1.getNumParticles())]).value_in_unit(amu)
         totalMass2 = sum([system2.getParticleMass(i) for i in range(system2.getNumParticles())]).value_in_unit(amu)
         self.assertAlmostEqual(totalMass1, totalMass2)
-    
+
     def test_Forces(self):
         """Compute forces and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""
-        
+
         pdb = PDBFile('systems/lysozyme-implicit.pdb')
         system = self.forcefield2.createSystem(pdb.topology)
         integrator = VerletIntegrator(0.001)
@@ -156,10 +156,10 @@ class TestForceField(unittest.TestCase):
             if diff > 0.1 and diff/norm(f1) > 1e-3:
                 numDifferences += 1
         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})
@@ -181,11 +181,11 @@ class TestForceField(unittest.TestCase):
         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')
@@ -347,11 +347,60 @@ class TestForceField(unittest.TestCase):
             system = forcefield.createSystem(pdb.topology, nonbondedMethod=test['nonbondedMethod'])
             # TODO: Test energies are finite?
 
+    def test_getUnmatchedResidues(self):
+        """Test retrieval of list of residues for which no templates are available."""
+        #
+        # Test where we generate parameters for only a ligand.
+        #
+
+        # Load the PDB file.
+        pdb = PDBFile(os.path.join('systems', 'T4-lysozyme-L99A-p-xylene-implicit.pdb'))
+        # Create a ForceField object.
+        forcefield = ForceField('amber99sb.xml', 'tip3p.xml')
+        # Get list of unmatched residues.
+        unmatched_residues = forcefield.getUnmatchedResidues(pdb.topology)
+        # Check results.
+        self.assertEqual(len(unmatched_residues), 1)
+        self.assertEqual(unmatched_residues[0].name, 'TMP')
+        self.assertEqual(unmatched_residues[0].id, 163)
+
+        # Load the PDB file.
+        pdb = PDBFile(os.path.join('systems', 'ala_ala_ala.pdb'))
+        # Create a ForceField object.
+        forcefield = ForceField('tip3p.xml')
+        # Get list of unmatched residues.
+        unmatched_residues = forcefield.getUnmatchedResidues(pdb.topology)
+        # Check results.
+        self.assertEqual(len(unmatched_residues), 3)
+        self.assertEqual(unmatched_residues[0].name, 'ALA')
+        self.assertEqual(unmatched_residues[0].chain, 'X')
+        self.assertEqual(unmatched_residues[0].id, 1)
+
+    def test_getUniqueUnmatchedResidues(self):
+        """Test retrieval of list of residues for which no templates are available."""
+        #
+        # Test where we generate parameters for only a ligand.
+        #
+
+        # Load the PDB file.
+        pdb = PDBFile(os.path.join('systems', 'nacl-water.pdb'))
+        # Create a ForceField object.
+        forcefield = ForceField('tip3p.xml')
+        # Get list of unmatched residues.
+        unmatched_residues = forcefield.getUnmatchedResidues(pdb.topology)
+        unique_unmatched_residues = forcefield.getUniqueUnmatchedResidues(pdb.topology)
+        # Check results.
+        self.assertEqual(len(unmatched_residues), 14)
+        self.assertEqual(len(unique_unmatched_residues), 2)
+        unique_names = set([ residue.name for residue in unique_unmatched_residues ])
+        self.assertTrue('NA' in unique_names)
+        self.assertTrue('CL' in unique_names)
+
 class AmoebaTestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method with the AMOEBA forcefield."""
- 
+
     def setUp(self):
-        """Set up the tests by loading the input pdb files and force field 
+        """Set up the tests by loading the input pdb files and force field
         xml files.
 
         """
@@ -436,4 +485,3 @@ class AmoebaTestForceField(unittest.TestCase):
 
 if __name__ == '__main__':
     unittest.main()
-