Bug fixes and code simplification

platforms/cuda/src/CudaKernels.cpp

@@ -1749,7 +1749,6 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
                     pmeDefines["GRID_SIZE_X"] = cu.intToString(dispersionGridSizeX);
                     pmeDefines["GRID_SIZE_Y"] = cu.intToString(dispersionGridSizeY);
                     pmeDefines["GRID_SIZE_Z"] = cu.intToString(dispersionGridSizeZ);
-                    pmeDefines["EPSILON_FACTOR"] = "1";
                     pmeDefines["RECIP_EXP_FACTOR"] = cu.doubleToString(M_PI*M_PI/(dispersionAlpha*dispersionAlpha));
                     pmeDefines["USE_LJPME"] = "1";
                     double invRCut6 = pow(force.getCutoffDistance(), -6);
@@ -1778,7 +1777,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
                     roundedZSize = PmeOrder*(int) ceil(dispersionGridSizeZ/(double) PmeOrder);
                     gridElements = max(gridElements, dispersionGridSizeX*dispersionGridSizeY*roundedZSize);
                 }
-                pmeGrid1.initialize(cu, gridElements, 2*elementSize, "originalPmeGrid");
+                pmeGrid1.initialize(cu, gridElements, 2*elementSize, "pmeGrid1");
                 pmeGrid2.initialize(cu, gridElements, 2*elementSize, "pmeGrid2");
                 cu.addAutoclearBuffer(pmeGrid2);
                 pmeBsplineModuliX.initialize(cu, gridSizeX, elementSize, "pmeBsplineModuliX");
@@ -2196,14 +2195,14 @@ double CudaCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeF
                 sort->sort(pmeAtomGridIndex);
             }
 
-            cu.clearBuffer(pmeGrid1);
-            void* spreadArgs[] = {&cu.getPosq().getDevicePointer(), &pmeGrid1.getDevicePointer(), cu.getPeriodicBoxSizePointer(),
+            cu.clearBuffer(pmeGrid2);
+            void* spreadArgs[] = {&cu.getPosq().getDevicePointer(), &pmeGrid2.getDevicePointer(), cu.getPeriodicBoxSizePointer(),
                     cu.getInvPeriodicBoxSizePointer(), cu.getPeriodicBoxVecXPointer(), cu.getPeriodicBoxVecYPointer(), cu.getPeriodicBoxVecZPointer(),
                     recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2], &pmeAtomGridIndex.getDevicePointer(),
                     &sigmaEpsilon.getDevicePointer()};
             cu.executeKernel(pmeDispersionSpreadChargeKernel, spreadArgs, cu.getNumAtoms(), 128);
 
-            void* finishSpreadArgs[] = {&pmeGrid1.getDevicePointer()};
+            void* finishSpreadArgs[] = {&pmeGrid2.getDevicePointer(), &pmeGrid1.getDevicePointer()};
             cu.executeKernel(pmeDispersionFinishSpreadChargeKernel, finishSpreadArgs, dispersionGridSizeX*dispersionGridSizeY*dispersionGridSizeZ, 256);
 
             if (useCudaFFT) {

platforms/cuda/src/kernels/pme.cu

@@ -52,7 +52,7 @@ extern "C" __global__ void gridSpreadCharge(const real4* __restrict__ posq, real
         real4 pos = posq[atom];
 #ifdef USE_LJPME
         const float2 sigEps = sigmaEpsilon[atom];
-        const real charge = 8*sigEps.x*sigEps.x*sigEps.x*sigEps.y*EPSILON_FACTOR;
+        const real charge = 8*sigEps.x*sigEps.x*sigEps.x*sigEps.y;
 #else
         const real charge = (CHARGE)*EPSILON_FACTOR;
 #endif
@@ -140,11 +140,8 @@ extern "C" __global__ void finishSpreadCharge(
     const unsigned int gridSize = GRID_SIZE_X*GRID_SIZE_Y*GRID_SIZE_Z;
     real scale = 1/(real) 0x100000000;
     for (int index = blockIdx.x*blockDim.x+threadIdx.x; index < gridSize; index += blockDim.x*gridDim.x) {
-        int xindex = index/(GRID_SIZE_Y*GRID_SIZE_Z);
-        int remainder = index-xindex*GRID_SIZE_Y*GRID_SIZE_Z;
-        int yindex = remainder/GRID_SIZE_Z;
-        int zindex = remainder-yindex*GRID_SIZE_Z;
-        int loadIndex = zindexTable[zindex]+(xindex*GRID_SIZE_Y+yindex)*blockSize;
+        int zindex = index%GRID_SIZE_Z;
+        int loadIndex = zindexTable[zindex] + blockSize*(int) (index/GRID_SIZE_Z);
 #if defined(USE_DOUBLE_PRECISION) || defined(USE_DETERMINISTIC_FORCES)
         grid2[index] = scale*grid1[loadIndex];
 #else