Patches can define virtual sites

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

@@ -358,6 +358,9 @@ class ForceField(object):
                     for bond in patch.findall('RemoveExternalBond'):
                         atom = ForceField._PatchAtomData(bond.attrib['atomName'])
                         patchData.deletedExternalBonds.append(atom)
+                    atomIndices = dict((atom.name, i) for i, atom in enumerate(patchData.addedAtoms[atomDescription.residue]+patchData.changedAtoms[atomDescription.residue]))
+                    for site in patch.findall('VirtualSite'):
+                        patchData.virtualSites[atomDescription.residue].append(ForceField._VirtualSiteData(site, atomIndices))
                     for residue in patch.findall('ApplyToResidue'):
                         name = residue.attrib['name']
                         if ':' in name:
@@ -740,6 +743,7 @@ class ForceField(object):
             self.addedExternalBonds = []
             self.deletedExternalBonds = []
             self.allAtomNames = set()
+            self.virtualSites = [[] for i in range(numResidues)]
 
         def createPatchedTemplates(self, templates):
             """Apply this patch to a set of templates, creating new modified ones."""
@@ -801,6 +805,16 @@ class ForceField(object):
                         newTemplate.addExternalBondByName(atom2.name)
                 for atom in self.addedExternalBonds:
                     newTemplate.addExternalBondByName(atom.name)
+
+                # Add new virtual sites.
+
+                indexMap = dict((i, newAtomIndex[atom.name]) for i, atom in enumerate(self.addedAtoms[index]+self.changedAtoms[index]))
+                for site in self.virtualSites[index]:
+                    newSite = deepcopy(site)
+                    newSite.index = indexMap[site.index]
+                    newSite.atoms = [indexMap[i] for i in site.atoms]
+                    newTemplate.virtualSites = [site for site in newTemplate.virtualSites if site.index != newSite.index]
+                    newTemplate.virtualSites.append(newSite)
             return newTemplates
 
     class _PatchAtomData(object):

wrappers/python/tests/TestPatches.py

@@ -32,6 +32,7 @@ class TestPatches(unittest.TestCase):
     <AddExternalBond atomName="A"/>
     <RemoveExternalBond atomName="C"/>
     <ApplyToResidue name="RES"/>
+    <VirtualSite type="average2" siteName="A" atomName1="B" atomName2="B" weight1="0.8" weight2="0.2"/>
   </Patch>
  </Patches>
 </ForceField>"""
@@ -47,6 +48,7 @@ class TestPatches(unittest.TestCase):
         self.assertEqual(1, len(patch.deletedBonds))
         self.assertEqual(1, len(patch.addedExternalBonds))
         self.assertEqual(1, len(patch.deletedExternalBonds))
+        self.assertEqual(1, len(patch.virtualSites))
         self.assertEqual(1, len(ff._templatePatches))
         self.assertEqual(1, len(ff._templatePatches['RES']))
         self.assertEqual('A', patch.addedAtoms[0][0].name)
@@ -67,6 +69,9 @@ class TestPatches(unittest.TestCase):
         self.assertEqual(0, patch.addedExternalBonds[0].residue)
         self.assertEqual('C', patch.deletedExternalBonds[0].name)
         self.assertEqual(0, patch.deletedExternalBonds[0].residue)
+        self.assertEqual(0, patch.virtualSites[0][0].index)
+        self.assertEqual([1, 1], patch.virtualSites[0][0].atoms)
+        self.assertEqual([0.8, 0.2], patch.virtualSites[0][0].weights)
         patch = list(ff._templatePatches['RES'])[0]
         self.assertEqual('Test', patch[0])
         self.assertEqual(0, patch[1])