ForceField supports LocalCoordinatesSite depending on arbitrary particles

docs-source/usersguide/application.rst

@@ -2051,10 +2051,9 @@ be represented with a virtual site.  The :code:`type` attribute may equal
 :code:`"localCoords"`\ , which correspond to the TwoParticleAverageSite, ThreeParticleAverageSite,
 OutOfPlaneSite, and LocalCoordinatesSite classes respectively.  The :code:`siteName`
 attribute gives the name of the atom to represent with a virtual site.  The atoms
-it is calculated based on are specified by :code:`atomName1`\ , :code:`atomName2`\ ,
-and (for virtual site classes that involve three atoms) :code:`atomName3`\ .
+it is calculated based on are specified by :code:`atomName1`\ , :code:`atomName2`\ , etc.
 (Some old force fields use the deprecated tags :code:`index`, :code:`atom1`,
-:code:`atom2`, and :code:`atom3` to refer to them by index instead of name.)
+:code:`atom2`, etc. to refer to them by index instead of name.)
 
 The remaining attributes are specific to the virtual site class, and specify the
 parameters for calculating the site position.  For a TwoParticleAverageSite,
@@ -2062,9 +2061,10 @@ they are :code:`weight1` and :code:`weight2`\ .  For a
 ThreeParticleAverageSite, they are :code:`weight1`\ , :code:`weight2`\ , and
 \ :code:`weight3`\ . For an OutOfPlaneSite, they are :code:`weight12`\ ,
 :code:`weight13`\ , and :code:`weightCross`\ . For a LocalCoordinatesSite, they
-are :code:`wo1`\ , :code:`wo2`\ , :code:`wo3`\ , :code:`wx1`\ , :code:`wx2`\ ,
-:code:`wx3`\ , :code:`wy1`\ , :code:`wy2`\ , :code:`wy3`\ , :code:`p1`\ ,
-:code:`p2`\ , and :code:`p3`\ .
+are :code:`p1`\ , :code:`p2`\ , and :code:`p3` (giving the x, y, and z coordinates
+of the site position in the local coordinate system), and :code:`wo1`\ ,
+:code:`wx1`\ , :code:`wy1`\ , :code:`wo2`\ , :code:`wx2`\ , :code:`wy2`\ , ...
+(giving the weights for computing the origin, x axis, and y axis).
 
 <Patches>
 =========

docs-source/usersguide/theory.rst

@@ -1572,17 +1572,18 @@ specific types of rules.  They are:
   :math:`\mathbf{r}_{13} = \mathbf{r}_{3}-\mathbf{r}_{1}`\ .  This allows
   the virtual site to be located outside the plane of the three particles.
 
-* LocalCoordinatesSite: The locations of three other particles are used to compute a local
+* LocalCoordinatesSite: The locations of several other particles are used to compute a local
   coordinate system, and the virtual site is placed 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:
+  system.  The number of particles used to define the coordinate system is user defined.
+  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 other particles:
 
 .. math::
-   \mathbf{o}={w}^{o}_{1}\mathbf{r}_{1} + {w}^{o}_{2}\mathbf{r}_{2} + {w}^{o}_{3}\mathbf{r}_{3}
+   \mathbf{o}={w}^{o}_{1}\mathbf{r}_{1} + {w}^{o}_{2}\mathbf{r}_{2} + ...
 
-   \mathbf{dx}={w}^{x}_{1}\mathbf{r}_{1} + {w}^{x}_{2}\mathbf{r}_{2} + {w}^{x}_{3}\mathbf{r}_{3}
+   \mathbf{dx}={w}^{x}_{1}\mathbf{r}_{1} + {w}^{x}_{2}\mathbf{r}_{2} + ...
 
-   \mathbf{dy}={w}^{y}_{1}\mathbf{r}_{1} + {w}^{y}_{2}\mathbf{r}_{2} + {w}^{y}_{3}\mathbf{r}_{3}
+   \mathbf{dy}={w}^{y}_{1}\mathbf{r}_{1} + {w}^{y}_{2}\mathbf{r}_{2} + ...
 
    \mathbf{dz}=\mathbf{dx}\times \mathbf{dy}
 ..

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

@@ -647,10 +647,15 @@ class ForceField(object):
                 numAtoms = 3
                 self.weights = [float(attrib['weight12']), float(attrib['weight13']), float(attrib['weightCross'])]
             elif self.type == 'localCoords':
-                numAtoms = 3
-                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'])]
+                numAtoms = 0
+                self.originWeights = []
+                self.xWeights = []
+                self.yWeights = []
+                while ('wo%d' % (numAtoms+1)) in attrib:
+                    numAtoms += 1
+                    self.originWeights.append(float(attrib['wo%d' % numAtoms]))
+                    self.xWeights.append(float(attrib['wx%d' % numAtoms]))
+                    self.yWeights.append(float(attrib['wy%d' % numAtoms]))
                 self.localPos = [float(attrib['p1']), float(attrib['p2']), float(attrib['p3'])]
             else:
                 raise ValueError('Unknown virtual site type: %s' % self.type)
@@ -1255,11 +1260,7 @@ class ForceField(object):
             elif site.type == 'outOfPlane':
                 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])))
+                sys.setVirtualSite(index, mm.LocalCoordinatesSite(atoms, site.originWeights, site.xWeights, site.yWeights, site.localPos))
 
         # Add forces to the System