API and reference implementation of LocalCoordinatesSite depending on arbitrary particles

openmmapi/include/openmm/VirtualSite.h

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2012-2014 Stanford University and the Authors.      *
+ * Portions copyright (c) 2012-2017 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -168,21 +168,21 @@ private:
 };
 
 /**
- * This is a VirtualSite that uses the locations of three other particles to compute a local
+ * This is a VirtualSite that uses the locations of several other particles to compute a local
  * coordinate system, then places the virtual site at a fixed location in that coordinate
  * system.  The origin of the coordinate system and the directions of its x and y axes
- * are each specified as a weighted sum of the locations of the three particles:
+ * are each specified as a weighted sum of the locations of the other particles:
  * 
- * origin = w<sup>o</sup><sub>1</sub>r<sub>1</sub> +  w<sup>o</sup><sub>2</sub>r<sub>2</sub> +  w<sup>o</sup><sub>3</sub>r<sub>3</sub>
+ * origin = w<sup>o</sup><sub>1</sub>r<sub>1</sub> +  w<sup>o</sup><sub>2</sub>r<sub>2</sub> +  w<sup>o</sup><sub>3</sub>r<sub>3</sub> + ...
  * 
- * xdir = w<sup>x</sup><sub>1</sub>r<sub>1</sub> +  w<sup>x</sup><sub>2</sub>r<sub>2</sub> +  w<sup>x</sup><sub>3</sub>r<sub>3</sub>
+ * xdir = w<sup>x</sup><sub>1</sub>r<sub>1</sub> +  w<sup>x</sup><sub>2</sub>r<sub>2</sub> +  w<sup>x</sup><sub>3</sub>r<sub>3</sub> + ...
  * 
- * ydir = w<sup>y</sup><sub>1</sub>r<sub>1</sub> +  w<sup>y</sup><sub>2</sub>r<sub>2</sub> +  w<sup>y</sup><sub>3</sub>r<sub>3</sub>
+ * ydir = w<sup>y</sup><sub>1</sub>r<sub>1</sub> +  w<sup>y</sup><sub>2</sub>r<sub>2</sub> +  w<sup>y</sup><sub>3</sub>r<sub>3</sub> + ...
  * 
- * For the origin, the three weights must add to one.  For example if
+ * For the origin, the weights must add to one.  For example if
  * (w<sup>o</sup><sub>1</sub>, w<sup>o</sup><sub>2</sub>, w<sup>o</sup><sub>3</sub>) = (1.0, 0.0, 0.0),
  * the origin of the local coordinate system is at the location of particle 1.  For xdir and ydir,
- * the weights must add to zero.  For excample, if
+ * the weights must add to zero.  For example, if
  * (w<sup>x</sup><sub>1</sub>, w<sup>x</sup><sub>2</sub>, w<sup>x</sup><sub>3</sub>) = (-1.0, 0.5, 0.5),
  * the x axis points from particle 1 toward the midpoint between particles 2 and 3.
  * 
@@ -197,6 +197,17 @@ public:
     /**
      * Create a new LocalCoordinatesSite virtual site.
      * 
+     * @param particles      the indices of the particle the site depends on
+     * @param originWeights  the weight factors for the particles when computing the origin location
+     * @param xWeights       the weight factors for the particles when computing xdir
+     * @param yWeights       the weight factors for the particles when computing ydir
+     * @param localPosition  the position of the virtual site in the local coordinate system
+     */
+    LocalCoordinatesSite(const std::vector<int>& particles, const std::vector<double>& originWeights, const std::vector<double>& xWeights, const std::vector<double>& yWeights, const Vec3& localPosition);
+    /**
+     * Create a new LocalCoordinatesSite virtual site.  This constructor assumes the site depends on
+     * exactly three other particles.
+     * 
      * @param particle1      the index of the first particle
      * @param particle2      the index of the second particle
      * @param particle3      the index of the third particle
@@ -207,23 +218,42 @@ public:
      */
     LocalCoordinatesSite(int particle1, int particle2, int particle3, const Vec3& originWeights, const Vec3& xWeights, const Vec3& yWeights, const Vec3& localPosition);
     /**
-     * Get the weight factors for the three particles when computing the origin location.
+     * Get the weight factors for the particles when computing the origin location.
+     */
+    void getOriginWeights(std::vector<double>& weights) const;
+    /**
+     * Get the weight factors for the particles when computing the origin location.
+     * This version assumes the site depends on exactly three other particles, and
+     * throws an exception if that is not the case.
+     */
+    Vec3 getOriginWeights() const;
+    /**
+     * Get the weight factors for the particles when computing xdir.
+     */
+    void getXWeights(std::vector<double>& weights) const;
+    /**
+     * Get the weight factors for the particles when computing xdir.
+     * This version assumes the site depends on exactly three other particles, and
+     * throws an exception if that is not the case.
      */
-    const Vec3& getOriginWeights() const;
+    Vec3 getXWeights() const;
     /**
-     * Get the weight factors for the three particles when computing xdir.
+     * Get the weight factors for the particles when computing ydir.
      */
-    const Vec3& getXWeights() const;
+    void getYWeights(std::vector<double>& weights) const;
     /**
-     * Get the weight factors for the three particles when computing ydir.
+     * Get the weight factors for the particles when computing ydir.
+     * This version assumes the site depends on exactly three other particles, and
+     * throws an exception if that is not the case.
      */
-    const Vec3& getYWeights() const;
+    Vec3 getYWeights() const;
     /**
      * Get the position of the virtual site in the local coordinate system.
      */
     const Vec3& getLocalPosition() const;
 private:
-    Vec3 originWeights, xWeights, yWeights, localPosition;
+    std::vector<double> originWeights, xWeights, yWeights;
+    Vec3 localPosition;
 };
 
 } // namespace OpenMM

openmmapi/src/VirtualSite.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2012-2014 Stanford University and the Authors.      *
+ * Portions copyright (c) 2012-2017 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -105,8 +105,30 @@ double OutOfPlaneSite::getWeightCross() const {
     return weightCross;
 }
 
-LocalCoordinatesSite::LocalCoordinatesSite(int particle1, int particle2, int particle3, const Vec3& originWeights, const Vec3& xWeights, const Vec3& yWeights, const Vec3& localPosition) :
+LocalCoordinatesSite::LocalCoordinatesSite(const vector<int>& particles, const vector<double>& originWeights, const vector<double>& xWeights, const vector<double>& yWeights, const Vec3& localPosition) :
         originWeights(originWeights), xWeights(xWeights), yWeights(yWeights), localPosition(localPosition) {
+    int numParticles = particles.size();
+    if (numParticles < 2)
+        throw OpenMMException("LocalCoordinatesSite: Must depend on at least two other particles");
+    if (originWeights.size() != numParticles || xWeights.size() != numParticles || yWeights.size() != numParticles)
+        throw OpenMMException("LocalCoordinatesSite: Number of weights does not match number of particles");
+    double originSum = 0, xSum = 0, ySum = 0;
+    for (int i = 0; i < numParticles; i++) {
+        originSum += originWeights[i];
+        xSum += xWeights[i];
+        ySum += yWeights[i];
+    }
+    if (fabs(originSum-1.0) > 1e-6)
+        throw OpenMMException("LocalCoordinatesSite: Weights for computing origin must add to 1");
+    if (fabs(xSum) > 1e-6)
+        throw OpenMMException("LocalCoordinatesSite: Weights for computing x axis must add to 0");
+    if (fabs(ySum) > 1e-6)
+        throw OpenMMException("LocalCoordinatesSite: Weights for computing y axis must add to 0");
+    setParticles(particles);
+}
+
+LocalCoordinatesSite::LocalCoordinatesSite(int particle1, int particle2, int particle3, const Vec3& originWeights, const Vec3& xWeights, const Vec3& yWeights, const Vec3& localPosition) :
+        localPosition(localPosition) {
     if (fabs(originWeights[0]+originWeights[1]+originWeights[2]-1.0) > 1e-6)
         throw OpenMMException("LocalCoordinatesSite: Weights for computing origin must add to 1");
     if (fabs(xWeights[0]+xWeights[1]+xWeights[2]) > 1e-6)
@@ -118,18 +140,45 @@ LocalCoordinatesSite::LocalCoordinatesSite(int particle1, int particle2, int par
     particles[1] = particle2;
     particles[2] = particle3;
     setParticles(particles);
+    this->originWeights.push_back(originWeights[0]);
+    this->originWeights.push_back(originWeights[1]);
+    this->originWeights.push_back(originWeights[2]);
+    this->xWeights.push_back(xWeights[0]);
+    this->xWeights.push_back(xWeights[1]);
+    this->xWeights.push_back(xWeights[2]);
+    this->yWeights.push_back(yWeights[0]);
+    this->yWeights.push_back(yWeights[1]);
+    this->yWeights.push_back(yWeights[2]);
+}
+
+void LocalCoordinatesSite::getOriginWeights(vector<double>& weights) const {
+    weights = originWeights;
+}
+
+Vec3 LocalCoordinatesSite::getOriginWeights() const {
+    if (originWeights.size() != 3)
+        throw OpenMMException("LocalCoordinatesSite: This version of getOriginWeights() requires the site to depend on three particles");
+    return Vec3(originWeights[0], originWeights[1], originWeights[2]);
+}
+
+void LocalCoordinatesSite::getXWeights(vector<double>& weights) const {
+    weights = xWeights;
 }
 
-const Vec3& LocalCoordinatesSite::getOriginWeights() const {
-    return originWeights;
+Vec3 LocalCoordinatesSite::getXWeights() const {
+    if (xWeights.size() != 3)
+        throw OpenMMException("LocalCoordinatesSite: This version of getXWeights() requires the site to depend on three particles");
+    return Vec3(xWeights[0], xWeights[1], xWeights[2]);
 }
 
-const Vec3& LocalCoordinatesSite::getXWeights() const {
-    return xWeights;
+void LocalCoordinatesSite::getYWeights(vector<double>& weights) const {
+    weights = yWeights;
 }
 
-const Vec3& LocalCoordinatesSite::getYWeights() const {
-    return yWeights;
+Vec3 LocalCoordinatesSite::getYWeights() const {
+    if (yWeights.size() != 3)
+        throw OpenMMException("LocalCoordinatesSite: This version of getYWeights() requires the site to depend on three particles");
+    return Vec3(yWeights[0], yWeights[1], yWeights[2]);
 }
 
 const Vec3& LocalCoordinatesSite::getLocalPosition() const {

platforms/reference/src/SimTKReference/ReferenceVirtualSites.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2012-2014 Stanford University and the Authors.      *
+ * Portions copyright (c) 2012-2017 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -70,22 +70,30 @@ void ReferenceVirtualSites::computePositions(const OpenMM::System& system, vecto
                 // A local coordinates site.
                 
                 const LocalCoordinatesSite& site = dynamic_cast<const LocalCoordinatesSite&>(system.getVirtualSite(i));
-                int p1 = site.getParticle(0), p2 = site.getParticle(1), p3 = site.getParticle(2);
-                Vec3 originWeights = site.getOriginWeights();
-                Vec3 xWeights = site.getXWeights();
-                Vec3 yWeights = site.getYWeights();
+                int numParticles = site.getNumParticles();
+                vector<double> originWeights, xWeights, yWeights;
+                site.getOriginWeights(originWeights);
+                site.getXWeights(xWeights);
+                site.getYWeights(yWeights);
+                Vec3 origin, xdir, ydir;
+                for (int j = 0; j < numParticles; j++) {
+                    Vec3 pos = atomCoordinates[site.getParticle(j)];
+                    origin += pos*originWeights[j];
+                    xdir += pos*xWeights[j];
+                    ydir += pos*yWeights[j];
+                }
                 Vec3 localPosition = site.getLocalPosition();
-                Vec3 origin = atomCoordinates[p1]*originWeights[0] + atomCoordinates[p2]*originWeights[1] + atomCoordinates[p3]*originWeights[2];
-                Vec3 xdir = atomCoordinates[p1]*xWeights[0] + atomCoordinates[p2]*xWeights[1] + atomCoordinates[p3]*xWeights[2];
-                Vec3 ydir = atomCoordinates[p1]*yWeights[0] + atomCoordinates[p2]*yWeights[1] + atomCoordinates[p3]*yWeights[2];
                 Vec3 zdir = xdir.cross(ydir);
-                xdir /= sqrt(xdir.dot(xdir));
-                zdir /= sqrt(zdir.dot(zdir));
+                double normXdir = sqrt(xdir.dot(xdir));
+                double normZdir = sqrt(zdir.dot(zdir));
+                if (normXdir > 0.0)
+                    xdir /= normXdir;
+                if (normZdir > 0.0)
+                    zdir /= normZdir;
                 ydir = zdir.cross(xdir);
                 atomCoordinates[i] = origin + xdir*localPosition[0] + ydir*localPosition[1] + zdir*localPosition[2];
             }
         }
-
 }
 
 void ReferenceVirtualSites::distributeForces(const OpenMM::System& system, const vector<OpenMM::Vec3>& atomCoordinates, vector<OpenMM::Vec3>& forces) {
@@ -133,71 +141,53 @@ void ReferenceVirtualSites::distributeForces(const OpenMM::System& system, const
                 // A local coordinates site.
                 
                 const LocalCoordinatesSite& site = dynamic_cast<const LocalCoordinatesSite&>(system.getVirtualSite(i));
-                int p1 = site.getParticle(0), p2 = site.getParticle(1), p3 = site.getParticle(2);
-                Vec3 originWeights = site.getOriginWeights();
-                Vec3 wx = site.getXWeights();
-                Vec3 wy = site.getYWeights();
+                int numParticles = site.getNumParticles();
+                vector<double> originWeights, wx, wy;
+                site.getOriginWeights(originWeights);
+                site.getXWeights(wx);
+                site.getYWeights(wy);
+                Vec3 xdir, ydir;
+                for (int j = 0; j < numParticles; j++) {
+                    Vec3 pos = atomCoordinates[site.getParticle(j)];
+                    xdir += pos*wx[j];
+                    ydir += pos*wy[j];
+                }
                 Vec3 localPosition = site.getLocalPosition();
-                Vec3 xdir = atomCoordinates[p1]*wx[0] + atomCoordinates[p2]*wx[1] + atomCoordinates[p3]*wx[2];
-                Vec3 ydir = atomCoordinates[p1]*wy[0] + atomCoordinates[p2]*wy[1] + atomCoordinates[p3]*wy[2];
                 Vec3 zdir = xdir.cross(ydir);
-                double invNormXdir = 1.0/sqrt(xdir.dot(xdir));
-                double invNormZdir = 1.0/sqrt(zdir.dot(zdir));
+                double normXdir = sqrt(xdir.dot(xdir));
+                double normZdir = sqrt(zdir.dot(zdir));
+                double invNormXdir = (normXdir > 0.0 ? 1.0/normXdir : 0.0);
+                double invNormZdir = (normZdir > 0.0 ? 1.0/normZdir : 0.0);
                 Vec3 dx = xdir*invNormXdir;
                 Vec3 dz = zdir*invNormZdir;
                 Vec3 dy = dz.cross(dx);
                 
                 // The derivatives for this case are very complicated.  They were computed with SymPy then simplified by hand.
                 
-                double t11 = (wx[0]*ydir[0]-wy[0]*xdir[0])*invNormZdir;
-                double t12 = (wx[0]*ydir[1]-wy[0]*xdir[1])*invNormZdir;
-                double t13 = (wx[0]*ydir[2]-wy[0]*xdir[2])*invNormZdir;
-                double t21 = (wx[1]*ydir[0]-wy[1]*xdir[0])*invNormZdir;
-                double t22 = (wx[1]*ydir[1]-wy[1]*xdir[1])*invNormZdir;
-                double t23 = (wx[1]*ydir[2]-wy[1]*xdir[2])*invNormZdir;
-                double t31 = (wx[2]*ydir[0]-wy[2]*xdir[0])*invNormZdir;
-                double t32 = (wx[2]*ydir[1]-wy[2]*xdir[1])*invNormZdir;
-                double t33 = (wx[2]*ydir[2]-wy[2]*xdir[2])*invNormZdir;
-                double sx1 = t13*dz[1]-t12*dz[2];
-                double sy1 = t11*dz[2]-t13*dz[0];
-                double sz1 = t12*dz[0]-t11*dz[1];
-                double sx2 = t23*dz[1]-t22*dz[2];
-                double sy2 = t21*dz[2]-t23*dz[0];
-                double sz2 = t22*dz[0]-t21*dz[1];
-                double sx3 = t33*dz[1]-t32*dz[2];
-                double sy3 = t31*dz[2]-t33*dz[0];
-                double sz3 = t32*dz[0]-t31*dz[1];
-                Vec3 wxScaled = wx*invNormXdir;
+                vector<double> wxScaled(numParticles);
+                for (int j = 0; j < numParticles; j++)
+                    wxScaled[j] = wx[j]*invNormXdir;
                 Vec3 fp1 = localPosition*f[0];
                 Vec3 fp2 = localPosition*f[1];
                 Vec3 fp3 = localPosition*f[2];
-                forces[p1][0] += fp1[0]*wxScaled[0]*(1-dx[0]*dx[0]) + fp1[2]*(dz[0]*sx1    ) + fp1[1]*((-dx[0]*dy[0]      )*wxScaled[0] + dy[0]*sx1 - dx[1]*t12 - dx[2]*t13) + f[0]*originWeights[0];
-                forces[p1][1] += fp1[0]*wxScaled[0]*( -dx[0]*dx[1]) + fp1[2]*(dz[0]*sy1+t13) + fp1[1]*((-dx[1]*dy[0]-dz[2])*wxScaled[0] + dy[0]*sy1 + dx[1]*t11);
-                forces[p1][2] += fp1[0]*wxScaled[0]*( -dx[0]*dx[2]) + fp1[2]*(dz[0]*sz1-t12) + fp1[1]*((-dx[2]*dy[0]+dz[1])*wxScaled[0] + dy[0]*sz1 + dx[2]*t11);
-                forces[p2][0] += fp1[0]*wxScaled[1]*(1-dx[0]*dx[0]) + fp1[2]*(dz[0]*sx2    ) + fp1[1]*((-dx[0]*dy[0]      )*wxScaled[1] + dy[0]*sx2 - dx[1]*t22 - dx[2]*t23) + f[0]*originWeights[1];
-                forces[p2][1] += fp1[0]*wxScaled[1]*( -dx[0]*dx[1]) + fp1[2]*(dz[0]*sy2+t23) + fp1[1]*((-dx[1]*dy[0]-dz[2])*wxScaled[1] + dy[0]*sy2 + dx[1]*t21);
-                forces[p2][2] += fp1[0]*wxScaled[1]*( -dx[0]*dx[2]) + fp1[2]*(dz[0]*sz2-t22) + fp1[1]*((-dx[2]*dy[0]+dz[1])*wxScaled[1] + dy[0]*sz2 + dx[2]*t21);
-                forces[p3][0] += fp1[0]*wxScaled[2]*(1-dx[0]*dx[0]) + fp1[2]*(dz[0]*sx3    ) + fp1[1]*((-dx[0]*dy[0]      )*wxScaled[2] + dy[0]*sx3 - dx[1]*t32 - dx[2]*t33) + f[0]*originWeights[2];
-                forces[p3][1] += fp1[0]*wxScaled[2]*( -dx[0]*dx[1]) + fp1[2]*(dz[0]*sy3+t33) + fp1[1]*((-dx[1]*dy[0]-dz[2])*wxScaled[2] + dy[0]*sy3 + dx[1]*t31);
-                forces[p3][2] += fp1[0]*wxScaled[2]*( -dx[0]*dx[2]) + fp1[2]*(dz[0]*sz3-t32) + fp1[1]*((-dx[2]*dy[0]+dz[1])*wxScaled[2] + dy[0]*sz3 + dx[2]*t31);
-                forces[p1][0] += fp2[0]*wxScaled[0]*( -dx[1]*dx[0]) + fp2[2]*(dz[1]*sx1-t13) - fp2[1]*(( dx[0]*dy[1]-dz[2])*wxScaled[0] - dy[1]*sx1 - dx[0]*t12);
-                forces[p1][1] += fp2[0]*wxScaled[0]*(1-dx[1]*dx[1]) + fp2[2]*(dz[1]*sy1    ) - fp2[1]*(( dx[1]*dy[1]      )*wxScaled[0] - dy[1]*sy1 + dx[0]*t11 + dx[2]*t13) + f[1]*originWeights[0];
-                forces[p1][2] += fp2[0]*wxScaled[0]*( -dx[1]*dx[2]) + fp2[2]*(dz[1]*sz1+t11) - fp2[1]*(( dx[2]*dy[1]+dz[0])*wxScaled[0] - dy[1]*sz1 - dx[2]*t12);
-                forces[p2][0] += fp2[0]*wxScaled[1]*( -dx[1]*dx[0]) + fp2[2]*(dz[1]*sx2-t23) - fp2[1]*(( dx[0]*dy[1]-dz[2])*wxScaled[1] - dy[1]*sx2 - dx[0]*t22);
-                forces[p2][1] += fp2[0]*wxScaled[1]*(1-dx[1]*dx[1]) + fp2[2]*(dz[1]*sy2    ) - fp2[1]*(( dx[1]*dy[1]      )*wxScaled[1] - dy[1]*sy2 + dx[0]*t21 + dx[2]*t23) + f[1]*originWeights[1];
-                forces[p2][2] += fp2[0]*wxScaled[1]*( -dx[1]*dx[2]) + fp2[2]*(dz[1]*sz2+t21) - fp2[1]*(( dx[2]*dy[1]+dz[0])*wxScaled[1] - dy[1]*sz2 - dx[2]*t22);
-                forces[p3][0] += fp2[0]*wxScaled[2]*( -dx[1]*dx[0]) + fp2[2]*(dz[1]*sx3-t33) - fp2[1]*(( dx[0]*dy[1]-dz[2])*wxScaled[2] - dy[1]*sx3 - dx[0]*t32);
-                forces[p3][1] += fp2[0]*wxScaled[2]*(1-dx[1]*dx[1]) + fp2[2]*(dz[1]*sy3    ) - fp2[1]*(( dx[1]*dy[1]      )*wxScaled[2] - dy[1]*sy3 + dx[0]*t31 + dx[2]*t33) + f[1]*originWeights[2];
-                forces[p3][2] += fp2[0]*wxScaled[2]*( -dx[1]*dx[2]) + fp2[2]*(dz[1]*sz3+t31) - fp2[1]*(( dx[2]*dy[1]+dz[0])*wxScaled[2] - dy[1]*sz3 - dx[2]*t32);
-                forces[p1][0] += fp3[0]*wxScaled[0]*( -dx[2]*dx[0]) + fp3[2]*(dz[2]*sx1+t12) + fp3[1]*((-dx[0]*dy[2]-dz[1])*wxScaled[0] + dy[2]*sx1 + dx[0]*t13);
-                forces[p1][1] += fp3[0]*wxScaled[0]*( -dx[2]*dx[1]) + fp3[2]*(dz[2]*sy1-t11) + fp3[1]*((-dx[1]*dy[2]+dz[0])*wxScaled[0] + dy[2]*sy1 + dx[1]*t13);
-                forces[p1][2] += fp3[0]*wxScaled[0]*(1-dx[2]*dx[2]) + fp3[2]*(dz[2]*sz1    ) + fp3[1]*((-dx[2]*dy[2]      )*wxScaled[0] + dy[2]*sz1 - dx[0]*t11 - dx[1]*t12) + f[2]*originWeights[0];
-                forces[p2][0] += fp3[0]*wxScaled[1]*( -dx[2]*dx[0]) + fp3[2]*(dz[2]*sx2+t22) + fp3[1]*((-dx[0]*dy[2]-dz[1])*wxScaled[1] + dy[2]*sx2 + dx[0]*t23);
-                forces[p2][1] += fp3[0]*wxScaled[1]*( -dx[2]*dx[1]) + fp3[2]*(dz[2]*sy2-t21) + fp3[1]*((-dx[1]*dy[2]+dz[0])*wxScaled[1] + dy[2]*sy2 + dx[1]*t23);
-                forces[p2][2] += fp3[0]*wxScaled[1]*(1-dx[2]*dx[2]) + fp3[2]*(dz[2]*sz2    ) + fp3[1]*((-dx[2]*dy[2]      )*wxScaled[1] + dy[2]*sz2 - dx[0]*t21 - dx[1]*t22) + f[2]*originWeights[1];
-                forces[p3][0] += fp3[0]*wxScaled[2]*( -dx[2]*dx[0]) + fp3[2]*(dz[2]*sx3+t32) + fp3[1]*((-dx[0]*dy[2]-dz[1])*wxScaled[2] + dy[2]*sx3 + dx[0]*t33);
-                forces[p3][1] += fp3[0]*wxScaled[2]*( -dx[2]*dx[1]) + fp3[2]*(dz[2]*sy3-t31) + fp3[1]*((-dx[1]*dy[2]+dz[0])*wxScaled[2] + dy[2]*sy3 + dx[1]*t33);
-                forces[p3][2] += fp3[0]*wxScaled[2]*(1-dx[2]*dx[2]) + fp3[2]*(dz[2]*sz3    ) + fp3[1]*((-dx[2]*dy[2]      )*wxScaled[2] + dy[2]*sz3 - dx[0]*t31 - dx[1]*t32) + f[2]*originWeights[2];
+                for (int j = 0; j < numParticles; j++) {
+                    double t1 = (wx[j]*ydir[0]-wy[j]*xdir[0])*invNormZdir;
+                    double t2 = (wx[j]*ydir[1]-wy[j]*xdir[1])*invNormZdir;
+                    double t3 = (wx[j]*ydir[2]-wy[j]*xdir[2])*invNormZdir;
+                    double sx = t3*dz[1]-t2*dz[2];
+                    double sy = t1*dz[2]-t3*dz[0];
+                    double sz = t2*dz[0]-t1*dz[1];
+                    int p = site.getParticle(j);
+                    forces[p][0] += fp1[0]*wxScaled[j]*(1-dx[0]*dx[0]) + fp1[2]*(dz[0]*sx   ) + fp1[1]*((-dx[0]*dy[0]      )*wxScaled[j] + dy[0]*sx - dx[1]*t2 - dx[2]*t3) + f[0]*originWeights[j];
+                    forces[p][1] += fp1[0]*wxScaled[j]*( -dx[0]*dx[1]) + fp1[2]*(dz[0]*sy+t3) + fp1[1]*((-dx[1]*dy[0]-dz[2])*wxScaled[j] + dy[0]*sy + dx[1]*t1);
+                    forces[p][2] += fp1[0]*wxScaled[j]*( -dx[0]*dx[2]) + fp1[2]*(dz[0]*sz-t2) + fp1[1]*((-dx[2]*dy[0]+dz[1])*wxScaled[j] + dy[0]*sz + dx[2]*t1);
+                    forces[p][0] += fp2[0]*wxScaled[j]*( -dx[1]*dx[0]) + fp2[2]*(dz[1]*sx-t3) - fp2[1]*(( dx[0]*dy[1]-dz[2])*wxScaled[j] - dy[1]*sx - dx[0]*t2);
+                    forces[p][1] += fp2[0]*wxScaled[j]*(1-dx[1]*dx[1]) + fp2[2]*(dz[1]*sy   ) - fp2[1]*(( dx[1]*dy[1]      )*wxScaled[j] - dy[1]*sy + dx[0]*t1 + dx[2]*t3) + f[1]*originWeights[j];
+                    forces[p][2] += fp2[0]*wxScaled[j]*( -dx[1]*dx[2]) + fp2[2]*(dz[1]*sz+t1) - fp2[1]*(( dx[2]*dy[1]+dz[0])*wxScaled[j] - dy[1]*sz - dx[2]*t2);
+                    forces[p][0] += fp3[0]*wxScaled[j]*( -dx[2]*dx[0]) + fp3[2]*(dz[2]*sx+t2) + fp3[1]*((-dx[0]*dy[2]-dz[1])*wxScaled[j] + dy[2]*sx + dx[0]*t3);
+                    forces[p][1] += fp3[0]*wxScaled[j]*( -dx[2]*dx[1]) + fp3[2]*(dz[2]*sy-t1) + fp3[1]*((-dx[1]*dy[2]+dz[0])*wxScaled[j] + dy[2]*sy + dx[1]*t3);
+                    forces[p][2] += fp3[0]*wxScaled[j]*(1-dx[2]*dx[2]) + fp3[2]*(dz[2]*sz   ) + fp3[1]*((-dx[2]*dy[2]      )*wxScaled[j] + dy[2]*sz - dx[0]*t1 - dx[1]*t2) + f[2]*originWeights[j];
+                }
            }
         }
 }

serialization/src/SystemProxy.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010-2015 Stanford University and the Authors.      *
+ * Portions copyright (c) 2010-2017 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -71,15 +71,23 @@ void SystemProxy::serialize(const void* object, SerializationNode& node) const {
             }
             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();
+                int numParticles = site.getNumParticles();
+                vector<double> wo, wx, wy;
+                site.getOriginWeights(wo);
+                site.getXWeights(wx);
+                site.getYWeights(wy);
                 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& siteNode = particle.createChildNode("LocalCoordinatesSite");
+                siteNode.setDoubleProperty("pos1", p[0]).setDoubleProperty("pos2", p[1]).setDoubleProperty("pos3", p[2]);
+                for (int j = 0; j < numParticles; j++) {
+                    stringstream ss;
+                    ss << (j+1);
+                    string index = ss.str();
+                    siteNode.setIntProperty("p"+index, site.getParticle(j));
+                    siteNode.setDoubleProperty("wo"+index, wo[j]);
+                    siteNode.setDoubleProperty("wx"+index, wx[j]);
+                    siteNode.setDoubleProperty("wy"+index, wy[j]);
+                }
             }
         }
     }
@@ -120,11 +128,21 @@ void* SystemProxy::deserialize(const SerializationNode& node) const {
                 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"));
+                    vector<int> particles;
+                    vector<double> wo, wx, wy;
+                    for (int j = 0; ; j++) {
+                        stringstream ss;
+                        ss << (j+1);
+                        string index = ss.str();
+                        if (!vsite.hasProperty("p"+index))
+                            break;
+                        particles.push_back(vsite.getIntProperty("p"+index));
+                        wo.push_back(vsite.getDoubleProperty("wo"+index));
+                        wx.push_back(vsite.getDoubleProperty("wx"+index));
+                        wy.push_back(vsite.getDoubleProperty("wy"+index));
+                    }
                     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));
+                    system->setVirtualSite(i, new LocalCoordinatesSite(particles, wo, wx, wy, p));
                 }
             }
         }

serialization/tests/TestSerializeSystem.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2010 Stanford University and the Authors.           *
+ * Portions copyright (c) 2010-2017 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -82,6 +82,23 @@ void compareSystems(System& system, System& system2) {
     ASSERT_EQUAL(0.1, site7.getWeight12());
     ASSERT_EQUAL(0.2, site7.getWeight13());
     ASSERT_EQUAL(0.5, site7.getWeightCross());
+    const LocalCoordinatesSite& site8 = dynamic_cast<const LocalCoordinatesSite&>(system2.getVirtualSite(8));
+    ASSERT_EQUAL(4, site8.getNumParticles());
+    ASSERT_EQUAL(4, site8.getParticle(0));
+    ASSERT_EQUAL(3, site8.getParticle(1));
+    ASSERT_EQUAL(2, site8.getParticle(2));
+    ASSERT_EQUAL(1, site8.getParticle(3));
+    ASSERT_EQUAL(Vec3(-0.5, 1.0, 1.5), site8.getLocalPosition());
+    vector<double> wo, wx, wy;
+    site8.getOriginWeights(wo);
+    site8.getXWeights(wx);
+    site8.getYWeights(wy);
+    vector<double> woExpected = {0.1, 0.2, 0.3, 0.4};
+    vector<double> wxExpected = {-1.0, 0.4, 0.4, 0.2};
+    vector<double> wyExpected = {0.3, 0.7, 0.0, -1.0};
+    ASSERT_EQUAL_CONTAINERS(woExpected, wo);
+    ASSERT_EQUAL_CONTAINERS(wxExpected, wx);
+    ASSERT_EQUAL_CONTAINERS(wyExpected, wy);
     ASSERT_EQUAL(system.getNumForces(), system2.getNumForces());
     for (int i = 0; i < system.getNumForces(); i++)
         ASSERT(typeid(system.getForce(i)) == typeid(system2.getForce(i)))
@@ -93,7 +110,7 @@ void testSerialization() {
     System system;
     for (int i = 0; i < 5; i++)
         system.addParticle(0.1*i+1);
-    for (int i = 0; i < 4; i++)
+    for (int i = 0; i < 5; i++)
         system.addParticle(0.0);
     system.addConstraint(0, 1, 3.0);
     system.addConstraint(1, 2, 2.5);
@@ -102,6 +119,7 @@ void testSerialization() {
     system.setVirtualSite(5, new TwoParticleAverageSite(0, 1, 0.3, 0.7));
     system.setVirtualSite(6, new ThreeParticleAverageSite(2, 4, 3, 0.5, 0.2, 0.3));
     system.setVirtualSite(7, new OutOfPlaneSite(0, 3, 1, 0.1, 0.2, 0.5));
+    system.setVirtualSite(8, new LocalCoordinatesSite({4, 3, 2, 1}, {0.1, 0.2, 0.3, 0.4}, {-1.0, 0.4, 0.4, 0.2}, {0.3, 0.7, 0.0, -1.0}, Vec3(-0.5, 1.0, 1.5)));
     system.addForce(new HarmonicBondForce());
 
     // Serialize and then deserialize it, then make sure the systems are identical.

tests/TestVirtualSites.h

@@ -206,30 +206,51 @@ void testOutOfPlane() {
     }
 }
 
+Vec3 computeWeightedPosition(const vector<Vec3>& positions, const vector<double>& weights) {
+    Vec3 sum;
+    for (int i = 0; i < weights.size(); i++)
+        sum += positions[i]*weights[i];
+    return sum;
+}
+
 /**
  * Test a LocalCoordinatesSite virtual site.
  */
-void testLocalCoordinates() {
-    const Vec3 originWeights(0.2, 0.3, 0.5);
-    const Vec3 xWeights(-1.0, 0.5, 0.5);
-    const Vec3 yWeights(0.0, -1.0, 1.0);
+void testLocalCoordinates(int numSiteParticles) {
+    vector<int> particles;
+    vector<double> originWeights, xWeights, yWeights;
+    if (numSiteParticles == 2) {
+        particles = {0, 1};
+        originWeights = {0.4, 0.6};
+        xWeights = {-1.0, 1.0};
+        yWeights = {1.0, -1.0};
+    }
+    else if (numSiteParticles == 3) {
+        particles = {0, 1, 2};
+        originWeights = {0.2, 0.3, 0.5};
+        xWeights = {-1.0, 0.5, 0.5};
+        yWeights = {0.0, -1.0, 1.0};
+    }
+    else if (numSiteParticles == 4) {
+        particles = {0, 1, 2, 3};
+        originWeights = {0.2, 0.3, 0.1, 0.4};
+        xWeights = {-1.0, 0.3, 0.3, 0.4};
+        yWeights = {0.5, 0.5, -0.5, -0.5};
+    }
     const Vec3 localPosition(0.4, 0.3, 0.2);
     System system;
-    system.addParticle(1.0);
-    system.addParticle(1.0);
-    system.addParticle(1.0);
+    for (int i = 0; i < numSiteParticles; i++)
+        system.addParticle(1.0);
     system.addParticle(0.0);
-    system.setVirtualSite(3, new LocalCoordinatesSite(0, 1, 2, originWeights, xWeights, yWeights, localPosition));
+    system.setVirtualSite(numSiteParticles, new LocalCoordinatesSite(particles, originWeights, xWeights, yWeights, localPosition));
     CustomExternalForce* forceField = new CustomExternalForce("2*x^2+3*y^2+4*z^2");
     system.addForce(forceField);
     vector<double> params;
-    forceField->addParticle(0, params);
-    forceField->addParticle(1, params);
-    forceField->addParticle(2, params);
-    forceField->addParticle(3, params);
+    for (int i = 0; i < numSiteParticles+1; i++)
+        forceField->addParticle(0, params);
     LangevinIntegrator integrator(300.0, 0.1, 0.002);
     Context context(system, integrator, platform);
-    vector<Vec3> positions(4), positions2(4), positions3(4);
+    vector<Vec3> positions(numSiteParticles+1), positions2(numSiteParticles+1), positions3(numSiteParticles+1);
     OpenMM_SFMT::SFMT sfmt;
     init_gen_rand(0, sfmt);
     for (int i = 0; i < 100; i++) {
@@ -237,12 +258,12 @@ void testLocalCoordinates() {
         
         Vec3 xdir, ydir, zdir;
         do {
-            for (int j = 0; j < 3; j++)
+            for (int j = 0; j < numSiteParticles; j++)
                 positions[j] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
-            xdir = positions[0]*xWeights[0] + positions[1]*xWeights[1] + positions[2]*xWeights[2];
-            ydir = positions[0]*yWeights[0] + positions[1]*yWeights[1] + positions[2]*yWeights[2];
+            xdir = computeWeightedPosition(positions, xWeights);
+            ydir = computeWeightedPosition(positions, yWeights);;
             zdir = xdir.cross(ydir);
-            if (sqrt(xdir.dot(xdir)) > 0.1 && sqrt(ydir.dot(ydir)) > 0.1 && sqrt(zdir.dot(zdir)) > 0.1)
+            if (sqrt(xdir.dot(xdir)) > 0.1 && (numSiteParticles == 2 || (sqrt(ydir.dot(ydir)) > 0.1 && sqrt(zdir.dot(zdir)) > 0.1)))
                 break; // These positions give a reasonable coordinate system.
         } while (true);
         context.setPositions(positions);
@@ -252,23 +273,23 @@ void testLocalCoordinates() {
         
         State state = context.getState(State::Positions | State::Forces);
         const vector<Vec3>& pos = state.getPositions();
-        Vec3 origin = pos[0]*originWeights[0] + pos[1]*originWeights[1] + pos[2]*originWeights[2];
+        Vec3 origin = computeWeightedPosition(pos, originWeights);;
         xdir /= sqrt(xdir.dot(xdir));
         zdir /= sqrt(zdir.dot(zdir));
         ydir = zdir.cross(xdir);
-        ASSERT_EQUAL_VEC(origin+xdir*localPosition[0]+ydir*localPosition[1]+zdir*localPosition[2], pos[3], 1e-5);
+        ASSERT_EQUAL_VEC(origin+xdir*localPosition[0]+ydir*localPosition[1]+zdir*localPosition[2], pos[numSiteParticles], 1e-5);
 
         // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
 
         double norm = 0.0;
-        for (int i = 0; i < 3; ++i) {
+        for (int i = 0; i < numSiteParticles; ++i) {
             Vec3 f = state.getForces()[i];
             norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
         }
         norm = std::sqrt(norm);
         const double delta = 1e-2;
         double step = 0.5*delta/norm;
-        for (int i = 0; i < 3; ++i) {
+        for (int i = 0; i < numSiteParticles; ++i) {
             Vec3 p = positions[i];
             Vec3 f = state.getForces()[i];
             positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
@@ -469,7 +490,9 @@ int main(int argc, char* argv[]) {
         testTwoParticleAverage();
         testThreeParticleAverage();
         testOutOfPlane();
-        testLocalCoordinates();
+        testLocalCoordinates(2);
+        testLocalCoordinates(3);
+        testLocalCoordinates(4);
         testConservationLaws();
         testOverlappingSites();
         runPlatformTests();

wrappers/generateWrappers.py

@@ -68,7 +68,15 @@ class WrapperGenerator:
     
     def __init__(self, inputDirname, output):
         self.skipClasses = ['OpenMM::Vec3', 'OpenMM::XmlSerializer', 'OpenMM::Kernel', 'OpenMM::KernelImpl', 'OpenMM::KernelFactory', 'OpenMM::ContextImpl', 'OpenMM::SerializationNode', 'OpenMM::SerializationProxy']
-        self.skipMethods = ['OpenMM::Context::getState', 'OpenMM::Platform::loadPluginsFromDirectory', 'OpenMM::Platform::getPluginLoadFailures', 'OpenMM::Context::createCheckpoint', 'OpenMM::Context::loadCheckpoint', 'OpenMM::Context::getMolecules']
+        self.skipMethods = ['State OpenMM::Context::getState',
+                            'void OpenMM::Context::createCheckpoint',
+                            'void OpenMM::Context::loadCheckpoint',
+                            'const std::vector<std::vector<int> >& OpenMM::Context::getMolecules',
+                            'static std::vector<std::string> OpenMM::Platform::getPluginLoadFailures',
+                            'static std::vector<std::string> OpenMM::Platform::loadPluginsFromDirectory',
+                            'Vec3 OpenMM::LocalCoordinatesSite::getOriginWeights',
+                            'Vec3 OpenMM::LocalCoordinatesSite::getXWeights',
+                            'Vec3 OpenMM::LocalCoordinatesSite::getYWeights']
         self.hideClasses = ['Kernel', 'KernelImpl', 'KernelFactory', 'ContextImpl', 'SerializationNode', 'SerializationProxy']
         self.nodeByID={}
 
@@ -119,9 +127,7 @@ class WrapperGenerator:
         for section in findNodes(classNode, "sectiondef", kind="public-static-func")+findNodes(classNode, "sectiondef", kind="public-func"):
             for memberNode in findNodes(section, "memberdef", kind="function", prot="public"):
                 methodDefinition = getText("definition", memberNode)
-                shortMethodDefinition = stripOpenMMPrefix(methodDefinition)
-                methodName = shortMethodDefinition.split()[-1]
-                if className+'::'+methodName in self.skipMethods:
+                if methodDefinition in self.skipMethods:
                     continue
                 methodList.append(memberNode)
         return methodList

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -108,6 +108,9 @@ SKIP_METHODS = [('State', 'getPositions'),
                 ('IntegrateDrudeSCFStepKernel',),
                 ('XmlSerializer',  'serialize'),
                 ('XmlSerializer',  'deserialize'),
+                ('LocalCoordinatesSite',  'getOriginWeights', 0),
+                ('LocalCoordinatesSite',  'getXWeights', 0),
+                ('LocalCoordinatesSite',  'getYWeights', 0),
 ]
 
 # The build script assumes method args that are non-const references are