Templates can match whole molecules (#5181)

* Templates can match whole molecules

* addExtraParticles() supports molecule templates

* Documentation on molecule templates

* Bug fix

docs-source/usersguide/application/06_creating_ffs.rst

@@ -1156,6 +1156,28 @@ the :code:`<Include>` tag) or the OpenMM data directory (the one containing
 built in force fields).
 
 
+Molecule Templates
+==================
+
+Although templates are normally matched to individual residues, it is also possible
+to write templates that are intended to be matched to entire molecules. This is a
+less common feature that is useful for force fields that parameterize, for example,
+an entire multi-residue peptide as a single unit.  It is especially useful in
+combination with template generators, which are discussed below.
+
+A force field first tries to assign atom types by matching each residue to a template.
+If it is unable to match some residues this way, it then tries merging all the
+residues in a molecule together and looks for a template that matches the complete
+molecule.
+
+When defining templates, you do not need to distinguish between residue and molecule
+templates.  They are all defined with the :code:`<Residue>` tag.  A molecule template
+can be recognized by the lack of any :code:`<ExternalBond>` tags, since a whole
+molecule by definition never has any external bonds.  When a molecule consists of
+only a single residue, there is of course no difference between a residue template
+and a molecule template.  The distinction only matters for multi-residue molecules.
+
+
 Using Multiple Files
 ********************
 

wrappers/python/openmm/app/data/charmm36_2024.xml

@@ -259286,11 +259286,10 @@ def find_improper(improper_types, atom_classes):
 residue_data = []
-for residue_index, residue in enumerate(topology.residues()):
+for residue in topology.residues():
     # Extract the residue or patch names and atom names.
-    template_data = templateForResidue[residue_index]
     templates = set()
-    if template_data is None:
+    if residue not in templateForResidue:
         # Multi-residue patch; fall back to parsing atom names because OpenMM
         # leaves None in templateForResidue in this case.
         atom_names = []
@@ -259308,6 +259307,7 @@ for residue_index, residue in enumerate(topology.residues()):
             continue
     else:
         # Extract names from template name.
+        template_data = templateForResidue[residue]
         for template_index, template_name in enumerate(template_data.name.split("-")):
             if template_index:
                 # Patch name.

wrappers/python/openmm/app/data/charmm_polar_2023.xml

@@ -30887,9 +30887,9 @@ def find_improper(improper_types, atom_classes):
 residue_data = []
-for residue_index, residue in enumerate(topology.residues()):
+for residue in topology.residues():
     # Extract the residue or patch names and atom names.
-    template_data = templateForResidue[residue_index]
+    template_data = templateForResidue[residue]
     templates = set()
     if template_data is None:
         # Multi-residue patch; fall back to parsing atom names because OpenMM
@@ -31531,9 +31531,9 @@ def extract_atom_name(raw_atom_name, templates):
 residue_data = []
-for residue_index, residue in enumerate(topology.residues()):
+for residue in topology.residues():
     # Extract the residue or patch names and atom names.
-    template_data = templateForResidue[residue_index]
+    template_data = templateForResidue[residue]
     templates = set()
     if template_data is None:
         # Multi-residue patch; fall back to parsing atom names because OpenMM

wrappers/python/openmm/app/forcefield.py

@@ -4,7 +4,7 @@ forcefield.py: Constructs OpenMM System objects based on a Topology and an XML f
 This is part of the OpenMM molecular simulation toolkit.
 See https://openmm.org/development.
 
-Portions copyright (c) 2012-2025 Stanford University and the Authors.
+Portions copyright (c) 2012-2026 Stanford University and the Authors.
 Authors: Peter Eastman, Mark Friedrichs
 Contributors: Evan Pretti
 
@@ -43,6 +43,7 @@ from difflib import SequenceMatcher
 import openmm as mm
 import openmm.unit as unit
 from . import element as elem
+from openmm.app.topology import MergedResidue
 from openmm.app.internal.singleton import Singleton
 from openmm.app.internal import compiled, amoebaforces
 from openmm.app.internal.argtracker import ArgTracker
@@ -641,14 +642,22 @@ class ForceField(object):
 
             # Record which atoms are bonded to each other atom
 
-            for i in range(len(self.bonds)):
-                bond = self.bonds[i]
+            for i, bond in enumerate(self.bonds):
                 self.bondedToAtom[bond.atom1].add(bond.atom2)
                 self.bondedToAtom[bond.atom2].add(bond.atom1)
                 self.atomBonds[bond.atom1].append(i)
                 self.atomBonds[bond.atom2].append(i)
             self.bondedToAtom = [sorted(b) for b in self.bondedToAtom]
 
+            # Record which residues are bonded to each other atom
+
+            bondedResidues = defaultdict(set)
+            for atom1, atom2 in topology.bonds():
+                if atom1.residue != atom2.residue:
+                    bondedResidues[atom1.residue].add(atom2.residue)
+                    bondedResidues[atom2.residue].add(atom1.residue)
+            self.bondedResidues = {x:list(y) for x, y in bondedResidues.items()}
+
         def addConstraint(self, system, atom1, atom2, distance):
             """Add a constraint to the system, avoiding duplicate constraints."""
             key = (min(atom1, atom2), max(atom1, atom2))
@@ -1264,12 +1273,12 @@ class ForceField(object):
         # Find the template matching each residue and assign atom types.
 
         templateForResidue = self._matchAllResiduesToTemplates(data, topology, residueTemplates, ignoreExternalBonds)
-        for res in topology.residues():
+        for res, template in templateForResidue.items():
             if res.name == 'HOH':
                 # Determine whether this should be a rigid water.
 
                 if rigidWater is None:
-                    rigidResidue[res.index] = templateForResidue[res.index].rigidWater
+                    rigidResidue[res.index] = template.rigidWater
                 elif rigidWater:
                     rigidResidue[res.index] = True
 
@@ -1439,7 +1448,7 @@ class ForceField(object):
 
     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()
+        templateForResidue = {}
         unmatchedResidues = []
         for chain in topology.chains():
             for res in chain.residues():
@@ -1457,7 +1466,7 @@ class ForceField(object):
                 else:
                     if recordParameters:
                         data.recordMatchedAtomParameters(res, template, matches)
-                    templateForResidue[res.index] = template
+                    templateForResidue[res] = template
 
         # Try to apply patches to find matches for any unmatched residues.
 
@@ -1465,17 +1474,32 @@ class ForceField(object):
             unmatchedResidues = _applyPatchesToMatchResidues(self, data, unmatchedResidues, templateForResidue, data.bondedToAtom, ignoreExternalBonds, ignoreExtraParticles)
 
         # If we still haven't found a match for a residue, attempt to use residue template generators to create
-        # new templates (and potentially atom types/parameters).
+        # new templates (and potentially atom types/parameters).  Also try to match templates to the whole molecule.
 
-        for res in unmatchedResidues:
+        while len(unmatchedResidues) > 0:
+            res = unmatchedResidues[0]
+            matchedRes = res
+            mol = self._mergeMolecule(data, res)
             # A template might have been generated on an earlier iteration of this loop.
             [template, matches] = self._getResidueTemplateMatches(res, data.bondedToAtom, ignoreExternalBonds=ignoreExternalBonds, ignoreExtraParticles=ignoreExtraParticles)
+            if matches is None and mol is not None:
+                # Try matching to the whole molecule.
+                [template, matches] = self._getResidueTemplateMatches(mol, data.bondedToAtom, ignoreExternalBonds=ignoreExternalBonds, ignoreExtraParticles=ignoreExtraParticles)
+                if matches is not None:
+                    matchedRes = mol
+                    unmatchedResidues = [r for r in unmatchedResidues if r not in mol.residues]
             if matches is None:
                 # Try all generators.
                 for generator in self._templateGenerators:
                     if generator(self, res):
                         # This generator has registered a new residue template that should match.
                         [template, matches] = self._getResidueTemplateMatches(res, data.bondedToAtom, ignoreExternalBonds=ignoreExternalBonds, ignoreExtraParticles=ignoreExtraParticles)
+                        if matches is None and mol is not None:
+                            # Try matching to the whole molecule.
+                            [template, matches] = self._getResidueTemplateMatches(mol, data.bondedToAtom, ignoreExternalBonds=ignoreExternalBonds, ignoreExtraParticles=ignoreExtraParticles)
+                            if matches is not None:
+                                matchedRes = mol
+                                unmatchedResidues = [r for r in unmatchedResidues if r not in mol.residues]
                         if matches is None:
                             # Something went wrong because the generated template does not match the residue signature.
                             raise Exception('The residue handler %s indicated it had correctly parameterized residue %s, but the generated template did not match the residue signature.' % (generator.__class__.__name__, str(res)))
@@ -1485,12 +1509,39 @@ class ForceField(object):
             if matches is None:
                 raise ValueError('No template found for residue %d (%s).  %s  For more information, see https://github.com/openmm/openmm/wiki/Frequently-Asked-Questions#template' % (res.index, res.name, _findMatchErrors(self, res)))
             else:
+                if res in unmatchedResidues:
+                    unmatchedResidues.remove(res)
                 if recordParameters:
-                    data.recordMatchedAtomParameters(res, template, matches)
-                templateForResidue[res.index] = template
+                    data.recordMatchedAtomParameters(matchedRes, template, matches)
+                templateForResidue[matchedRes] = template
         return templateForResidue
 
 
+    def _mergeMolecule(self, data, residue):
+        """Starting from one residue, identify all others it is bonded to and create a MergedResidue for the whole molecule."""
+        if residue not in data.bondedResidues:
+            # This residue is not bonded to any other, so there is nothing to merge.
+            return None
+        visited = set()
+        molecule = set()
+        residueStack = [residue]
+        neighborStack = [0]
+        while len(residueStack) > 0:
+            currentRes = residueStack[-1]
+            bonded = data.bondedResidues[currentRes]
+            molecule.add(currentRes)
+            visited.add(currentRes)
+            while neighborStack[-1] < len(bonded) and bonded[neighborStack[-1]] in visited:
+                neighborStack[-1] +=1
+            if neighborStack[-1] < len(bonded):
+                residueStack.append(bonded[neighborStack[-1]])
+                neighborStack.append(0)
+            else:
+                residueStack.pop()
+                neighborStack.pop()
+        return MergedResidue(list(molecule))
+
+
 def _findBondsForExclusions(data, sys):
     """Create a list of bonds to use when identifying exclusions."""
     bondIndices = []
@@ -1643,7 +1694,7 @@ def _applyPatchesToMatchResidues(forcefield, data, residues, templateForResidue,
             unmatchedResidues.append(res)
         else:
             data.recordMatchedAtomParameters(res, template, matches)
-            templateForResidue[res.index] = template
+            templateForResidue[res] = template
     if len(unmatchedResidues) == 0:
         return []
 

wrappers/python/openmm/app/modeller.py

@@ -4,7 +4,7 @@ modeller.py: Provides tools for editing molecular models
 This is part of the OpenMM molecular simulation toolkit.
 See https://openmm.org/development.
 
-Portions copyright (c) 2012-2025 Stanford University and the Authors.
+Portions copyright (c) 2012-2026 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors: 
 
@@ -47,7 +47,7 @@ import sys
 import xml.etree.ElementTree as etree
 from copy import deepcopy
 from math import ceil, floor, sqrt
-from collections import defaultdict
+from collections import defaultdict, namedtuple
 
 class Modeller(object):
     """Modeller provides tools for editing molecular models, such as adding water or missing hydrogens.
@@ -1156,6 +1156,50 @@ class Modeller(object):
 
         templates = forcefield._matchAllResiduesToTemplates(ForceField._SystemData(self.topology), self.topology, residueTemplates, False, True, False)
 
+        # Identify extra points to add.  Record information about each one we will need when adding it.
+
+        ExtraPoint = namedtuple('ExtraPoint', ['matchedResidue', 'template', 'templateAtom', 'matchingAtoms', 'templateAtomPositions'])
+        extraPoints = defaultdict(list)
+        from openmm.app.topology import MergedResidue
+        for matchedResidue, template in templates.items():
+            if len(template.atoms) != len(list(matchedResidue.atoms())):
+                matches = compiled.matchResidueToTemplate(matchedResidue, template, bondedToAtom, ignoreExternalBonds, True)
+                atomsNoEP = [a for a in matchedResidue.atoms() if a.element is not None]
+                templateAtomsNoEP = [a for a in template.atoms if a.element is not None]
+                matchingAtoms = {}
+                for atom, match in zip(atomsNoEP, matches):
+                    templateAtomName = templateAtomsNoEP[match].name
+                    for templateAtom in template.atoms:
+                        if templateAtom.name == templateAtomName:
+                            matchingAtoms[templateAtom] = atom
+                templateAtomPositions = len(template.atoms)*[None]
+                for index, atom in enumerate(template.atoms):
+                    if atom in matchingAtoms:
+                        templateAtomPositions[index] = self.positions[matchingAtoms[atom].index].value_in_unit(nanometer)
+                for index, atom in enumerate(template.atoms):
+                    if atom.element is None:
+                        # This is an extra particle.  Figure out what residue it belongs to.
+                        residue = None
+                        if not isinstance(matchedResidue, MergedResidue):
+                            residue = matchedResidue
+                        else:
+                            if atom.type in drudeTypeMap:
+                                # It's a Drude particle.
+                                for atom2 in template.atoms:
+                                    if atom2.type in drudeTypeMap[atom.type]:
+                                        residue = matchingAtoms[atom2].residue
+                                        break
+                            if residue is None:
+                                for site in template.virtualSites:
+                                    if site.index == index:
+                                        # It's a virtual site.
+                                        residue = matchingAtoms[template.atoms[site.atoms[0]]].residue
+                                        break
+                            if residue is None:
+                                # Neither a virtual site nor a Drude particle???  Just guess at the residue.
+                                residue = matchedResidue.residues[0]
+                        extraPoints[residue].append(ExtraPoint(matchedResidue, template, atom, matchingAtoms, templateAtomPositions))
+
         # Create the new Topology.
 
         newTopology = Topology()
@@ -1170,8 +1214,7 @@ class Modeller(object):
                 newResidue = newTopology.addResidue(residue.name, newChain, residue.id, residue.insertionCode)
                 if residue in residueTemplates:
                     newResidueTemplates[newResidue] = residueTemplates[residue]
-                template = templates[residue.index]
-                if len(template.atoms) == len(list(residue.atoms())):
+                if residue not in extraPoints:
                     # Just copy the residue over.
 
                     for atom in residue.atoms():
@@ -1179,18 +1222,6 @@ class Modeller(object):
                         newAtoms[atom] = newAtom
                         newPositions.append(deepcopy(self.positions[atom.index]))
                 else:
-                    # Record the corresponding atoms.
-
-                    matches = compiled.matchResidueToTemplate(residue, template, bondedToAtom, ignoreExternalBonds, True)
-                    atomsNoEP = [a for a in residue.atoms() if a.element is not None]
-                    templateAtomsNoEP = [a for a in template.atoms if a.element is not None]
-                    matchingAtoms = {}
-                    for atom, match in zip(atomsNoEP, matches):
-                        templateAtomName = templateAtomsNoEP[match].name
-                        for templateAtom in template.atoms:
-                            if templateAtom.name == templateAtomName:
-                                matchingAtoms[templateAtom] = atom
-
                     # Add the regular atoms.
 
                     for atom in residue.atoms():
@@ -1200,61 +1231,61 @@ class Modeller(object):
 
                     # Add the extra points.
 
-                    templateAtomPositions = len(template.atoms)*[None]
-                    for index, atom in enumerate(template.atoms):
-                        if atom in matchingAtoms:
-                            templateAtomPositions[index] = self.positions[matchingAtoms[atom].index].value_in_unit(nanometer)
                     newExtraPoints = {}
-                    for index, atom in enumerate(template.atoms):
-                        if atom.element is None:
-                            newExtraPoints[atom] = newTopology.addAtom(atom.name, None, newResidue)
-                            position = None
-                            for site in template.virtualSites:
-                                if site.index == index:
-                                    # This is a virtual site.  Compute its position by the correct rule.
-
-                                    try:
-                                        if site.type == 'average2':
-                                            position = site.weights[0]*templateAtomPositions[site.atoms[0]] + site.weights[1]*templateAtomPositions[site.atoms[1]]
-                                        elif site.type == 'average3':
-                                            position = site.weights[0]*templateAtomPositions[site.atoms[0]] + site.weights[1]*templateAtomPositions[site.atoms[1]] + site.weights[2]*templateAtomPositions[site.atoms[2]]
-                                        elif site.type == 'outOfPlane':
-                                            v1 = templateAtomPositions[site.atoms[1]] - templateAtomPositions[site.atoms[0]]
-                                            v2 = templateAtomPositions[site.atoms[2]] - templateAtomPositions[site.atoms[0]]
-                                            cross = Vec3(v1[1]*v2[2]-v1[2]*v2[1], v1[2]*v2[0]-v1[0]*v2[2], v1[0]*v2[1]-v1[1]*v2[0])
-                                            position = templateAtomPositions[site.atoms[0]] + site.weights[0]*v1 + site.weights[1]*v2 + site.weights[2]*cross
-                                        elif site.type == 'localCoords':
-                                            origin = unit.sum([templateAtomPositions[atom]*weight for atom, weight in zip(site.atoms, site.originWeights)])
-                                            xdir = unit.sum([templateAtomPositions[atom]*weight for atom, weight in zip(site.atoms, site.xWeights)])
-                                            ydir = unit.sum([templateAtomPositions[atom]*weight for atom, weight in zip(site.atoms, site.yWeights)])
-                                            zdir = Vec3(xdir[1]*ydir[2]-xdir[2]*ydir[1], xdir[2]*ydir[0]-xdir[0]*ydir[2], xdir[0]*ydir[1]-xdir[1]*ydir[0])
-                                            xdir /= norm(xdir)
-                                            zdir /= norm(zdir)
-                                            ydir = Vec3(zdir[1]*xdir[2]-zdir[2]*xdir[1], zdir[2]*xdir[0]-zdir[0]*xdir[2], zdir[0]*xdir[1]-zdir[1]*xdir[0])
-                                            position = origin + xdir*site.localPos[0] + ydir*site.localPos[1] + zdir*site.localPos[2]
-                                    except:
-                                        # This can happen if the virtual site depends on another virtual site whose position
-                                        # hasn't been set yet.  Ignore the error.  We'll put it at a random position (see below),
-                                        # which will get replaced with the correct position at the start of the simulation.
-
-                                        pass
-                            if position is None and atom.type in drudeTypeMap:
-                                # This is a Drude particle.  Put it on top of its parent atom.
-
-                                for atom2, pos in zip(template.atoms, templateAtomPositions):
-                                    if atom2.type in drudeTypeMap[atom.type]:
-                                        position = deepcopy(pos)
-                            if position is None:
-                                # We couldn't figure out the correct position.  Put it at a random position near the center of the residue,
-                                # and we'll try to fix it later based on bonds.
-
-                                knownPositions = [x for x in templateAtomPositions if x is not None]
-                                position = Vec3(random.gauss(0, 1), random.gauss(0, 1), random.gauss(0, 1))+(unit.sum(knownPositions)/len(knownPositions))
-                                missingPositions.add(len(newPositions))
-                            newPositions.append(position*nanometer)
+                    for extraPoint in extraPoints[residue]:
+                        atom = extraPoint.templateAtom
+                        template = extraPoint.template
+                        templateAtomPositions = extraPoint.templateAtomPositions
+                        newExtraPoints[atom] = newTopology.addAtom(atom.name, None, newResidue)
+                        position = None
+                        for site in extraPoint.template.virtualSites:
+                            if template.atoms[site.index] == atom:
+                                # This is a virtual site.  Compute its position by the correct rule.
+
+                                try:
+                                    if site.type == 'average2':
+                                        position = site.weights[0]*templateAtomPositions[site.atoms[0]] + site.weights[1]*templateAtomPositions[site.atoms[1]]
+                                    elif site.type == 'average3':
+                                        position = site.weights[0]*templateAtomPositions[site.atoms[0]] + site.weights[1]*templateAtomPositions[site.atoms[1]] + site.weights[2]*templateAtomPositions[site.atoms[2]]
+                                    elif site.type == 'outOfPlane':
+                                        v1 = templateAtomPositions[site.atoms[1]] - templateAtomPositions[site.atoms[0]]
+                                        v2 = templateAtomPositions[site.atoms[2]] - templateAtomPositions[site.atoms[0]]
+                                        cross = Vec3(v1[1]*v2[2]-v1[2]*v2[1], v1[2]*v2[0]-v1[0]*v2[2], v1[0]*v2[1]-v1[1]*v2[0])
+                                        position = templateAtomPositions[site.atoms[0]] + site.weights[0]*v1 + site.weights[1]*v2 + site.weights[2]*cross
+                                    elif site.type == 'localCoords':
+                                        origin = unit.sum([templateAtomPositions[atom]*weight for atom, weight in zip(site.atoms, site.originWeights)])
+                                        xdir = unit.sum([templateAtomPositions[atom]*weight for atom, weight in zip(site.atoms, site.xWeights)])
+                                        ydir = unit.sum([templateAtomPositions[atom]*weight for atom, weight in zip(site.atoms, site.yWeights)])
+                                        zdir = Vec3(xdir[1]*ydir[2]-xdir[2]*ydir[1], xdir[2]*ydir[0]-xdir[0]*ydir[2], xdir[0]*ydir[1]-xdir[1]*ydir[0])
+                                        xdir /= norm(xdir)
+                                        zdir /= norm(zdir)
+                                        ydir = Vec3(zdir[1]*xdir[2]-zdir[2]*xdir[1], zdir[2]*xdir[0]-zdir[0]*xdir[2], zdir[0]*xdir[1]-zdir[1]*xdir[0])
+                                        position = origin + xdir*site.localPos[0] + ydir*site.localPos[1] + zdir*site.localPos[2]
+                                except:
+                                    # This can happen if the virtual site depends on another virtual site whose position
+                                    # hasn't been set yet.  Ignore the error.  We'll put it at a random position (see below),
+                                    # which will get replaced with the correct position at the start of the simulation.
+
+                                    pass
+                        if position is None and atom.type in drudeTypeMap:
+                            # This is a Drude particle.  Put it on top of its parent atom.
+
+                            for atom2, pos in zip(template.atoms, templateAtomPositions):
+                                if atom2.type in drudeTypeMap[atom.type]:
+                                    position = deepcopy(pos)
+                        if position is None:
+                            # We couldn't figure out the correct position.  Put it at a random position near the center of the residue,
+                            # and we'll try to fix it later based on bonds.
+
+                            knownPositions = [x for x in templateAtomPositions if x is not None]
+                            position = Vec3(random.gauss(0, 1), random.gauss(0, 1), random.gauss(0, 1))+(unit.sum(knownPositions)/len(knownPositions))
+                            missingPositions.add(len(newPositions))
+                        newPositions.append(position*nanometer)
 
                     # Add bonds involving the extra points.
 
+                    template = extraPoints[residue][0].template
+                    matchingAtoms = extraPoints[residue][0].matchingAtoms
                     for atom1, atom2 in template.bonds:
                         atom1 = template.atoms[atom1]
                         atom2 = template.atoms[atom2]

wrappers/python/openmm/app/topology.py

@@ -4,7 +4,7 @@ topology.py: Used for storing topological information about a system.
 This is part of the OpenMM molecular simulation toolkit.
 See https://openmm.org/development.
 
-Portions copyright (c) 2012-2025 Stanford University and the Authors.
+Portions copyright (c) 2012-2026 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors:
 
@@ -515,3 +515,32 @@ class Bond(namedtuple('Bond', ['atom1', 'atom2'])):
             s = "%s, order=%d" % (s, self.order)
         s += ")"
         return s
+
+class MergedResidue(Residue):
+    """A MergedResidue is a pseudo-residue created by merging multiple Residues into a single object.  It is never
+    contained in a Topology, but is sometimes created for use in parameterization and modelling.  A MergedResidue
+    differs from an ordinary Residue in the follow ways.
+
+    1. It contains a list of the Residue objects that were merged to create it.
+    2. It does not own its Atoms.  They are the same objects found in the original Residues.
+    3. Its chain field refers to the Chain containing the first Residue that was merged.  Because a MergedResidue might
+       span multiple chains, it may not be contained entirely in that Chain.
+    4. Its index field is set to -1, since it has no index within the Topology.
+    """
+    def __init__(self, residues: list[Residue]):
+        """Create a MergedResidue by combining a list of Residues."""
+        if len(residues) == 0:
+            raise ValueError('A MergedResidue must contain at least one Residue')
+        for res in residues[1:]:
+            if res.chain.topology != residues[0].chain.topology:
+                raise ValueError('All Residues in a MergedResidue must belong to the same Topology')
+        ## The list of Residues that were merged.
+        self.residues = residues
+        self.name = 'Merged'
+        self.index = -1
+        self.chain = residues[0].chain
+        self.id = None
+        self.insertionCode = ''
+        self._atoms = []
+        for res in residues:
+            self._atoms += res._atoms

wrappers/python/tests/TestForceField.py

@@ -1870,6 +1870,111 @@ self.scriptExecuted = True
         self.assertTrue(abs(energy1 - energy_amber) < energy_tolerance)
         self.assertTrue(abs(energy1 - energy2) < energy_tolerance)
 
+    def testWholeMolecule(self):
+        """Test matching a template to a whole molecule."""
+        xml = """
+<ForceField>
+  <AtomTypes>
+    <Type class="C" element="C" mass="12.01" name="C" />
+    <Type class="CT" element="C" mass="12.01" name="CT" />
+    <Type class="CX" element="C" mass="12.01" name="CX"/>
+    <Type class="H" element="H" mass="1.008" name="H"/>
+    <Type class="HC" element="H" mass="1.008" name="HC" />
+    <Type class="H1" element="H" mass="1.008" name="H1" />
+    <Type class="N" element="N" mass="14.01" name="N" />
+    <Type class="O" element="O" mass="16.0" name="O" />
+  </AtomTypes>
+  <Residues>
+    <Residue name="Alanine-Dipeptide">
+      <Atom charge="0.1123" name="ACE-H1" type="HC" />
+      <Atom charge="-0.3662" name="ACE-CH3" type="CT" />
+      <Atom charge="0.1123" name="ACE-H2" type="HC" />
+      <Atom charge="0.1123" name="ACE-H3" type="HC" />
+      <Atom charge="0.5972" name="ACE-C" type="C" />
+      <Atom charge="-0.5679" name="ACE-O" type="O" />
+      <Bond atomName1="ACE-H1" atomName2="ACE-CH3" />
+      <Bond atomName1="ACE-CH3" atomName2="ACE-H2" />
+      <Bond atomName1="ACE-CH3" atomName2="ACE-H3" />
+      <Bond atomName1="ACE-CH3" atomName2="ACE-C" />
+      <Bond atomName1="ACE-C" atomName2="ACE-O" />
+      <Atom charge="-0.4157" name="ALA-N" type="N" />
+      <Atom charge="0.2719" name="ALA-H" type="H" />
+      <Atom charge="0.0337" name="ALA-CA" type="CX" />
+      <Atom charge="0.0823" name="ALA-HA" type="H1" />
+      <Atom charge="-0.1825" name="ALA-CB" type="CT" />
+      <Atom charge="0.0603" name="ALA-HB1" type="HC" />
+      <Atom charge="0.0603" name="ALA-HB2" type="HC" />
+      <Atom charge="0.0603" name="ALA-HB3" type="HC" />
+      <Atom charge="0.5973" name="ALA-C" type="C" />
+      <Atom charge="-0.5679" name="ALA-O" type="O" />
+      <Bond atomName1="ALA-N" atomName2="ALA-H" />
+      <Bond atomName1="ALA-N" atomName2="ALA-CA" />
+      <Bond atomName1="ALA-CA" atomName2="ALA-HA" />
+      <Bond atomName1="ALA-CA" atomName2="ALA-CB" />
+      <Bond atomName1="ALA-CA" atomName2="ALA-C" />
+      <Bond atomName1="ALA-CB" atomName2="ALA-HB1" />
+      <Bond atomName1="ALA-CB" atomName2="ALA-HB2" />
+      <Bond atomName1="ALA-CB" atomName2="ALA-HB3" />
+      <Bond atomName1="ALA-C" atomName2="ALA-O" />
+      <Atom charge="-0.4157" name="NME-N" type="N" />
+      <Atom charge="0.2719" name="NME-H" type="H" />
+      <Atom charge="-0.149" name="NME-C" type="CT" />
+      <Atom charge="0.0976" name="NME-H1" type="H1" />
+      <Atom charge="0.0976" name="NME-H2" type="H1" />
+      <Atom charge="0.0976" name="NME-H3" type="H1" />
+      <Bond atomName1="NME-N" atomName2="NME-H" />
+      <Bond atomName1="NME-N" atomName2="NME-C" />
+      <Bond atomName1="NME-C" atomName2="NME-H1" />
+      <Bond atomName1="NME-C" atomName2="NME-H2" />
+      <Bond atomName1="NME-C" atomName2="NME-H3" />
+      <Bond atomName1="ACE-C" atomName2="ALA-N" />
+      <Bond atomName1="ALA-C" atomName2="NME-N" />
+    </Residue>
+    <!-- A template that matches just the ACE with different parameters -->
+    <Residue name="ACE">
+      <Atom charge="0.1123" name="ACE-H1" type="HC" />
+      <Atom charge="-0.3662" name="ACE-CH3" type="CT" />
+      <Atom charge="-10.0" name="ACE-H2" type="HC" />
+      <Atom charge="0.1123" name="ACE-H3" type="HC" />
+      <Atom charge="10.0" name="ACE-C" type="C" />
+      <Atom charge="-0.5679" name="ACE-O" type="O" />
+      <Bond atomName1="ACE-H1" atomName2="ACE-CH3" />
+      <Bond atomName1="ACE-CH3" atomName2="ACE-H2" />
+      <Bond atomName1="ACE-CH3" atomName2="ACE-H3" />
+      <Bond atomName1="ACE-CH3" atomName2="ACE-C" />
+      <Bond atomName1="ACE-C" atomName2="ACE-O" />
+      <ExternalBond atomName="ACE-C" />
+    </Residue>
+  </Residues>
+  <NonbondedForce coulomb14scale="0.8333333333333334" lj14scale="0.5">
+    <UseAttributeFromResidue name="charge"/>
+    <Atom epsilon="0.359824" sigma="0.3399669508423535" type="C"/>
+    <Atom epsilon="0.4577296" sigma="0.3399669508423535" type="CT"/>
+    <Atom epsilon="0.4577296" sigma="0.3399669508423535" type="CX"/>
+    <Atom epsilon="0.06568879999999999" sigma="0.2649532787749369" type="HC"/>
+    <Atom epsilon="0.06568879999999999" sigma="0.10690784617684071" type="H"/>
+    <Atom epsilon="0.06568879999999999" sigma="0.2471353044121301" type="H1"/>
+    <Atom epsilon="0.7112800000000001" sigma="0.3249998523775958" type="N"/>
+    <Atom epsilon="0.87864" sigma="0.2959921901149463" type="O"/>
+  </NonbondedForce>
+</ForceField>"""
+        pdb = PDBFile('systems/alanine-dipeptide-implicit.pdb')
+        ff = ForceField(StringIO(xml))
+        system = ff.createSystem(pdb.topology)
+        nonbonded = next(f for f in system.getForces() if isinstance(f, NonbondedForce))
+
+        def checkAtom(resName, atomName, expected):
+            for atom in pdb.topology.atoms():
+                if atom.name == atomName and atom.residue.name == resName:
+                    params = nonbonded.getParticleParameters(atom.index)
+                    self.assertEqual(expected, params)
+                    return
+            raise ValueError(f'{resName} {atomName} not found')
+
+        checkAtom('ACE', 'C', [0.5972*elementary_charge, 0.3399669508423535*nanometers, 0.359824*kilojoules_per_mole])
+        checkAtom('ACE', 'H2', [0.1123*elementary_charge, 0.2649532787749369*nanometers, 0.06568879999999999*kilojoules_per_mole])
+        checkAtom('ALA', 'CA', [0.0337*elementary_charge, 0.3399669508423535*nanometers, 0.4577296*kilojoules_per_mole])
+        checkAtom('NME', 'N', [-0.4157*elementary_charge, 0.3249998523775958*nanometers, 0.7112800000000001*kilojoules_per_mole])
 
 class AmoebaTestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method with the AMOEBA forcefield."""

wrappers/python/tests/TestModeller.py

@@ -1230,6 +1230,46 @@ class TestModeller(unittest.TestCase):
                 self.assertTrue(dist > (expectedDist-0.01)*nanometers and dist < (expectedDist+0.01)*nanometers)
 
 
+    def testWholeMoleculeExtraParticles(self):
+        """Test adding extra particles based on matching a template to a whole molecule."""
+        pdb = PDBFile(StringIO("""
+ATOM      1  C   AAA     1A    0.000     0.000   0.000  1.00  0.00           C
+ATOM      2  C   BBB     2A    1.000     0.000   0.000  1.00  0.00           C
+ATOM      3  O   BBB     2A    2.000     0.000   0.000  1.00  0.00           O
+CONECT    1    2
+CONECT    2    1    3
+CONECT    3    2
+END"""))
+        ff = ForceField(StringIO("""
+<ForceField>
+  <AtomTypes>
+    <Type class="C" element="C" mass="12.01" name="C" />
+    <Type class="O" element="O" mass="16.0" name="O" />
+    <Type class="V" mass="0.0" name="V" />
+  </AtomTypes>
+  <Residues>
+    <Residue name="Molecule">
+      <Atom name="C1" type="C" />
+      <Atom name="C2" type="C" />
+      <Atom name="O" type="O" />
+      <Atom name="V" type="V" />
+      <Bond atomName1="C1" atomName2="C2"/>
+      <Bond atomName1="C2" atomName2="O"/>
+      <VirtualSite type="average2" siteName="V" atomName1="C2" atomName2="O" weight1="0.25" weight2="0.75"/>
+    </Residue>
+  </Residues>
+</ForceField>"""))
+        modeller = Modeller(pdb.topology, pdb.positions)
+        modeller.addExtraParticles(ff)
+        residues = list(modeller.topology.residues())
+        self.assertEqual(2, len(residues))
+        self.assertEqual(1, len(residues[0]))
+        self.assertEqual(3, len(residues[1]))
+        self.assertVecAlmostEqual(Vec3(1.75, 0.0, 0.0), modeller.positions[3].value_in_unit(angstrom))
+        system = ff.createSystem(modeller.topology)
+        self.assertTrue(system.isVirtualSite(3))
+
+
     def test_addMembrane(self):
         """Test adding a membrane to a realistic system."""