Merge pull request #426 from peastman/vsite

Additional changes to support LocalCoordinatesSite

serialization/src/SystemProxy.cpp

@@ -67,6 +67,18 @@ void SystemProxy::serialize(const void* object, SerializationNode& node) const {
                 const OutOfPlaneSite& site = dynamic_cast<const OutOfPlaneSite&>(system.getVirtualSite(i));
                 particle.createChildNode("OutOfPlaneSite").setIntProperty("p1", site.getParticle(0)).setIntProperty("p2", site.getParticle(1)).setIntProperty("p3", site.getParticle(2)).setDoubleProperty("w12", site.getWeight12()).setDoubleProperty("w13", site.getWeight13()).setDoubleProperty("wc", site.getWeightCross());
             }
+            else if (typeid(system.getVirtualSite(i)) == typeid(LocalCoordinatesSite)) {
+                const LocalCoordinatesSite& site = dynamic_cast<const LocalCoordinatesSite&>(system.getVirtualSite(i));
+                Vec3 wo = site.getOriginWeights();
+                Vec3 wx = site.getXWeights();
+                Vec3 wy = site.getYWeights();
+                Vec3 p = site.getLocalPosition();
+                particle.createChildNode("LocalCoordinatesSite").setIntProperty("p1", site.getParticle(0)).setIntProperty("p2", site.getParticle(1)).setIntProperty("p3", site.getParticle(2)).
+                        setDoubleProperty("wo1", wo[0]).setDoubleProperty("wo2", wo[1]).setDoubleProperty("wo3", wo[2]).
+                        setDoubleProperty("wx1", wx[0]).setDoubleProperty("wx2", wx[1]).setDoubleProperty("wx3", wx[2]).
+                        setDoubleProperty("wy1", wy[0]).setDoubleProperty("wy2", wy[1]).setDoubleProperty("wy3", wy[2]).
+                        setDoubleProperty("pos1", p[0]).setDoubleProperty("pos2", p[1]).setDoubleProperty("pos3", p[2]);
+            }
         }
     }
     SerializationNode& constraints = node.createChildNode("Constraints");
@@ -105,6 +117,13 @@ void* SystemProxy::deserialize(const SerializationNode& node) const {
                     system->setVirtualSite(i, new ThreeParticleAverageSite(vsite.getIntProperty("p1"), vsite.getIntProperty("p2"), vsite.getIntProperty("p3"), vsite.getDoubleProperty("w1"), vsite.getDoubleProperty("w2"), vsite.getDoubleProperty("w3")));
                 else if (vsite.getName() == "OutOfPlaneSite")
                     system->setVirtualSite(i, new OutOfPlaneSite(vsite.getIntProperty("p1"), vsite.getIntProperty("p2"), vsite.getIntProperty("p3"), vsite.getDoubleProperty("w12"), vsite.getDoubleProperty("w13"), vsite.getDoubleProperty("wc")));
+                else if (vsite.getName() == "LocalCoordinatesSite") {
+                    Vec3 wo(vsite.getDoubleProperty("wo1"), vsite.getDoubleProperty("wo2"), vsite.getDoubleProperty("wo3"));
+                    Vec3 wx(vsite.getDoubleProperty("wx1"), vsite.getDoubleProperty("wx2"), vsite.getDoubleProperty("wx3"));
+                    Vec3 wy(vsite.getDoubleProperty("wy1"), vsite.getDoubleProperty("wy2"), vsite.getDoubleProperty("wy3"));
+                    Vec3 p(vsite.getDoubleProperty("pos1"), vsite.getDoubleProperty("pos2"), vsite.getDoubleProperty("pos3"));
+                    system->setVirtualSite(i, new LocalCoordinatesSite(vsite.getIntProperty("p1"), vsite.getIntProperty("p2"), vsite.getIntProperty("p3"), wo, wx, wy, p));
+                }
             }
         }
         const SerializationNode& constraints = node.getChildNode("Constraints");

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

@@ -267,9 +267,14 @@ class ForceField(object):
             elif self.type == 'outOfPlane':
                 self.atoms = [int(attrib['atom1']), int(attrib['atom2']), int(attrib['atom3'])]
                 self.weights = [float(attrib['weight12']), float(attrib['weight13']), float(attrib['weightCross'])]
+            elif self.type == 'localCoords':
+                self.atoms = [int(attrib['atom1']), int(attrib['atom2']), int(attrib['atom3'])]
+                self.originWeights = [float(attrib['wo1']), float(attrib['wo2']), float(attrib['wo3'])]
+                self.xWeights = [float(attrib['wx1']), float(attrib['wx2']), float(attrib['wx3'])]
+                self.yWeights = [float(attrib['wy1']), float(attrib['wy2']), float(attrib['wy3'])]
+                self.localPos = [float(attrib['p1']), float(attrib['p2']), float(attrib['p3'])]
             else:
                 raise ValueError('Unknown virtual site type: %s' % self.type)
-            self.atoms = [x-self.index for x in self.atoms]
 
     def createSystem(self, topology, nonbondedMethod=NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
                      constraints=None, rigidWater=True, removeCMMotion=True, hydrogenMass=None, **args):
@@ -326,11 +331,12 @@ class ForceField(object):
                             break
                 if matches is None:
                     raise ValueError('No template found for residue %d (%s).  %s' % (res.index+1, res.name, _findMatchErrors(self, res)))
+                matchAtoms = dict(zip(matches, res.atoms()))
                 for atom, match in zip(res.atoms(), matches):
                     data.atomType[atom] = template.atoms[match].type
                     for site in template.virtualSites:
                         if match == site.index:
-                            data.virtualSites[atom] = site
+                            data.virtualSites[atom] = (site, [matchAtoms[i].index for i in site.atoms])
 
         # Create the System and add atoms
 
@@ -445,14 +451,20 @@ class ForceField(object):
         # Add virtual sites
 
         for atom in data.virtualSites:
-            site = data.virtualSites[atom]
+            (site, atoms) = data.virtualSites[atom]
             index = atom.index
             if site.type == 'average2':
-                sys.setVirtualSite(index, mm.TwoParticleAverageSite(index+site.atoms[0], index+site.atoms[1], site.weights[0], site.weights[1]))
+                sys.setVirtualSite(index, mm.TwoParticleAverageSite(atoms[0], atoms[1], site.weights[0], site.weights[1]))
             elif site.type == 'average3':
-                sys.setVirtualSite(index, mm.ThreeParticleAverageSite(index+site.atoms[0], index+site.atoms[1], index+site.atoms[2], site.weights[0], site.weights[1], site.weights[2]))
+                sys.setVirtualSite(index, mm.ThreeParticleAverageSite(atoms[0], atoms[1], atoms[2], site.weights[0], site.weights[1], site.weights[2]))
             elif site.type == 'outOfPlane':
-                sys.setVirtualSite(index, mm.OutOfPlaneSite(index+site.atoms[0], index+site.atoms[1], index+site.atoms[2], site.weights[0], site.weights[1], site.weights[2]))
+                sys.setVirtualSite(index, mm.OutOfPlaneSite(atoms[0], atoms[1], atoms[2], site.weights[0], site.weights[1], site.weights[2]))
+            elif site.type == 'localCoords':
+                sys.setVirtualSite(index, mm.LocalCoordinatesSite(atoms[0], atoms[1], atoms[2],
+                                                                  mm.Vec3(site.originWeights[0], site.originWeights[1], site.originWeights[2]),
+                                                                  mm.Vec3(site.xWeights[0], site.xWeights[1], site.xWeights[2]),
+                                                                  mm.Vec3(site.yWeights[0], site.yWeights[1], site.yWeights[2]),
+                                                                  mm.Vec3(site.localPos[0], site.localPos[1], site.localPos[2])))
 
         # Add forces to the System
 
@@ -563,7 +575,7 @@ def _findAtomMatches(atoms, template, bondedTo, externalBonds, matches, hasMatch
     name = atoms[position].name
     for i in range(len(atoms)):
         atom = template.atoms[i]
-        if (atom.element == elem or (atom.element is None and atom.name == name)) and not hasMatch[i] and len(atom.bondedTo) == len(bondedTo[position]) and atom.externalBonds == externalBonds[position]:
+        if ((atom.element is not None and atom.element == elem) or (atom.element is None and atom.name == name)) and not hasMatch[i] and len(atom.bondedTo) == len(bondedTo[position]) and atom.externalBonds == externalBonds[position]:
             # See if the bonds for this identification are consistent
 
             allBondsMatch = all((bonded > position or matches[bonded] in atom.bondedTo for bonded in bondedTo[position]))