OpenCL implementation of virtual sites

platforms/opencl/src/OpenCLContext.cpp

@@ -38,10 +38,12 @@
 #include "hilbert.h"
 #include "openmm/Platform.h"
 #include "openmm/System.h"
+#include "openmm/VirtualSite.h"
 #include <algorithm>
 #include <fstream>
 #include <iostream>
 #include <sstream>
+#include <typeinfo>
 
 using namespace OpenMM;
 using namespace std;
@@ -230,8 +232,10 @@ OpenCLContext::~OpenCLContext() {
 }
 
 void OpenCLContext::initialize(const System& system) {
-    for (int i = 0; i < numAtoms; i++)
-        (*velm)[i].w = (float) (1.0/system.getParticleMass(i));
+    for (int i = 0; i < numAtoms; i++) {
+        double mass = system.getParticleMass(i);
+        (*velm)[i].w = (float) (mass == 0.0 ? 0.0 : 1.0/mass);
+    }
     velm->upload();
     bonded->initialize(system);
     numForceBuffers = platformData.contexts.size();
@@ -462,7 +466,77 @@ struct OpenCLContext::Molecule {
     vector<vector<int> > groups;
 };
 
+/**
+ * This class ensures that atom reordering doesn't break virtual sites.
+ */
+class OpenCLContext::VirtualSiteInfo : public OpenCLForceInfo {
+public:
+    VirtualSiteInfo(const System& system) : OpenCLForceInfo(0) {
+        for (int i = 0; i < system.getNumParticles(); i++) {
+            if (system.isVirtualSite(i)) {
+                siteTypes.push_back(&typeid(system.getVirtualSite(i)));
+                vector<int> particles;
+                particles.push_back(i);
+                for (int j = 0; j < system.getVirtualSite(i).getNumParticles(); j++)
+                    particles.push_back(system.getVirtualSite(i).getParticle(j));
+                siteParticles.push_back(particles);
+                vector<double> weights;
+                if (dynamic_cast<const TwoParticleAverageSite*>(&system.getVirtualSite(i)) != NULL) {
+                    // A two particle average.
+
+                    const TwoParticleAverageSite& site = dynamic_cast<const TwoParticleAverageSite&>(system.getVirtualSite(i));
+                    weights.push_back(site.getWeight(0));
+                    weights.push_back(site.getWeight(1));
+                }
+                else if (dynamic_cast<const ThreeParticleAverageSite*>(&system.getVirtualSite(i)) != NULL) {
+                    // A three particle average.
+
+                    const ThreeParticleAverageSite& site = dynamic_cast<const ThreeParticleAverageSite&>(system.getVirtualSite(i));
+                    weights.push_back(site.getWeight(0));
+                    weights.push_back(site.getWeight(1));
+                    weights.push_back(site.getWeight(2));
+                }
+                else if (dynamic_cast<const OutOfPlaneSite*>(&system.getVirtualSite(i)) != NULL) {
+                    // An out of plane site.
+
+                    const OutOfPlaneSite& site = dynamic_cast<const OutOfPlaneSite&>(system.getVirtualSite(i));
+                    weights.push_back(site.getWeight12());
+                    weights.push_back(site.getWeight13());
+                    weights.push_back(site.getWeightCross());
+                }
+                siteWeights.push_back(weights);
+            }
+        }
+    }
+    int getNumParticleGroups() {
+        return siteTypes.size();
+    }
+    void getParticlesInGroup(int index, std::vector<int>& particles) {
+        particles = siteParticles[index];
+    }
+    bool areGroupsIdentical(int group1, int group2) {
+        if (siteTypes[group1] != siteTypes[group2])
+            return false;
+        int numParticles = siteWeights[group1].size();
+        if (siteWeights[group2].size() != numParticles)
+            return false;
+        for (int i = 0; i < numParticles; i++)
+            if (siteWeights[group1][i] != siteWeights[group2][i])
+                return false;
+        return true;
+    }
+private:
+    vector<const type_info*> siteTypes;
+    vector<vector<int> > siteParticles;
+    vector<vector<double> > siteWeights;
+};
+
+
 void OpenCLContext::findMoleculeGroups(const System& system) {
+    // Add a ForceInfo that makes sure reordering doesn't break virtual sites.
+    
+    addForce(new VirtualSiteInfo(system));
+    
     // First make a list of every other atom to which each atom is connect by a constraint or force group.
 
     vector<vector<int> > atomBonds(system.getNumParticles());

platforms/opencl/src/OpenCLContext.h

@@ -476,6 +476,7 @@ public:
 private:
     struct Molecule;
     struct MoleculeGroup;
+    class VirtualSiteInfo;
     void findMoleculeGroups(const System& system);
     static void tagAtomsInMolecule(int atom, int molecule, std::vector<int>& atomMolecule, std::vector<std::vector<int> >& atomBonds);
     double time;

platforms/opencl/src/OpenCLIntegrationUtilities.cpp

@@ -28,6 +28,7 @@
 #include "OpenCLArray.h"
 #include "OpenCLKernelSources.h"
 #include "openmm/HarmonicAngleForce.h"
+#include "openmm/VirtualSite.h"
 #include "quern.h"
 #include "OpenCLExpressionUtilities.h"
 #include <algorithm>
@@ -91,10 +92,13 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
         random(NULL), randomSeed(NULL), randomPos(0), stepSize(NULL), ccmaAtoms(NULL), ccmaDistance(NULL),
         ccmaReducedMass(NULL), ccmaAtomConstraints(NULL), ccmaNumAtomConstraints(NULL), ccmaConstraintMatrixColumn(NULL),
         ccmaConstraintMatrixValue(NULL), ccmaDelta1(NULL), ccmaDelta2(NULL), ccmaConverged(NULL),
-        ccmaConvergedBuffer(NULL), hasInitializedPosConstraintKernels(false), hasInitializedVelConstraintKernels(false) {
+        ccmaConvergedBuffer(NULL), vsite2AvgAtoms(NULL), vsite2AvgWeights(NULL), vsite3AvgAtoms(NULL), vsite3AvgWeights(NULL),
+        vsiteOutOfPlaneAtoms(NULL), vsiteOutOfPlaneWeights(NULL), hasInitializedPosConstraintKernels(false), hasInitializedVelConstraintKernels(false) {
     // Create workspace arrays.
 
     posDelta = new OpenCLArray<mm_float4>(context, context.getPaddedNumAtoms(), "posDelta");
+    vector<mm_float4> deltas(posDelta->getSize(), mm_float4(0.0, 0.0, 0.0, 0.0));
+    posDelta->upload(deltas);
     stepSize = new OpenCLArray<mm_float2>(context, 1, "stepSize", true);
     stepSize->set(0, mm_float2(0.0f, 0.0f));
     stepSize->upload();
@@ -515,6 +519,87 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
         ccmaVelForceKernel = cl::Kernel(ccmaProgram, "computeConstraintForce");
         ccmaVelUpdateKernel = cl::Kernel(ccmaProgram, "updateAtomPositions");
     }
+    
+    // Build the list of virtual sites.
+    
+    vector<mm_int4> vsite2AvgAtomVec;
+    vector<mm_float2> vsite2AvgWeightVec;
+    vector<mm_int4> vsite3AvgAtomVec;
+    vector<mm_float4> vsite3AvgWeightVec;
+    vector<mm_int4> vsiteOutOfPlaneAtomVec;
+    vector<mm_float4> vsiteOutOfPlaneWeightVec;
+    for (int i = 0; i < numAtoms; i++) {
+        if (system.isVirtualSite(i)) {
+            if (dynamic_cast<const TwoParticleAverageSite*>(&system.getVirtualSite(i)) != NULL) {
+                // A two particle average.
+                
+                const TwoParticleAverageSite& site = dynamic_cast<const TwoParticleAverageSite&>(system.getVirtualSite(i));
+                vsite2AvgAtomVec.push_back(mm_int4(i, site.getParticle(0), site.getParticle(1), 0));
+                vsite2AvgWeightVec.push_back(mm_float2((float) site.getWeight(0), (float) site.getWeight(1)));
+            }
+            else if (dynamic_cast<const ThreeParticleAverageSite*>(&system.getVirtualSite(i)) != NULL) {
+                // A three particle average.
+                
+                const ThreeParticleAverageSite& site = dynamic_cast<const ThreeParticleAverageSite&>(system.getVirtualSite(i));
+                vsite3AvgAtomVec.push_back(mm_int4(i, site.getParticle(0), site.getParticle(1), site.getParticle(2)));
+                vsite3AvgWeightVec.push_back(mm_float4((float) site.getWeight(0), (float) site.getWeight(1), (float) site.getWeight(2), 0.0f));
+            }
+            else if (dynamic_cast<const OutOfPlaneSite*>(&system.getVirtualSite(i)) != NULL) {
+                // An out of plane site.
+                
+                const OutOfPlaneSite& site = dynamic_cast<const OutOfPlaneSite&>(system.getVirtualSite(i));
+                vsiteOutOfPlaneAtomVec.push_back(mm_int4(i, site.getParticle(0), site.getParticle(1), site.getParticle(2)));
+                vsiteOutOfPlaneWeightVec.push_back(mm_float4((float) site.getWeight12(), (float) site.getWeight13(), (float) site.getWeightCross(), 0.0f));
+            }
+        }
+    }
+    int num2Avg = vsite2AvgAtomVec.size();
+    int num3Avg = vsite3AvgAtomVec.size();
+    int numOutOfPlane = vsiteOutOfPlaneAtomVec.size();
+    vsite2AvgAtoms = new OpenCLArray<mm_int4>(context, max(1, num2Avg), "vsite2AvgAtoms");
+    vsite2AvgWeights = new OpenCLArray<mm_float2>(context, max(1, num2Avg), "vsite2AvgWeights");
+    vsite3AvgAtoms = new OpenCLArray<mm_int4>(context, max(1, num3Avg), "vsite3AvgAtoms");
+    vsite3AvgWeights = new OpenCLArray<mm_float4>(context, max(1, num3Avg), "vsite3AvgWeights");
+    vsiteOutOfPlaneAtoms = new OpenCLArray<mm_int4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneAtoms");
+    vsiteOutOfPlaneWeights = new OpenCLArray<mm_float4>(context, max(1, numOutOfPlane), "vsiteOutOfPlaneWeights");
+    if (num2Avg > 0) {
+        vsite2AvgAtoms->upload(vsite2AvgAtomVec);
+        vsite2AvgWeights->upload(vsite2AvgWeightVec);
+    }
+    if (num3Avg > 0) {
+        vsite3AvgAtoms->upload(vsite3AvgAtomVec);
+        vsite3AvgWeights->upload(vsite3AvgWeightVec);
+    }
+    if (numOutOfPlane > 0) {
+        vsiteOutOfPlaneAtoms->upload(vsiteOutOfPlaneAtomVec);
+        vsiteOutOfPlaneWeights->upload(vsiteOutOfPlaneWeightVec);
+    }
+    
+    // Create the kernels for virtual sites.
+
+    map<string, string> defines;
+    defines["NUM_2_AVERAGE"] = OpenCLExpressionUtilities::intToString(num2Avg);
+    defines["NUM_3_AVERAGE"] = OpenCLExpressionUtilities::intToString(num3Avg);
+    defines["NUM_OUT_OF_PLANE"] = OpenCLExpressionUtilities::intToString(numOutOfPlane);
+    cl::Program ccmaProgram = context.createProgram(OpenCLKernelSources::virtualSites, defines);
+    vsitePositionKernel = cl::Kernel(ccmaProgram, "computeVirtualSites");
+    vsitePositionKernel.setArg<cl::Buffer>(0, context.getPosq().getDeviceBuffer());
+    vsitePositionKernel.setArg<cl::Buffer>(1, vsite2AvgAtoms->getDeviceBuffer());
+    vsitePositionKernel.setArg<cl::Buffer>(2, vsite2AvgWeights->getDeviceBuffer());
+    vsitePositionKernel.setArg<cl::Buffer>(3, vsite3AvgAtoms->getDeviceBuffer());
+    vsitePositionKernel.setArg<cl::Buffer>(4, vsite3AvgWeights->getDeviceBuffer());
+    vsitePositionKernel.setArg<cl::Buffer>(5, vsiteOutOfPlaneAtoms->getDeviceBuffer());
+    vsitePositionKernel.setArg<cl::Buffer>(6, vsiteOutOfPlaneWeights->getDeviceBuffer());
+    vsiteForceKernel = cl::Kernel(ccmaProgram, "distributeForces");
+    vsiteForceKernel.setArg<cl::Buffer>(0, context.getPosq().getDeviceBuffer());
+    vsiteForceKernel.setArg<cl::Buffer>(1, context.getForce().getDeviceBuffer());
+    vsiteForceKernel.setArg<cl::Buffer>(2, vsite2AvgAtoms->getDeviceBuffer());
+    vsiteForceKernel.setArg<cl::Buffer>(3, vsite2AvgWeights->getDeviceBuffer());
+    vsiteForceKernel.setArg<cl::Buffer>(4, vsite3AvgAtoms->getDeviceBuffer());
+    vsiteForceKernel.setArg<cl::Buffer>(5, vsite3AvgWeights->getDeviceBuffer());
+    vsiteForceKernel.setArg<cl::Buffer>(6, vsiteOutOfPlaneAtoms->getDeviceBuffer());
+    vsiteForceKernel.setArg<cl::Buffer>(7, vsiteOutOfPlaneWeights->getDeviceBuffer());
+    numVsites = num2Avg+num3Avg+numOutOfPlane;
 }
 
 OpenCLIntegrationUtilities::~OpenCLIntegrationUtilities() {
@@ -556,6 +641,18 @@ OpenCLIntegrationUtilities::~OpenCLIntegrationUtilities() {
         delete ccmaConverged;
     if (ccmaConvergedBuffer != NULL)
         delete ccmaConvergedBuffer;
+    if (vsite2AvgAtoms != NULL)
+        delete vsite2AvgAtoms;
+    if (vsite2AvgWeights != NULL)
+        delete vsite2AvgWeights;
+    if (vsite3AvgAtoms != NULL)
+        delete vsite3AvgAtoms;
+    if (vsite3AvgWeights != NULL)
+        delete vsite3AvgWeights;
+    if (vsiteOutOfPlaneAtoms != NULL)
+        delete vsiteOutOfPlaneAtoms;
+    if (vsiteOutOfPlaneWeights != NULL)
+        delete vsiteOutOfPlaneWeights;
 }
 
 void OpenCLIntegrationUtilities::applyConstraints(double tol) {
@@ -660,6 +757,16 @@ void OpenCLIntegrationUtilities::applyConstraints(bool constrainVelocities, doub
     }
 }
 
+void OpenCLIntegrationUtilities::computeVirtualSites() {
+    if (numVsites > 0)
+        context.executeKernel(vsitePositionKernel, numVsites);
+}
+
+void OpenCLIntegrationUtilities::distributeForcesFromVirtualSites() {
+    if (numVsites > 0)
+        context.executeKernel(vsiteForceKernel, numVsites);
+}
+
 void OpenCLIntegrationUtilities::initRandomNumberGenerator(unsigned int randomNumberSeed) {
     if (random != NULL) {
         if (randomNumberSeed != lastSeed)

platforms/opencl/src/OpenCLIntegrationUtilities.h

@@ -84,6 +84,14 @@ public:
      * @return the index in the array at which to start reading
      */
     int prepareRandomNumbers(int numValues);
+    /**
+     * Compute the positions of virtual sites.
+     */
+    void computeVirtualSites();
+    /**
+     * Distribute forces from virtual sites to the atoms they are based on.
+     */
+    void distributeForcesFromVirtualSites();
 private:
     void applyConstraints(bool constrainVelocities, double tol);
     OpenCLContext& context;
@@ -93,6 +101,7 @@ private:
     cl::Kernel ccmaPosForceKernel, ccmaVelForceKernel;
     cl::Kernel ccmaMultiplyKernel;
     cl::Kernel ccmaPosUpdateKernel, ccmaVelUpdateKernel;
+    cl::Kernel vsitePositionKernel, vsiteForceKernel;
     cl::Kernel randomKernel;
     OpenCLArray<mm_float4>* posDelta;
     OpenCLArray<mm_int4>* settleAtoms;
@@ -114,8 +123,14 @@ private:
     OpenCLArray<cl_int>* ccmaConverged;
     cl::Buffer* ccmaConvergedBuffer;
     cl_int* ccmaConvergedMemory;
+    OpenCLArray<mm_int4>* vsite2AvgAtoms;
+    OpenCLArray<mm_float2>* vsite2AvgWeights;
+    OpenCLArray<mm_int4>* vsite3AvgAtoms;
+    OpenCLArray<mm_float4>* vsite3AvgWeights;
+    OpenCLArray<mm_int4>* vsiteOutOfPlaneAtoms;
+    OpenCLArray<mm_float4>* vsiteOutOfPlaneWeights;
     int randomPos;
-    int lastSeed;
+    int lastSeed, numVsites;
     bool hasInitializedPosConstraintKernels, hasInitializedVelConstraintKernels;
     struct ShakeCluster;
     struct ConstraintOrderer;

platforms/opencl/src/OpenCLKernels.cpp

@@ -108,6 +108,7 @@ double OpenCLCalcForcesAndEnergyKernel::finishComputation(ContextImpl& context,
     cl.getBondedUtilities().computeInteractions();
     cl.getNonbondedUtilities().computeInteractions();
     cl.reduceForces();
+    cl.getIntegrationUtilities().distributeForcesFromVirtualSites();
     double sum = 0.0f;
     if (includeEnergy) {
         OpenCLArray<cl_float>& energy = cl.getEnergyBuffer();
@@ -221,6 +222,7 @@ void OpenCLApplyConstraintsKernel::initialize(const System& system) {
 
 void OpenCLApplyConstraintsKernel::apply(ContextImpl& context, double tol) {
     cl.getIntegrationUtilities().applyConstraints(tol);
+    cl.getIntegrationUtilities().computeVirtualSites();
 }
 
 class OpenCLBondForceInfo : public OpenCLForceInfo {
@@ -3216,6 +3218,7 @@ void OpenCLIntegrateVerletStepKernel::execute(ContextImpl& context, const Verlet
     // Call the second integration kernel.
 
     cl.executeKernel(kernel2, numAtoms);
+    integration.computeVirtualSites();
 
     // Update the time and step count.
 
@@ -3292,6 +3295,7 @@ void OpenCLIntegrateLangevinStepKernel::execute(ContextImpl& context, const Lang
     // Call the second integration kernel.
 
     cl.executeKernel(kernel2, numAtoms);
+    integration.computeVirtualSites();
 
     // Update the time and step count.
 
@@ -3351,6 +3355,7 @@ void OpenCLIntegrateBrownianStepKernel::execute(ContextImpl& context, const Brow
     // Call the second integration kernel.
 
     cl.executeKernel(kernel2, numAtoms);
+    integration.computeVirtualSites();
 
     // Update the time and step count.
 
@@ -3411,6 +3416,7 @@ void OpenCLIntegrateVariableVerletStepKernel::execute(ContextImpl& context, cons
     // Call the second integration kernel.
 
     cl.executeKernel(kernel2, numAtoms);
+    integration.computeVirtualSites();
 
     // Update the time and step count.
 
@@ -3488,6 +3494,7 @@ void OpenCLIntegrateVariableLangevinStepKernel::execute(ContextImpl& context, co
     // Call the second integration kernel.
 
     cl.executeKernel(kernel2, numAtoms);
+    integration.computeVirtualSites();
 
     // Update the time and step count.
 
@@ -3958,6 +3965,7 @@ void OpenCLIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegr
         else if (stepType[i] == CustomIntegrator::ConstrainPositions) {
             cl.getIntegrationUtilities().applyConstraints(integrator.getConstraintTolerance());
             cl.executeKernel(kernels[i][0], numAtoms);
+            cl.getIntegrationUtilities().computeVirtualSites();
         }
         else if (stepType[i] == CustomIntegrator::ConstrainVelocities) {
             cl.getIntegrationUtilities().applyVelocityConstraints(integrator.getConstraintTolerance());

platforms/opencl/src/kernels/brownian.cl

@@ -7,7 +7,8 @@ __kernel void integrateBrownianPart1(float tauDeltaT, float noiseAmplitude, __gl
     randomIndex += get_global_id(0);
     for (int index = get_global_id(0); index < NUM_ATOMS; index += get_global_size(0)) {
         float invMass = velm[index].w;
-        posDelta[index] = (float4) (tauDeltaT*invMass*force[index].xyz + noiseAmplitude*sqrt(invMass)*random[randomIndex].xyz, 0.0f);
+        if (invMass != 0.0)
+            posDelta[index] = (float4) (tauDeltaT*invMass*force[index].xyz + noiseAmplitude*sqrt(invMass)*random[randomIndex].xyz, 0.0f);
         randomIndex += get_global_size(0);
     }
 }
@@ -18,8 +19,10 @@ __kernel void integrateBrownianPart1(float tauDeltaT, float noiseAmplitude, __gl
 
 __kernel void integrateBrownianPart2(float oneOverDeltaT, __global float4* posq, __global float4* velm, __global const float4* restrict posDelta) {
     for (int index = get_global_id(0); index < NUM_ATOMS; index += get_global_size(0)) {
-        float4 delta = posDelta[index];
-        velm[index].xyz = oneOverDeltaT*delta.xyz;
-        posq[index].xyz = posq[index].xyz + delta.xyz;
+        if (velm[index].w != 0.0) {
+            float4 delta = posDelta[index];
+            velm[index].xyz = oneOverDeltaT*delta.xyz;
+            posq[index].xyz = posq[index].xyz + delta.xyz;
+        }
     }
 }

platforms/opencl/src/kernels/customIntegratorPerDof.cl

@@ -30,9 +30,11 @@ __kernel void computePerDof(__global float4* restrict posq, __global float4* res
         float4 f = force[index];
         float mass = 1.0f/velocity.w;
 #endif
-        float4 gaussian = gaussianValues[randomIndex];
-        float4 uniform = uniformValues[index];
-        COMPUTE_STEP
+        if (velocity.w != 0.0) {
+            float4 gaussian = gaussianValues[randomIndex];
+            float4 uniform = uniformValues[index];
+            COMPUTE_STEP
+        }
         randomIndex += get_global_size(0);
         index += get_global_size(0);
     }

platforms/opencl/src/kernels/langevin.cl

@@ -18,10 +18,12 @@ __kernel void integrateLangevinPart1(__global float4* restrict velm, __global co
     randomIndex += index;
     while (index < NUM_ATOMS) {
         float4 velocity = velm[index];
-        float sqrtInvMass = sqrt(velocity.w);
-        velocity.xyz = vscale*velocity.xyz + fscale*velocity.w*force[index].xyz + noisescale*sqrtInvMass*random[randomIndex].xyz;
-        velm[index] = velocity;
-        posDelta[index] = stepSize*velocity;
+        if (velocity.w != 0.0) {
+            float sqrtInvMass = sqrt(velocity.w);
+            velocity.xyz = vscale*velocity.xyz + fscale*velocity.w*force[index].xyz + noisescale*sqrtInvMass*random[randomIndex].xyz;
+            velm[index] = velocity;
+            posDelta[index] = stepSize*velocity;
+        }
         randomIndex += get_global_size(0);
         index += get_global_size(0);
     }
@@ -39,17 +41,19 @@ __kernel void integrateLangevinPart2(__global float4* restrict posq, __global co
 #endif
     int index = get_global_id(0);
     while (index < NUM_ATOMS) {
-        float4 pos = posq[index];
-        float4 delta = posDelta[index];
         float4 vel = velm[index];
-        pos.xyz += delta.xyz;
+        if (vel.w != 0.0) {
+            float4 pos = posq[index];
+            float4 delta = posDelta[index];
+            pos.xyz += delta.xyz;
 #ifdef SUPPORTS_DOUBLE_PRECISION
-        vel.xyz = convert_float4(invStepSize*convert_double4(delta)).xyz;
+            vel.xyz = convert_float4(invStepSize*convert_double4(delta)).xyz;
 #else
-        vel.xyz = invStepSize*delta.xyz;
+            vel.xyz = invStepSize*delta.xyz;
 #endif
-        posq[index] = pos;
-        velm[index] = vel;
+            posq[index] = pos;
+            velm[index] = vel;
+        }
         index += get_global_size(0);
     }
 }

platforms/opencl/src/kernels/verlet.cl

@@ -12,12 +12,14 @@ __kernel void integrateVerletPart1(int numAtoms, __global const float2* restrict
     float dtVel = 0.5f*(stepSize.x+stepSize.y);
     int index = get_global_id(0);
     while (index < numAtoms) {
-        float4 pos = posq[index];
         float4 velocity = velm[index];
-        velocity.xyz += force[index].xyz*dtVel*velocity.w;
-        pos.xyz = velocity.xyz*dtPos;
-        posDelta[index] = pos;
-        velm[index] = velocity;
+        if (velocity.w != 0.0) {
+            float4 pos = posq[index];
+            velocity.xyz += force[index].xyz*dtVel*velocity.w;
+            pos.xyz = velocity.xyz*dtPos;
+            posDelta[index] = pos;
+            velm[index] = velocity;
+        }
         index += get_global_size(0);
     }
 }
@@ -38,17 +40,19 @@ __kernel void integrateVerletPart2(int numAtoms, __global float2* restrict dt, _
     barrier(CLK_LOCAL_MEM_FENCE);
     int index = get_global_id(0);
     while (index < numAtoms) {
-        float4 pos = posq[index];
-        float4 delta = posDelta[index];
         float4 velocity = velm[index];
-        pos.xyz += delta.xyz;
+        if (velocity.w != 0.0) {
+            float4 pos = posq[index];
+            float4 delta = posDelta[index];
+            pos.xyz += delta.xyz;
 #ifdef SUPPORTS_DOUBLE_PRECISION
-        velocity.xyz = convert_float4(convert_double4(delta)*oneOverDt).xyz;
+            velocity.xyz = convert_float4(convert_double4(delta)*oneOverDt).xyz;
 #else
-        velocity.xyz = delta.xyz*oneOverDt;
+            velocity.xyz = delta.xyz*oneOverDt;
 #endif
-        posq[index] = pos;
-        velm[index] = velocity;
+            posq[index] = pos;
+            velm[index] = velocity;
+        }
         index += get_global_size(0);
     }
 }

platforms/opencl/src/kernels/virtualSites.cl

+/**
+ * Compute the positions of virtual sites
+ */
+__kernel void computeVirtualSites(__global float4* restrict posq, __global const int4* restrict avg2Atoms, __global const float2* restrict avg2Weights,
+        __global const int4* restrict avg3Atoms, __global const float4* restrict avg3Weights,
+        __global const int4* restrict outOfPlaneAtoms, __global const float4* restrict outOfPlaneWeights) {
+    
+    // Two particle average sites.
+    
+    for (int index = get_global_id(0); index < NUM_2_AVERAGE; index += get_global_size(0)) {
+        int4 atoms = avg2Atoms[index];
+        float2 weights = avg2Weights[index];
+        float4 pos = posq[atoms.x];
+        float4 pos1 = posq[atoms.y];
+        float4 pos2 = posq[atoms.z];
+        pos.xyz = pos1.xyz*weights.x + pos2.xyz*weights.y;
+        posq[atoms.x] = pos;
+    }
+    
+    // Three particle average sites.
+    
+    for (int index = get_global_id(0); index < NUM_3_AVERAGE; index += get_global_size(0)) {
+        int4 atoms = avg3Atoms[index];
+        float4 weights = avg3Weights[index];
+        float4 pos = posq[atoms.x];
+        float4 pos1 = posq[atoms.y];
+        float4 pos2 = posq[atoms.z];
+        float4 pos3 = posq[atoms.w];
+        pos.xyz = pos1.xyz*weights.x + pos2.xyz*weights.y + pos3.xyz*weights.z;
+        posq[atoms.x] = pos;
+    }
+    
+    // Out of plane sites.
+    
+    for (int index = get_global_id(0); index < NUM_OUT_OF_PLANE; index += get_global_size(0)) {
+        int4 atoms = outOfPlaneAtoms[index];
+        float4 weights = outOfPlaneWeights[index];
+        float4 pos = posq[atoms.x];
+        float4 pos1 = posq[atoms.y];
+        float4 pos2 = posq[atoms.z];
+        float4 pos3 = posq[atoms.w];
+        float4 v12 = pos2-pos1;
+        float4 v13 = pos3-pos1;
+        pos.xyz = pos1.xyz + v12.xyz*weights.x + v13.xyz*weights.y + cross(v12, v13).xyz*weights.z;
+        posq[atoms.x] = pos;
+    }
+}
+
+/**
+ * Distribute forces from virtual sites to the atoms they are based on.
+ */
+__kernel void distributeForces(__global const float4* restrict posq, __global float4* restrict force,
+        __global const int4* restrict avg2Atoms, __global const float2* restrict avg2Weights,
+        __global const int4* restrict avg3Atoms, __global const float4* restrict avg3Weights,
+        __global const int4* restrict outOfPlaneAtoms, __global const float4* restrict outOfPlaneWeights) {
+    
+    // Two particle average sites.
+    
+    for (int index = get_global_id(0); index < NUM_2_AVERAGE; index += get_global_size(0)) {
+        int4 atoms = avg2Atoms[index];
+        float2 weights = avg2Weights[index];
+        float4 f = force[atoms.x];
+        float4 f1 = force[atoms.y];
+        float4 f2 = force[atoms.z];
+        f1.xyz += f.xyz*weights.x;
+        f2.xyz += f.xyz*weights.y;
+        force[atoms.y] = f1;
+        force[atoms.z] = f2;
+    }
+    
+    // Three particle average sites.
+    
+    for (int index = get_global_id(0); index < NUM_3_AVERAGE; index += get_global_size(0)) {
+        int4 atoms = avg3Atoms[index];
+        float4 weights = avg3Weights[index];
+        float4 f = force[atoms.x];
+        float4 f1 = force[atoms.y];
+        float4 f2 = force[atoms.z];
+        float4 f3 = force[atoms.w];
+        f1.xyz += f.xyz*weights.x;
+        f2.xyz += f.xyz*weights.y;
+        f3.xyz += f.xyz*weights.z;
+        force[atoms.y] = f1;
+        force[atoms.z] = f2;
+        force[atoms.w] = f3;
+    }
+    
+    // Out of plane sites.
+    
+    for (int index = get_global_id(0); index < NUM_OUT_OF_PLANE; index += get_global_size(0)) {
+        int4 atoms = outOfPlaneAtoms[index];
+        float4 weights = outOfPlaneWeights[index];
+        float4 pos1 = posq[atoms.y];
+        float4 pos2 = posq[atoms.z];
+        float4 pos3 = posq[atoms.w];
+        float4 v12 = pos2-pos1;
+        float4 v13 = pos3-pos1;
+        float4 f = force[atoms.x];
+        float4 f1 = force[atoms.y];
+        float4 f2 = force[atoms.z];
+        float4 f3 = force[atoms.w];
+        float4 fp2 = (float4) (weights.x*f.x - weights.z*v13.z*f.y + weights.z*v13.y*f[2],
+                   weights.z*v13.z*f.x + weights.x*f.y - weights.z*v13.x*f[2],
+                  -weights.z*v13.y*f.x + weights.z*v13.x*f.y + weights.x*f.z, 0.0f);
+        float4 fp3 = (float4) (weights.y*f.x + weights.z*v12[2]*f.y - weights.z*v12.y*f.z,
+                  -weights.z*v12[2]*f.x + weights.y*f.y + weights.z*v12.x*f.z,
+                   weights.z*v12.y*f.x - weights.z*v12.x*f.y + weights.y*f[2], 0.0f);
+        f1.xyz += f.xyz-fp2.xyz-fp3.xyz;
+        f2.xyz += fp2.xyz;
+        f3.xyz += fp3.xyz;
+        force[atoms.y] = f1;
+        force[atoms.z] = f2;
+        force[atoms.w] = f3;
+    }
+}

platforms/opencl/tests/TestOpenCLVirtualSites.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.                                     *
+ * -------------------------------------------------------------------------- */
+
+/**
+ * This tests the OpenCL implementation of virtual sites.
+ */
+
+#include "../../../tests/AssertionUtilities.h"
+#include "openmm/Context.h"
+#include "OpenCLPlatform.h"
+#include "openmm/CustomBondForce.h"
+#include "openmm/CustomExternalForce.h"
+#include "openmm/LangevinIntegrator.h"
+#include "openmm/NonbondedForce.h"
+#include "openmm/System.h"
+#include "openmm/VerletIntegrator.h"
+#include "openmm/VirtualSite.h"
+#include "sfmt/SFMT.h"
+#include <iostream>
+#include <vector>
+
+using namespace OpenMM;
+using namespace std;
+
+/**
+ * Check that massless particles are handled correctly.
+ */
+void testMasslessParticle() {
+    OpenCLPlatform platform;
+    System system;
+    system.addParticle(0.0);
+    system.addParticle(1.0);
+    CustomBondForce* bonds = new CustomBondForce("-1/r");
+    system.addForce(bonds);
+    vector<double> params;
+    bonds->addBond(0, 1, params);
+    VerletIntegrator integrator(0.002);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(2);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    context.setPositions(positions);
+    vector<Vec3> velocities(2);
+    velocities[0] = Vec3(0, 0, 0);
+    velocities[1] = Vec3(0, 1, 0);
+    context.setVelocities(velocities);
+    
+    // The second particle should move in a circular orbit around the first one.
+    // Compare it to the analytical solution.
+    
+    for (int i = 0; i < 1000; ++i) {
+        State state = context.getState(State::Positions | State::Velocities | State::Forces);
+        double time = state.getTime();
+        ASSERT_EQUAL_VEC(Vec3(), state.getPositions()[0], 0.0);
+        ASSERT_EQUAL_VEC(Vec3(), state.getVelocities()[0], 0.0);
+        ASSERT_EQUAL_VEC(Vec3(cos(time), sin(time), 0), state.getPositions()[1], 0.01);
+        ASSERT_EQUAL_VEC(Vec3(-sin(time), cos(time), 0), state.getVelocities()[1], 0.01);
+        integrator.step(1);
+    }
+}
+
+/**
+ * Test a TwoParticleAverageSite virtual site.
+ */
+void testTwoParticleAverage() {
+    OpenCLPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(0.0);
+    system.setVirtualSite(2, new TwoParticleAverageSite(0, 1, 0.8, 0.2));
+    CustomExternalForce* forceField = new CustomExternalForce("-a*x");
+    system.addForce(forceField);
+    forceField->addPerParticleParameter("a");
+    vector<double> params(1);
+    params[0] = 0.1;
+    forceField->addParticle(0, params);
+    params[0] = 0.2;
+    forceField->addParticle(1, params);
+    params[0] = 0.3;
+    forceField->addParticle(2, params);
+    LangevinIntegrator integrator(300.0, 0.1, 0.002);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    context.setPositions(positions);
+    context.applyConstraints(0.0001);
+    for (int i = 0; i < 1000; i++) {
+        State state = context.getState(State::Positions | State::Forces);
+        const vector<Vec3>& pos = state.getPositions();
+        ASSERT_EQUAL_VEC(pos[0]*0.8+pos[1]*0.2, pos[2], 1e-5);
+        ASSERT_EQUAL_VEC(Vec3(0.1+0.3*0.8, 0, 0), state.getForces()[0], 1e-4);
+        ASSERT_EQUAL_VEC(Vec3(0.2+0.3*0.2, 0, 0), state.getForces()[1], 1e-4);
+        integrator.step(1);
+    }
+}
+
+/**
+ * Test a ThreeParticleAverageSite virtual site.
+ */
+void testThreeParticleAverage() {
+    OpenCLPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(0.0);
+    system.setVirtualSite(3, new ThreeParticleAverageSite(0, 1, 2, 0.2, 0.3, 0.5));
+    CustomExternalForce* forceField = new CustomExternalForce("-a*x");
+    system.addForce(forceField);
+    forceField->addPerParticleParameter("a");
+    vector<double> params(1);
+    params[0] = 0.1;
+    forceField->addParticle(0, params);
+    params[0] = 0.2;
+    forceField->addParticle(1, params);
+    params[0] = 0.3;
+    forceField->addParticle(2, params);
+    params[0] = 0.4;
+    forceField->addParticle(3, params);
+    LangevinIntegrator integrator(300.0, 0.1, 0.002);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(4);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    positions[2] = Vec3(0, 1, 0);
+    context.setPositions(positions);
+    context.applyConstraints(0.0001);
+    for (int i = 0; i < 1000; i++) {
+        State state = context.getState(State::Positions | State::Forces);
+        const vector<Vec3>& pos = state.getPositions();
+        ASSERT_EQUAL_VEC(pos[0]*0.2+pos[1]*0.3+pos[2]*0.5, pos[3], 1e-5);
+        ASSERT_EQUAL_VEC(Vec3(0.1+0.4*0.2, 0, 0), state.getForces()[0], 1e-5);
+        ASSERT_EQUAL_VEC(Vec3(0.2+0.4*0.3, 0, 0), state.getForces()[1], 1e-5);
+        ASSERT_EQUAL_VEC(Vec3(0.3+0.4*0.5, 0, 0), state.getForces()[2], 1e-5);
+        integrator.step(1);
+    }
+}
+
+/**
+ * Test an OutOfPlaneSite virtual site.
+ */
+void testOutOfPlane() {
+    OpenCLPlatform platform;
+    System system;
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(0.0);
+    system.setVirtualSite(3, new OutOfPlaneSite(0, 1, 2, 0.3, 0.4, 0.5));
+    CustomExternalForce* forceField = new CustomExternalForce("-a*x");
+    system.addForce(forceField);
+    forceField->addPerParticleParameter("a");
+    vector<double> params(1);
+    params[0] = 0.1;
+    forceField->addParticle(0, params);
+    params[0] = 0.2;
+    forceField->addParticle(1, params);
+    params[0] = 0.3;
+    forceField->addParticle(2, params);
+    params[0] = 0.4;
+    forceField->addParticle(3, params);
+    LangevinIntegrator integrator(300.0, 0.1, 0.002);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(4);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(1, 0, 0);
+    positions[2] = Vec3(0, 1, 0);
+    context.setPositions(positions);
+    context.applyConstraints(0.0001);
+    for (int i = 0; i < 1000; i++) {
+        State state = context.getState(State::Positions | State::Forces);
+        const vector<Vec3>& pos = state.getPositions();
+        Vec3 v12 = pos[1]-pos[0];
+        Vec3 v13 = pos[2]-pos[0];
+        Vec3 cross = v12.cross(v13);
+        ASSERT_EQUAL_VEC(pos[0]+v12*0.3+v13*0.4+cross*0.5, pos[3], 1e-5);
+        const vector<Vec3>& f = state.getForces();
+        ASSERT_EQUAL_VEC(Vec3(0.1+0.2+0.3+0.4, 0, 0), f[0]+f[1]+f[2], 1e-5);
+        Vec3 f2(0.4*0.3, 0.4*0.5*v13[2], -0.4*0.5*v13[1]);
+        Vec3 f3(0.4*0.4, -0.4*0.5*v12[2], 0.4*0.5*v12[1]);
+        ASSERT_EQUAL_VEC(Vec3(0.1+0.4, 0, 0)-f2-f3, f[0], 1e-5);
+        ASSERT_EQUAL_VEC(Vec3(0.2, 0, 0)+f2, f[1], 1e-5);
+        ASSERT_EQUAL_VEC(Vec3(0.3, 0, 0)+f3, f[2], 1e-5);
+        integrator.step(1);
+    }
+}
+
+/**
+ * Make sure that energy, linear momentum, and angular momentum are all conserved
+ * when using virtual sites.
+ */
+void testConservationLaws() {
+    OpenCLPlatform platform;
+    System system;
+    NonbondedForce* forceField = new NonbondedForce();
+    system.addForce(forceField);
+    vector<Vec3> positions;
+    
+    // Create a linear molecule with a TwoParticleAverage virtual site.
+    
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(0.0);
+    system.setVirtualSite(2, new TwoParticleAverageSite(0, 1, 0.4, 0.6));
+    system.addConstraint(0, 1, 2.0);
+    for (int i = 0; i < 3; i++) {
+        forceField->addParticle(0, 1, 10);
+        for (int j = 0; j < i; j++)
+            forceField->addException(i, j, 0, 1, 0);
+    }
+    positions.push_back(Vec3(0, 0, 0));
+    positions.push_back(Vec3(2, 0, 0));
+    positions.push_back(Vec3());
+    
+    // Create a planar molecule with a ThreeParticleAverage virtual site.
+    
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(0.0);
+    system.setVirtualSite(6, new ThreeParticleAverageSite(3, 4, 5, 0.3, 0.5, 0.2));
+    system.addConstraint(3, 4, 1.0);
+    system.addConstraint(3, 5, 1.0);
+    system.addConstraint(4, 5, sqrt(2.0));
+    for (int i = 0; i < 4; i++) {
+        forceField->addParticle(0, 1, 10);
+        for (int j = 0; j < i; j++)
+            forceField->addException(i+3, j+3, 0, 1, 0);
+    }
+    positions.push_back(Vec3(0, 0, 1));
+    positions.push_back(Vec3(1, 0, 1));
+    positions.push_back(Vec3(0, 1, 1));
+    positions.push_back(Vec3());
+    
+    // Create a tetrahedral molecule with an OutOfPlane virtual site.
+    
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(1.0);
+    system.addParticle(0.0);
+    system.setVirtualSite(10, new OutOfPlaneSite(7, 8, 9, 0.3, 0.5, 0.2));
+    system.addConstraint(7, 8, 1.0);
+    system.addConstraint(7, 9, 1.0);
+    system.addConstraint(8, 9, sqrt(2.0));
+    for (int i = 0; i < 4; i++) {
+        forceField->addParticle(0, 1, 10);
+        for (int j = 0; j < i; j++)
+            forceField->addException(i+7, j+7, 0, 1, 0);
+    }
+    positions.push_back(Vec3(1, 0, -1));
+    positions.push_back(Vec3(2, 0, -1));
+    positions.push_back(Vec3(1, 1, -1));
+    positions.push_back(Vec3());
+
+    // Simulate it and check conservation laws.
+    
+    VerletIntegrator integrator(0.002);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.applyConstraints(0.0001);
+    int numParticles = system.getNumParticles();
+    double initialEnergy;
+    Vec3 initialMomentum, initialAngularMomentum;
+    for (int i = 0; i < 1000; i++) {
+        State state = context.getState(State::Positions | State::Velocities | State::Forces | State::Energy);
+        const vector<Vec3>& pos = state.getPositions();
+        const vector<Vec3>& vel = state.getVelocities();
+        const vector<Vec3>& f = state.getForces();
+        double energy = state.getPotentialEnergy();
+        for (int j = 0; j < numParticles; j++) {
+            Vec3 v = vel[j] + f[j]*0.5*integrator.getStepSize();
+            energy += 0.5*system.getParticleMass(j)*v.dot(v);
+        }
+        if (i == 0)
+            initialEnergy = energy;
+        else
+            ASSERT_EQUAL_TOL(initialEnergy, energy, 0.01);
+        Vec3 momentum;
+        for (int j = 0; j < numParticles; j++)
+            momentum += vel[j]*system.getParticleMass(j);
+        if (i == 0)
+            initialMomentum = momentum;
+        else
+            ASSERT_EQUAL_VEC(initialMomentum, momentum, 0.01);
+        Vec3 angularMomentum;
+        for (int j = 0; j < numParticles; j++)
+            angularMomentum += pos[j].cross(vel[j])*system.getParticleMass(j);
+        if (i == 0)
+            initialAngularMomentum = angularMomentum;
+        else
+            ASSERT_EQUAL_VEC(initialAngularMomentum, angularMomentum, 0.01);
+        integrator.step(1);
+    }
+}
+
+/**
+ * Make sure that atom reordering respects virtual sites.
+ */
+void testReordering() {
+    const double cutoff = 2.0;
+    const double boxSize = 20.0;
+    OpenCLPlatform platform;
+    System system;
+    NonbondedForce* nonbonded = new NonbondedForce();
+    system.addForce(nonbonded);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffNonPeriodic);
+    nonbonded->setCutoffDistance(cutoff);
+    vector<Vec3> positions;
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    
+    // Create linear molecules with TwoParticleAverage virtual sites.
+    
+    for (int i = 0; i < 50; i++) {
+        int start = system.getNumParticles();
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        system.addParticle(0.0);
+        system.setVirtualSite(start+2, new TwoParticleAverageSite(start, start+1, 0.4, 0.6));
+        system.addConstraint(start, start+1, 2.0);
+        for (int i = 0; i < 3; i++) {
+            nonbonded->addParticle(0, 0.2, 1);
+            for (int j = 0; j < i; j++)
+                nonbonded->addException(start+i, start+j, 0, 1, 0);
+        }
+        Vec3 pos(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
+        positions.push_back(pos);
+        positions.push_back(pos+Vec3(2, 0, 0));
+        positions.push_back(Vec3());
+    }
+    
+    // Create planar molecules with ThreeParticleAverage virtual sites.
+    
+    for (int i = 0; i < 50; i++) {
+        int start = system.getNumParticles();
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        system.addParticle(0.0);
+        system.setVirtualSite(start+3, new ThreeParticleAverageSite(start, start+1, start+2, 0.3, 0.5, 0.2));
+        system.addConstraint(start, start+1, 1.0);
+        system.addConstraint(start, start+2, 1.0);
+        system.addConstraint(start+1, start+2, sqrt(2.0));
+        for (int i = 0; i < 4; i++) {
+            nonbonded->addParticle(0, 0.2, 1);
+            for (int j = 0; j < i; j++)
+                nonbonded->addException(start+i, start+j, 0, 1, 0);
+        }
+        Vec3 pos(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
+        positions.push_back(pos);
+        positions.push_back(pos+Vec3(1, 0, 0));
+        positions.push_back(pos+Vec3(0, 1, 0));
+        positions.push_back(Vec3());
+    }
+    
+    // Create tetrahedral molecules with OutOfPlane virtual sites.
+    
+    for (int i = 0; i < 50; i++) {
+        int start = system.getNumParticles();
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        system.addParticle(0.0);
+        system.setVirtualSite(start+3, new OutOfPlaneSite(start, start+1, start+2, 0.3, 0.5, 0.2));
+        system.addConstraint(start, start+1, 1.0);
+        system.addConstraint(start, start+2, 1.0);
+        system.addConstraint(start+1, start+2, sqrt(2.0));
+        for (int i = 0; i < 4; i++) {
+            nonbonded->addParticle(0, 0.2, 1);
+            for (int j = 0; j < i; j++)
+                nonbonded->addException(start+i, start+j, 0, 1, 0);
+        }
+        Vec3 pos(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
+        positions.push_back(pos);
+        positions.push_back(pos+Vec3(1, 0, 0));
+        positions.push_back(pos+Vec3(0, 1, 0));
+        positions.push_back(Vec3());
+    }
+
+    // Simulate it and check conservation laws.
+    
+    LangevinIntegrator integrator(300.0, 0.1, 0.002);
+    Context context(system, integrator, platform);
+    context.setPositions(positions);
+    context.applyConstraints(0.0001);
+    for (int i = 0; i < 1000; i++) {
+        State state = context.getState(State::Positions);
+        const vector<Vec3>& pos = state.getPositions();
+        for (int j = 0; j < 150; j += 3)
+            ASSERT_EQUAL_VEC(pos[j]*0.4+pos[j+1]*0.6, pos[j+2], 1e-5);
+        for (int j = 150; j < 350; j += 4)
+            ASSERT_EQUAL_VEC(pos[j]*0.3+pos[j+1]*0.5+pos[j+2]*0.2, pos[j+3], 1e-5);
+        for (int j = 350; j < 550; j += 4) {
+            Vec3 v12 = pos[j+1]-pos[j];
+            Vec3 v13 = pos[j+2]-pos[j];
+            Vec3 cross = v12.cross(v13);
+            ASSERT_EQUAL_VEC(pos[j]+v12*0.3+v13*0.5+cross*0.2, pos[j+3], 1e-5);
+        }
+        integrator.step(1);
+    }
+}
+
+int main() {
+    try {
+        testMasslessParticle();
+        testTwoParticleAverage();
+        testThreeParticleAverage();
+        testOutOfPlane();
+        testConservationLaws();
+        testReordering();
+    }
+    catch(const exception& e) {
+        cout << "exception: " << e.what() << endl;
+        return 1;
+    }
+    cout << "Done" << endl;
+    return 0;
+}

platforms/reference/tests/TestReferenceVirtualSites.cpp

@@ -120,6 +120,7 @@ void testTwoParticleAverage() {
         ASSERT_EQUAL_VEC(pos[0]*0.8+pos[1]*0.2, pos[2], 1e-10);
         ASSERT_EQUAL_VEC(Vec3(0.1+0.3*0.8, 0, 0), state.getForces()[0], 1e-10);
         ASSERT_EQUAL_VEC(Vec3(0.2+0.3*0.2, 0, 0), state.getForces()[1], 1e-10);
+        integrator.step(1);
     }
 }
 
@@ -161,6 +162,7 @@ void testThreeParticleAverage() {
         ASSERT_EQUAL_VEC(Vec3(0.1+0.4*0.2, 0, 0), state.getForces()[0], 1e-10);
         ASSERT_EQUAL_VEC(Vec3(0.2+0.4*0.3, 0, 0), state.getForces()[1], 1e-10);
         ASSERT_EQUAL_VEC(Vec3(0.3+0.4*0.5, 0, 0), state.getForces()[2], 1e-10);
+        integrator.step(1);
     }
 }
 
@@ -209,6 +211,7 @@ void testOutOfPlane() {
         ASSERT_EQUAL_VEC(Vec3(0.1+0.4, 0, 0)-f2-f3, f[0], 1e-10);
         ASSERT_EQUAL_VEC(Vec3(0.2, 0, 0)+f2, f[1], 1e-10);
         ASSERT_EQUAL_VEC(Vec3(0.3, 0, 0)+f3, f[2], 1e-10);
+        integrator.step(1);
     }
 }