RPMD supports virtual sites

plugins/rpmd/openmmapi/src/RPMDIntegrator.cpp

@@ -54,12 +54,8 @@ RPMDIntegrator::RPMDIntegrator(int numCopies, double temperature, double frictio
 void RPMDIntegrator::initialize(ContextImpl& contextRef) {
     if (owner != NULL && &contextRef.getOwner() != owner)
         throw OpenMMException("This Integrator is already bound to a context");
-    const System& system = contextRef.getSystem();
-    if (system.getNumConstraints() > 0)
+    if (contextRef.getSystem().getNumConstraints() > 0)
         throw OpenMMException("RPMDIntegrator cannot be used with Systems that include constraints");
-    for (int i = 0; i < system.getNumParticles(); i++)
-        if (system.isVirtualSite(i))
-            throw OpenMMException("RPMDIntegrator cannot be used with Systems that include virtual sites");
     context = &contextRef;
     owner = &contextRef.getOwner();
     kernel = context->getPlatform().createKernel(IntegrateRPMDStepKernel::Name(), contextRef);

plugins/rpmd/platforms/cuda/src/CudaRpmdKernels.cpp

@@ -187,6 +187,7 @@ void CudaIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
         void* copyToContextArgs[] = {&velocities->getDevicePointer(), &cu.getVelm().getDevicePointer(), &positions->getDevicePointer(),
                 &cu.getPosq().getDevicePointer(), &cu.getAtomIndexArray().getDevicePointer(), &i};
         cu.executeKernel(copyToContextKernel, copyToContextArgs, cu.getNumAtoms());
+        context.computeVirtualSites();
         context.updateContextState();
         context.calcForcesAndEnergy(true, false);
         if (cu.getAtomsWereReordered() && cu.getNonbondedUtilities().getUsePeriodic()) {

plugins/rpmd/platforms/cuda/src/kernels/rpmd.cu

@@ -38,6 +38,8 @@ extern "C" __global__ void applyPileThermostat(mixed4* velm, float4* random, uns
     for (int particle = (blockIdx.x*blockDim.x+threadIdx.x)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) {
         mixed4 particleVelm = velm[particle+indexInBlock*PADDED_NUM_ATOMS];
         mixed invMass = particleVelm.w;
+        if (invMass == 0)
+            continue;
         mixed c3_0 = c2_0*SQRT(nkT*invMass);
         
         // Forward FFT.
@@ -100,6 +102,8 @@ extern "C" __global__ void integrateStep(mixed4* posq, mixed4* velm, long long*
         int index = particle+indexInBlock*PADDED_NUM_ATOMS;
         int forceIndex = particle+indexInBlock*PADDED_NUM_ATOMS*3;
         mixed4 particleVelm = velm[index];
+        if (particleVelm.w == 0)
+            continue;
         particleVelm.x += forceScale*force[forceIndex]*(0.5f*dt*particleVelm.w);
         particleVelm.y += forceScale*force[forceIndex+PADDED_NUM_ATOMS]*(0.5f*dt*particleVelm.w);
         particleVelm.z += forceScale*force[forceIndex+PADDED_NUM_ATOMS*2]*(0.5f*dt*particleVelm.w);
@@ -118,6 +122,8 @@ extern "C" __global__ void integrateStep(mixed4* posq, mixed4* velm, long long*
     for (int particle = (blockIdx.x*blockDim.x+threadIdx.x)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) {
         mixed4 particlePosq = posq[particle+indexInBlock*PADDED_NUM_ATOMS];
         mixed4 particleVelm = velm[particle+indexInBlock*PADDED_NUM_ATOMS];
+        if (particleVelm.w == 0)
+            continue;
         
         // Forward FFT.
         
@@ -173,6 +179,8 @@ extern "C" __global__ void advanceVelocities(mixed4* velm, long long* force, mix
         int index = particle+indexInBlock*PADDED_NUM_ATOMS;
         int forceIndex = particle+indexInBlock*PADDED_NUM_ATOMS*3;
         mixed4 particleVelm = velm[index];
+        if (particleVelm.w == 0)
+            continue;
         particleVelm.x += forceScale*force[forceIndex]*(0.5f*dt*particleVelm.w);
         particleVelm.y += forceScale*force[forceIndex+PADDED_NUM_ATOMS]*(0.5f*dt*particleVelm.w);
         particleVelm.z += forceScale*force[forceIndex+PADDED_NUM_ATOMS*2]*(0.5f*dt*particleVelm.w);

plugins/rpmd/platforms/cuda/tests/TestCudaRpmd.cpp

@@ -38,9 +38,11 @@
 #include "openmm/Context.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/HarmonicBondForce.h"
+#include "openmm/NonbondedForce.h"
 #include "openmm/Platform.h"
 #include "openmm/System.h"
 #include "openmm/RPMDIntegrator.h"
+#include "openmm/VirtualSite.h"
 #include "SimTKUtilities/SimTKOpenMMUtilities.h"
 #include "sfmt/SFMT.h"
 #include <iostream>
@@ -272,6 +274,87 @@ void testCMMotionRemoval() {
     }
 }
 
+void testVirtualSites() {
+    const int gridSize = 3;
+    const int numMolecules = gridSize*gridSize*gridSize;
+    const int numParticles = numMolecules*3;
+    const int numCopies = 10;
+    const double spacing = 2.0;
+    const double cutoff = 3.0;
+    const double boxSize = spacing*(gridSize+1);
+    const double temperature = 300.0;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    system.addForce(bonds);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    system.addForce(nonbonded);
+
+    // Create a cloud of molecules.
+
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numMolecules; i++) {
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        system.addParticle(0.0);
+        nonbonded->addParticle(-0.2, 0.2, 0.2);
+        nonbonded->addParticle(0.1, 0.2, 0.2);
+        nonbonded->addParticle(0.1, 0.2, 0.2);
+        nonbonded->addException(3*i, 3*i+1, 0, 1, 0);
+        nonbonded->addException(3*i, 3*i+2, 0, 1, 0);
+        nonbonded->addException(3*i+1, 3*i+2, 0, 1, 0);
+        bonds->addBond(3*i, 3*i+1, 1.0, 10000.0);
+        system.setVirtualSite(3*i+2, new TwoParticleAverageSite(3*i, 3*i+1, 0.5, 0.5));
+    }
+    RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
+    Platform& platform = Platform::getPlatformByName("CUDA");
+    Context context(system, integ, platform);
+    for (int copy = 0; copy < numCopies; copy++) {
+        for (int i = 0; i < gridSize; i++)
+            for (int j = 0; j < gridSize; j++)
+                for (int k = 0; k < gridSize; k++) {
+                    Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
+                    int index = k+gridSize*(j+gridSize*i);
+                    positions[3*index] = pos;
+                    positions[3*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
+                    positions[3*index+2] = Vec3();
+                }
+        integ.setPositions(copy, positions);
+    }
+
+    // Check the temperature and virtual site locations.
+    
+    const int numSteps = 1000;
+    integ.step(1000);
+    vector<double> ke(numCopies, 0.0);
+    for (int i = 0; i < numSteps; i++) {
+        integ.step(1);
+        vector<State> state(numCopies);
+        for (int j = 0; j < numCopies; j++) {
+            state[j] = integ.getState(j, State::Positions | State::Velocities | State::Forces, true);
+            const vector<Vec3>& pos = state[j].getPositions();
+            for (int k = 0; k < numMolecules; k++)
+                ASSERT_EQUAL_VEC((pos[3*k]+pos[3*k+1])*0.5, pos[3*k+2], 1e-5);
+        }
+        for (int j = 0; j < numParticles; j++) {
+            for (int k = 0; k < numCopies; k++) {
+                Vec3 v = state[k].getVelocities()[j];
+                ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
+            }
+        }
+    }
+    double meanKE = 0.0;
+    for (int i = 0; i < numCopies; i++)
+        meanKE += ke[i];
+    meanKE /= numSteps*numCopies;
+    double expectedKE = 0.5*numCopies*(2*numMolecules)*3*BOLTZ*temperature;
+    ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
+}
+
 int main(int argc, char* argv[]) {
     try {
         registerRPMDCudaKernelFactories();
@@ -280,6 +363,7 @@ int main(int argc, char* argv[]) {
         testFreeParticles();
         testParaHydrogen();
         testCMMotionRemoval();
+        testVirtualSites();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

plugins/rpmd/platforms/opencl/src/OpenCLRpmdKernels.cpp

@@ -198,6 +198,7 @@ void OpenCLIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
     for (int i = 0; i < numCopies; i++) {
         copyToContextKernel.setArg<cl_int>(5, i);
         cl.executeKernel(copyToContextKernel, cl.getNumAtoms());
+        context.computeVirtualSites();
         context.updateContextState();
         context.calcForcesAndEnergy(true, false);
         if (cl.getAtomsWereReordered() && cl.getNonbondedUtilities().getUsePeriodic()) {

plugins/rpmd/platforms/opencl/src/kernels/rpmd.cl

@@ -38,6 +38,8 @@ __kernel void applyPileThermostat(__global mixed4* velm, __global float4* random
     for (int particle = get_global_id(0)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) {
         mixed4 particleVelm = velm[particle+indexInBlock*PADDED_NUM_ATOMS];
         mixed invMass = particleVelm.w;
+        if (invMass == 0)
+            continue;
         mixed c3_0 = c2_0*sqrt(nkT*invMass);
         
         // Forward FFT.
@@ -97,6 +99,8 @@ __kernel void integrateStep(__global mixed4* posq, __global mixed4* velm, __glob
     for (int particle = get_global_id(0)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) {
         int index = particle+indexInBlock*PADDED_NUM_ATOMS;
         mixed4 particleVelm = velm[index];
+        if (particleVelm.w == 0)
+            continue;
         particleVelm.xyz += convert_mixed4(force[index]).xyz*(0.5f*dt*particleVelm.w);
         velm[index] = particleVelm;
     }
@@ -113,6 +117,8 @@ __kernel void integrateStep(__global mixed4* posq, __global mixed4* velm, __glob
     for (int particle = get_global_id(0)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) {
         mixed4 particlePosq = posq[particle+indexInBlock*PADDED_NUM_ATOMS];
         mixed4 particleVelm = velm[particle+indexInBlock*PADDED_NUM_ATOMS];
+        if (particleVelm.w == 0)
+            continue;
         
         // Forward FFT.
         
@@ -166,6 +172,8 @@ __kernel void advanceVelocities(__global mixed4* velm, __global real4* force, mi
     for (int particle = get_global_id(0)/NUM_COPIES; particle < NUM_ATOMS; particle += numBlocks) {
         int index = particle+indexInBlock*PADDED_NUM_ATOMS;
         mixed4 particleVelm = velm[index];
+        if (particleVelm.w == 0)
+            continue;
         particleVelm.xyz += convert_mixed4(force[index]).xyz*(0.5f*dt*particleVelm.w);
         velm[index] = particleVelm;
     }

plugins/rpmd/platforms/opencl/tests/TestOpenCLRpmd.cpp

@@ -38,9 +38,11 @@
 #include "openmm/Context.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/HarmonicBondForce.h"
+#include "openmm/NonbondedForce.h"
 #include "openmm/Platform.h"
 #include "openmm/System.h"
 #include "openmm/RPMDIntegrator.h"
+#include "openmm/VirtualSite.h"
 #include "SimTKUtilities/SimTKOpenMMUtilities.h"
 #include "sfmt/SFMT.h"
 #include <iostream>
@@ -273,6 +275,87 @@ void testCMMotionRemoval() {
     }
 }
 
+void testVirtualSites() {
+    const int gridSize = 3;
+    const int numMolecules = gridSize*gridSize*gridSize;
+    const int numParticles = numMolecules*3;
+    const int numCopies = 10;
+    const double spacing = 2.0;
+    const double cutoff = 3.0;
+    const double boxSize = spacing*(gridSize+1);
+    const double temperature = 300.0;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    system.addForce(bonds);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    system.addForce(nonbonded);
+
+    // Create a cloud of molecules.
+
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numMolecules; i++) {
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        system.addParticle(0.0);
+        nonbonded->addParticle(-0.2, 0.2, 0.2);
+        nonbonded->addParticle(0.1, 0.2, 0.2);
+        nonbonded->addParticle(0.1, 0.2, 0.2);
+        nonbonded->addException(3*i, 3*i+1, 0, 1, 0);
+        nonbonded->addException(3*i, 3*i+2, 0, 1, 0);
+        nonbonded->addException(3*i+1, 3*i+2, 0, 1, 0);
+        bonds->addBond(3*i, 3*i+1, 1.0, 10000.0);
+        system.setVirtualSite(3*i+2, new TwoParticleAverageSite(3*i, 3*i+1, 0.5, 0.5));
+    }
+    RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
+    Platform& platform = Platform::getPlatformByName("OpenCL");
+    Context context(system, integ, platform);
+    for (int copy = 0; copy < numCopies; copy++) {
+        for (int i = 0; i < gridSize; i++)
+            for (int j = 0; j < gridSize; j++)
+                for (int k = 0; k < gridSize; k++) {
+                    Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
+                    int index = k+gridSize*(j+gridSize*i);
+                    positions[3*index] = pos;
+                    positions[3*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
+                    positions[3*index+2] = Vec3();
+                }
+        integ.setPositions(copy, positions);
+    }
+
+    // Check the temperature and virtual site locations.
+    
+    const int numSteps = 1000;
+    integ.step(1000);
+    vector<double> ke(numCopies, 0.0);
+    for (int i = 0; i < numSteps; i++) {
+        integ.step(1);
+        vector<State> state(numCopies);
+        for (int j = 0; j < numCopies; j++) {
+            state[j] = integ.getState(j, State::Positions | State::Velocities | State::Forces, true);
+            const vector<Vec3>& pos = state[j].getPositions();
+            for (int k = 0; k < numMolecules; k++)
+                ASSERT_EQUAL_VEC((pos[3*k]+pos[3*k+1])*0.5, pos[3*k+2], 1e-5);
+        }
+        for (int j = 0; j < numParticles; j++) {
+            for (int k = 0; k < numCopies; k++) {
+                Vec3 v = state[k].getVelocities()[j];
+                ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
+            }
+        }
+    }
+    double meanKE = 0.0;
+    for (int i = 0; i < numCopies; i++)
+        meanKE += ke[i];
+    meanKE /= numSteps*numCopies;
+    double expectedKE = 0.5*numCopies*(2*numMolecules)*3*BOLTZ*temperature;
+    ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
+}
+
 int main(int argc, char* argv[]) {
     try {
         registerRPMDOpenCLKernelFactories();
@@ -281,6 +364,7 @@ int main(int argc, char* argv[]) {
         testFreeParticles();
         testParaHydrogen();
         testCMMotionRemoval();
+        testVirtualSites();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

plugins/rpmd/platforms/reference/src/ReferenceRpmdKernels.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) 2011 Stanford University and the Authors.           *
+ * Portions copyright (c) 2011-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -86,6 +86,7 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
         for (int i = 0; i < numCopies; i++) {
             pos = positions[i];
             vel = velocities[i];
+            context.computeVirtualSites();
             context.calcForcesAndEnergy(true, false);
             forces[i] = f;
         }
@@ -102,6 +103,8 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
     const RealOpenMM c1_0 = exp(-halfdt*integrator.getFriction());
     const RealOpenMM c2_0 = sqrt(1.0-c1_0*c1_0);
     for (int particle = 0; particle < numParticles; particle++) {
+        if (system.getParticleMass(particle) == 0.0)
+            continue;
         const RealOpenMM c3_0 = c2_0*sqrt(nkT/system.getParticleMass(particle));
         for (int component = 0; component < 3; component++) {
             for (int k = 0; k < numCopies; k++)
@@ -136,11 +139,14 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
     
     for (int i = 0; i < numCopies; i++)
         for (int j = 0; j < numParticles; j++)
-            velocities[i][j] += forces[i][j]*(halfdt/system.getParticleMass(j));
+            if (system.getParticleMass(j) != 0.0)
+                velocities[i][j] += forces[i][j]*(halfdt/system.getParticleMass(j));
     
     // Evolve the free ring polymer by transforming to the frequency domain.
 
     for (int particle = 0; particle < numParticles; particle++) {
+        if (system.getParticleMass(particle) == 0.0)
+            continue;
         for (int component = 0; component < 3; component++) {
             for (int k = 0; k < numCopies; k++) {
                 q[k] = t_complex(scale*positions[k][particle][component], 0.0);
@@ -173,6 +179,7 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
     for (int i = 0; i < numCopies; i++) {
         pos = positions[i];
         vel = velocities[i];
+        context.computeVirtualSites();
         context.updateContextState();
         positions[i] = pos;
         velocities[i] = vel;
@@ -184,11 +191,14 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
     
     for (int i = 0; i < numCopies; i++)
         for (int j = 0; j < numParticles; j++)
-            velocities[i][j] += forces[i][j]*(halfdt/system.getParticleMass(j));
+            if (system.getParticleMass(j) != 0.0)
+                velocities[i][j] += forces[i][j]*(halfdt/system.getParticleMass(j));
 
     // Apply the PILE-L thermostat again.
     
     for (int particle = 0; particle < numParticles; particle++) {
+        if (system.getParticleMass(particle) == 0.0)
+            continue;
         const RealOpenMM c3_0 = c2_0*sqrt(nkT/system.getParticleMass(particle));
         for (int component = 0; component < 3; component++) {
             for (int k = 0; k < numCopies; k++)
@@ -224,101 +234,6 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
     context.setTime(context.getTime()+dt);
 }
 
-//void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDIntegrator& integrator, bool forcesAreValid) {
-//    int numCopies = positions.size();
-//    int numParticles = positions[0].size();
-//    vector<RealVec>& pos = extractPositions(context);
-//    vector<RealVec>& vel = extractVelocities(context);
-//    vector<RealVec>& f = extractForces(context);
-//    
-//    // Loop over copies and compute the force on each one.
-//    
-//    for (int i = 0; i < numCopies; i++) {
-//        pos = positions[i];
-//        vel = velocities[i];
-//        context.calcForcesAndEnergy(true, false);
-//        forces[i] = f;
-//    }
-//
-//    // Update velocities.
-//    
-//    const RealOpenMM dt = integrator.getStepSize();
-//    const System& system = context.getSystem();
-//    for (int i = 0; i < numCopies; i++)
-//        for (int j = 0; j < numParticles; j++)
-//            velocities[i][j] += forces[i][j]*(dt/system.getParticleMass(j));
-//    
-//    // Evolve the free ring polymer by transforming to the frequency domain.
-//
-//    vector<t_complex> p(numCopies);
-//    vector<t_complex> q(numCopies);
-//    const RealOpenMM scale = 1.0/sqrt(numCopies);
-//    const RealOpenMM hbar = 1.054571628e-34*AVOGADRO/(1000*1e-12);
-//    const RealOpenMM nkT = numCopies*BOLTZ*integrator.getTemperature();
-//    const RealOpenMM twown = 2.0*nkT/hbar;
-//    for (int particle = 0; particle < numParticles; particle++) {
-//        for (int component = 0; component < 3; component++) {
-//            for (int k = 0; k < numCopies; k++) {
-//                q[k] = t_complex(scale*positions[k][particle][component], 0.0);
-//                p[k] = t_complex(scale*velocities[k][particle][component]*system.getParticleMass(particle), 0.0);
-//            }
-//            fftpack_exec_1d(fft, FFTPACK_FORWARD, &q[0], &q[0]);
-//            fftpack_exec_1d(fft, FFTPACK_FORWARD, &p[0], &p[0]);
-//            q[0] += p[0]*(dt/system.getParticleMass(particle));
-//            for (int k = 1; k < numCopies; k++) {
-//                const RealOpenMM wk = twown*sin(k*M_PI/numCopies);
-//                const RealOpenMM wt = wk*dt;
-//                const RealOpenMM sinwt2 = sin(wt/2);
-//                const RealOpenMM wm = wk*system.getParticleMass(particle);
-//                const t_complex pprime = p[k] - q[k]*(2.0*wm*sinwt2); // Advance momentum from t-dt/2 to t+dt/2
-//                q[k] = p[k]*(2.0*sinwt2/wm) + q[k]*(1.0-4.0*sinwt2*sinwt2); // Advance position from t to t+dt
-//                p[k] = pprime;
-//            }
-//            fftpack_exec_1d(fft, FFTPACK_BACKWARD, &q[0], &q[0]);
-//            fftpack_exec_1d(fft, FFTPACK_BACKWARD, &p[0], &p[0]);
-//            for (int k = 0; k < numCopies; k++) {
-//                positions[k][particle][component] = scale*q[k].re;
-//                velocities[k][particle][component] = scale*p[k].re/system.getParticleMass(particle);
-//            }
-//        }
-//    }
-//
-//    // Apply the PILE-L thermostat.
-//    
-//    const RealOpenMM c1_0 = exp(-dt*integrator.getFriction());
-//    const RealOpenMM c2_0 = sqrt(1.0-c1_0*c1_0);
-//    for (int particle = 0; particle < numParticles; particle++) {
-//        const RealOpenMM c3_0 = c2_0*sqrt(system.getParticleMass(particle)*nkT);
-//        for (int component = 0; component < 3; component++) {
-//            for (int k = 0; k < numCopies; k++)
-//                p[k] = t_complex(scale*velocities[k][particle][component]*system.getParticleMass(particle), 0.0);
-//            fftpack_exec_1d(fft, FFTPACK_FORWARD, &p[0], &p[0]);
-//            
-//            // Apply a local Langevin thermostat to the centroid mode.
-//
-//            p[0].re = p[0].re*c1_0 + c3_0*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
-//
-//            // Use critical damping white noise for the remaining modes.
-//            
-//            for (int k = 1; k <= numCopies/2; k++) {
-//                const bool isCenter = (numCopies%2 == 0 && k == numCopies/2);
-//                const RealOpenMM wk = twown*sin(k*M_PI/numCopies);
-//                const RealOpenMM c1 = exp(-2.0*wk*dt);
-//                const RealOpenMM c2 = sqrt((1.0-c1*c1)/2) * (isCenter ? sqrt(2.0) : 1.0);
-//                const RealOpenMM c3 = c2*sqrt(system.getParticleMass(particle)*nkT);
-//                RealOpenMM rand1 = c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber();
-//                RealOpenMM rand2 = (isCenter ? 0.0 : c3*SimTKOpenMMUtilities::getNormallyDistributedRandomNumber());
-//                p[k] = p[k]*c1 + t_complex(rand1, rand2);
-//                if (k < numCopies-k)
-//                    p[numCopies-k] = p[numCopies-k]*c1 + t_complex(rand1, -rand2);                
-//            }
-//            fftpack_exec_1d(fft, FFTPACK_BACKWARD, &p[0], &p[0]);
-//            for (int k = 0; k < numCopies; k++)
-//                velocities[k][particle][component] = scale*p[k].re/system.getParticleMass(particle);
-//        }
-//    }
-//}
-
 double ReferenceIntegrateRPMDStepKernel::computeKineticEnergy(ContextImpl& context, const RPMDIntegrator& integrator) {
     System& system = context.getSystem();
     int numParticles = system.getNumParticles();

plugins/rpmd/platforms/reference/tests/TestReferenceRpmd.cpp

@@ -37,9 +37,11 @@
 #include "openmm/CMMotionRemover.h"
 #include "openmm/Context.h"
 #include "openmm/HarmonicBondForce.h"
+#include "openmm/NonbondedForce.h"
 #include "openmm/Platform.h"
 #include "openmm/System.h"
 #include "openmm/RPMDIntegrator.h"
+#include "openmm/VirtualSite.h"
 #include "SimTKUtilities/SimTKOpenMMUtilities.h"
 #include "sfmt/SFMT.h"
 #include <iostream>
@@ -154,10 +156,92 @@ void testCMMotionRemoval() {
     }
 }
 
+void testVirtualSites() {
+    const int gridSize = 3;
+    const int numMolecules = gridSize*gridSize*gridSize;
+    const int numParticles = numMolecules*3;
+    const int numCopies = 10;
+    const double spacing = 2.0;
+    const double cutoff = 3.0;
+    const double boxSize = spacing*(gridSize+1);
+    const double temperature = 300.0;
+    System system;
+    system.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
+    HarmonicBondForce* bonds = new HarmonicBondForce();
+    system.addForce(bonds);
+    NonbondedForce* nonbonded = new NonbondedForce();
+    nonbonded->setCutoffDistance(cutoff);
+    nonbonded->setNonbondedMethod(NonbondedForce::CutoffPeriodic);
+    system.addForce(nonbonded);
+
+    // Create a cloud of molecules.
+
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numMolecules; i++) {
+        system.addParticle(1.0);
+        system.addParticle(1.0);
+        system.addParticle(0.0);
+        nonbonded->addParticle(-0.2, 0.2, 0.2);
+        nonbonded->addParticle(0.1, 0.2, 0.2);
+        nonbonded->addParticle(0.1, 0.2, 0.2);
+        nonbonded->addException(3*i, 3*i+1, 0, 1, 0);
+        nonbonded->addException(3*i, 3*i+2, 0, 1, 0);
+        nonbonded->addException(3*i+1, 3*i+2, 0, 1, 0);
+        bonds->addBond(3*i, 3*i+1, 1.0, 10000.0);
+        system.setVirtualSite(3*i+2, new TwoParticleAverageSite(3*i, 3*i+1, 0.5, 0.5));
+    }
+    RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
+    Platform& platform = Platform::getPlatformByName("Reference");
+    Context context(system, integ, platform);
+    for (int copy = 0; copy < numCopies; copy++) {
+        for (int i = 0; i < gridSize; i++)
+            for (int j = 0; j < gridSize; j++)
+                for (int k = 0; k < gridSize; k++) {
+                    Vec3 pos = Vec3(spacing*(i+0.02*genrand_real2(sfmt)), spacing*(j+0.02*genrand_real2(sfmt)), spacing*(k+0.02*genrand_real2(sfmt)));
+                    int index = k+gridSize*(j+gridSize*i);
+                    positions[3*index] = pos;
+                    positions[3*index+1] = Vec3(pos[0]+1.0, pos[1], pos[2]);
+                    positions[3*index+2] = Vec3();
+                }
+        integ.setPositions(copy, positions);
+    }
+
+    // Check the temperature and virtual site locations.
+    
+    const int numSteps = 1000;
+    integ.step(1000);
+    vector<double> ke(numCopies, 0.0);
+    for (int i = 0; i < numSteps; i++) {
+        integ.step(1);
+        vector<State> state(numCopies);
+        for (int j = 0; j < numCopies; j++) {
+            state[j] = integ.getState(j, State::Positions | State::Velocities | State::Forces, true);
+            const vector<Vec3>& pos = state[j].getPositions();
+            for (int k = 0; k < numMolecules; k++)
+                ASSERT_EQUAL_VEC((pos[3*k]+pos[3*k+1])*0.5, pos[3*k+2], 1e-5);
+        }
+        for (int j = 0; j < numParticles; j++) {
+            for (int k = 0; k < numCopies; k++) {
+                Vec3 v = state[k].getVelocities()[j];
+                ke[k] += 0.5*system.getParticleMass(j)*v.dot(v);
+            }
+        }
+    }
+    double meanKE = 0.0;
+    for (int i = 0; i < numCopies; i++)
+        meanKE += ke[i];
+    meanKE /= numSteps*numCopies;
+    double expectedKE = 0.5*numCopies*(2*numMolecules)*3*BOLTZ*temperature;
+    ASSERT_USUALLY_EQUAL_TOL(expectedKE, meanKE, 1e-2);
+}
+
 int main() {
     try {
         testFreeParticles();
         testCMMotionRemoval();
+        testVirtualSites();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;