Merge pull request #2257 from peastman/selfenergy

Fixed inconsistency in computing Ewald self energy

platforms/cuda/src/CudaKernels.cpp

@@ -1677,7 +1677,7 @@ void CudaCalcNonbondedForceKernel::initialize(const System& system, const Nonbon
             paramsDefines["INCLUDE_EWALD"] = "1";
             paramsDefines["EWALD_SELF_ENERGY_SCALE"] = cu.doubleToString(ONE_4PI_EPS0*alpha/sqrt(M_PI));
             for (int i = 0; i < numParticles; i++)
-                ewaldSelfEnergy += baseParticleParamVec[i].x*baseParticleParamVec[i].x*ONE_4PI_EPS0*alpha/sqrt(M_PI);
+                ewaldSelfEnergy -= baseParticleParamVec[i].x*baseParticleParamVec[i].x*ONE_4PI_EPS0*alpha/sqrt(M_PI);
 
             // Create the reciprocal space kernels.
 
@@ -2117,7 +2117,8 @@ double CudaCalcNonbondedForceKernel::execute(ContextImpl& context, bool includeF
             cuEventRecord(paramsSyncEvent, cu.getCurrentStream());
             cuStreamWaitEvent(pmeStream, paramsSyncEvent, 0);
         }
-        energy = 0.0; // The Ewald self energy was computed in the kernel.
+        if (hasOffsets)
+            energy = 0.0; // The Ewald self energy was computed in the kernel.
         recomputeParams = false;
     }
     
@@ -2348,10 +2349,12 @@ void CudaCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& context,
     
     ewaldSelfEnergy = 0.0;
     if (nonbondedMethod == Ewald || nonbondedMethod == PME || nonbondedMethod == LJPME) {
-        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 (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))

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

@@ -1667,7 +1667,7 @@ void OpenCLCalcNonbondedForceKernel::initialize(const System& system, const Nonb
             paramsDefines["INCLUDE_EWALD"] = "1";
             paramsDefines["EWALD_SELF_ENERGY_SCALE"] = cl.doubleToString(ONE_4PI_EPS0*alpha/sqrt(M_PI));
             for (int i = 0; i < numParticles; i++)
-                ewaldSelfEnergy += baseParticleParamVec[i].x*baseParticleParamVec[i].x*ONE_4PI_EPS0*alpha/sqrt(M_PI);
+                ewaldSelfEnergy -= baseParticleParamVec[i].x*baseParticleParamVec[i].x*ONE_4PI_EPS0*alpha/sqrt(M_PI);
 
             // Create the reciprocal space kernels.
 
@@ -2181,7 +2181,8 @@ double OpenCLCalcNonbondedForceKernel::execute(ContextImpl& context, bool includ
             cl.getQueue().enqueueMarker(&events[0]);
             pmeQueue.enqueueWaitForEvents(events);
         }
-        energy = 0.0; // The Ewald self energy was computed in the kernel.
+        if (hasOffsets)
+            energy = 0.0; // The Ewald self energy was computed in the kernel.
         recomputeParams = false;
     }
     
@@ -2487,10 +2488,12 @@ void OpenCLCalcNonbondedForceKernel::copyParametersToContext(ContextImpl& contex
     
     ewaldSelfEnergy = 0.0;
     if (nonbondedMethod == Ewald || nonbondedMethod == PME || nonbondedMethod == LJPME) {
-        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 (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))

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
     }