Created API for virtual sites

openmmapi/include/openmm/CustomIntegrator.h

@@ -79,6 +79,10 @@ namespace OpenMM {
  * velocity along any constrained distance is 0.</li>
  * </ul>
  * 
+ * Like all integrators, CustomIntegrator ignores any particle whose mass is 0.
+ * It is skipped when doing per-DOF computations, and is not included when
+ * computing sums over degrees of freedom.
+ * 
  * In addition to the variables you define by calling addGlobalVariable() and
  * addPerDofVariable(), the integrator provides the following pre-defined
  * variables:

openmmapi/include/openmm/System.h

@@ -39,6 +39,7 @@
 namespace OpenMM {
 
 class OPENMM_EXPORT Force;
+class OPENMM_EXPORT VirtualSite;
 
 /**
  * This class represents a molecular system.  The definition of a System involves
@@ -55,6 +56,11 @@ class OPENMM_EXPORT Force;
  * forces are defined by objects that extend the Force class.  After creating a
  * System, call addParticle() once for each particle, addConstraint() for each constraint,
  * and addForce() for each Force.
+ * 
+ * In addition, particles may be designated as "virtual sites".  These are particles
+ * whose positions are computed automatically based on the positions of other particles.
+ * To define a virtual site, call addVirtualSite(), passing in a VirtualSite object
+ * that defines the rules for computing its position.
  */
 
 class OPENMM_EXPORT System {
@@ -71,7 +77,10 @@ public:
         return masses.size();
     }
     /**
-     * Add a particle to the System.
+     * Add a particle to the System.  If the mass is 0, Integrators will ignore
+     * the particle and not modify its position or velocity.  This is most often
+     * used for virtual sites, but can also be used as a way to prevent a particle
+     * from moving.
      *
      * @param mass   the mass of the particle (in atomic mass units)
      * @return the index of the particle that was added
@@ -81,7 +90,10 @@ public:
         return masses.size()-1;
     }
     /**
-     * Get the mass (in atomic mass units) of a particle.
+     * Get the mass (in atomic mass units) of a particle.  If the mass is 0, Integrators will ignore
+     * the particle and not modify its position or velocity.  This is most often
+     * used for virtual sites, but can also be used as a way to prevent a particle
+     * from moving.
      * 
      * @param index the index of the particle for which to get the mass
      */
@@ -89,7 +101,10 @@ public:
         return masses[index];
     }
     /**
-     * Set the mass (in atomic mass units) of a particle.
+     * Set the mass (in atomic mass units) of a particle.  If the mass is 0, Integrators will ignore
+     * the particle and not modify its position or velocity.  This is most often
+     * used for virtual sites, but can also be used as a way to prevent a particle
+     * from moving.
      * 
      * @param index  the index of the particle for which to set the mass
      * @param mass   the mass of the particle
@@ -104,7 +119,8 @@ public:
         return constraints.size();
     }
     /**
-     * Add a constraint to the System.
+     * Add a constraint to the System.  Particles whose mass is 0 cannot participate
+     * in constraints.
      *
      * @param particle1 the index of the first particle involved in the constraint
      * @param particle2 the index of the second particle involved in the constraint
@@ -122,7 +138,8 @@ public:
      */
     void getConstraintParameters(int index, int& particle1, int& particle2, double& distance) const;
     /**
-     * Set the parameters defining a distance constraint.
+     * Set the parameters defining a distance constraint.  Particles whose mass is 0 cannot participate
+     * in constraints.
      * 
      * @param index     the index of the constraint for which to set parameters
      * @param particle1 the index of the first particle involved in the constraint
@@ -164,6 +181,32 @@ public:
     Force& getForce(int index) {
         return *forces[index];
     }
+    /**
+     * Add a VirtualSite to the System.  The VirtualSite should have been created on the heap with the
+     * "new" operator.  The System takes over ownership of it, and deletes the VirtualSite when the
+     * System itself is deleted.
+     *
+     * @param virtualSite   a pointer to the VirtualSite object to be added
+     * @return the index within the System of the VirtualSite that was added
+     */
+    int addVirtualSite(VirtualSite* virtualSite) {
+        virtualSites.push_back(virtualSite);
+        return virtualSites.size()-1;
+    }
+    /**
+     * Get the number of VirtualSite objects that have been added to the System.
+     */
+    int getNumVirtualSites() const {
+        return virtualSites.size();
+    }
+    /**
+     * Get a const reference to one of the VirtualSites in this System.
+     * 
+     * @param index  the index of the VirtualSite to get
+     */
+    const VirtualSite& getVirtualSite(int index) const {
+        return *virtualSites[index];
+    }
     /**
      * Get the default values of the vectors defining the axes of the periodic box (measured in nm).  Any newly
      * created Context will have its box vectors set to these.  They will affect
@@ -196,6 +239,7 @@ private:
     std::vector<double> masses;
     std::vector<ConstraintInfo> constraints;
     std::vector<Force*> forces;
+    std::vector<VirtualSite*> virtualSites;
 };
 
 /**

openmmapi/include/openmm/VirtualSite.h

+#ifndef OPENMM_VIRTUALSITE_H_
+#define OPENMM_VIRTUALSITE_H_
+
+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2012 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "internal/windowsExport.h"
+#include <vector>
+
+namespace OpenMM {
+
+/**
+ * A VirtualSite describes the rules for computing a particle's position based on
+ * other particles.  This is an abstract class.  Subclasses define particular rules.
+ * To define a virtual site, create an instance of a VirtualSite subclass and then
+ * call addVirtualSite() on the System.
+ */
+
+class OPENMM_EXPORT VirtualSite {
+public:
+    class TwoParticleAverage;
+    class ThreeParticleAverage;
+    class OutOfPlane;
+    virtual ~VirtualSite() {
+    }
+    /**
+     * Get the number of particles this virtual site depends on.
+     */
+    int getNumParticles() const;
+    /**
+     * Get the index of a particle this virtual site depends on.
+     * 
+     * @param particle    the particle to get (between 0 and getNumParticles())
+     * @return the index of the particle in the System
+     */
+    int getParticle(int particle) const;
+protected:
+    VirtualSite() {
+    }
+    void setParticles(const std::vector<int>& particleIndices);
+private:
+    std::vector<int> particles;
+};
+
+/**
+ * This is a VirtualSite that computes the particle location as a weighted average
+ * of two other particle's locations.  Assuming the weights add up to 1, this means
+ * the virtual site is on the line passing through the two particles.
+ */
+class VirtualSite::TwoParticleAverage : public VirtualSite {
+public:
+    /**
+     * Create a new TwoParticleAverage virtual site.  Normally weight1 and weight2
+     * should add up to 1, although this is not strictly required.
+     * 
+     * @param particle1    the index of the first particle
+     * @param particle2    the index of the second particle
+     * @param weight1      the weight factor (between 0 and 1) for the first particle
+     * @param weight2      the weight factor (between 0 and 1) for the second particle
+     */
+    TwoParticleAverage(int particle1, int particle2, double weight1, double weight2);
+    /**
+     * Get the weight factor used for a particle this virtual site depends on.
+     * 
+     * @param particle    the particle to get (between 0 and getNumParticles())
+     * @return the weight factor used for that particle
+     */
+    double getWeight(int particle) const;
+private:
+    double weight1, weight2;
+};
+
+/**
+ * This is a VirtualSite that computes the particle location as a weighted average
+ * of three other particle's locations.  Assuming the weights add up to 1, this means
+ * the virtual site is in the plane of the three particles.
+ */
+class VirtualSite::ThreeParticleAverage : public VirtualSite {
+public:
+    /**
+     * Create a new ThreeParticleAverage virtual site.  Normally the weights
+     * should add up to 1, although this is not strictly required.
+     * 
+     * @param particle1    the index of the first particle
+     * @param particle2    the index of the second particle
+     * @param particle3    the index of the third particle
+     * @param weight1      the weight factor (between 0 and 1) for the first particle
+     * @param weight2      the weight factor (between 0 and 1) for the second particle
+     * @param weight2      the weight factor (between 0 and 1) for the third particle
+     */
+    ThreeParticleAverage(int particle1, int particle2, int particle3, double weight1, double weight2, double weight3);
+    /**
+     * Get the weight factor used for a particle this virtual site depends on.
+     * 
+     * @param particle    the particle to get (between 0 and getNumParticles())
+     * @return the weight factor used for that particle
+     */
+    double getWeight(int particle) const;
+private:
+    double weight1, weight2, weight3;
+};
+
+/**
+ * This is a VirtualSite that computes the particle location based on three other
+ * particles' locations.  If r<sub>1</sub> is the location of particle 1,
+ * r<sub>12</sub> is the vector from particle 1 to particle 2, and
+ * r<sub>13</sub> is the vector from particle 1 to particle 3, then the virtual
+ * site location is given by
+ * 
+ * r<sub>1</sub> + w<sub>12</sub>r<sub>12</sub> + w<sub>13</sub>r<sub>13</sub> + w<sub>cross</sub>(r<sub>12</sub>&times;r<sub>13</sub>)
+ * 
+ * The three weight factors are user-specified.  This allows the virtual site location
+ * to be out of the plane of the three particles.
+ */
+class VirtualSite::OutOfPlane : public VirtualSite {
+public:
+    /**
+     * Create a new OutOfPlane virtual site.
+     * 
+     * @param particle1    the index of the first particle
+     * @param particle2    the index of the second particle
+     * @param particle3    the index of the third particle
+     * @param weight12     the weight factor for the vector from particle1 to particle2
+     * @param weight13     the weight factor for the vector from particle1 to particle3
+     * @param weightCross  the weight factor for the cross product
+     */
+    OutOfPlane(int particle1, int particle2, int particle3, double weight12, double weight13, double weightCross);
+    /**
+     * Get the weight factor for the vector from particle1 to particle2.
+     */
+    double getWeight12() const;
+    /**
+     * Get the weight factor for the vector from particle1 to particle3.
+     */
+    double getWeight13() const;
+    /**
+     * Get the weight factor for the cross product.
+     */
+    double getWeightCross() const;
+private:
+    double weight12, weight13, weightCross;
+};
+
+} // namespace OpenMM
+
+#endif /*OPENMM_VIRTUALSITE_H_*/

openmmapi/src/System.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) 2008-2009 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -44,6 +44,8 @@ System::System() {
 System::~System() {
     for (int i = 0; i < (int) forces.size(); ++i)
         delete forces[i];
+    for (int i = 0; i < (int) virtualSites.size(); ++i)
+        delete virtualSites[i];
 }
 
 int System::addConstraint(int particle1, int particle2, double distance) {

openmmapi/src/VirtualSite.cpp

+/* -------------------------------------------------------------------------- *
+ *                                   OpenMM                                   *
+ * -------------------------------------------------------------------------- *
+ * This is part of the OpenMM molecular simulation toolkit originating from   *
+ * Simbios, the NIH National Center for Physics-Based Simulation of           *
+ * Biological Structures at Stanford, funded under the NIH Roadmap for        *
+ * Medical Research, grant U54 GM072970. See https://simtk.org.               *
+ *                                                                            *
+ * Portions copyright (c) 2012 Stanford University and the Authors.           *
+ * Authors: Peter Eastman                                                     *
+ * Contributors:                                                              *
+ *                                                                            *
+ * Permission is hereby granted, free of charge, to any person obtaining a    *
+ * copy of this software and associated documentation files (the "Software"), *
+ * to deal in the Software without restriction, including without limitation  *
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
+ * and/or sell copies of the Software, and to permit persons to whom the      *
+ * Software is furnished to do so, subject to the following conditions:       *
+ *                                                                            *
+ * The above copyright notice and this permission notice shall be included in *
+ * all copies or substantial portions of the Software.                        *
+ *                                                                            *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
+ * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
+ * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
+ * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
+ * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
+ * -------------------------------------------------------------------------- */
+
+#include "openmm/VirtualSite.h"
+#include "openmm/OpenMMException.h"
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+void VirtualSite::setParticles(const vector<int>& particleIndices) {
+    particles = particleIndices;
+}
+
+int VirtualSite::getNumParticles() const {
+    return particles.size();
+}
+
+int VirtualSite::getParticle(int particle) const {
+    return particles[particle];
+}
+
+VirtualSite::TwoParticleAverage::TwoParticleAverage(int particle1, int particle2, double weight1, double weight2) :
+        weight1(weight1), weight2(weight2) {
+    vector<int> particles(2);
+    particles[0] = particle1;
+    particles[1] = particle2;
+    setParticles(particles);
+}
+
+double VirtualSite::TwoParticleAverage::getWeight(int particle) const {
+    if (particle == 0)
+        return weight1;
+    if (particle == 1)
+        return weight2;
+    throw OpenMMException("Illegal index for particle");
+}
+
+VirtualSite::ThreeParticleAverage::ThreeParticleAverage(int particle1, int particle2, int particle3, double weight1, double weight2, double weight3) :
+        weight1(weight1), weight2(weight2), weight3(weight3) {
+    vector<int> particles(3);
+    particles[0] = particle1;
+    particles[1] = particle2;
+    particles[2] = particle3;
+    setParticles(particles);
+}
+
+double VirtualSite::ThreeParticleAverage::getWeight(int particle) const {
+    if (particle == 0)
+        return weight1;
+    if (particle == 1)
+        return weight2;
+    if (particle == 2)
+        return weight3;
+    throw OpenMMException("Illegal index for particle");
+}
+
+VirtualSite::OutOfPlane::OutOfPlane(int particle1, int particle2, int particle3, double weight12, double weight13, double weightCross) :
+        weight12(weight12), weight13(weight13), weightCross(weightCross) {
+    vector<int> particles(3);
+    particles[0] = particle1;
+    particles[1] = particle2;
+    particles[2] = particle3;
+    setParticles(particles);
+}
+
+double VirtualSite::OutOfPlane::getWeight12() const {
+    return weight12;
+}
+
+double VirtualSite::OutOfPlane::getWeight13() const {
+    return weight13;
+}
+
+double VirtualSite::OutOfPlane::getWeightCross() const {
+    return weightCross;
+}