Bug fixes to AmoebaGeneralizedKirkwoodForce

plugins/amoeba/platforms/cuda2/src/AmoebaCudaKernels.cpp

@@ -1778,7 +1778,7 @@ void CudaCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System& syst
     defines["GK_FD"] = cu.doubleToString(2*(1-solventDielectric)/(1+2*solventDielectric));
     defines["GK_FQ"] = cu.doubleToString(3*(1-solventDielectric)/(2+3*solventDielectric));
     defines["EPSILON_FACTOR"] = cu.doubleToString(138.9354558456);
-    defines["ENERGY_SCALE_FACTOR"] = cu.doubleToString(138.9354558456/solventDielectric);
+    defines["ENERGY_SCALE_FACTOR"] = cu.doubleToString(138.9354558456/force.getSoluteDielectric());
     if (multipoles->getPolarizationType() == AmoebaMultipoleForce::Direct)
         defines["DIRECT_POLARIZATION"] = "";
     stringstream forceSource;
@@ -1843,11 +1843,6 @@ void CudaCalcAmoebaGeneralizedKirkwoodForceKernel::finishComputation(CudaArray&
         &labFrameQuadrupoles.getDevicePointer(), &inducedDipoleS->getDevicePointer(), &inducedDipolePolarS->getDevicePointer(),
         &bornRadii->getDevicePointer(), &bornForce->getDevicePointer()};
     cu.executeKernel(gkForceKernel, gkForceArgs, numForceThreadBlocks*forceThreadBlockSize, forceThreadBlockSize);
-    printf("bornForce\n");
-    vector<long long> f;
-    bornForce->download(f);
-    for (int i = 0; i < cu.getNumAtoms(); i++)
-        printf("%d %g\n", i, f[i]/(double) 0xFFFFFFFF);
 
     // Compute cavity term...
     

plugins/amoeba/platforms/cuda2/src/kernels/amoebaGk.cu

@@ -191,6 +191,7 @@ inline __device__ void zeroAtomData(AtomData2& data) {
     data.force = make_real3(0);
     data.bornForce = 0;
 }
+
 /**
  * Compute electrostatic interactions.
  */
@@ -245,7 +246,7 @@ extern "C" __global__ void computeGKForces(
                 
                 for (unsigned int j = 0; j < TILE_SIZE; j++) {
                     int atom2 = y*TILE_SIZE+j;
-                    if (atom1 != atom2 && atom1 < NUM_ATOMS && atom2 < NUM_ATOMS) {
+                    if (atom1 < NUM_ATOMS && atom2 < NUM_ATOMS) {
                         real3 tempForce;
                         real tempEnergy;
                         computeOneInteractionF1(data, localData[tbx+j], tempEnergy, tempForce);
@@ -264,7 +265,7 @@ extern "C" __global__ void computeGKForces(
                 zeroAtomData(data);
                 for (unsigned int j = 0; j < TILE_SIZE; j++) {
                     int atom2 = y*TILE_SIZE+j;
-                    if (atom1 != atom2 && atom1 < NUM_ATOMS && atom2 < NUM_ATOMS) {
+                    if (atom1 < NUM_ATOMS && atom2 < NUM_ATOMS) {
                         real3 tempTorque;
                         computeOneInteractionT1(data, localData[tbx+j], tempTorque);
                         computeOneInteractionT2(data, localData[tbx+j], tempTorque);
@@ -434,7 +435,7 @@ __device__ void computeBornChainRuleInteraction(AtomData3& atom1, AtomData3& ato
 }
 
 /**
- * Compute electrostatic interactions.
+ * Compute chain rule terms.
  */
 extern "C" __global__ void computeChainRuleForce(
         unsigned long long* __restrict__ forceBuffers, const real4* __restrict__ posq, unsigned int startTileIndex, unsigned int numTileIndices,
@@ -472,10 +473,11 @@ extern "C" __global__ void computeChainRuleForce(
                 localData[threadIdx.x].scaledRadius = data.scaledRadius;
                 localData[threadIdx.x].bornRadius = data.bornRadius;
                 localData[threadIdx.x].bornForce = data.bornForce;
+                localData[threadIdx.x].force = make_real3(0);
                 
                 // Compute forces.
                 
-                for (unsigned int j = 0; j < TILE_SIZE; j++) {
+                for (unsigned int j = (tgx+1)&(TILE_SIZE-1); j != tgx; j = (j+1)&(TILE_SIZE-1)) {
                     int atom2 = y*TILE_SIZE+j;
                     if (atom1 != atom2 && atom1 < NUM_ATOMS && atom2 < NUM_ATOMS) {
                         real3 tempForce;
@@ -484,9 +486,9 @@ extern "C" __global__ void computeChainRuleForce(
                         localData[tbx+j].force += tempForce;
                     }
                 }
-                atomicAdd(&forceBuffers[atom1], static_cast<unsigned long long>((long long) (data.force.x*0xFFFFFFFF)));
-                atomicAdd(&forceBuffers[atom1+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (data.force.y*0xFFFFFFFF)));
-                atomicAdd(&forceBuffers[atom1+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) (data.force.z*0xFFFFFFFF)));
+                atomicAdd(&forceBuffers[atom1], static_cast<unsigned long long>((long long) ((data.force.x+localData[threadIdx.x].force.x)*0xFFFFFFFF)));
+                atomicAdd(&forceBuffers[atom1+PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) ((data.force.y+localData[threadIdx.x].force.y)*0xFFFFFFFF)));
+                atomicAdd(&forceBuffers[atom1+2*PADDED_NUM_ATOMS], static_cast<unsigned long long>((long long) ((data.force.z+localData[threadIdx.x].force.z)*0xFFFFFFFF)));
             }
             else {
                 // This is an off-diagonal tile.

plugins/amoeba/platforms/cuda2/src/kernels/gkEDiffPairForce.cu

+#define APPLY_SCALE
+
 #if defined F1
 __device__ void computeOneEDiffInteractionF1(AtomData4& atom1, volatile AtomData4& atom2, float dScale, float pScale, real& outputEnergy, real3& outputForce) {
 #elif defined T1

plugins/amoeba/platforms/cuda2/src/kernels/gkPairForce.cu

@@ -1080,7 +1080,7 @@ __device__ void computeOneInteractionB1B2(AtomData2& atom1, volatile AtomData2&
 
 #if defined F1
 
-    outputEnergy += energy;
+    outputEnergy = energy;
 
     if ((xr != 0 || yr != 0 || zr != 0)) {
         force.x = dedx;