Merge pull request #2452 from peastman/template

Allow multiple matching templates if they assign identical parameters

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

@@ -632,6 +632,29 @@ class ForceField(object):
             atom = self.getAtomIndexByName(atom_name)
             self.addExternalBond(atom)
 
+        def areParametersIdentical(self, template2, matchingAtoms, matchingAtoms2):
+            """Get whether this template and another one both assign identical atom types and parameters to all atoms.
+            
+            Parameters
+            ----------
+            template2: _TemplateData
+                the template to compare this one to
+            matchingAtoms: list
+                the indices of atoms in this template that atoms of the residue are matched to
+            matchingAtoms2: list
+                the indices of atoms in template2 that atoms of the residue are matched to
+            """
+            atoms1 = [self.atoms[m] for m in matchingAtoms]
+            atoms2 = [template2.atoms[m] for m in matchingAtoms2]
+            if any(a1.type != a2.type or a1.parameters != a2.parameters for a1,a2 in zip(atoms1, atoms2)):
+                return False
+            # Properly comparing virtual sites really needs a much more complicated analysis.  This simple check
+            # could easily fail for templates containing vsites, even if they're actually identical.  Since we
+            # currently have no force fields that include both patches and vsites, I'm not going to worry about it now.
+            if self.virtualSites != template2.virtualSites:
+                return False
+            return True
+
     class _TemplateAtomData(object):
         """Inner class used to encapsulate data about an atom in a residue template definition."""
         def __init__(self, name, type, element, parameters={}):
@@ -866,7 +889,7 @@ class ForceField(object):
 
         Returns
         -------
-        template : _ForceFieldTemplate
+        template : _TemplateData
             The matching forcefield residue template, or None if no matches are found.
         matches : list
             a list specifying which atom of the template each atom of the residue
@@ -888,7 +911,13 @@ class ForceField(object):
                 template = allMatches[0][0]
                 matches = allMatches[0][1]
             elif len(allMatches) > 1:
-                raise Exception('Multiple matching templates found for residue %d (%s): %s.' % (res.index+1, res.name, ', '.join(match[0].name for match in allMatches)))
+                # We found multiple matches.  This is OK if and only if they assign identical types and parameters to all atoms.
+                t1, m1 = allMatches[0]
+                for t2, m2 in allMatches[1:]:
+                    if not t1.areParametersIdentical(t2, m1, m2):
+                        raise Exception('Multiple non-identical matching templates found for residue %d (%s): %s.' % (res.index+1, res.name, ', '.join(match[0].name for match in allMatches)))
+                template = allMatches[0][0]
+                matches = allMatches[0][1]
         return [template, matches]
 
     def _buildBondedToAtomList(self, topology):

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

@@ -69,8 +69,8 @@ class PDBFile(object):
 
         Parameters
         ----------
-        file : string
-            the name of the file to load
+        file : string or file
+            the name of the file to load.  Alternatively you can pass an open file object.
         extraParticleIdentifier : string='EP'
             if this value appears in the element column for an ATOM record, the Atom's element will be set to None to mark it as an extra particle
         """

wrappers/python/tests/TestForceField.py

@@ -477,7 +477,7 @@ class TestForceField(unittest.TestCase):
         self.assertEqual(unmatched_residues[0].chain.id, 'X')
         self.assertEqual(unmatched_residues[0].id, '1')
 
-    def test_ggenerateTemplatesForUnmatchedResidues(self):
+    def test_generateTemplatesForUnmatchedResidues(self):
         """Test generation of blank forcefield residue templates for unmatched residues."""
         #
         # Test where we generate parameters for only a ligand.
@@ -914,6 +914,33 @@ class TestForceField(unittest.TestCase):
             forcefield = ForceField(ff)
             system = forcefield.createSystem(pdb.topology)
 
+    def test_IdenticalTemplates(self):
+        """Test a case where patches produce two identical templates."""
+        ff = ForceField('charmm36.xml')
+        pdb = PDBFile(StringIO("""
+ATOM      1  N   HIS     1A   -2.670    -0.476   0.475  1.00  0.00           N
+ATOM      2  HT1 HIS     1A   -2.645    -1.336   1.036  1.00  0.00           H
+ATOM      3  HT2 HIS     1A   -2.859    -0.751  -0.532  1.00  0.00           H
+ATOM      4  HT3 HIS     1A   -3.415     0.201   0.731  1.00  0.00           H
+ATOM      5  CA  HIS     1A   -1.347     0.163   0.471  1.00  0.00           C
+ATOM      6  HA  HIS     1A   -1.111     0.506   1.479  1.00  0.00           H
+ATOM      7  CB  HIS     1A   -0.352    -0.857  -0.040  1.00  0.00           C
+ATOM      8  HB1 HIS     1A   -0.360    -1.741   0.636  1.00  0.00           H
+ATOM      9  HB2 HIS     1A   -0.640    -1.175  -1.046  1.00  0.00           H
+ATOM     10  CG  HIS     1A    1.003    -0.275  -0.063  1.00  0.00           C
+ATOM     11  CD2 HIS     1A    2.143    -0.931  -0.476  1.00  0.00           C
+ATOM     12  HD2 HIS     1A    2.217    -1.952  -0.840  1.00  0.00           H
+ATOM     13  NE2 HIS     1A    3.137    -0.024  -0.328  1.00  0.00           N
+ATOM     14  HE2 HIS     1A    4.132    -0.238  -0.565  1.00  0.00           H
+ATOM     15  CE1 HIS     1A    2.649     1.130   0.150  1.00  0.00           C
+ATOM     16  HE1 HIS     1A    3.233     2.020   0.360  1.00  0.00           H
+ATOM     17  ND1 HIS     1A    1.323     0.973   0.314  1.00  0.00           N
+ATOM     18  C   HIS     1A   -1.465     1.282  -0.497  1.00  0.00           C
+ATOM     19  OT1 HIS     1A   -2.108     2.309  -0.180  1.00  0.00           O
+ATOM     20  OT2 HIS     1A   -0.864     1.172  -1.737  1.00  0.00           O
+END"""))
+        # If the check is not done correctly, this will throw an exception.
+        ff.createSystem(pdb.topology)
 
 class AmoebaTestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method with the AMOEBA forcefield."""