Disable reaction field when using implicit solvent (see bug 1881)

openmmapi/include/openmm/GBSAOBCForce.h

@@ -40,13 +40,19 @@ namespace OpenMM {
 
 /**
  * This class implements an implicit solvation force using the GBSA-OBC model.
- * <p>
+ * 
  * To use this class, create a GBSAOBCForce object, then call addParticle() once for each particle in the
  * System to define its parameters.  The number of particles for which you define GBSA parameters must
  * be exactly equal to the number of particles in the System, or else an exception will be thrown when you
  * try to create a Context.  After a particle has been added, you can modify its force field parameters
  * by calling setParticleParameters().  This will have no effect on Contexts that already exist unless you
  * call updateParametersInContext().
+ * 
+ * When using this Force, the System should also include a NonbondedForce, and both objects must specify
+ * identical charges for all particles.  Otherwise, the results will not be correct.  Furthermore, if the
+ * nonbonded method is set to CutoffNonPeriodic or CutoffPeriodic, you should call setReactionFieldDielectric(1.0)
+ * on the NonbondedForce to turn off the reaction field approximation, which does not produce correct results
+ * when combined with GBSA.
  */
 
 class OPENMM_EXPORT GBSAOBCForce : public Force {
@@ -70,7 +76,7 @@ public:
          */
         CutoffPeriodic = 2,
     };
-    /*
+    /**
      * Create a GBSAOBCForce.
      */
     GBSAOBCForce();

wrappers/python/simtk/openmm/app/forcefield.py

@@ -1046,6 +1046,13 @@ class GBSAOBCGenerator:
             force.setSolventDielectric(float(args['solventDielectric']))
         sys.addForce(force)
 
+    def postprocessSystem(self, sys, data, args):
+        # Disable the reaction field approximation, since it produces bad results when combined with GB.
+        
+        for force in sys.getForces():
+            if isinstance(force, mm.NonbondedForce):
+                force.setReactionFieldDielectric(1.0)
+
 parsers["GBSAOBCForce"] = GBSAOBCGenerator.parseElement
 
 

wrappers/python/simtk/openmm/app/gromacstopfile.py

@@ -448,6 +448,7 @@ class GromacsTopFile(object):
             gb.setSoluteDielectric(soluteDielectric)
             gb.setSolventDielectric(solventDielectric)
             sys.addForce(gb)
+            nb.setReactionFieldDielectric(1.0)
         elif implicitSolvent is not None:
             raise ValueError('Illegal value for implicitSolvent')
         bonds = None

wrappers/python/simtk/openmm/app/internal/amber_file_parser.py

@@ -824,6 +824,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
             gb.setCutoffDistance(nonbondedCutoff)
         else:
             raise Exception("Illegal nonbonded method for use with GBSA")
+        force.setReactionFieldDielectric(1.0)
 
     # TODO: Add GBVI terms?