Allow multiple matching templates if they assign identical parameters

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

@@ -888,7 +888,23 @@ class ForceField(object):
                 template = allMatches[0][0]
                 matches = allMatches[0][1]
             elif len(allMatches) > 1:
+                # We found multiple matches.  This is OK if and only if they assign identical types and parameters to all atoms.
+                t1, m1 = allMatches[0]
+                atoms1 = [t1.atoms[m] for m in m1]
+                allIdentical = True
+                for t2, m2 in allMatches[1:]:
+                    atoms2 = [t2.atoms[m] for m in m2]
+                    if any(a1.type != a2.type or a1.parameters != a2.parameters for a1,a2 in zip(atoms1, atoms2)):
+                        allIdentical = 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 t1.virtualSites != t2.virtualSites:
+                        allIdentical = False
+                if not allIdentical:
                     raise Exception('Multiple 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

@@ -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."""