Revised the virtual site API to fix some problems with it

openmmapi/include/openmm/System.h

@@ -59,7 +59,7 @@ class OPENMM_EXPORT VirtualSite;
  * 
  * 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
+ * To define a virtual site, call setVirtualSite(), passing in a VirtualSite object
  * that defines the rules for computing its position.
  */
 
@@ -112,6 +112,31 @@ public:
     void setParticleMass(int index, double mass) {
         masses[index] = mass;
     }
+    /**
+     * Set a particle to be a virtual site.  The VirtualSite object should have
+     * been created on the heap with the "new" operator.  The System takes over
+     * ownership of it, and deletes it when the System itself is deleted.
+     *
+     * @param index        the index of the particle that should be treated as a
+     *                     virtual site
+     * @param virtualSite  a pointer to the VirtualSite object describing it
+     */
+    void setVirtualSite(int index, VirtualSite* virtualSite);
+    /**
+     * Get whether a particle is a VirtualSite.
+     * 
+     * @param index  the index of the particle to check
+     */
+    bool isVirtualSite(int index) const {
+        return (index < virtualSites.size() && virtualSites[index] != NULL);
+    }
+    /**
+     * Get VirtualSite object for a particle.  If the particle is not a virtual
+     * site, this throws an exception.
+     * 
+     * @param index  the index of the particle to get
+     */
+    const VirtualSite& getVirtualSite(int index) const;
     /**
      * Get the number of distance constraints in this System.
      */
@@ -181,32 +206,6 @@ 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

openmmapi/include/openmm/VirtualSite.h

@@ -41,7 +41,7 @@ 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.
+ * call setVirtualSite() on the System.
  */
 
 class OPENMM_EXPORT VirtualSite {

openmmapi/src/ContextImpl.cpp

@@ -37,6 +37,7 @@
 #include "openmm/internal/ForceImpl.h"
 #include "openmm/internal/ContextImpl.h"
 #include "openmm/State.h"
+#include "openmm/VirtualSite.h"
 #include <map>
 #include <utility>
 #include <vector>
@@ -51,6 +52,29 @@ ContextImpl::ContextImpl(Context& owner, System& system, Integrator& integrator,
          owner(owner), system(system), integrator(integrator), hasInitializedForces(false), platform(platform), platformData(NULL) {
     if (system.getNumParticles() == 0)
         throw OpenMMException("Cannot create a Context for a System with no particles");
+    
+    // Check for errors in virtual sites and massless particles.
+    
+    for (int i = 0; i < system.getNumParticles(); i++) {
+        if (system.isVirtualSite(i)) {
+            if (system.getParticleMass(i) != 0.0)
+                throw OpenMMException("Virtual site has nonzero mass");
+            const VirtualSite& site = system.getVirtualSite(i);
+            for (int j = 0; j < site.getNumParticles(); j++)
+                if (system.isVirtualSite(site.getParticle(j)))
+                    throw OpenMMException("A virtual site cannot depend on another virtual site");
+        }
+    }
+    for (int i = 0; i < system.getNumConstraints(); i++) {
+        int particle1, particle2;
+        double distance;
+        system.getConstraintParameters(i, particle1, particle2, distance);
+        if (system.getParticleMass(particle1) == 0.0 || system.getParticleMass(particle2) == 0.0)
+            throw OpenMMException("A constraint cannot involve a massless particle");
+    }
+    
+    // Find the list of kernels required.
+    
     vector<string> kernelNames;
     kernelNames.push_back(CalcKineticEnergyKernel::Name());
     kernelNames.push_back(CalcForcesAndEnergyKernel::Name());
@@ -69,6 +93,9 @@ ContextImpl::ContextImpl(Context& owner, System& system, Integrator& integrator,
         this->platform = platform = &Platform::findPlatform(kernelNames);
     else if (!platform->supportsKernels(kernelNames))
         throw OpenMMException("Specified a Platform for a Context which does not support all required kernels");
+    
+    // Create and initialize kernels and other objects.
+    
     platform->contextCreated(*this, properties);
     initializeForcesKernel = platform->createKernel(CalcForcesAndEnergyKernel::Name(), *this);
     dynamic_cast<CalcForcesAndEnergyKernel&>(initializeForcesKernel.getImpl()).initialize(system);

openmmapi/src/System.cpp

@@ -48,6 +48,18 @@ System::~System() {
         delete virtualSites[i];
 }
 
+void System::setVirtualSite(int index, VirtualSite* virtualSite) {
+    if (index >= virtualSites.size())
+        virtualSites.resize(getNumParticles(), NULL);
+    virtualSites[index] = virtualSite;
+}
+
+const VirtualSite& System::getVirtualSite(int index) const {
+    if (index >= virtualSites.size() || virtualSites[index] == NULL)
+        throw OpenMMException("This particle is not a virtual site");
+    return *virtualSites[index];
+}
+
 int System::addConstraint(int particle1, int particle2, double distance) {
     constraints.push_back(ConstraintInfo(particle1, particle2, distance));
     return constraints.size()-1;