Fix for PME energy

plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaCudaPME.cu

@@ -901,11 +901,14 @@ void kComputeInducedDipoleForceAndEnergy_kernel()
         float* phid = &cAmoebaSim.pPhid[10*i];
         float4 f = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
 
+        energy += cSim.pmeGridSize.x*cSim.invPeriodicBoxSizeX*inducedDipole[0]*phi[1];
+        energy += cSim.pmeGridSize.y*cSim.invPeriodicBoxSizeY*inducedDipole[1]*phi[2];
+        energy += cSim.pmeGridSize.z*cSim.invPeriodicBoxSizeZ*inducedDipole[2]*phi[3];
+
         for (int k = 0; k < 3; k++) {
             int j1 = deriv1[k+1];
             int j2 = deriv2[k+1];
             int j3 = deriv3[k+1];
-            energy += inducedDipole[k]*phi[k+1];
 
             f.x += (inducedDipole[k]+inducedDipolePolar[k])*phi[j1] + inducedDipole[k]*phip[j1] + inducedDipolePolar[k]*phid[j1];
             f.y += (inducedDipole[k]+inducedDipolePolar[k])*phi[j2] + inducedDipole[k]*phip[j2] + inducedDipolePolar[k]*phid[j2];
@@ -923,17 +926,10 @@ void kComputeInducedDipoleForceAndEnergy_kernel()
             f.z += multipole[k]*phidp[deriv3[k]];
         }
 
-
         f.x *= 0.5f*cAmoebaSim.electric*cSim.pmeGridSize.x*cSim.invPeriodicBoxSizeX;
         f.y *= 0.5f*cAmoebaSim.electric*cSim.pmeGridSize.y*cSim.invPeriodicBoxSizeY;
         f.z *= 0.5f*cAmoebaSim.electric*cSim.pmeGridSize.z*cSim.invPeriodicBoxSizeZ;
 
-/*
-        f.x *= 0.5f*cAmoebaSim.electric;
-        f.y *= 0.5f*cAmoebaSim.electric;
-        f.z *= 0.5f*cAmoebaSim.electric;
-*/
-
         float4 force = cSim.pForce4[i];
         force.x -= f.x;
         force.y -= f.y;

plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaCudaPmeDirectElectrostatic.cu

@@ -1399,6 +1399,7 @@ void cudaComputeAmoebaPmeDirectElectrostatic( amoebaGpuContext amoebaGpu )
    @param amoebaGpu        amoebaGpu context
 
    --------------------------------------------------------------------------------------- */
+extern double kReduceEnergy(  gpuContext gpu);
 
 void cudaComputeAmoebaPmeElectrostatic( amoebaGpuContext amoebaGpu )
 {
@@ -1423,7 +1424,10 @@ void cudaComputeAmoebaPmeElectrostatic( amoebaGpuContext amoebaGpu )
             index                          += 3;
         }
         //zeroForce( amoebaGpu );
+        double dem = kReduceEnergy( gpu );
         kCalculateAmoebaPMEInducedDipoleForces( amoebaGpu );
+        double dep = kReduceEnergy( gpu );
+        fprintf( stderr, "Recip Em=%15.7e ep=%15.7e  ttl=%15.7e", dem/4.184, (dep-dem)/4.184, dep/4.184 );
         copyForce( amoebaGpu, conversion );
         VectorOfDoubleVectors outputVector1;
         cudaLoadCudaFloatArray( gpu->natoms,  3, amoebaGpu->psForce,      outputVector1, gpu->psAtomIndex->_pSysData );