Implemented surface area term for GK

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

@@ -1807,6 +1807,12 @@ void CudaCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System& syst
     defines["ENERGY_SCALE_FACTOR"] = cu.doubleToString(138.9354558456/force.getSoluteDielectric());
     if (multipoles->getPolarizationType() == AmoebaMultipoleForce::Direct)
         defines["DIRECT_POLARIZATION"] = "";
+    includeSurfaceArea = force.getIncludeCavityTerm();
+    if (includeSurfaceArea) {
+        defines["SURFACE_AREA_FACTOR"] = cu.doubleToString(force.getSurfaceAreaFactor());
+        defines["PROBE_RADIUS"] = cu.doubleToString(force.getProbeRadius());
+        defines["DIELECTRIC_OFFSET"] = cu.doubleToString(0.009);
+    }
     stringstream forceSource;
     forceSource << CudaKernelSources::vectorOps;
     forceSource << CudaAmoebaKernelSources::amoebaGk;
@@ -1836,6 +1842,8 @@ void CudaCalcAmoebaGeneralizedKirkwoodForceKernel::initialize(const System& syst
     gkForceKernel = cu.getKernel(module, "computeGKForces");
     chainRuleKernel = cu.getKernel(module, "computeChainRuleForce");
     ediffKernel = cu.getKernel(module, "computeEDiffForce");
+    if (includeSurfaceArea)
+        surfaceAreaKernel = cu.getKernel(module, "computeSurfaceAreaForce");
     cu.addForce(new ForceInfo(force));
 }
 
@@ -1864,13 +1872,23 @@ void CudaCalcAmoebaGeneralizedKirkwoodForceKernel::finishComputation(CudaArray&
     int numTileIndices = nb.getNumTiles();
     int numForceThreadBlocks = nb.getNumForceThreadBlocks();
     int forceThreadBlockSize = nb.getForceThreadBlockSize();
+    
+    // Compute the GK force.
+    
     void* gkForceArgs[] = {&cu.getForce().getDevicePointer(), &torque.getDevicePointer(), &cu.getEnergyBuffer().getDevicePointer(),
         &cu.getPosq().getDevicePointer(), &startTileIndex, &numTileIndices, &labFrameDipoles.getDevicePointer(),
         &labFrameQuadrupoles.getDevicePointer(), &inducedDipoleS->getDevicePointer(), &inducedDipolePolarS->getDevicePointer(),
         &bornRadii->getDevicePointer(), &bornForce->getDevicePointer()};
     cu.executeKernel(gkForceKernel, gkForceArgs, numForceThreadBlocks*forceThreadBlockSize, forceThreadBlockSize);
 
-    // Compute cavity term...
+    // Compute the surface area force.
+    
+    if (includeSurfaceArea) {
+        void* surfaceAreaArgs[] = {&bornForce->getDevicePointer(), &cu.getEnergyBuffer().getDevicePointer(), &params->getDevicePointer(), &bornRadii->getDevicePointer()};
+        cu.executeKernel(surfaceAreaKernel, surfaceAreaArgs, cu.getNumAtoms());
+    }
+    
+    // Apply the remaining terms.
     
     void* chainRuleArgs[] = {&cu.getForce().getDevicePointer(), &cu.getPosq().getDevicePointer(), &startTileIndex, &numTileIndices,
         &params->getDevicePointer(), &bornRadii->getDevicePointer(), &bornForce->getDevicePointer()};

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

@@ -486,6 +486,7 @@ private:
     class ForceInfo;
     CudaContext& cu;
     System& system;
+    bool includeSurfaceArea;
     CudaArray* params;
     CudaArray* bornSum;
     CudaArray* bornRadii;
@@ -495,7 +496,7 @@ private:
     CudaArray* inducedFieldPolar;
     CudaArray* inducedDipoleS;
     CudaArray* inducedDipolePolarS;
-    CUfunction computeBornSumKernel, reduceBornSumKernel, gkForceKernel, chainRuleKernel, ediffKernel;
+    CUfunction computeBornSumKernel, reduceBornSumKernel, surfaceAreaKernel, gkForceKernel, chainRuleKernel, ediffKernel;
 };
 
 /**

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

@@ -20,6 +20,27 @@ extern "C" __global__ void reduceBornSum(const long long* __restrict__ bornSum,
     }
 }
 
+#ifdef SURFACE_AREA_FACTOR
+/**
+ * Apply the surface area term to the force and energy.
+ */
+extern "C" __global__ void computeSurfaceAreaForce(long long* __restrict__ bornForce, real* __restrict__ energyBuffer, const float2* __restrict__ params, const real* __restrict__ bornRadii) {
+    real energy = 0;
+    for (unsigned int index = blockIdx.x*blockDim.x+threadIdx.x; index < NUM_ATOMS; index += blockDim.x*gridDim.x) {
+        real bornRadius = bornRadii[index];
+        float radius = params[index].x;
+        real r = radius + DIELECTRIC_OFFSET + PROBE_RADIUS;
+        real ratio6 = (radius+DIELECTRIC_OFFSET)/bornRadius;
+        ratio6 = ratio6*ratio6*ratio6;
+        ratio6 = ratio6*ratio6;
+        real saTerm = SURFACE_AREA_FACTOR * r * r * ratio6;
+        bornForce[index] += (long long) (saTerm*0xFFFFFFFF/bornRadius);
+        energy += saTerm;
+    }
+    energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] -= energy/6;
+}
+#endif
+
 /**
  * Data structure used by computeBornSum().
  */

plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaGeneralizedKirkwoodForce.cpp

@@ -358,9 +358,6 @@ static void compareForcesEnergy( std::string& testName, double expectedEnergy, d
     }
 #endif
 
-    for( unsigned int ii = 0; ii < forces.size(); ii++ )
-        std::cout << forces[ii]<<" "<<expectedForces[ii]<< std::endl;
-    std::cout << energy<<" "<<expectedEnergy<< std::endl;
     for( unsigned int ii = 0; ii < forces.size(); ii++ ){
         ASSERT_EQUAL_VEC_MOD( expectedForces[ii], forces[ii], tolerance, testName );
     }