Fixed inconsistency in computing Ewald self energy

platforms/cuda/src/CudaKernels.cpp

@@ -2117,6 +2117,7 @@ double CudaCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeF
             cuEventRecord(paramsSyncEvent, cu.getCurrentStream());
             cuStreamWaitEvent(pmeStream, paramsSyncEvent, 0);
         }
+        if (hasOffsets)
             energy = 0.0; // The Ewald self energy was computed in the kernel.
         recomputeParams = false;
     }
@@ -2348,12 +2349,14 @@ void CudaCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context,
     
     ewaldSelfEnergy = 0.0;
     if (nonbondedMethod == Ewald || nonbondedMethod == PME || nonbondedMethod == LJPME) {
+        if (cu.getContextIndex() == 0) {
             for (int i = 0; i < force.getNumParticles(); i++) {
                 ewaldSelfEnergy -= baseParticleParamVec[i].x*baseParticleParamVec[i].x*ONE_4PI_EPS0*alpha/sqrt(M_PI);
                 if (doLJPME)
                     ewaldSelfEnergy += baseParticleParamVec[i].z*pow(baseParticleParamVec[i].y*dispersionAlpha, 6)/3.0;
             }
         }
+    }
     if (force.getUseDispersionCorrection() && cu.getContextIndex() == 0 && (nonbondedMethod == CutoffPeriodic || nonbondedMethod == Ewald || nonbondedMethod == PME))
         dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(context.getSystem(), force);
     cu.invalidateMolecules();

platforms/cuda/src/kernels/nonbondedParameters.cu

@@ -31,12 +31,14 @@ extern "C" __global__ void computeParameters(mixed* __restrict__ energyBuffer, b
         charge[i] = params.x;
 #endif
         sigmaEpsilon[i] = make_float2(0.5f*params.y, 2*SQRT(params.z));
-#ifdef INCLUDE_EWALD
+#ifdef HAS_OFFSETS
+    #ifdef INCLUDE_EWALD
         energy -= EWALD_SELF_ENERGY_SCALE*params.x*params.x;
-#endif
-#ifdef INCLUDE_LJPME
+    #endif
+    #ifdef INCLUDE_LJPME
         real sig3 = params.y*params.y*params.y;
         energy += LJPME_SELF_ENERGY_SCALE*sig3*sig3*params.z;
+    #endif
 #endif
     }
 

platforms/opencl/src/OpenCLKernels.cpp

@@ -2181,6 +2181,7 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ
             cl.getQueue().enqueueMarker(&events[0]);
             pmeQueue.enqueueWaitForEvents(events);
         }
+        if (hasOffsets)
             energy = 0.0; // The Ewald self energy was computed in the kernel.
         recomputeParams = false;
     }
@@ -2487,12 +2488,14 @@ void OpenCLCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& contex
     
     ewaldSelfEnergy = 0.0;
     if (nonbondedMethod == Ewald || nonbondedMethod == PME || nonbondedMethod == LJPME) {
+        if (cl.getContextIndex() == 0) {
             for (int i = 0; i < force.getNumParticles(); i++) {
                 ewaldSelfEnergy -= baseParticleParamVec[i].x*baseParticleParamVec[i].x*ONE_4PI_EPS0*alpha/sqrt(M_PI);
                 if (doLJPME)
                     ewaldSelfEnergy += baseParticleParamVec[i].z*pow(baseParticleParamVec[i].y*dispersionAlpha, 6)/3.0;
             }
         }
+    }
     if (force.getUseDispersionCorrection() && cl.getContextIndex() == 0 && (nonbondedMethod == CutoffPeriodic || nonbondedMethod == Ewald || nonbondedMethod == PME))
         dispersionCoefficient = NonbondedForceImpl::calcDispersionCorrection(context.getSystem(), force);
     cl.invalidateMolecules(info);

platforms/opencl/src/kernels/nonbondedParameters.cl

@@ -31,12 +31,14 @@ __kernel void computeParameters(__global mixed* restrict energyBuffer, int inclu
         charge[i] = params.x;
 #endif
         sigmaEpsilon[i] = (float2) (0.5f*params.y, 2*SQRT(params.z));
-#ifdef INCLUDE_EWALD
+#ifdef HAS_OFFSETS
+    #ifdef INCLUDE_EWALD
         energy -= EWALD_SELF_ENERGY_SCALE*params.x*params.x;
-#endif
-#ifdef INCLUDE_LJPME
+    #endif
+    #ifdef INCLUDE_LJPME
         real sig3 = params.y*params.y*params.y;
         energy += LJPME_SELF_ENERGY_SCALE*sig3*sig3*params.z;
+    #endif
 #endif
     }