Added missing lambda factor in second OBC force loop of free energy plugin

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

@@ -153,6 +153,9 @@ void SetCalculateObcGbsaSoftcoreBornSumSim( gpuContext gpu );
 extern "C"
 void SetCalculateObcGbsaSoftcoreNonPolarScalingFactorsSim( float* nonPolarScalingFactors );
 
+extern "C"
+void SetCalculateObcGbsaSoftcoreNonPolarScalingFactorsObc2Sim( float* nonPolarScalingFactors );
+
 // this method and kCalculateObcGbsaSoftcoreForces2() are being
 // used until changes in OpenMM version are made 
 extern "C"

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

@@ -59,6 +59,7 @@ int GpuObcGbsaSoftcore::setNonPolarScalingFactors( unsigned int particleIndex, f
 int GpuObcGbsaSoftcore::upload( gpuContext gpu ){
     _psNonPolarScalingFactors->Upload();
     SetCalculateObcGbsaSoftcoreNonPolarScalingFactorsSim( getGpuNonPolarScalingFactors() );
+    SetCalculateObcGbsaSoftcoreNonPolarScalingFactorsObc2Sim( getGpuNonPolarScalingFactors() );
     return 0;
 }
 

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

@@ -232,6 +232,7 @@ __global__ void kReduceObcGbsaSoftcoreBornSum_kernel()
         
         // Now calculate Born radius and OBC term.
         sum                    *= 0.5f * atom.x;
+        sum                    *= gbsaSimDev.pNonPolarScalingFactors[pos];
         float sum2              = sum * sum;
         float sum3              = sum * sum2;
         float tanhSum           = tanh(cSim.alphaOBC * sum - cSim.betaOBC * sum2 + cSim.gammaOBC * sum3);

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

@@ -34,6 +34,7 @@
 using namespace std;
 
 #include "gputypes.h"
+#include "GpuObcGbsaSoftcore.h"
 
 struct Atom {
     float x;
@@ -41,14 +42,20 @@ struct Atom {
     float z;
     float r;
     float sr;
+    float npScale;
     float fx;
     float fy;
     float fz;
     float fb;
 };
 
+struct cudaFreeEnergySimulationObcGbsaSoftcore {
+    float* pNonPolarScalingFactors;
+};
+struct cudaFreeEnergySimulationObcGbsaSoftcore gbsaSimObc2;
 
 static __constant__ cudaGmxSimulation cSim;
+static __constant__ cudaFreeEnergySimulationObcGbsaSoftcore gbsaSimDev;
 
 extern "C"
 void SetCalculateObcGbsaSoftcoreForces2Sim(gpuContext gpu)
@@ -58,6 +65,17 @@ void SetCalculateObcGbsaSoftcoreForces2Sim(gpuContext gpu)
     RTERROR(status, "cudaMemcpyToSymbol: SetSim copy to cSim failed");
 }
 
+extern "C" 
+void SetCalculateObcGbsaSoftcoreNonPolarScalingFactorsObc2Sim( float* nonPolarScalingFactors )
+{
+    cudaError_t status;
+    gbsaSimObc2.pNonPolarScalingFactors = nonPolarScalingFactors;
+    status                              = cudaMemcpyToSymbol(gbsaSimDev, &gbsaSimObc2, sizeof(cudaFreeEnergySimulationObcGbsaSoftcore));
+    RTERROR(status, "cudaMemcpyToSymbol: SetCalculateObcGbsaSoftcoreNonPolarScalingFactorsObc2Sim");
+
+    //(void) fprintf( stderr, "In SetCalculateObcGbsaSoftcoreNonPolarScalingFactorsObc2Sim\n" );
+}
+
 void GetCalculateObcGbsaSoftcoreForces2Sim(gpuContext gpu)
 {
     cudaError_t status;

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

@@ -55,6 +55,7 @@ __global__ void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned
         float4 apos                     = cSim.pPosq[i];
         float2 a                        = cSim.pObcData[i];
         float fb                        = cSim.pBornForce[i];
+        float  nonPolarScaleDataI       = gbsaSimDev.pNonPolarScalingFactors[i];
         unsigned int tbx                = threadIdx.x - tgx;
         unsigned int tj                 = tgx;
         Atom* psA                       = &sA[tbx];
@@ -72,6 +73,7 @@ __global__ void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned
             sA[threadIdx.x].r           = a.x;
             sA[threadIdx.x].sr          = a.y;
             sA[threadIdx.x].fb          = fb;
+            sA[threadIdx.x].npScale     = nonPolarScaleDataI;
 
             for (unsigned int j = (tgx+1)&(GRID-1); j != tgx; j = (j+1)&(GRID-1))
             {
@@ -105,6 +107,7 @@ __global__ void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned
                 float term          =  0.125f *
                                       (1.000f + psA[j].sr * psA[j].sr * r2Inverse) * t3 +
                                        0.250f * t1 * r2Inverse;
+                term               *= psA[j].npScale*nonPolarScaleDataI;
                 float dE            = fb * term;
 
 #if defined USE_PERIODIC
@@ -151,6 +154,7 @@ __global__ void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned
                 float4 temp                 = cSim.pPosq[j];
                 float2 temp1                = cSim.pObcData[j];
                 sA[threadIdx.x].fb          = cSim.pBornForce[j];
+                sA[threadIdx.x].npScale     = gbsaSimDev.pNonPolarScalingFactors[j];
                 sA[threadIdx.x].x           = temp.x;
                 sA[threadIdx.x].y           = temp.y;
                 sA[threadIdx.x].z           = temp.z;
@@ -208,6 +212,7 @@ __global__ void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned
                     float term              =  0.125f *
                                               (1.000f + psA[tj].sr * psA[tj].sr * r2Inverse) * t3J +
                                                0.250f * t1J * r2Inverse;
+                    term                   *= psA[tj].npScale*nonPolarScaleDataI;
                     float dE                = fb * term;
 
 #if defined USE_PERIODIC
@@ -235,6 +240,7 @@ __global__ void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned
                     term                    =  0.125f *
                                               (1.000f + sr2 * r2Inverse) * t3I +
                                                0.250f * t1I * r2Inverse;
+                    term                   *= psA[tj].npScale*nonPolarScaleDataI;
                     dE                      = psA[tj].fb * term;
 
                     float rj = psA[tj].r;
@@ -306,6 +312,7 @@ __global__ void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned
                         float term              =  0.125f *
                                                   (1.000f + psA[j].sr * psA[j].sr * r2Inverse) * t3J +
                                                    0.250f * t1J * r2Inverse;
+                        term                   *= psA[j].npScale*nonPolarScaleDataI;
                         float dE                = fb * term;
 
     #if defined USE_PERIODIC
@@ -333,6 +340,7 @@ __global__ void METHOD_NAME(kCalculateObcGbsaSoftcore, Forces2_kernel)(unsigned
                         term                    =  0.125f *
                                                   (1.000f + sr2 * r2Inverse) * t3I +
                                                    0.250f * t1I * r2Inverse;
+                        term                   *= psA[j].npScale*nonPolarScaleDataI;
                         dE                      = psA[j].fb * term;
 
                         float rj = psA[j].r;

plugins/freeEnergy/platforms/reference/src/gbsa/CpuObcSoftcore.cpp

@@ -241,7 +241,7 @@ void CpuObcSoftcore::computeBornRadii( vector<RealVec>& atomCoordinates,  std::v
 
 //FILE* logFile = SimTKOpenMMLog::getSimTKOpenMMLogFile( );
 //FILE* logFile = NULL;
-//FILE* logFile = fopen( "bR", "w" );
+//FILE* logFile = fopen( "bRSoft", "w" );
 
    for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
      
@@ -290,22 +290,13 @@ void CpuObcSoftcore::computeBornRadii( vector<RealVec>& atomCoordinates,  std::v
                   term += two*( radiusIInverse - l_ij);
                }
                sum += nonPolarScaleFactors[atomJ]*term;
-
-/*
-if( logFile && atomI == 0 ){
-   (void) fprintf( logFile, "\nRR %d %d r=%.4f rads[%.6f %.6f] scl=[%.3f %.3f] sum=%12.6e %12.6e %12.6e %12.6e",
-                   atomI, atomJ, r, offsetRadiusI, offsetRadiusJ, scaledRadiusFactor[atomI], scaledRadiusFactor[atomJ], 0.5f*sum,
-                   l_ij, u_ij, term );
-}
-*/
-
             }
          }
       }
  
       // OBC-specific code (Eqs. 6-8 in paper)
 
-      sum                  *= half*offsetRadiusI;
+      sum                  *= nonPolarScaleFactors[atomI]*half*offsetRadiusI;
       RealOpenMM sum2       = sum*sum;
       RealOpenMM sum3       = sum*sum2;
       RealOpenMM tanhSum    = TANH( alphaObc*sum - betaObc*sum2 + gammaObc*sum3 );
@@ -317,7 +308,7 @@ if( logFile && atomI == 0 ){
 
 #if 0
 if( logFile && atomI >= 0 ){
-   (void) fprintf( logFile, "\nRRQ %d sum %12.6e tanhS %12.6e radI %.5f %.5f born %18.10e obc %12.6e",
+   (void) fprintf( logFile, "RRQ %d sum %12.6e tanhS %12.6e radI %.5f %.5f born %18.10e obc %12.6e\n",
                    atomI, sum, tanhSum, radiusI, offsetRadiusI, bornRadii[atomI], obcChain[atomI] );
 }
 #endif
@@ -327,8 +318,15 @@ if( logFile && atomI >= 0 ){
 #if 0
 if( logFile ){
       (void) fclose( logFile );
+      logFile = fopen( "bRSoftJ", "w" );
+      for( int atomI = 0; atomI < numberOfAtoms; atomI++ ){
+          (void) fprintf( logFile, "%6d %18.10e %18.10e\n",
+                          atomI, bornRadii[atomI], obcChain[atomI] );
+      }
+      (void) fclose( logFile );
 }
 #endif
+
 }
 
 /**---------------------------------------------------------------------------------------
@@ -476,6 +474,8 @@ void CpuObcSoftcore::computeBornEnergyForces( vector<RealOpenMM>& bornRadii, vec
       computeAceNonPolarForce( obcSoftcoreParameters, bornRadii, &obcEnergy, bornForces );
    }
  
+   const RealOpenMM* nonPolarScaleFactors   = obcSoftcoreParameters->getNonPolarScaleFactors();
+
    // ---------------------------------------------------------------------------------------
 
    // first main loop
@@ -497,27 +497,23 @@ void CpuObcSoftcore::computeBornEnergyForces( vector<RealOpenMM>& bornRadii, vec
          RealOpenMM deltaY                 = deltaR[ReferenceForce::YIndex];
          RealOpenMM deltaZ                 = deltaR[ReferenceForce::ZIndex];
 
-         // 3 FLOP
-
          RealOpenMM alpha2_ij              = bornRadii[atomI]*bornRadii[atomJ];
          RealOpenMM D_ij                   = r2/(four*alpha2_ij);
 
-         // exp + 2 + sqrt FLOP 
-
          RealOpenMM expTerm                = EXP( -D_ij );
          RealOpenMM denominator2           = r2 + alpha2_ij*expTerm; 
          RealOpenMM denominator            = SQRT( denominator2 ); 
 
-         // 6 FLOP
+         // charges are assumed to be scaled on input so nonPolarScaleFactor is not needed
 
-         RealOpenMM Gpol               = (partialChargeI*partialCharges[atomJ])/denominator; 
-         RealOpenMM dGpol_dr           = -Gpol*( one - fourth*expTerm )/denominator2;  
+         //RealOpenMM nonPolarScaleFactor    = atomI != atomJ ? nonPolarScaleFactors[atomI]*nonPolarScaleFactors[atomJ] : nonPolarScaleFactors[atomI];
 
-         // 5 FLOP
+         RealOpenMM Gpol                   = (partialChargeI*partialCharges[atomJ])/denominator; 
+         //           Gpol                  *= nonPolarScaleFactor;
 
          RealOpenMM dGpol_dalpha2_ij       = -half*Gpol*expTerm*( one + D_ij )/denominator2;
 
-         // 11 FLOP
+         RealOpenMM dGpol_dr               = -Gpol*( one - fourth*expTerm )/denominator2;  
 
          if( atomI != atomJ ){
 
@@ -547,8 +543,6 @@ void CpuObcSoftcore::computeBornEnergyForces( vector<RealOpenMM>& bornRadii, vec
       }
    }
 
-   //obcEnergy *= getEnergyConversionFactor();
-
    // ---------------------------------------------------------------------------------------
 
    // second main loop
@@ -630,6 +624,7 @@ void CpuObcSoftcore::computeBornEnergyForces( vector<RealOpenMM>& bornRadii, vec
                RealOpenMM r2Inverse     = rInverse*rInverse;
 
                RealOpenMM t3            = eighth*(one + scaledRadiusJ2*r2Inverse)*(l_ij2 - u_ij2) + fourth*LN( u_ij/l_ij )*r2Inverse;
+                          t3           *= nonPolarScaleFactors[atomI]*nonPolarScaleFactors[atomJ];
 
                RealOpenMM de            = bornForces[atomI]*t3*rInverse;