Fixed bug in CustomIntegrator's handling of massless particles (see bug 1871)

platforms/cuda/src/CudaKernels.cpp

@@ -4554,7 +4554,7 @@ void CudaIntegrateCustomStepKernel::initialize(const System& system, const Custo
     numGlobalVariables = integrator.getNumGlobalVariables();
     int elementSize = (cu.getUseDoublePrecision() || cu.getUseMixedPrecision() ? sizeof(double) : sizeof(float));
     globalValues = new CudaArray(cu, max(1, numGlobalVariables), elementSize, "globalVariables");
-    sumBuffer = new CudaArray(cu, 3*system.getNumParticles(), elementSize, "sumBuffer");
+    sumBuffer = new CudaArray(cu, ((3*system.getNumParticles()+3)/4)*4, elementSize, "sumBuffer");
     potentialEnergy = new CudaArray(cu, 1, cu.getEnergyBuffer().getElementSize(), "potentialEnergy");
     kineticEnergy = new CudaArray(cu, 1, elementSize, "kineticEnergy");
     perDofValues = new CudaParameterSet(cu, integrator.getNumPerDofVariables(), 3*system.getNumParticles(), "perDofVariables", false, cu.getUseDoublePrecision() || cu.getUseMixedPrecision());
@@ -5103,6 +5103,7 @@ void CudaIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegrat
             kernelArgs[i][0][10] = &randomIndex;
             if (requiredUniform[i] > 0)
                 cu.executeKernel(randomKernel, &randomArgs[0], numAtoms);
+            cu.clearBuffer(*sumBuffer);
             cu.executeKernel(kernels[i][0], &kernelArgs[i][0][0], numAtoms);
             cu.executeKernel(kernels[i][1], &kernelArgs[i][1][0], CudaContext::ThreadBlockSize, CudaContext::ThreadBlockSize);
         }
@@ -5150,6 +5151,7 @@ double CudaIntegrateCustomStepKernel::computeKineticEnergy(ContextImpl& context,
     int randomIndex = 0;
     kineticEnergyArgs[1] = &posCorrection;
     kineticEnergyArgs[10] = &randomIndex;
+    cu.clearBuffer(*sumBuffer);
     cu.executeKernel(kineticEnergyKernel, &kineticEnergyArgs[0], cu.getNumAtoms());
     void* args[] = {&sumBuffer->getDevicePointer(), &kineticEnergy->getDevicePointer()};
     cu.executeKernel(sumKineticEnergyKernel, args, CudaContext::ThreadBlockSize, CudaContext::ThreadBlockSize);

platforms/cuda/tests/TestCudaCustomIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2012 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -345,20 +345,20 @@ void testSum() {
     OpenMM_SFMT::SFMT sfmt;
     init_gen_rand(0, sfmt);
     for (int i = 0; i < numParticles; i++) {
-        system.addParticle(1.5);
-        nb->addParticle(i%2 == 0 ? 1 : -1, 0.1, 1);
+        system.addParticle(i%10 == 0 ? 0.0 : 1.5);
+        nb->addParticle(i%2 == 0 ? 0.1 : -0.1, 0.1, 1);
         bool close = true;
         while (close) {
             positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
             close = false;
             for (int j = 0; j < i; ++j) {
                 Vec3 delta = positions[i]-positions[j];
-                if (delta.dot(delta) < 0.1)
+                if (delta.dot(delta) < 1)
                     close = true;
             }
         }
     }
-    CustomIntegrator integrator(0.01);
+    CustomIntegrator integrator(0.005);
     integrator.addGlobalVariable("ke", 0);
     integrator.addComputePerDof("v", "v+dt*f/m");
     integrator.addComputePerDof("x", "x+dt*v");
@@ -368,10 +368,8 @@ void testSum() {
     
     // See if the sum is being computed correctly.
     
-    State state = context.getState(State::Energy);
-    const double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
     for (int i = 0; i < 100; ++i) {
-        state = context.getState(State::Energy);
+        State state = context.getState(State::Energy);
         ASSERT_EQUAL_TOL(state.getKineticEnergy(), integrator.getGlobalVariable(0), 1e-5);
         integrator.step(1);
     }

platforms/opencl/src/OpenCLKernels.cpp

@@ -4786,7 +4786,7 @@ void OpenCLIntegrateCustomStepKernel::initialize(const System& system, const Cus
     numGlobalVariables = integrator.getNumGlobalVariables();
     int elementSize = (cl.getUseDoublePrecision() || cl.getUseMixedPrecision() ? sizeof(double) : sizeof(float));
     globalValues = new OpenCLArray(cl, max(1, numGlobalVariables), elementSize, "globalVariables");
-    sumBuffer = new OpenCLArray(cl, 3*system.getNumParticles(), elementSize, "sumBuffer");
+    sumBuffer = new OpenCLArray(cl, ((3*system.getNumParticles()+3)/4)*4, elementSize, "sumBuffer");
     potentialEnergy = new OpenCLArray(cl, 1, cl.getEnergyBuffer().getElementSize(), "potentialEnergy");
     kineticEnergy = new OpenCLArray(cl, 1, elementSize, "kineticEnergy");
     perDofValues = new OpenCLParameterSet(cl, integrator.getNumPerDofVariables(), 3*system.getNumParticles(), "perDofVariables", false, cl.getUseDoublePrecision() || cl.getUseMixedPrecision());
@@ -5329,6 +5329,7 @@ void OpenCLIntegrateCustomStepKernel::execute(ContextImpl& context, CustomIntegr
             kernels[i][0].setArg<cl_uint>(10, integration.prepareRandomNumbers(requiredGaussian[i]));
             if (requiredUniform[i] > 0)
                 cl.executeKernel(randomKernel, numAtoms);
+            cl.clearBuffer(*sumBuffer);
             cl.executeKernel(kernels[i][0], numAtoms);
             cl.executeKernel(kernels[i][1], OpenCLContext::ThreadBlockSize, OpenCLContext::ThreadBlockSize);
         }
@@ -5375,6 +5376,7 @@ double OpenCLIntegrateCustomStepKernel::computeKineticEnergy(ContextImpl& contex
             cl.executeKernel(sumPotentialEnergyKernel, OpenCLContext::ThreadBlockSize, OpenCLContext::ThreadBlockSize);
         forcesAreValid = true;
     }
+    cl.clearBuffer(*sumBuffer);
     cl.executeKernel(kineticEnergyKernel, cl.getNumAtoms());
     cl.executeKernel(sumKineticEnergyKernel, OpenCLContext::ThreadBlockSize, OpenCLContext::ThreadBlockSize);
     if (cl.getUseDoublePrecision() || cl.getUseMixedPrecision()) {

platforms/opencl/tests/TestOpenCLCustomIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2011 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -345,20 +345,20 @@ void testSum() {
     OpenMM_SFMT::SFMT sfmt;
     init_gen_rand(0, sfmt);
     for (int i = 0; i < numParticles; i++) {
-        system.addParticle(1.5);
-        nb->addParticle(i%2 == 0 ? 1 : -1, 0.1, 1);
+        system.addParticle(i%10 == 0 ? 0.0 : 1.5);
+        nb->addParticle(i%2 == 0 ? 0.1 : -0.1, 0.1, 1);
         bool close = true;
         while (close) {
             positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
             close = false;
             for (int j = 0; j < i; ++j) {
                 Vec3 delta = positions[i]-positions[j];
-                if (delta.dot(delta) < 0.1)
+                if (delta.dot(delta) < 1)
                     close = true;
             }
         }
     }
-    CustomIntegrator integrator(0.01);
+    CustomIntegrator integrator(0.005);
     integrator.addGlobalVariable("ke", 0);
     integrator.addComputePerDof("v", "v+dt*f/m");
     integrator.addComputePerDof("x", "x+dt*v");
@@ -368,10 +368,8 @@ void testSum() {
     
     // See if the sum is being computed correctly.
     
-    State state = context.getState(State::Energy);
-    const double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
     for (int i = 0; i < 100; ++i) {
-        state = context.getState(State::Energy);
+        State state = context.getState(State::Energy);
         ASSERT_EQUAL_TOL(state.getKineticEnergy(), integrator.getGlobalVariable(0), 1e-5);
         integrator.step(1);
     }

platforms/reference/tests/TestReferenceCustomIntegrator.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2011 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -336,20 +336,20 @@ void testSum() {
     OpenMM_SFMT::SFMT sfmt;
     init_gen_rand(0, sfmt);
     for (int i = 0; i < numParticles; i++) {
-        system.addParticle(1.5);
-        nb->addParticle(i%2 == 0 ? 1 : -1, 0.1, 1);
+        system.addParticle(i%10 == 0 ? 0.0 : 1.5);
+        nb->addParticle(i%2 == 0 ? 0.1 : -0.1, 0.1, 1);
         bool close = true;
         while (close) {
             positions[i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
             close = false;
             for (int j = 0; j < i; ++j) {
                 Vec3 delta = positions[i]-positions[j];
-                if (delta.dot(delta) < 0.1)
+                if (delta.dot(delta) < 1)
                     close = true;
             }
         }
     }
-    CustomIntegrator integrator(0.01);
+    CustomIntegrator integrator(0.005);
     integrator.addGlobalVariable("ke", 0);
     integrator.addComputePerDof("v", "v+dt*f/m");
     integrator.addComputePerDof("x", "x+dt*v");
@@ -359,10 +359,8 @@ void testSum() {
     
     // See if the sum is being computed correctly.
     
-    State state = context.getState(State::Energy);
-    const double initialEnergy = state.getKineticEnergy()+state.getPotentialEnergy();
     for (int i = 0; i < 100; ++i) {
-        state = context.getState(State::Energy);
+        State state = context.getState(State::Energy);
         ASSERT_EQUAL_TOL(state.getKineticEnergy(), integrator.getGlobalVariable(0), 1e-5);
         integrator.step(1);
     }