Merge pull request #654 from peastman/sasa

GBSAOBCForce allows the surface area energy to be changed

platforms/reference/tests/TestReferenceGBSAOBCForce.cpp

@@ -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() {
+    ReferencePlatform 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);
 }
 
@@ -190,6 +217,7 @@ void testForce() {
 int main() {
     try {
         testSingleParticle();
+        testGlobalSettings();
         testCutoffAndPeriodic();
         testForce();
     }

serialization/src/GBSAOBCForceProxy.cpp

@@ -42,13 +42,14 @@ GBSAOBCForceProxy::GBSAOBCForceProxy() : SerializationProxy("GBSAOBCForce") {
 }
 
 void GBSAOBCForceProxy::serialize(const void* object, SerializationNode& node) const {
-    node.setIntProperty("version", 1);
+    node.setIntProperty("version", 2);
     const GBSAOBCForce& force = *reinterpret_cast<const GBSAOBCForce*>(object);
     node.setIntProperty("forceGroup", force.getForceGroup());
     node.setIntProperty("method", (int) force.getNonbondedMethod());
     node.setDoubleProperty("cutoff", force.getCutoffDistance());
     node.setDoubleProperty("soluteDielectric", force.getSoluteDielectric());
     node.setDoubleProperty("solventDielectric", force.getSolventDielectric());
+    node.setDoubleProperty("surfaceAreaEnergy", force.getSurfaceAreaEnergy());
     SerializationNode& particles = node.createChildNode("Particles");
     for (int i = 0; i < force.getNumParticles(); i++) {
         double charge, radius, scale;
@@ -58,7 +59,8 @@ void GBSAOBCForceProxy::serialize(const void* object, SerializationNode& node) c
 }
 
 void* GBSAOBCForceProxy::deserialize(const SerializationNode& node) const {
-    if (node.getIntProperty("version") != 1)
+    int version = node.getIntProperty("version");
+    if (version < 1 || version > 2)
         throw OpenMMException("Unsupported version number");
     GBSAOBCForce* force = new GBSAOBCForce();
     try {
@@ -67,6 +69,8 @@ void* GBSAOBCForceProxy::deserialize(const SerializationNode& node) const {
         force->setCutoffDistance(node.getDoubleProperty("cutoff"));
         force->setSoluteDielectric(node.getDoubleProperty("soluteDielectric"));
         force->setSolventDielectric(node.getDoubleProperty("solventDielectric"));
+        if (version > 1)
+            force->setSurfaceAreaEnergy(node.getDoubleProperty("surfaceAreaEnergy"));
         const SerializationNode& particles = node.getChildNode("Particles");
         for (int i = 0; i < (int) particles.getChildren().size(); i++) {
             const SerializationNode& particle = particles.getChildren()[i];

serialization/tests/TestSerializeGBSAOBCForce.cpp

@@ -47,6 +47,7 @@ void testSerialization() {
     force.setCutoffDistance(2.0);
     force.setSoluteDielectric(5.1);
     force.setSolventDielectric(50.0);
+    force.setSurfaceAreaEnergy(1.7);
     force.addParticle(1, 0.1, 0.01);
     force.addParticle(0.5, 0.2, 0.02);
     force.addParticle(-0.5, 0.3, 0.03);
@@ -65,6 +66,7 @@ void testSerialization() {
     ASSERT_EQUAL(force.getCutoffDistance(), force2.getCutoffDistance());
     ASSERT_EQUAL(force.getSoluteDielectric(), force2.getSoluteDielectric());
     ASSERT_EQUAL(force.getSolventDielectric(), force2.getSolventDielectric());
+    ASSERT_EQUAL(force.getSurfaceAreaEnergy(), force2.getSurfaceAreaEnergy());
     ASSERT_EQUAL(force.getNumParticles(), force2.getNumParticles());
     for (int i = 0; i < force.getNumParticles(); i++) {
         double charge1, radius1, scale1;

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -378,6 +378,7 @@ UNITS = {
 ("GBSAOBCForce", "getParticleParameters")
  : (None, ('unit.elementary_charge',
            'unit.nanometer', None)),
+("GBSAOBCForce", "getSurfaceAreaEnergy") : ('unit.kilojoule_per_mole/unit.nanometer/unit.nanometer', ()),
 ("GBVIForce", "getBornRadiusScalingMethod") : (None, ()),
 ("GBVIForce", "getQuinticLowerLimitFactor") : (None, ()),
 ("GBVIForce", "getQuinticUpperBornRadiusLimit") : ('unit.nanometer', ()),