RPMDIntegrator invokes updateContextState() (see bug 1798)

plugins/rpmd/platforms/cuda/src/CudaRpmdKernels.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-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -128,9 +128,8 @@ void CudaIntegrateRPMDStepKernel::initialize(const System& system, const RPMDInt
     pileKernel = cu.getKernel(module, "applyPileThermostat");
     stepKernel = cu.getKernel(module, "integrateStep");
     velocitiesKernel = cu.getKernel(module, "advanceVelocities");
-    copyPositionsToContextKernel = cu.getKernel(module, "copyPositionsToContext");
-    copyVelocitiesToContextKernel = cu.getKernel(module, "copyVelocitiesToContext");
-    copyForcesFromContextKernel = cu.getKernel(module, "copyForcesFromContext");
+    copyToContextKernel = cu.getKernel(module, "copyDataToContext");
+    copyFromContextKernel = cu.getKernel(module, "copyDataFromContext");
     translateKernel = cu.getKernel(module, "applyCellTranslations");
 }
 
@@ -185,11 +184,11 @@ void CudaIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDIntegr
 
 void CudaIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
     for (int i = 0; i < numCopies; i++) {
-        void* copyToContextArgs[] = {&positions->getDevicePointer(), &cu.getPosq().getDevicePointer(), &cu.getAtomIndexArray().getDevicePointer(), &i};
-        cu.executeKernel(copyPositionsToContextKernel, copyToContextArgs, cu.getNumAtoms());
+        void* copyToContextArgs[] = {&velocities->getDevicePointer(), &cu.getVelm().getDevicePointer(), &positions->getDevicePointer(),
+                &cu.getPosq().getDevicePointer(), &cu.getAtomIndexArray().getDevicePointer(), &i};
+        cu.executeKernel(copyToContextKernel, copyToContextArgs, cu.getNumAtoms());
+        context.updateContextState();
         context.calcForcesAndEnergy(true, false);
-        void* copyFromContextArgs[] = {&cu.getForce().getDevicePointer(), &forces->getDevicePointer(), &cu.getAtomIndexArray().getDevicePointer(), &i};
-        cu.executeKernel(copyForcesFromContextKernel, copyFromContextArgs, cu.getNumAtoms());
         if (cu.getAtomsWereReordered() && cu.getNonbondedUtilities().getUsePeriodic()) {
             // Atoms may have been translated into a different periodic box, so apply
             // the same translation to all the beads.
@@ -197,6 +196,9 @@ void CudaIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
             void* args[] = {&positions->getDevicePointer(), &cu.getPosq().getDevicePointer(), &cu.getAtomIndexArray().getDevicePointer(), &i};
             cu.executeKernel(translateKernel, args, cu.getNumAtoms());
         }
+        void* copyFromContextArgs[] = {&cu.getForce().getDevicePointer(), &forces->getDevicePointer(), &cu.getVelm().getDevicePointer(),
+                &velocities->getDevicePointer(), &cu.getPosq().getDevicePointer(), &positions->getDevicePointer(), &cu.getAtomIndexArray().getDevicePointer(), &i};
+        cu.executeKernel(copyFromContextKernel, copyFromContextArgs, cu.getNumAtoms());
     }
 }
 
@@ -267,10 +269,9 @@ void CudaIntegrateRPMDStepKernel::setVelocities(int copy, const vector<Vec3>& ve
 }
 
 void CudaIntegrateRPMDStepKernel::copyToContext(int copy, ContextImpl& context) {
-    void* copyPositionsArgs[] = {&positions->getDevicePointer(), &cu.getPosq().getDevicePointer(), &cu.getAtomIndexArray().getDevicePointer(), &copy};
-    cu.executeKernel(copyPositionsToContextKernel, copyPositionsArgs, cu.getNumAtoms());
-    void* copyVelocitiesArgs[] = {&velocities->getDevicePointer(), &cu.getVelm().getDevicePointer(), &cu.getAtomIndexArray().getDevicePointer(), &copy};
-    cu.executeKernel(copyVelocitiesToContextKernel, copyVelocitiesArgs, cu.getNumAtoms());
+    void* copyArgs[] = {&velocities->getDevicePointer(), &cu.getVelm().getDevicePointer(), &positions->getDevicePointer(),
+            &cu.getPosq().getDevicePointer(), &cu.getAtomIndexArray().getDevicePointer(), &copy};
+    cu.executeKernel(copyToContextKernel, copyArgs, cu.getNumAtoms());
 }
 
 string CudaIntegrateRPMDStepKernel::createFFT(int size, const string& variable, bool forward) {

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

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2011-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -91,7 +91,7 @@ private:
     CudaArray* forces;
     CudaArray* positions;
     CudaArray* velocities;
-    CUfunction pileKernel, stepKernel, velocitiesKernel, copyPositionsToContextKernel, copyVelocitiesToContextKernel, copyForcesFromContextKernel, translateKernel;
+    CUfunction pileKernel, stepKernel, velocitiesKernel, copyToContextKernel, copyFromContextKernel, translateKernel;
 };
 
 } // namespace OpenMM

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

@@ -181,35 +181,33 @@ extern "C" __global__ void advanceVelocities(mixed4* velm, long long* force, mix
 }
 
 /**
- * Copy a set of positions from the integrator's arrays to the context.
+ * Copy a set of positions and velocities from the integrator's arrays to the context.
  */
-extern "C" __global__ void copyPositionsToContext(mixed4* src, real4* dst, int* order, int copy) {
+extern "C" __global__ void copyDataToContext(mixed4* srcVel, mixed4* dstVel, mixed4* srcPos, real4* dstPos, int* order, int copy) {
     const int base = copy*PADDED_NUM_ATOMS;
     for (int particle = blockIdx.x*blockDim.x+threadIdx.x; particle < NUM_ATOMS; particle += blockDim.x*gridDim.x) {
-        mixed4 posq = src[base+order[particle]];
-        dst[particle] = make_real4(posq.x, posq.y, posq.z, posq.w);
+        int index = base+order[particle];
+        dstVel[particle] = srcVel[index];
+        mixed4 posq = srcPos[index];
+        dstPos[particle] = make_real4(posq.x, posq.y, posq.z, posq.w);
     }
 }
 
 /**
- * Copy a set of velocities from the integrator's arrays to the context.
+ * Copy a set of positions, velocities, and forces from the context to the integrator's arrays.
  */
-extern "C" __global__ void copyVelocitiesToContext(mixed4* src, mixed4* dst, int* order, int copy) {
+extern "C" __global__ void copyDataFromContext(long long* srcForce, long long* dstForce, mixed4* srcVel, mixed4* dstVel,
+        real4* srcPos, mixed4* dstPos, int* order, int copy) {
     const int base = copy*PADDED_NUM_ATOMS;
     for (int particle = blockIdx.x*blockDim.x+threadIdx.x; particle < NUM_ATOMS; particle += blockDim.x*gridDim.x) {
-        dst[particle] = src[base+order[particle]];
-    }
-}
+        int index = order[particle];
+        dstForce[base*3+index] = srcForce[particle];
+        dstForce[base*3+index+PADDED_NUM_ATOMS] = srcForce[particle+PADDED_NUM_ATOMS];
+        dstForce[base*3+index+PADDED_NUM_ATOMS*2] = srcForce[particle+PADDED_NUM_ATOMS*2];
+        dstVel[base+index] = srcVel[particle];
+        real4 posq = srcPos[particle];
+        dstPos[base+index] = make_mixed4(posq.x, posq.y, posq.z, posq.w);
 
-/**
- * Copy a set of forces from the context to the integrator's arrays.
- */
-extern "C" __global__ void copyForcesFromContext(long long* src, long long* dst, int* order, int copy) {
-    const int base = copy*PADDED_NUM_ATOMS*3;
-    for (int particle = blockIdx.x*blockDim.x+threadIdx.x; particle < NUM_ATOMS; particle += blockDim.x*gridDim.x) {
-        dst[base+order[particle]] = src[particle];
-        dst[base+order[particle]+PADDED_NUM_ATOMS] = src[particle+PADDED_NUM_ATOMS];
-        dst[base+order[particle]+PADDED_NUM_ATOMS*2] = src[particle+PADDED_NUM_ATOMS*2];
     }
 }
 

plugins/rpmd/platforms/cuda/tests/TestCudaRpmd.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:                                                              *
  *                                                                            *
@@ -34,6 +34,7 @@
  */
 
 #include "openmm/internal/AssertionUtilities.h"
+#include "openmm/CMMotionRemover.h"
 #include "openmm/Context.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/HarmonicBondForce.h"
@@ -222,6 +223,55 @@ void testParaHydrogen() {
     ASSERT_USUALLY_EQUAL_TOL(60.0, 1.5*temperature+meanKE, 0.02);
 }
 
+Vec3 calcCM(const vector<Vec3>& values, System& system) {
+    Vec3 cm;
+    for (int j = 0; j < system.getNumParticles(); ++j) {
+        cm[0] += values[j][0]*system.getParticleMass(j);
+        cm[1] += values[j][1]*system.getParticleMass(j);
+        cm[2] += values[j][2]*system.getParticleMass(j);
+    }
+    return cm;
+}
+
+void testCMMotionRemoval() {
+    const int numParticles = 100;
+    const int numCopies = 30;
+    const double temperature = 300.0;
+    const double mass = 1.0;
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(mass);
+    system.addForce(new CMMotionRemover());
+    RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
+    Platform& platform = Platform::getPlatformByName("CUDA");
+    Context context(system, integ, platform);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numCopies; i++)
+    {
+        for (int j = 0; j < numParticles; j++)
+            positions[j] = Vec3(0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt));
+        Vec3 cmPos = calcCM(positions, system);
+        for (int j = 0; j < numParticles; j++)
+            positions[j] -= cmPos*(1/(mass*numParticles));
+        integ.setPositions(i, positions);
+    }
+    
+    // Make sure the CMMotionRemover is getting applied.
+    
+    for (int i = 0; i < 200; ++i) {
+        integ.step(1);
+        Vec3 pos;
+        for (int j = 0; j < numCopies; j++) {
+            State state = integ.getState(0, State::Positions | State::Velocities);
+            pos += calcCM(state.getPositions(), system);
+        }
+        pos *= 1.0/numCopies;
+        ASSERT_EQUAL_VEC(Vec3(), pos, 0.5);
+    }
+}
+
 int main(int argc, char* argv[]) {
     try {
         registerRPMDCudaKernelFactories();
@@ -229,6 +279,7 @@ int main(int argc, char* argv[]) {
             Platform::getPlatformByName("CUDA").setPropertyDefaultValue("CudaPrecision", string(argv[1]));
         testFreeParticles();
         testParaHydrogen();
+        testCMMotionRemoval();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

plugins/rpmd/platforms/opencl/src/OpenCLRpmdKernels.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-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -108,15 +108,12 @@ void OpenCLIntegrateRPMDStepKernel::initialize(const System& system, const RPMDI
     pileKernel = cl::Kernel(program, "applyPileThermostat");
     stepKernel = cl::Kernel(program, "integrateStep");
     velocitiesKernel = cl::Kernel(program, "advanceVelocities");
-    copyPositionsToContextKernel = cl::Kernel(program, "copyPositionsToContext");
-    copyVelocitiesToContextKernel = cl::Kernel(program, "copyVelocitiesToContext");
-    copyForcesFromContextKernel = cl::Kernel(program, "copyForcesFromContext");
+    copyToContextKernel = cl::Kernel(program, "copyDataToContext");
+    copyFromContextKernel = cl::Kernel(program, "copyDataFromContext");
     translateKernel = cl::Kernel(program, "applyCellTranslations");
 }
 
-void OpenCLIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDIntegrator& integrator, bool forcesAreValid) {
-    OpenCLIntegrationUtilities& integration = cl.getIntegrationUtilities();
-    if (!hasInitializedKernel) {
+void OpenCLIntegrateRPMDStepKernel::initializeKernels(ContextImpl& context) {
     hasInitializedKernel = true;
     pileKernel.setArg<cl::Buffer>(0, velocities->getDeviceBuffer());
     stepKernel.setArg<cl::Buffer>(0, positions->getDeviceBuffer());
@@ -127,16 +124,24 @@ void OpenCLIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDInte
     translateKernel.setArg<cl::Buffer>(0, positions->getDeviceBuffer());
     translateKernel.setArg<cl::Buffer>(1, cl.getPosq().getDeviceBuffer());
     translateKernel.setArg<cl::Buffer>(2, cl.getAtomIndexArray().getDeviceBuffer());
-        copyPositionsToContextKernel.setArg<cl::Buffer>(0, positions->getDeviceBuffer());
-        copyPositionsToContextKernel.setArg<cl::Buffer>(1, cl.getPosq().getDeviceBuffer());
-        copyPositionsToContextKernel.setArg<cl::Buffer>(2, cl.getAtomIndexArray().getDeviceBuffer());
-        copyVelocitiesToContextKernel.setArg<cl::Buffer>(0, velocities->getDeviceBuffer());
-        copyVelocitiesToContextKernel.setArg<cl::Buffer>(1, cl.getVelm().getDeviceBuffer());
-        copyVelocitiesToContextKernel.setArg<cl::Buffer>(2, cl.getAtomIndexArray().getDeviceBuffer());
-        copyForcesFromContextKernel.setArg<cl::Buffer>(0, cl.getForce().getDeviceBuffer());
-        copyForcesFromContextKernel.setArg<cl::Buffer>(1, forces->getDeviceBuffer());
-        copyForcesFromContextKernel.setArg<cl::Buffer>(2, cl.getAtomIndexArray().getDeviceBuffer());
-    }
+    copyToContextKernel.setArg<cl::Buffer>(0, velocities->getDeviceBuffer());
+    copyToContextKernel.setArg<cl::Buffer>(1, cl.getVelm().getDeviceBuffer());
+    copyToContextKernel.setArg<cl::Buffer>(2, positions->getDeviceBuffer());
+    copyToContextKernel.setArg<cl::Buffer>(3, cl.getPosq().getDeviceBuffer());
+    copyToContextKernel.setArg<cl::Buffer>(4, cl.getAtomIndexArray().getDeviceBuffer());
+    copyFromContextKernel.setArg<cl::Buffer>(0, cl.getForce().getDeviceBuffer());
+    copyFromContextKernel.setArg<cl::Buffer>(1, forces->getDeviceBuffer());
+    copyFromContextKernel.setArg<cl::Buffer>(2, cl.getVelm().getDeviceBuffer());
+    copyFromContextKernel.setArg<cl::Buffer>(3, velocities->getDeviceBuffer());
+    copyFromContextKernel.setArg<cl::Buffer>(4, cl.getPosq().getDeviceBuffer());
+    copyFromContextKernel.setArg<cl::Buffer>(5, positions->getDeviceBuffer());
+    copyFromContextKernel.setArg<cl::Buffer>(6, cl.getAtomIndexArray().getDeviceBuffer());
+}
+
+void OpenCLIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDIntegrator& integrator, bool forcesAreValid) {
+    OpenCLIntegrationUtilities& integration = cl.getIntegrationUtilities();
+    if (!hasInitializedKernel)
+        initializeKernels(context);
     
     // Loop over copies and compute the force on each one.
     
@@ -191,11 +196,10 @@ void OpenCLIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDInte
 
 void OpenCLIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
     for (int i = 0; i < numCopies; i++) {
-        copyPositionsToContextKernel.setArg<cl_int>(3, i);
-        cl.executeKernel(copyPositionsToContextKernel, cl.getNumAtoms());
+        copyToContextKernel.setArg<cl_int>(5, i);
+        cl.executeKernel(copyToContextKernel, cl.getNumAtoms());
+        context.updateContextState();
         context.calcForcesAndEnergy(true, false);
-        copyForcesFromContextKernel.setArg<cl_int>(3, i);
-        cl.executeKernel(copyForcesFromContextKernel, cl.getNumAtoms());
         if (cl.getAtomsWereReordered() && cl.getNonbondedUtilities().getUsePeriodic()) {
             // Atoms may have been translated into a different periodic box, so apply
             // the same translation to all the beads.
@@ -203,6 +207,8 @@ void OpenCLIntegrateRPMDStepKernel::computeForces(ContextImpl& context) {
             translateKernel.setArg<cl_int>(3, i);
             cl.executeKernel(translateKernel, cl.getNumAtoms());
         }
+        copyFromContextKernel.setArg<cl_int>(7, i);
+        cl.executeKernel(copyFromContextKernel, cl.getNumAtoms());
     }
 }
 
@@ -261,10 +267,10 @@ void OpenCLIntegrateRPMDStepKernel::setVelocities(int copy, const vector<Vec3>&
 }
 
 void OpenCLIntegrateRPMDStepKernel::copyToContext(int copy, ContextImpl& context) {
-    copyPositionsToContextKernel.setArg<cl_int>(3, copy);
-    cl.executeKernel(copyPositionsToContextKernel, cl.getNumAtoms());
-    copyVelocitiesToContextKernel.setArg<cl_int>(3, copy);
-    cl.executeKernel(copyVelocitiesToContextKernel, cl.getNumAtoms());
+    if (!hasInitializedKernel)
+        initializeKernels(context);
+    copyToContextKernel.setArg<cl_int>(5, copy);
+    cl.executeKernel(copyToContextKernel, cl.getNumAtoms());
 }
 
 string OpenCLIntegrateRPMDStepKernel::createFFT(int size, const string& variable, bool forward) {

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

@@ -9,7 +9,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2011-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -84,6 +84,7 @@ public:
      */
     void copyToContext(int copy, ContextImpl& context);
 private:
+    void initializeKernels(ContextImpl& context);
     void computeForces(ContextImpl& context);
     std::string createFFT(int size, const std::string& variable, bool forward);
     OpenCLContext& cl;
@@ -92,7 +93,7 @@ private:
     OpenCLArray* forces;
     OpenCLArray* positions;
     OpenCLArray* velocities;
-    cl::Kernel pileKernel, stepKernel, velocitiesKernel, copyPositionsToContextKernel, copyVelocitiesToContextKernel, copyForcesFromContextKernel, translateKernel;
+    cl::Kernel pileKernel, stepKernel, velocitiesKernel, copyToContextKernel, copyFromContextKernel, translateKernel;
 };
 
 } // namespace OpenMM

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

@@ -172,32 +172,29 @@ __kernel void advanceVelocities(__global mixed4* velm, __global real4* force, mi
 }
 
 /**
- * Copy a set of positions from the integrator's arrays to the context.
+ * Copy a set of positions and velocities from the integrator's arrays to the context.
  */
-__kernel void copyPositionsToContext(__global mixed4* src, __global real4* dst, __global int* order, int copy) {
+__kernel void copyDataToContext(__global mixed4* srcVel, __global mixed4* dstVel, __global mixed4* srcPos,
+        __global real4* dstPos, __global int* order, int copy) {
     const int base = copy*PADDED_NUM_ATOMS;
     for (int particle = get_global_id(0); particle < NUM_ATOMS; particle += get_global_size(0)) {
-        dst[particle] = convert_real4(src[base+order[particle]]);
+        int index = base+order[particle];
+        dstVel[particle] = srcVel[index];
+        dstPos[particle] = convert_real4(srcPos[index]);
     }
 }
 
 /**
- * Copy a set of velocities from the integrator's arrays to the context.
+ * Copy a set of positions, velocities, and forces from the context to the integrator's arrays.
  */
-__kernel void copyVelocitiesToContext(__global mixed4* src, __global mixed4* dst, __global int* order, int copy) {
+__kernel void copyDataFromContext(__global real4* srcForce, __global real4* dstForce, __global mixed4* srcVel,
+        __global mixed4* dstVel, __global real4* srcPos, __global mixed4* dstPos, __global int* order, int copy) {
     const int base = copy*PADDED_NUM_ATOMS;
     for (int particle = get_global_id(0); particle < NUM_ATOMS; particle += get_global_size(0)) {
-        dst[particle] = src[base+order[particle]];
-    }
-}
-
-/**
- * Copy a set forces from the context to the integrator's arrays.
- */
-__kernel void copyForcesFromContext(__global real4* src, __global real4* dst, __global int* order, int copy) {
-    const int base = copy*PADDED_NUM_ATOMS;
-    for (int particle = get_global_id(0); particle < NUM_ATOMS; particle += get_global_size(0)) {
-        dst[base+order[particle]] = src[particle];
+        int index = base+order[particle];
+        dstForce[index] = srcForce[particle];
+        dstVel[index] = srcVel[particle];
+        dstPos[index] = convert_mixed4(srcPos[particle]);
     }
 }
 

plugins/rpmd/platforms/opencl/tests/TestOpenCLRpmd.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-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -34,6 +34,7 @@
  */
 
 #include "openmm/internal/AssertionUtilities.h"
+#include "openmm/CMMotionRemover.h"
 #include "openmm/Context.h"
 #include "openmm/CustomNonbondedForce.h"
 #include "openmm/HarmonicBondForce.h"
@@ -223,6 +224,55 @@ void testParaHydrogen() {
     ASSERT_USUALLY_EQUAL_TOL(60.0, 1.5*temperature+meanKE, 0.02);
 }
 
+Vec3 calcCM(const vector<Vec3>& values, System& system) {
+    Vec3 cm;
+    for (int j = 0; j < system.getNumParticles(); ++j) {
+        cm[0] += values[j][0]*system.getParticleMass(j);
+        cm[1] += values[j][1]*system.getParticleMass(j);
+        cm[2] += values[j][2]*system.getParticleMass(j);
+    }
+    return cm;
+}
+
+void testCMMotionRemoval() {
+    const int numParticles = 100;
+    const int numCopies = 30;
+    const double temperature = 300.0;
+    const double mass = 1.0;
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(mass);
+    system.addForce(new CMMotionRemover());
+    RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
+    Platform& platform = Platform::getPlatformByName("OpenCL");
+    Context context(system, integ, platform);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numCopies; i++)
+    {
+        for (int j = 0; j < numParticles; j++)
+            positions[j] = Vec3(0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt));
+        Vec3 cmPos = calcCM(positions, system);
+        for (int j = 0; j < numParticles; j++)
+            positions[j] -= cmPos*(1/(mass*numParticles));
+        integ.setPositions(i, positions);
+    }
+    
+    // Make sure the CMMotionRemover is getting applied.
+    
+    for (int i = 0; i < 200; ++i) {
+        integ.step(1);
+        Vec3 pos;
+        for (int j = 0; j < numCopies; j++) {
+            State state = integ.getState(0, State::Positions | State::Velocities);
+            pos += calcCM(state.getPositions(), system);
+        }
+        pos *= 1.0/numCopies;
+        ASSERT_EQUAL_VEC(Vec3(), pos, 0.5);
+    }
+}
+
 int main(int argc, char* argv[]) {
     try {
         registerRPMDOpenCLKernelFactories();
@@ -230,6 +280,7 @@ int main(int argc, char* argv[]) {
             Platform::getPlatformByName("OpenCL").setPropertyDefaultValue("OpenCLPrecision", string(argv[1]));
         testFreeParticles();
         testParaHydrogen();
+        testCMMotionRemoval();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

plugins/rpmd/platforms/reference/src/ReferenceRpmdKernels.cpp

@@ -186,7 +186,6 @@ void ReferenceIntegrateRPMDStepKernel::execute(ContextImpl& context, const RPMDI
         for (int j = 0; j < numParticles; j++)
             velocities[i][j] += forces[i][j]*(halfdt/system.getParticleMass(j));
 
-
     // Apply the PILE-L thermostat again.
     
     for (int particle = 0; particle < numParticles; particle++) {

plugins/rpmd/platforms/reference/tests/TestReferenceRpmd.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:                                                              *
  *                                                                            *
@@ -34,6 +34,7 @@
  */
 
 #include "openmm/internal/AssertionUtilities.h"
+#include "openmm/CMMotionRemover.h"
 #include "openmm/Context.h"
 #include "openmm/HarmonicBondForce.h"
 #include "openmm/Platform.h"
@@ -104,9 +105,59 @@ void testFreeParticles() {
     ASSERT_USUALLY_EQUAL_TOL(expectedRg, sqrt(meanRg2), 1e-3);
 }
 
+Vec3 calcCM(const vector<Vec3>& values, System& system) {
+    Vec3 cm;
+    for (int j = 0; j < system.getNumParticles(); ++j) {
+        cm[0] += values[j][0]*system.getParticleMass(j);
+        cm[1] += values[j][1]*system.getParticleMass(j);
+        cm[2] += values[j][2]*system.getParticleMass(j);
+    }
+    return cm;
+}
+
+void testCMMotionRemoval() {
+    const int numParticles = 100;
+    const int numCopies = 30;
+    const double temperature = 300.0;
+    const double mass = 1.0;
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(mass);
+    system.addForce(new CMMotionRemover());
+    RPMDIntegrator integ(numCopies, temperature, 10.0, 0.001);
+    Platform& platform = Platform::getPlatformByName("Reference");
+    Context context(system, integ, platform);
+    OpenMM_SFMT::SFMT sfmt;
+    init_gen_rand(0, sfmt);
+    vector<Vec3> positions(numParticles);
+    for (int i = 0; i < numCopies; i++)
+    {
+        for (int j = 0; j < numParticles; j++)
+            positions[j] = Vec3(0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt), 0.02*genrand_real2(sfmt));
+        Vec3 cmPos = calcCM(positions, system);
+        for (int j = 0; j < numParticles; j++)
+            positions[j] -= cmPos*(1/(mass*numParticles));
+        integ.setPositions(i, positions);
+    }
+    
+    // Make sure the CMMotionRemover is getting applied.
+    
+    for (int i = 0; i < 200; ++i) {
+        integ.step(1);
+        Vec3 pos;
+        for (int j = 0; j < numCopies; j++) {
+            State state = integ.getState(0, State::Positions | State::Velocities);
+            pos += calcCM(state.getPositions(), system);
+        }
+        pos *= 1.0/numCopies;
+        ASSERT_EQUAL_VEC(Vec3(), pos, 0.5);
+    }
+}
+
 int main() {
     try {
         testFreeParticles();
+        testCMMotionRemoval();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;