NonPolarPrefactor for GBSA OBC was missing for CudaPlatform

plugins/freeEnergy/platforms/cuda/src/CudaFreeEnergyKernels.cpp

@@ -61,7 +61,7 @@ const std::string GBSA_OBC_SOFTCORE_FORCE         = "ObcSoftcore";
 const std::string GBVI_FORCE                      = "GBVI";
 const std::string GBVI_SOFTCORE_FORCE             = "GBVISoftcore";
 
-static void getForceMap(const System& system, MapStringInt& forceMap) { 
+static void getForceMap(const System& system, MapStringInt& forceMap, FILE* log) { 
 
     // check forces and relevant parameters
 
@@ -69,6 +69,7 @@ static void getForceMap(const System& system, MapStringInt& forceMap) {
 
         int hit                 = 0;
         const Force& force      = system.getForce(i);
+         std::string forceName  = "NA";
 
         // bond
 
@@ -77,6 +78,7 @@ static void getForceMap(const System& system, MapStringInt& forceMap) {
             try {
                const HarmonicBondForce& harmonicBondForce = dynamic_cast<const HarmonicBondForce&>(force);
                forceMap[HARMONIC_BOND_FORCE]              = 1;
+               forceName                                  = HARMONIC_BOND_FORCE;
                hit++;
             } catch( std::bad_cast ){
             }
@@ -89,6 +91,7 @@ static void getForceMap(const System& system, MapStringInt& forceMap) {
             try {
                const HarmonicAngleForce& harmonicAngleForce = dynamic_cast<const HarmonicAngleForce&>(force);
                forceMap[HARMONIC_ANGLE_FORCE]               = 1;
+               forceName                                    = HARMONIC_ANGLE_FORCE;
                hit++;
             } catch( std::bad_cast ){
             }
@@ -101,6 +104,7 @@ static void getForceMap(const System& system, MapStringInt& forceMap) {
             try {
                const PeriodicTorsionForce & periodicTorsionForce = dynamic_cast<const PeriodicTorsionForce&>(force);
                forceMap[PERIODIC_TORSION_FORCE]                  = 1;
+               forceName                                         = PERIODIC_TORSION_FORCE;
                hit++;
             } catch( std::bad_cast ){
             }
@@ -112,6 +116,7 @@ static void getForceMap(const System& system, MapStringInt& forceMap) {
             try {
                const RBTorsionForce& rBTorsionForce = dynamic_cast<const RBTorsionForce&>(force);
                forceMap[RB_TORSION_FORCE]           = 1;
+               forceName                            = RB_TORSION_FORCE;
                hit++;
             } catch( std::bad_cast ){
             }
@@ -123,6 +128,7 @@ static void getForceMap(const System& system, MapStringInt& forceMap) {
             try {
                const NonbondedForce& nbForce = dynamic_cast<const NonbondedForce&>(force);
                forceMap[NB_FORCE]            = 1;
+               forceName                     = NB_FORCE;
             } catch( std::bad_cast ){
             }
         }
@@ -133,6 +139,7 @@ static void getForceMap(const System& system, MapStringInt& forceMap) {
             try {
                const NonbondedSoftcoreForce& nbForce = dynamic_cast<const NonbondedSoftcoreForce&>(force);
                forceMap[NB_SOFTCORE_FORCE]           = 1;
+               forceName                             = NB_SOFTCORE_FORCE;
             } catch( std::bad_cast ){
             }
         }
@@ -143,6 +150,7 @@ static void getForceMap(const System& system, MapStringInt& forceMap) {
             try {
                const GBSAOBCForce& obcForce       = dynamic_cast<const GBSAOBCForce&>(force);
                forceMap[GBSA_OBC_FORCE]           = 1;
+               forceName                          = GBSA_OBC_FORCE;
                hit++;
             } catch( std::bad_cast ){
             }
@@ -154,6 +162,7 @@ static void getForceMap(const System& system, MapStringInt& forceMap) {
             try {
                const GBSAOBCSoftcoreForce& obcForce = dynamic_cast<const GBSAOBCSoftcoreForce&>(force);
                forceMap[GBSA_OBC_SOFTCORE_FORCE]    = 1;
+               forceName                            = GBSA_OBC_SOFTCORE_FORCE;
                hit++;
             } catch( std::bad_cast ){
             }
@@ -164,7 +173,8 @@ static void getForceMap(const System& system, MapStringInt& forceMap) {
         if( !hit ){
             try {
                const GBVIForce& obcForce  = dynamic_cast<const GBVIForce&>(force);
-               forceMap[GBVI_FORCE] = 1;
+               forceMap[GBVI_FORCE]       = 1;
+               forceName                  = GBVI_FORCE;
                hit++;
             } catch( std::bad_cast ){
             }
@@ -176,10 +186,15 @@ static void getForceMap(const System& system, MapStringInt& forceMap) {
             try {
                const GBVISoftcoreForce& gbviForce = dynamic_cast<const GBVISoftcoreForce&>(force);
                forceMap[GBVI_SOFTCORE_FORCE]      = 1;
+               forceName                          = GBVI_SOFTCORE_FORCE;
                hit++;
             } catch( std::bad_cast ){
             }
         }
+
+        if( log ){
+            (void) fprintf( stderr, "Map: Force %d %s\n", i, forceName.c_str() );
+        }
      }
 }
 
@@ -208,7 +223,7 @@ void CudaFreeEnergyCalcNonbondedSoftcoreForceKernel::initialize(const System& sy
     // check forces and relevant parameters
 
     MapStringInt forceMap;
-    getForceMap( system, forceMap);
+    getForceMap( system, forceMap, log);
 
     int softcore        = 0;
     if( forceMap.find( GBSA_OBC_FORCE ) != forceMap.end() ){
@@ -470,13 +485,16 @@ void CudaFreeEnergyCalcGBSAOBCSoftcoreForceKernel::initialize(const System& syst
     _gpuContext* gpu = data.gpu;
 
     MapStringInt forceMap;
-    getForceMap( system, forceMap);
+    getForceMap( system, forceMap, log);
 
     // check that nonbonded (non-softcore is not active)
 
     if( forceMap.find( NB_FORCE ) != forceMap.end() ){ 
         throw OpenMMException( "Mixing NonbondedForce and GBSAOBCSoftoreForce not allowed -- use NonbondedSoftcoreForce " );
     }
+    if( forceMap.find( NB_SOFTCORE_FORCE ) == forceMap.end() ){ 
+        throw OpenMMException( "NonbondedSoftcore force must be included w/ GBSAOBCSoftcore force." );
+    }
 
     int numParticles = system.getNumParticles();
 
@@ -495,6 +513,7 @@ void CudaFreeEnergyCalcGBSAOBCSoftcoreForceKernel::initialize(const System& syst
     }
     gpuObcGbsaSoftcore = gpuSetObcSoftcoreParameters(gpu, static_cast<float>( force.getSoluteDielectric()),
                                                      static_cast<float>( force.getSolventDielectric()),
+                                                     static_cast<float>( force.getNonPolarPrefactor()),
                                                      radius, scale, charge, nonPolarScalingFactors  );
 }
 
@@ -585,7 +604,7 @@ void CudaFreeEnergyCalcGBVISoftcoreForceKernel::initialize(const System& system,
     // check forces and relevant parameters
 
     MapStringInt forceMap;
-    getForceMap( system, forceMap);
+    getForceMap( system, forceMap, log);
 
     // check that nonbonded (non-softcore is not active)
 

plugins/freeEnergy/platforms/cuda/src/kernels/GpuFreeEnergyCudaKernels.h

@@ -127,7 +127,8 @@ extern void kReduceGBVIBornForcesQuinticScaling( gpuContext gpu );
  */
 
 extern "C" 
-GpuObcGbsaSoftcore* gpuSetObcSoftcoreParameters(gpuContext gpu, float innerDielectric, float solventDielectric, const std::vector<float>& radius, const std::vector<float>& scale,
+GpuObcGbsaSoftcore* gpuSetObcSoftcoreParameters(gpuContext gpu, float innerDielectric, float solventDielectric, float nonPolarPrefactor,
+                                                const std::vector<float>& radius, const std::vector<float>& scale,
                                                 const std::vector<float>& charge, const std::vector<float>& nonPolarScalingFactors);
 
 // delete supplemtentary objects, ...

plugins/freeEnergy/platforms/cuda/src/kernels/GpuNonbondedSoftcore.cpp

@@ -71,7 +71,24 @@ int GpuNonbondedSoftcore::setParticleSoftCoreLJLambda( unsigned int particleInde
 // upload SoftCoreLJLambda array
 
 int GpuNonbondedSoftcore::upload( gpuContext gpu ){
+
+// ---------------------------------------------------------------------------------------
+
+   static const std::string methodName    = "GpuNonbondedSoftcore::upload";
+
+// ---------------------------------------------------------------------------------------
+
     _psSoftcoreLJLambda->Upload();
+
+#define DUMP_PARAMETERS 0
+#if (DUMP_PARAMETERS == 1)
+    (void) fprintf( stderr, "%s %u %u\n", methodName.c_str(), gpu->natoms, gpu->sim.paddedNumberOfAtoms );
+    for (unsigned int ii = 0; ii < gpu->natoms; ii++)
+    {   
+       (void) fprintf( stderr, "%6u %13.6e\n", ii,  (*_psSoftcoreLJLambda)[ii] );
+    }
+#endif
+
     SetCalculateCDLJSoftcoreSupplementarySim( getGpuParticleSoftCoreLJLambda() );
     SetCalculateCDLJObcGbsaSoftcoreSupplementary1Sim( getGpuParticleSoftCoreLJLambda() );
     return 0;

plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateCDLJObcGbsaSoftcoreForces1.cu

@@ -121,6 +121,7 @@ __device__ float fastErfc(float r)
 #define METHOD_NAME(a, b) a##N2ByWarp##b
 #include "kCalculateCDLJObcGbsaSoftcoreForces1.h"
 #undef USE_SOFTCORE_LJ
+#undef USE_OUTPUT_BUFFER_PER_WARP
 
 // Include versions of the kernel with cutoffs.
 
@@ -196,7 +197,7 @@ void kCalculateCDLJObcGbsaSoftcoreForces1(gpuContext gpu )
                    kCalculateCDLJObcGbsaSoftcoreN2Forces1_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
                            sizeof(Atom)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pWorkUnit);
    
-            LAUNCHERROR("kCalculateCDLJObcGbsaN2Forces1");
+            LAUNCHERROR("kCalculateCDLJObcGbsaSoftcoreForces1");
             break;
 #if 0
         case CUTOFF:

plugins/freeEnergy/platforms/cuda/src/kernels/kCalculateObcGbsaSoftcoreBornSum.cu

@@ -176,9 +176,11 @@ void kReduceObcGbsaBornSum(gpuContext gpu)
  */
 
 extern "C"
-void gpuSetObcSoftcoreParameters(gpuContext gpu, float innerDielectric, float solventDielectric, const std::vector<float>& radius, const std::vector<float>& scale,
+void gpuSetObcSoftcoreParameters(gpuContext gpu, float innerDielectric, float solventDielectric, float nonPolarPrefactor,
+                                 const std::vector<float>& radius, const std::vector<float>& scale,
                                  const std::vector<float>& charge, const std::vector<float>& nonPolarScalingFactors)
 {
+
 // ---------------------------------------------------------------------------------------
 
    static const float dielectricOffset    =    0.009f;
@@ -187,7 +189,6 @@ void gpuSetObcSoftcoreParameters(gpuContext gpu, float innerDielectric, float so
 
 // ---------------------------------------------------------------------------------------
 
-
     unsigned int atoms                     = radius.size();
 
     // initialize parameters
@@ -195,6 +196,7 @@ void gpuSetObcSoftcoreParameters(gpuContext gpu, float innerDielectric, float so
 //    gpu->bIncludeGBSA = true;
     GpuObcGbsaSoftcore* gpuObcGbsaSoftcore = new GpuObcGbsaSoftcore();
     gpuObcGbsaSoftcore->initializeNonPolarScalingFactors( gpu->sim.paddedNumberOfAtoms );
+    gpu->sim.surfaceAreaFactor             =  -6.0f*PI*4.0f*nonPolarPrefactor*1000.0f*0.4184f;
     for (unsigned int i = 0; i < atoms; i++)
     {
             (*gpu->psObcData)[i].x = radius[i] - dielectricOffset;
@@ -205,12 +207,13 @@ void gpuSetObcSoftcoreParameters(gpuContext gpu, float innerDielectric, float so
     }
 
     // diagnostics
-
+#define DUMP_PARAMETERS 0
 #if (DUMP_PARAMETERS == 1)
     (void) fprintf( stderr, "%s %u %u\n", methodName.c_str(), gpu->natoms, gpu->sim.paddedNumberOfAtoms );
     for (unsigned int i = 0; i < atoms; i++)
     {
        (void) fprintf( stderr, "%6u %13.6e %13.6e %8.3f %8.3f\n", i,  (*gpu->psObcData)[i].x, (*gpu->psObcData)[i].y, (*gpu->psPosq4)[i].w , nonPolarScalingFactors[i] );
+    }
 #endif
 
     // dummy out extra atom data