GBSAOBCForce allows the surface area energy to be changed

docs-source/usersguide/theory.rst

@@ -489,12 +489,12 @@ The surface area term is given by\ :cite:`Schaefer1998`\ :cite:`Ponder`
 
 
 .. math::
-   E=4\pi \cdot 2\text{.}\text{26}\sum _{i}{\left({r}_{i}+{r}_{\mathit{solvent}}\right)}^{2}{\left(\frac{{r}_{i}}{{R}_{i}}\right)}^{6}
+   E=E_{SA} \cdot 4\pi \sum _{i}{\left({r}_{i}+{r}_{\mathit{solvent}}\right)}^{2}{\left(\frac{{r}_{i}}{{R}_{i}}\right)}^{6}
 
 
 where :math:`r_i` is the atomic radius of particle *i*\ , :math:`r_i` is
-its Born radius, and :math:`r_\mathit{solvent}` is the solvent radius, which is taken
-to be 0.14 nm.
+its atomic radius, and :math:`r_\mathit{solvent}` is the solvent radius, which is taken
+to be 0.14 nm.  The default value for the energy scale :math:`E_{SA}` is 2.25936 kJ/mol/nm\ :sup:`2`\ .
 
 
 AndersenThermostat

openmmapi/include/openmm/GBSAOBCForce.h

@@ -138,6 +138,18 @@ public:
     void setSoluteDielectric(double dielectric) {
         soluteDielectric = dielectric;
     }
+    /**
+     * Get the energy scale for the surface energy term, measured in kJ/mol/nm^2.
+     */
+    double getSurfaceAreaEnergy() const {
+        return surfaceAreaEnergy;
+    }
+    /**
+     * Set the energy scale for the surface energy term, measured in kJ/mol/nm^2.
+     */
+    void setSurfaceAreaEnergy(double energy) {
+        surfaceAreaEnergy = energy;
+    }
     /**
      * Get the method used for handling long range nonbonded interactions.
      */
@@ -177,7 +189,7 @@ protected:
 private:
     class ParticleInfo;
     NonbondedMethod nonbondedMethod;
-    double cutoffDistance, solventDielectric, soluteDielectric;
+    double cutoffDistance, solventDielectric, soluteDielectric, surfaceAreaEnergy;
     std::vector<ParticleInfo> particles;
 };
 

openmmapi/src/GBSAOBCForce.cpp

@@ -37,7 +37,7 @@
 
 using namespace OpenMM;
 
-GBSAOBCForce::GBSAOBCForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), solventDielectric(78.3), soluteDielectric(1.0) {
+GBSAOBCForce::GBSAOBCForce() : nonbondedMethod(NoCutoff), cutoffDistance(1.0), solventDielectric(78.3), soluteDielectric(1.0), surfaceAreaEnergy(2.25936) {
 }
 
 int GBSAOBCForce::addParticle(double charge, double radius, double scalingFactor) {

platforms/cpu/include/CpuGBSAOBCForce.h

@@ -66,6 +66,11 @@ public:
      */
     void setSolventDielectric(float dielectric);
     
+    /**
+     * Set the surface area energy.
+     */
+    void setSurfaceAreaEnergy(float energy);
+    
     /**
      * Get the per-particle parameters (offset radius, scaled radius).
      */
@@ -96,7 +101,7 @@ private:
     bool cutoff;
     bool periodic;
     float periodicBoxSize[3];
-    float cutoffDistance, soluteDielectric, solventDielectric;
+    float cutoffDistance, soluteDielectric, solventDielectric, surfaceAreaFactor;
     std::vector<std::pair<float, float> > particleParams;        
     AlignedArray<float> bornRadii;
     std::vector<AlignedArray<float> > threadBornForces;

platforms/cpu/src/CpuGBSAOBCForce.cpp

@@ -77,6 +77,10 @@ void CpuGBSAOBCForce::setSolventDielectric(float dielectric) {
     solventDielectric = dielectric;
 }
 
+void CpuGBSAOBCForce::setSurfaceAreaEnergy(float energy) {
+    surfaceAreaFactor = 4*M_PI*energy;
+}
+
 const std::vector<std::pair<float, float> >& CpuGBSAOBCForce::getParticleParameters() const {
     return particleParams;
 }
@@ -211,7 +215,6 @@ void CpuGBSAOBCForce::threadComputeForce(ThreadPool& threads, int threadIndex) {
     // Calculate ACE surface area term.
 
     const float probeRadius = 0.14f;
-    const float surfaceAreaFactor = 28.3919551;
     double energy = 0.0;
     AlignedArray<float>& bornForces = threadBornForces[threadIndex];
     for (int i = 0; i < numParticles; i++)

platforms/cpu/src/CpuKernels.cpp

@@ -802,6 +802,7 @@ void CpuCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCFo
     obc.setParticleParameters(particleParams);
     obc.setSolventDielectric((float) force.getSolventDielectric());
     obc.setSoluteDielectric((float) force.getSoluteDielectric());
+    obc.setSurfaceAreaEnergy((float) force.getSurfaceAreaEnergy());
     if (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff)
         obc.setUseCutoff((float) force.getCutoffDistance());
     data.isPeriodic = (force.getNonbondedMethod() == GBSAOBCForce::CutoffPeriodic);

platforms/cpu/tests/TestCpuGBSAOBCForce.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -67,7 +67,7 @@ void testSingleParticle() {
     double eps0 = EPSILON0;
     double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/forceField->getSoluteDielectric()-1.0/forceField->getSolventDielectric())/bornRadius;
     double extendedRadius = 0.15+0.14; // probe radius
-    double nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.15/bornRadius, 6.0); // Where did this formula come from?  Just copied it from CpuImplicitSolvent.cpp
+    double nonpolarEnergy = 4*PI_M*2.25936*extendedRadius*extendedRadius*std::pow(0.15/bornRadius, 6.0);
     ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
     
     // Change the parameters and see if it is still correct.
@@ -77,8 +77,35 @@ void testSingleParticle() {
     state = context.getState(State::Energy);
     bornRadius = 0.25-0.009; // dielectric offset
     bornEnergy = (-0.4*0.4/(8*PI_M*eps0))*(1.0/forceField->getSoluteDielectric()-1.0/forceField->getSolventDielectric())/bornRadius;
-    extendedRadius = bornRadius+0.14;
-    nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.25/bornRadius, 6.0);
+    extendedRadius = 0.25+0.14;
+    nonpolarEnergy = 4*PI_M*2.25936*extendedRadius*extendedRadius*std::pow(0.25/bornRadius, 6.0);
+    ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
+}
+
+void testGlobalSettings() {
+    CpuPlatform platform;
+    System system;
+    system.addParticle(2.0);
+    LangevinIntegrator integrator(0, 0.1, 0.01);
+    GBSAOBCForce* forceField = new GBSAOBCForce();
+    forceField->addParticle(0.5, 0.15, 1);
+    const double soluteDielectric = 2.1;
+    const double solventDielectric = 35.0;
+    const double surfaceAreaEnergy = 0.75;
+    forceField->setSoluteDielectric(soluteDielectric);
+    forceField->setSolventDielectric(solventDielectric);
+    forceField->setSurfaceAreaEnergy(surfaceAreaEnergy);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(1);
+    positions[0] = Vec3(0, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Energy);
+    double bornRadius = 0.15-0.009; // dielectric offset
+    double eps0 = EPSILON0;
+    double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/soluteDielectric-1.0/solventDielectric)/bornRadius;
+    double extendedRadius = 0.15+0.14; // probe radius
+    double nonpolarEnergy = 4*PI_M*surfaceAreaEnergy*extendedRadius*extendedRadius*std::pow(0.15/bornRadius, 6.0);
     ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
 }
 
@@ -228,6 +255,7 @@ int main() {
             return 0;
         }
         testSingleParticle();
+        testGlobalSettings();
         testCutoffAndPeriodic();
         for (int i = 5; i < 11; i++) {
             testForce(i*i*i, NonbondedForce::NoCutoff, GBSAOBCForce::NoCutoff);

platforms/cuda/include/CudaKernels.h

@@ -719,7 +719,7 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force);
 private:
-    double prefactor;
+    double prefactor, surfaceAreaFactor;
     bool hasCreatedKernels;
     int maxTiles;
     CudaContext& cu;

platforms/cuda/src/CudaKernels.cpp

@@ -2480,6 +2480,7 @@ void CudaCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOBCF
     posq.upload(&temp[0]);
     params->upload(paramsVector);
     prefactor = -ONE_4PI_EPS0*((1.0/force.getSoluteDielectric())-(1.0/force.getSolventDielectric()));
+    surfaceAreaFactor = -6.0*4*M_PI*force.getSurfaceAreaEnergy();
     bool useCutoff = (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff);
     bool usePeriodic = (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff && force.getNonbondedMethod() != GBSAOBCForce::CutoffNonPeriodic);
     string source = CudaKernelSources::gbsaObc2;
@@ -2506,6 +2507,7 @@ double CudaCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool includeFor
         defines["CUTOFF_SQUARED"] = cu.doubleToString(nb.getCutoffDistance()*nb.getCutoffDistance());
         defines["CUTOFF"] = cu.doubleToString(nb.getCutoffDistance());
         defines["PREFACTOR"] = cu.doubleToString(prefactor);
+        defines["SURFACE_AREA_FACTOR"] = cu.doubleToString(surfaceAreaFactor);
         defines["NUM_ATOMS"] = cu.intToString(cu.getNumAtoms());
         defines["PADDED_NUM_ATOMS"] = cu.intToString(cu.getPaddedNumAtoms());
         defines["NUM_BLOCKS"] = cu.intToString(cu.getNumAtomBlocks());

platforms/cuda/src/kernels/gbsaObc1.cu

 #define DIELECTRIC_OFFSET 0.009f
 #define PROBE_RADIUS 0.14f
-#define SURFACE_AREA_FACTOR -170.351730667551f //-6.0f*3.14159265358979323846f*0.0216f*1000.0f*0.4184f;
 #define WARPS_PER_GROUP (FORCE_WORK_GROUP_SIZE/TILE_SIZE)
 
 /**

platforms/cuda/src/kernels/nonbonded.cu

@@ -398,6 +398,15 @@ extern "C" __global__ void computeNonbonded(
 #endif                
                 LOAD_LOCAL_PARAMETERS_FROM_GLOBAL
             }
+            else {
+#ifdef ENABLE_SHUFFLE
+                shflPosq = make_real4(0, 0, 0, 0);
+#else
+                localData[threadIdx.x].x = 0;
+                localData[threadIdx.x].y = 0;
+                localData[threadIdx.x].z = 0;
+#endif
+            }
 #ifdef USE_PERIODIC
             if (singlePeriodicCopy) {
                 // The box is small enough that we can just translate all the atoms into a single periodic

platforms/cuda/tests/TestCudaGBSAOBCForce.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -72,7 +72,7 @@ void testSingleParticle() {
     double eps0 = EPSILON0;
     double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/gbsa->getSoluteDielectric()-1.0/gbsa->getSolventDielectric())/bornRadius;
     double extendedRadius = 0.15+0.14; // probe radius
-    double nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.15/bornRadius, 6.0); // Where did this formula come from?  Just copied it from CpuImplicitSolvent.cpp
+    double nonpolarEnergy = 4*PI_M*2.25936*extendedRadius*extendedRadius*std::pow(0.15/bornRadius, 6.0);
     ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
     
     // Change the parameters and see if it is still correct.
@@ -82,8 +82,34 @@ void testSingleParticle() {
     state = context.getState(State::Energy);
     bornRadius = 0.25-0.009; // dielectric offset
     bornEnergy = (-0.4*0.4/(8*PI_M*eps0))*(1.0/gbsa->getSoluteDielectric()-1.0/gbsa->getSolventDielectric())/bornRadius;
-    extendedRadius = bornRadius+0.14;
-    nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.25/bornRadius, 6.0);
+    extendedRadius = 0.25+0.14;
+    nonpolarEnergy = 4*PI_M*2.25936*extendedRadius*extendedRadius*std::pow(0.25/bornRadius, 6.0);
+    ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
+}
+
+void testGlobalSettings() {
+    System system;
+    system.addParticle(2.0);
+    LangevinIntegrator integrator(0, 0.1, 0.01);
+    GBSAOBCForce* forceField = new GBSAOBCForce();
+    forceField->addParticle(0.5, 0.15, 1);
+    const double soluteDielectric = 2.1;
+    const double solventDielectric = 35.0;
+    const double surfaceAreaEnergy = 0.75;
+    forceField->setSoluteDielectric(soluteDielectric);
+    forceField->setSolventDielectric(solventDielectric);
+    forceField->setSurfaceAreaEnergy(surfaceAreaEnergy);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(1);
+    positions[0] = Vec3(0, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Energy);
+    double bornRadius = 0.15-0.009; // dielectric offset
+    double eps0 = EPSILON0;
+    double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/soluteDielectric-1.0/solventDielectric)/bornRadius;
+    double extendedRadius = 0.15+0.14; // probe radius
+    double nonpolarEnergy = 4*PI_M*surfaceAreaEnergy*extendedRadius*extendedRadius*std::pow(0.15/bornRadius, 6.0);
     ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
 }
 
@@ -229,6 +255,7 @@ int main(int argc, char* argv[]) {
         if (argc > 1)
             platform.setPropertyDefaultValue("CudaPrecision", string(argv[1]));
         testSingleParticle();
+        testGlobalSettings();
         testCutoffAndPeriodic();
         for (int i = 5; i < 11; i++) {
             testForce(i*i*i, NonbondedForce::NoCutoff, GBSAOBCForce::NoCutoff);

platforms/opencl/include/OpenCLKernels.h

@@ -721,7 +721,7 @@ public:
      */
     void copyParametersToContext(ContextImpl& context, const GBSAOBCForce& force);
 private:
-    double prefactor;
+    double prefactor, surfaceAreaFactor;
     bool hasCreatedKernels;
     int maxTiles;
     OpenCLContext& cl;

platforms/opencl/src/OpenCLKernels.cpp

@@ -2505,6 +2505,7 @@ void OpenCLCalcGBSAOBCForceKernel::initialize(const System& system, const GBSAOB
         cl.getPosq().upload(posqf);
     params->upload(paramsVector);
     prefactor = -ONE_4PI_EPS0*((1.0/force.getSoluteDielectric())-(1.0/force.getSolventDielectric()));
+    surfaceAreaFactor = -6.0*4*M_PI*force.getSurfaceAreaEnergy();
     bool useCutoff = (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff);
     bool usePeriodic = (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff && force.getNonbondedMethod() != GBSAOBCForce::CutoffNonPeriodic);
     string source = OpenCLKernelSources::gbsaObc2;
@@ -2530,6 +2531,7 @@ double OpenCLCalcGBSAOBCForceKernel::execute(ContextImpl& context, bool includeF
         defines["CUTOFF_SQUARED"] = cl.doubleToString(nb.getCutoffDistance()*nb.getCutoffDistance());
         defines["CUTOFF"] = cl.doubleToString(nb.getCutoffDistance());
         defines["PREFACTOR"] = cl.doubleToString(prefactor);
+        defines["SURFACE_AREA_FACTOR"] = cl.doubleToString(surfaceAreaFactor);
         defines["NUM_ATOMS"] = cl.intToString(cl.getNumAtoms());
         defines["PADDED_NUM_ATOMS"] = cl.intToString(cl.getPaddedNumAtoms());
         defines["NUM_BLOCKS"] = cl.intToString(cl.getNumAtomBlocks());

platforms/opencl/src/kernels/gbsaObcReductions.cl

 #define DIELECTRIC_OFFSET 0.009f
 #define PROBE_RADIUS 0.14f
-#define SURFACE_AREA_FACTOR -170.351730667551f //-6.0f*3.14159265358979323846f*0.0216f*1000.0f*0.4184f;
 
 /**
  * Reduce the Born sums to compute the Born radii.

platforms/opencl/src/kernels/nonbonded.cl

@@ -277,6 +277,11 @@ __kernel void computeNonbonded(
                 localData[localAtomIndex].fy = 0;
                 localData[localAtomIndex].fz = 0;
             }
+            else {
+                localData[localAtomIndex].x = 0;
+                localData[localAtomIndex].y = 0;
+                localData[localAtomIndex].z = 0;
+            }
             SYNC_WARPS;
 #ifdef USE_PERIODIC
             if (singlePeriodicCopy) {

platforms/opencl/tests/TestOpenCLGBSAOBCForce.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2008-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2008-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -72,7 +72,7 @@ void testSingleParticle() {
     double eps0 = EPSILON0;
     double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/gbsa->getSoluteDielectric()-1.0/gbsa->getSolventDielectric())/bornRadius;
     double extendedRadius = 0.15+0.14; // probe radius
-    double nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.15/bornRadius, 6.0); // Where did this formula come from?  Just copied it from CpuImplicitSolvent.cpp
+    double nonpolarEnergy = 4*PI_M*2.25936*extendedRadius*extendedRadius*std::pow(0.15/bornRadius, 6.0);
     ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
     
     // Change the parameters and see if it is still correct.
@@ -82,8 +82,34 @@ void testSingleParticle() {
     state = context.getState(State::Energy);
     bornRadius = 0.25-0.009; // dielectric offset
     bornEnergy = (-0.4*0.4/(8*PI_M*eps0))*(1.0/gbsa->getSoluteDielectric()-1.0/gbsa->getSolventDielectric())/bornRadius;
-    extendedRadius = bornRadius+0.14;
-    nonpolarEnergy = CAL2JOULE*PI_M*0.0216*(10*extendedRadius)*(10*extendedRadius)*std::pow(0.25/bornRadius, 6.0);
+    extendedRadius = 0.25+0.14;
+    nonpolarEnergy = 4*PI_M*2.25936*extendedRadius*extendedRadius*std::pow(0.25/bornRadius, 6.0);
+    ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
+}
+
+void testGlobalSettings() {
+    System system;
+    system.addParticle(2.0);
+    LangevinIntegrator integrator(0, 0.1, 0.01);
+    GBSAOBCForce* forceField = new GBSAOBCForce();
+    forceField->addParticle(0.5, 0.15, 1);
+    const double soluteDielectric = 2.1;
+    const double solventDielectric = 35.0;
+    const double surfaceAreaEnergy = 0.75;
+    forceField->setSoluteDielectric(soluteDielectric);
+    forceField->setSolventDielectric(solventDielectric);
+    forceField->setSurfaceAreaEnergy(surfaceAreaEnergy);
+    system.addForce(forceField);
+    Context context(system, integrator, platform);
+    vector<Vec3> positions(1);
+    positions[0] = Vec3(0, 0, 0);
+    context.setPositions(positions);
+    State state = context.getState(State::Energy);
+    double bornRadius = 0.15-0.009; // dielectric offset
+    double eps0 = EPSILON0;
+    double bornEnergy = (-0.5*0.5/(8*PI_M*eps0))*(1.0/soluteDielectric-1.0/solventDielectric)/bornRadius;
+    double extendedRadius = 0.15+0.14; // probe radius
+    double nonpolarEnergy = 4*PI_M*surfaceAreaEnergy*extendedRadius*extendedRadius*std::pow(0.15/bornRadius, 6.0);
     ASSERT_EQUAL_TOL((bornEnergy+nonpolarEnergy), state.getPotentialEnergy(), 0.01);
 }
 
@@ -227,6 +253,7 @@ void testForce(int numParticles, NonbondedForce::NonbondedMethod method, GBSAOBC
 int main() {
     try {
         testSingleParticle();
+        testGlobalSettings();
         testCutoffAndPeriodic();
         for (int i = 5; i < 11; i++) {
             testForce(i*i*i, NonbondedForce::NoCutoff, GBSAOBCForce::NoCutoff);

platforms/reference/include/ObcParameters.h

@@ -148,6 +148,14 @@ class ObcParameters {
 
       RealOpenMM getPi4Asolv( void ) const;
 
+      /**---------------------------------------------------------------------------------------
+      
+         Set pi4Asolv
+      
+         --------------------------------------------------------------------------------------- */
+
+      void setPi4Asolv( RealOpenMM pi4Asolv );
+
       /**---------------------------------------------------------------------------------------
       
          Get solvent dielectric

platforms/reference/src/ReferenceKernels.cpp

@@ -1098,6 +1098,7 @@ void ReferenceCalcGBSAOBCForceKernel::initialize(const System& system, const GBS
     obcParameters->setScaledRadiusFactors(scaleFactors);
     obcParameters->setSolventDielectric( static_cast<RealOpenMM>(force.getSolventDielectric()) );
     obcParameters->setSoluteDielectric( static_cast<RealOpenMM>(force.getSoluteDielectric()) );
+    obcParameters->setPi4Asolv(4*M_PI*force.getSurfaceAreaEnergy());
     if (force.getNonbondedMethod() != GBSAOBCForce::NoCutoff)
         obcParameters->setUseCutoff(static_cast<RealOpenMM>(force.getCutoffDistance()));
     isPeriodic = (force.getNonbondedMethod() == GBSAOBCForce::CutoffPeriodic);

platforms/reference/src/gbsa/ObcParameters.cpp

@@ -258,6 +258,9 @@ RealOpenMM ObcParameters::getPi4Asolv( void ) const {
     return _pi4Asolv;
 }
 
+void ObcParameters::setPi4Asolv(RealOpenMM pi4Asolv) {
+    _pi4Asolv = pi4Asolv;
+}
 
 /**---------------------------------------------------------------------------------------