Merge pull request #2174 from andysim/mpid

Fix bug in CUDA LJPME energy without Coulomb terms

platforms/cuda/src/CudaKernels.cpp

@@ -1801,7 +1801,6 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
                 pmeAtomGridIndex.initialize<int2>(cu, numParticles, "pmeAtomGridIndex");
                 int energyElementSize = (cu.getUseDoublePrecision() || cu.getUseMixedPrecision() ? sizeof(double) : sizeof(float));
                 pmeEnergyBuffer.initialize(cu, cu.getNumThreadBlocks()*CudaContext::ThreadBlockSize, energyElementSize, "pmeEnergyBuffer");
-                cu.clearBuffer(pmeEnergyBuffer);
                 sort = new CudaSort(cu, new SortTrait(), cu.getNumAtoms());
                 int cufftVersion;
                 cufftGetVersion(&cufftVersion);
@@ -2192,8 +2191,6 @@ double CudaCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeF
             cu.executeKernel(pmeInterpolateForceKernel, interpolateArgs, cu.getNumAtoms(), 128);
         }
 
-        // As written, we check only the Electrostatic grid pointer to get here.  We could separate them out, but for
-        // now we assume that LJPME can only be used if electrostatic PME is also active.
         if (doLJPME && hasLJ) {
             if (!hasCoulomb) {
                 void* gridIndexArgs[] = {&cu.getPosq().getDevicePointer(), &pmeAtomGridIndex.getDevicePointer(), cu.getPeriodicBoxSizePointer(),
@@ -2202,6 +2199,7 @@ double CudaCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeF
                 cu.executeKernel(pmeDispersionGridIndexKernel, gridIndexArgs, cu.getNumAtoms());
 
                 sort->sort(pmeAtomGridIndex);
+                cu.clearBuffer(pmeEnergyBuffer);
             }
 
             cu.clearBuffer(pmeGrid2);

platforms/opencl/src/OpenCLKernels.cpp

@@ -1760,7 +1760,6 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
                 pmeAtomGridIndex.initialize<mm_int2>(cl, numParticles, "pmeAtomGridIndex");
                 int energyElementSize = (cl.getUseDoublePrecision() || cl.getUseMixedPrecision() ? sizeof(double) : sizeof(float));
                 pmeEnergyBuffer.initialize(cl, cl.getNumThreadBlocks()*OpenCLContext::ThreadBlockSize, energyElementSize, "pmeEnergyBuffer");
-                cl.clearBuffer(pmeEnergyBuffer);
                 sort = new OpenCLSort(cl, new SortTrait(), cl.getNumAtoms());
                 fft = new OpenCLFFT3D(cl, gridSizeX, gridSizeY, gridSizeZ, true);
                 if (doLJPME)
@@ -2377,6 +2376,7 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ
                 pmeDispersionEvalEnergyKernel.setArg<mm_float4>(6, recipBoxVectorsFloat[1]);
                 pmeDispersionEvalEnergyKernel.setArg<mm_float4>(7, recipBoxVectorsFloat[2]);
             }
+            if (!hasCoulomb) cl.clearBuffer(pmeEnergyBuffer);
             if (includeEnergy)
                 cl.executeKernel(pmeDispersionEvalEnergyKernel, gridSizeX*gridSizeY*gridSizeZ);
             cl.executeKernel(pmeDispersionConvolutionKernel, gridSizeX*gridSizeY*gridSizeZ);

tests/TestDispersionPME.h

@@ -1307,6 +1307,11 @@ void testWater125DpmeVsLongCutoffNoExclusions() {
     State state = context.getState(State::Forces | State::Energy);
     double energy = state.getPotentialEnergy();
 
+    // Make another call to the get the energy to test that the call is
+    // idempotent when Coulomb terms are absent.
+    double secondEnergy = context.getState(State::Energy).getPotentialEnergy();
+    ASSERT_EQUAL_TOL(secondEnergy, energy, 5E-5);
+
 //Gromacs reference values.  See comments in testWater2DpmeEnergiesForcesNoExclusions() for details.
 //Coordinates are from make_waterbox, and the .gro file looks like
 //