Bug fixes

platforms/cuda/src/CudaKernels.cpp

@@ -1668,7 +1668,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
             cosSinSums.initialize(cu, (2*kmaxx-1)*(2*kmaxy-1)*(2*kmaxz-1), elementSize, "cosSinSums");
         }
     }
-    else if ((nonbondedMethod == PME && hasCoulomb) || doLJPME) {
+    else if (((nonbondedMethod == PME || nonbondedMethod == LJPME) && hasCoulomb) || doLJPME) {
         // Compute the PME parameters.
 
         NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSizeX, gridSizeY, gridSizeZ, false);
@@ -1934,7 +1934,6 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
 
     string source = cu.replaceStrings(CudaKernelSources::coulombLennardJones, defines);
     charges.initialize(cu, cu.getPaddedNumAtoms(), cu.getUseDoublePrecision() ? sizeof(double) : sizeof(float), "charges");
-    if (hasCoulomb) {
     map<string, string> replacements;
     if (usePosqCharges) {
         cu.setCharges(chargeVec);
@@ -1953,18 +1952,16 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
         replacements["CHARGE1"] = prefix+"charge1";
         replacements["CHARGE2"] = prefix+"charge2";
     }
-        source = cu.replaceStrings(source, replacements);
+    if (hasCoulomb)
         cu.getNonbondedUtilities().addParameter(CudaNonbondedUtilities::ParameterInfo(prefix+"charge", "real", 1, charges.getElementSize(), charges.getDevicePointer()));
-    }
     if (hasLJ) {
         sigmaEpsilon.initialize<float2>(cu, cu.getPaddedNumAtoms(), "sigmaEpsilon");
         sigmaEpsilon.upload(sigmaEpsilonVector);
-        map<string, string> replacements;
         replacements["SIGMA_EPSILON1"] = prefix+"sigmaEpsilon1";
         replacements["SIGMA_EPSILON2"] = prefix+"sigmaEpsilon2";
-        source = cu.replaceStrings(source, replacements);
         cu.getNonbondedUtilities().addParameter(CudaNonbondedUtilities::ParameterInfo(prefix+"sigmaEpsilon", "float", 2, sizeof(float2), sigmaEpsilon.getDevicePointer()));
     }
+    source = cu.replaceStrings(source, replacements);
     cu.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source, force.getForceGroup(), true);
 
     // Initialize the exceptions.
@@ -2032,6 +2029,7 @@ double CudaCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeF
 
         // Execute the reciprocal space kernels.
 
+        if (hasCoulomb) {
             void* gridIndexArgs[] = {&cu.getPosq().getDevicePointer(), &pmeAtomGridIndex.getDevicePointer(), cu.getPeriodicBoxSizePointer(),
                     cu.getInvPeriodicBoxSizePointer(), cu.getPeriodicBoxVecXPointer(), cu.getPeriodicBoxVecYPointer(), cu.getPeriodicBoxVecZPointer(),
                     recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2]};
@@ -2087,10 +2085,11 @@ double CudaCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeF
                     recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2], &pmeAtomGridIndex.getDevicePointer(),
                     &charges.getDevicePointer()};
             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) {
+        if (doLJPME && hasLJ) {
             void* gridIndexArgs[] = {&cu.getPosq().getDevicePointer(), &pmeAtomGridIndex.getDevicePointer(), cu.getPeriodicBoxSizePointer(),
                     cu.getInvPeriodicBoxSizePointer(), cu.getPeriodicBoxVecXPointer(), cu.getPeriodicBoxVecYPointer(), cu.getPeriodicBoxVecZPointer(),
                     recipBoxVectorPointer[0], recipBoxVectorPointer[1], recipBoxVectorPointer[2]};

platforms/opencl/src/OpenCLKernels.cpp

@@ -1657,7 +1657,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
             cosSinSums.initialize(cl, (2*kmaxx-1)*(2*kmaxy-1)*(2*kmaxz-1), elementSize, "cosSinSums");
         }
     }
-    else if ((nonbondedMethod == PME && hasCoulomb) || doLJPME) {
+    else if (((nonbondedMethod == PME || nonbondedMethod == LJPME) && hasCoulomb) || doLJPME) {
         // Compute the PME parameters.
 
         NonbondedForceImpl::calcPMEParameters(system, force, alpha, gridSizeX, gridSizeY, gridSizeZ, false);
@@ -1865,7 +1865,6 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
     
     string source = cl.replaceStrings(OpenCLKernelSources::coulombLennardJones, defines);
     charges.initialize(cl, cl.getPaddedNumAtoms(), cl.getUseDoublePrecision() ? sizeof(double) : sizeof(float), "charges");
-    if (hasCoulomb) {
     map<string, string> replacements;
     if (usePosqCharges) {
         cl.setCharges(chargeVec);
@@ -1884,18 +1883,16 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
         replacements["CHARGE1"] = prefix+"charge1";
         replacements["CHARGE2"] = prefix+"charge2";
     }
-        source = cl.replaceStrings(source, replacements);
+    if (hasCoulomb)
         cl.getNonbondedUtilities().addParameter(OpenCLNonbondedUtilities::ParameterInfo(prefix+"charge", "real", 1, charges.getElementSize(), charges.getDeviceBuffer()));
-    }
     if (hasLJ) {
         sigmaEpsilon.initialize<mm_float2>(cl, cl.getPaddedNumAtoms(), "sigmaEpsilon");
         sigmaEpsilon.upload(sigmaEpsilonVector);
-        map<string, string> replacements;
         replacements["SIGMA_EPSILON1"] = prefix+"sigmaEpsilon1";
         replacements["SIGMA_EPSILON2"] = prefix+"sigmaEpsilon2";
-        source = cl.replaceStrings(source, replacements);
         cl.getNonbondedUtilities().addParameter(OpenCLNonbondedUtilities::ParameterInfo(prefix+"sigmaEpsilon", "float", 2, sizeof(cl_float2), sigmaEpsilon.getDeviceBuffer()));
     }
+    source = cl.replaceStrings(source, replacements);
     cl.getNonbondedUtilities().addInteraction(useCutoff, usePeriodic, true, force.getCutoffDistance(), exclusionList, source, force.getForceGroup());
 
     // Initialize the exceptions.
@@ -2096,6 +2093,7 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ
         
         // Execute the reciprocal space kernels.
 
+        if (hasCoulomb) {
             setPeriodicBoxArgs(cl, pmeUpdateBsplinesKernel, 4);
             if (cl.getUseDoublePrecision()) {
                 pmeUpdateBsplinesKernel.setArg<mm_double4>(9, recipBoxVectors[0]);
@@ -2187,8 +2185,9 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ
                 cl.executeKernel(pmeInterpolateForceKernel, 2*cl.getDevice().getInfo<CL_DEVICE_MAX_COMPUTE_UNITS>(), 1);
             else
                 cl.executeKernel(pmeInterpolateForceKernel, cl.getNumAtoms());
+        }
         
-        if (doLJPME) {
+        if (doLJPME && hasLJ) {
             setPeriodicBoxArgs(cl, pmeDispersionUpdateBsplinesKernel, 4);
             if (cl.getUseDoublePrecision()) {
                 pmeDispersionUpdateBsplinesKernel.setArg<mm_double4>(9, recipBoxVectors[0]);
@@ -2267,7 +2266,7 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ
             if (includeEnergy)
                 cl.executeKernel(pmeDispersionEvalEnergyKernel, gridSizeX*gridSizeY*gridSizeZ);
             cl.executeKernel(pmeDispersionConvolutionKernel, gridSizeX*gridSizeY*gridSizeZ);
-            fft->execFFT(pmeGrid2, pmeGrid, false);
+            dispersionFft->execFFT(pmeGrid2, pmeGrid, false);
             setPeriodicBoxArgs(cl, pmeDispersionInterpolateForceKernel, 3);
             if (cl.getUseDoublePrecision()) {
                 pmeDispersionInterpolateForceKernel.setArg<mm_double4>(8, recipBoxVectors[0]);

platforms/opencl/src/kernels/gbsaObc_cpu.cl

@@ -18,7 +18,7 @@ __kernel void computeBornSum(
 #else
         __global real* restrict global_bornSum,
 #endif
-        __global const real4* restrict posq, , __global const real* restrict charge, __global const float2* restrict global_params,
+        __global const real4* restrict posq, __global const real* restrict charge, __global const float2* restrict global_params,
 #ifdef USE_CUTOFF
         __global const int* restrict tiles, __global const unsigned int* restrict interactionCount, real4 periodicBoxSize, real4 invPeriodicBoxSize,
         real4 periodicBoxVecX, real4 periodicBoxVecY, real4 periodicBoxVecZ, unsigned int maxTiles, __global const real4* restrict blockCenter,

tests/TestDispersionPME.h

@@ -1269,7 +1269,7 @@ void testWater2DpmeEnergiesForcesWithExclusions() {
     const vector<Vec3>& forces = state.getForces();
 
 
-    ASSERT_EQUAL_TOL(refenergy, energy, 1E-4);
+    ASSERT_EQUAL_TOL(refenergy, energy, 5E-4);
     for (int n = 0; n < numAtoms; ++n)
         ASSERT_EQUAL_VEC(refforces[n], forces[n], 5E-4);
 }