Fixed errors in adding extra particles with patches

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

@@ -1295,7 +1295,7 @@ class ForceField(object):
         return sys
 
 
-    def _matchAllResiduesToTemplates(self, data, topology, residueTemplates, ignoreExternalBonds, ignoreExtraParticles=False):
+    def _matchAllResiduesToTemplates(self, data, topology, residueTemplates, ignoreExternalBonds, ignoreExtraParticles=False, recordParameters=True):
         """Return a list of which template matches each residue in the topology, and assign atom types."""
         templateForResidue = [None]*topology.getNumResidues()
         unmatchedResidues = []
@@ -1313,6 +1313,7 @@ class ForceField(object):
                 if matches is None:
                     unmatchedResidues.append(res)
                 else:
+                    if recordParameters:
                         data.recordMatchedAtomParameters(res, template, matches)
                     templateForResidue[res.index] = template
 
@@ -1342,6 +1343,7 @@ class ForceField(object):
             if matches is None:
                 raise ValueError('No template found for residue %d (%s).  %s' % (res.index+1, res.name, _findMatchErrors(self, res)))
             else:
+                if recordParameters:
                     data.recordMatchedAtomParameters(res, template, matches)
                 templateForResidue[res.index] = template
         return templateForResidue

wrappers/python/simtk/openmm/app/internal/compiled.pyx

@@ -99,7 +99,7 @@ def matchResidueToTemplate(res, template, bondedToAtom, bint ignoreExternalBonds
         templateBondedTo = {}
         for i, atom in enumerate(template.atoms):
             if atom.element is not None:
-                templateBondedTo[atom] = [j for j in atom.bondedTo if template.atoms[j].element is not None]
+                templateBondedTo[atom] = [templateAtoms.index(template.atoms[j]) for j in atom.bondedTo if template.atoms[j].element is not None]
     else:
         templateAtoms = template.atoms
         templateBondedTo = dict((atom, atom.bondedTo) for atom in template.atoms)
@@ -191,13 +191,13 @@ def matchResidueToTemplate(res, template, bondedToAtom, bint ignoreExternalBonds
     return None
 
 
-def _getAtomMatchCandidates(templateAtoms, bondedTo, matches, candidates, position):
+def _getAtomMatchCandidates(templateAtoms, bondedTo, templateBondedTo, matches, candidates, position):
     """Get a list of template atoms that are potential matches for the next atom."""
     for bonded in bondedTo[position]:
         if bonded < position:
             # This atom is bonded to another one for which we already have a match, so only consider
             # template atoms that *that* one is bonded to.
-            return templateAtoms[matches[bonded]].bondedTo
+            return templateBondedTo[templateAtoms[matches[bonded]]]
     return candidates[position]
 
 
@@ -206,7 +206,7 @@ def _findAtomMatches(templateAtoms, bondedTo, templateBondedTo, matches, hasMatc
     if position == len(matches):
         return True
     cdef int i
-    for i in _getAtomMatchCandidates(templateAtoms, bondedTo, matches, candidates, position):
+    for i in _getAtomMatchCandidates(templateAtoms, bondedTo, templateBondedTo, matches, candidates, position):
         atom = templateAtoms[i]
         if not hasMatch[i] and i in candidates[position]:
             # See if the bonds for this identification are consistent

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

@@ -1048,7 +1048,7 @@ class Modeller(object):
 
         # Identify the template to use for each residue.
 
-        templates = forcefield._matchAllResiduesToTemplates(ForceField._SystemData(self.topology), self.topology, {}, False, True)
+        templates = forcefield._matchAllResiduesToTemplates(ForceField._SystemData(self.topology), self.topology, {}, False, True, False)
 
         # Create the new Topology.