Implicit solvent for modern force fields (#3214)

* Created OBC2 that works with current force fields

* Created HCT, OBC1, GBn, and GBn2 that works with current force fields

* Added documentation for GB models

* Updates to documentation and tests based on comments

* Added formula for screening parameter

docs-source/usersguide/application/02_running_sims.rst

@@ -593,6 +593,59 @@ such as :file:`charmm36/water.xml`, which specifies the default CHARMM water mod
 The converted parameter sets come from the `CHARMM36 July 2017 update <http://mackerell.umaryland.edu/charmm_ff.shtml>`_
 and were converted using the `openmm-forcefields <https://github.com/choderalab/openmm-forcefields>`_ package and `parmed <https://github.com/parmed/parmed>`_.
 
+Implicit Solvent
+----------------
+
+The Amber and CHARMM force fields described above can be used with any of the Generalized
+Born implicit solvent models from AMBER.  To use them, include an extra file when
+creating the ForceField.  For example,
+::
+
+    forcefield = ForceField('amber14-all.xml', 'implicit/gbn2.xml')
+
+.. tabularcolumns:: |l|L|
+
+==========================  ==================================================================================================================================
+File                        Implicit Solvent Model
+==========================  ==================================================================================================================================
+:file:`implicit/hct.xml`    Hawkins-Cramer-Truhlar GBSA model\ :cite:`Hawkins1995` (corresponds to igb=1 in AMBER)
+:file:`implicit/obc1.xml`   Onufriev-Bashford-Case GBSA model\ :cite:`Onufriev2004` using the GB\ :sup:`OBC`\ I parameters (corresponds to igb=2 in AMBER).
+:file:`implicit/obc2.xml`   Onufriev-Bashford-Case GBSA model\ :cite:`Onufriev2004` using the GB\ :sup:`OBC`\ II parameters (corresponds to igb=5 in AMBER).
+:file:`implicit/gbn.xml`    GBn solvation model\ :cite:`Mongan2007` (corresponds to igb=7 in AMBER).
+:file:`implicit/gbn2.xml`   GBn2 solvation model\ :cite:`Nguyen2013` (corresponds to igb=8 in AMBER).
+==========================  ==================================================================================================================================
+
+The only nonbonded methods that are supported with implicit solvent are :code:`NoCutoff` (the default),
+:code:`CutoffNonPeriodic`, and :code:`CutoffPeriodic.`  If you choose to use a nonbonded cutoff with
+implicit solvent, it is usually best to set the cutoff distance larger than is typical with explicit solvent.
+A cutoff of 2 nm gives good results in most cases.  Periodic boundary conditions are not usually used
+with implicit solvent.  In fact, the lack of need for periodicity and the artifacts it creates is one
+of the advantages of implicit solvent.  The option is still offered, since it could be useful in some
+unusual situations.
+
+You can further control the solvation model in a few ways.  First, you can
+specify the dielectric constants to use for the solute and solvent:
+::
+
+    system = forcefield.createSystem(topology, soluteDielectric=1.0, solventDielectric=80.0)
+
+If they are not specified, the solute and solvent dielectric constants default to 1.0 and
+78.5, respectively.
+
+You also can model the effect of a non-zero salt concentration by specifying the
+Debye-Huckel screening parameter\ :cite:`Srinivasan1999`:
+::
+
+    system = forcefield.createSystem(topology, implicitSolventKappa=1.0/nanometer)
+
+The screening parameter can be calculated as
+
+.. math::
+  \kappa = 367.434915 \sqrt{\frac{I}{\epsilon T}}
+
+where :math:`I` is the ionic strength in moles/liter, :math:`\epsilon` is the solvent
+dielectric constant, and :math:`T` is the temperature in Kelvin.
+
 AMOEBA
 ------
 

wrappers/python/openmm/app/data/implicit/gbn.xml

+<ForceField>
+ <Script>
+import openmm
+import openmm.app as app
+import openmm.unit as unit
+
+# Find the NonbondedForce.  We need it to look up charges, and also to set the reaction field dielectric to 1.
+
+nonbonded = [f for f in sys.getForces() if isinstance(f, openmm.NonbondedForce)]
+if len(nonbonded) != 1:
+    raise ValueError('Implicit solvent requires the System to contain a single NonbondedForce')
+nonbonded = nonbonded[0]
+nonbonded.setReactionFieldDielectric(1)
+
+# Construct the CustomGBForce.
+
+from openmm.app.internal.customgbforces import GBSAGBnForce
+argMap = {'soluteDielectric':'soluteDielectric', 'solventDielectric':'solventDielectric', 'implicitSolventKappa':'kappa'}
+solventArgs = {'SA':'ACE'}
+for key in argMap:
+    if key in args:
+        solventArgs[argMap[key]] = args[key]
+if nonbondedMethod != app.NoCutoff:
+    solventArgs['cutoff'] = nonbondedCutoff.value_in_unit(unit.nanometers)
+force = GBSAGBnForce(**solventArgs)
+params = GBSAGBnForce.getStandardParameters(topology)
+for i, p in enumerate(params):
+    charge, sigma, epsilon = nonbonded.getParticleParameters(i)
+    force.addParticle([charge, p[0], p[1]])
+force.finalize()
+
+# Set the nonbonded method and cutoff distance.
+
+if nonbondedMethod == app.NoCutoff:
+    force.setNonbondedMethod(openmm.CustomNonbondedForce.NoCutoff)
+elif nonbondedMethod == app.CutoffNonPeriodic:
+    force.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffNonPeriodic)
+    force.setCutoffDistance(nonbondedCutoff)
+elif nonbondedMethod == app.CutoffPeriodic:
+    force.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffPeriodic)
+    force.setCutoffDistance(nonbondedCutoff)
+else:
+    raise ValueError("Illegal nonbonded method for use with implicit solvent")
+sys.addForce(force)
+ </Script>
+</ForceField>
\ No newline at end of file

wrappers/python/openmm/app/data/implicit/gbn2.xml

+<ForceField>
+ <Script>
+import openmm
+import openmm.app as app
+import openmm.unit as unit
+
+# Find the NonbondedForce.  We need it to look up charges, and also to set the reaction field dielectric to 1.
+
+nonbonded = [f for f in sys.getForces() if isinstance(f, openmm.NonbondedForce)]
+if len(nonbonded) != 1:
+    raise ValueError('Implicit solvent requires the System to contain a single NonbondedForce')
+nonbonded = nonbonded[0]
+nonbonded.setReactionFieldDielectric(1)
+
+# Construct the CustomGBForce.
+
+from openmm.app.internal.customgbforces import GBSAGBn2Force
+argMap = {'soluteDielectric':'soluteDielectric', 'solventDielectric':'solventDielectric', 'implicitSolventKappa':'kappa'}
+solventArgs = {'SA':'ACE'}
+for key in argMap:
+    if key in args:
+        solventArgs[argMap[key]] = args[key]
+if nonbondedMethod != app.NoCutoff:
+    solventArgs['cutoff'] = nonbondedCutoff.value_in_unit(unit.nanometers)
+force = GBSAGBn2Force(**solventArgs)
+params = GBSAGBn2Force.getStandardParameters(topology)
+for i, p in enumerate(params):
+    charge, sigma, epsilon = nonbonded.getParticleParameters(i)
+    force.addParticle([charge, p[0], p[1], p[2], p[3], p[4]])
+force.finalize()
+
+# Set the nonbonded method and cutoff distance.
+
+if nonbondedMethod == app.NoCutoff:
+    force.setNonbondedMethod(openmm.CustomNonbondedForce.NoCutoff)
+elif nonbondedMethod == app.CutoffNonPeriodic:
+    force.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffNonPeriodic)
+    force.setCutoffDistance(nonbondedCutoff)
+elif nonbondedMethod == app.CutoffPeriodic:
+    force.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffPeriodic)
+    force.setCutoffDistance(nonbondedCutoff)
+else:
+    raise ValueError("Illegal nonbonded method for use with implicit solvent")
+sys.addForce(force)
+ </Script>
+</ForceField>
\ No newline at end of file

wrappers/python/openmm/app/data/implicit/hct.xml

+<ForceField>
+ <Script>
+import openmm
+import openmm.app as app
+import openmm.unit as unit
+
+# Find the NonbondedForce.  We need it to look up charges, and also to set the reaction field dielectric to 1.
+
+nonbonded = [f for f in sys.getForces() if isinstance(f, openmm.NonbondedForce)]
+if len(nonbonded) != 1:
+    raise ValueError('Implicit solvent requires the System to contain a single NonbondedForce')
+nonbonded = nonbonded[0]
+nonbonded.setReactionFieldDielectric(1)
+
+# Construct the CustomGBForce.
+
+from openmm.app.internal.customgbforces import GBSAHCTForce
+argMap = {'soluteDielectric':'soluteDielectric', 'solventDielectric':'solventDielectric', 'implicitSolventKappa':'kappa'}
+solventArgs = {'SA':'ACE'}
+for key in argMap:
+    if key in args:
+        solventArgs[argMap[key]] = args[key]
+if nonbondedMethod != app.NoCutoff:
+    solventArgs['cutoff'] = nonbondedCutoff.value_in_unit(unit.nanometers)
+force = GBSAHCTForce(**solventArgs)
+params = GBSAHCTForce.getStandardParameters(topology)
+for i, p in enumerate(params):
+    charge, sigma, epsilon = nonbonded.getParticleParameters(i)
+    force.addParticle([charge, p[0], p[1]])
+force.finalize()
+
+# Set the nonbonded method and cutoff distance.
+
+if nonbondedMethod == app.NoCutoff:
+    force.setNonbondedMethod(openmm.CustomNonbondedForce.NoCutoff)
+elif nonbondedMethod == app.CutoffNonPeriodic:
+    force.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffNonPeriodic)
+    force.setCutoffDistance(nonbondedCutoff)
+elif nonbondedMethod == app.CutoffPeriodic:
+    force.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffPeriodic)
+    force.setCutoffDistance(nonbondedCutoff)
+else:
+    raise ValueError("Illegal nonbonded method for use with implicit solvent")
+sys.addForce(force)
+ </Script>
+</ForceField>
\ No newline at end of file

wrappers/python/openmm/app/data/implicit/obc1.xml

+<ForceField>
+ <Script>
+import openmm
+import openmm.app as app
+import openmm.unit as unit
+
+# Find the NonbondedForce.  We need it to look up charges, and also to set the reaction field dielectric to 1.
+
+nonbonded = [f for f in sys.getForces() if isinstance(f, openmm.NonbondedForce)]
+if len(nonbonded) != 1:
+    raise ValueError('Implicit solvent requires the System to contain a single NonbondedForce')
+nonbonded = nonbonded[0]
+nonbonded.setReactionFieldDielectric(1)
+
+# Construct the CustomGBForce.
+
+from openmm.app.internal.customgbforces import GBSAOBC1Force
+argMap = {'soluteDielectric':'soluteDielectric', 'solventDielectric':'solventDielectric', 'implicitSolventKappa':'kappa'}
+solventArgs = {'SA':'ACE'}
+for key in argMap:
+    if key in args:
+        solventArgs[argMap[key]] = args[key]
+if nonbondedMethod != app.NoCutoff:
+    solventArgs['cutoff'] = nonbondedCutoff.value_in_unit(unit.nanometers)
+force = GBSAOBC1Force(**solventArgs)
+params = GBSAOBC1Force.getStandardParameters(topology)
+for i, p in enumerate(params):
+    charge, sigma, epsilon = nonbonded.getParticleParameters(i)
+    force.addParticle([charge, p[0], p[1]])
+force.finalize()
+
+# Set the nonbonded method and cutoff distance.
+
+if nonbondedMethod == app.NoCutoff:
+    force.setNonbondedMethod(openmm.CustomNonbondedForce.NoCutoff)
+elif nonbondedMethod == app.CutoffNonPeriodic:
+    force.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffNonPeriodic)
+    force.setCutoffDistance(nonbondedCutoff)
+elif nonbondedMethod == app.CutoffPeriodic:
+    force.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffPeriodic)
+    force.setCutoffDistance(nonbondedCutoff)
+else:
+    raise ValueError("Illegal nonbonded method for use with implicit solvent")
+sys.addForce(force)
+ </Script>
+</ForceField>
\ No newline at end of file

wrappers/python/openmm/app/data/implicit/obc2.xml

+<ForceField>
+ <Script>
+import openmm
+import openmm.app as app
+import openmm.unit as unit
+
+# Find the NonbondedForce.  We need it to look up charges, and also to set the reaction field dielectric to 1.
+
+nonbonded = [f for f in sys.getForces() if isinstance(f, openmm.NonbondedForce)]
+if len(nonbonded) != 1:
+    raise ValueError('Implicit solvent requires the System to contain a single NonbondedForce')
+nonbonded = nonbonded[0]
+nonbonded.setReactionFieldDielectric(1)
+
+# Construct the CustomGBForce.
+
+from openmm.app.internal.customgbforces import GBSAOBC2Force
+argMap = {'soluteDielectric':'soluteDielectric', 'solventDielectric':'solventDielectric', 'implicitSolventKappa':'kappa'}
+solventArgs = {'SA':'ACE'}
+for key in argMap:
+    if key in args:
+        solventArgs[argMap[key]] = args[key]
+if nonbondedMethod != app.NoCutoff:
+    solventArgs['cutoff'] = nonbondedCutoff.value_in_unit(unit.nanometers)
+force = GBSAOBC2Force(**solventArgs)
+params = GBSAOBC2Force.getStandardParameters(topology)
+for i, p in enumerate(params):
+    charge, sigma, epsilon = nonbonded.getParticleParameters(i)
+    force.addParticle([charge, p[0], p[1]])
+force.finalize()
+
+# Set the nonbonded method and cutoff distance.
+
+if nonbondedMethod == app.NoCutoff:
+    force.setNonbondedMethod(openmm.CustomNonbondedForce.NoCutoff)
+elif nonbondedMethod == app.CutoffNonPeriodic:
+    force.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffNonPeriodic)
+    force.setCutoffDistance(nonbondedCutoff)
+elif nonbondedMethod == app.CutoffPeriodic:
+    force.setNonbondedMethod(openmm.CustomNonbondedForce.CutoffPeriodic)
+    force.setCutoffDistance(nonbondedCutoff)
+else:
+    raise ValueError("Illegal nonbonded method for use with implicit solvent")
+sys.addForce(force)
+ </Script>
+</ForceField>
\ No newline at end of file

wrappers/python/setup.py

@@ -152,7 +152,7 @@ def buildKeywordDictionary(major_version_num=MAJOR_VERSION_NUM,
                                           "openmm.app.internal.pdbx.writer"]
     setupKeywords["data_files"]        = []
     setupKeywords["package_data"]      = {"openmm" : [],
-                                          "openmm.app" : ['data/*.xml', 'data/*.pdb', 'data/amber14/*.xml', 'data/charmm36/*.xml'],
+                                          "openmm.app" : ['data/*.xml', 'data/*.pdb', 'data/amber14/*.xml', 'data/charmm36/*.xml', 'data/implicit/*.xml'],
                                           "openmm.app.internal" : []}
     setupKeywords["platforms"]         = ["Linux", "Mac OS X", "Windows"]
     setupKeywords["description"]       = \

wrappers/python/tests/TestForceField.py

@@ -300,6 +300,27 @@ class TestForceField(unittest.TestCase):
                 numDifferences += 1
         self.assertTrue(numDifferences < system.getNumParticles()/20) # Tolerate occasional differences from numerical error
 
+    def test_ImplicitSolventForces(self):
+        """Compute forces for different implicit solvent types, and compare them to ones generated with AmberPrmtopFile."""
+
+        solventType = ['hct', 'obc1', 'obc2', 'gbn', 'gbn2']
+        nonbondedMethod = [NoCutoff, CutoffNonPeriodic, CutoffNonPeriodic, NoCutoff, NoCutoff]
+        kappa = [0.0, 0.0, 1.698295227342757, 1.698295227342757, 0.0]
+        file = [None, 'OBC1_NonPeriodic', 'OBC2_NonPeriodic_Salt', None, 'GBn2_NoCutoff']
+        for i in range(len(file)):
+            forcefield = ForceField('amber96.xml', f'implicit/{solventType[i]}.xml')
+            system = forcefield.createSystem(self.pdb2.topology, nonbondedMethod=nonbondedMethod[i], implicitSolventKappa=kappa[i])
+            integrator = VerletIntegrator(0.001)
+            context = Context(system, integrator, Platform.getPlatformByName("Reference"))
+            context.setPositions(self.pdb2.positions)
+            state1 = context.getState(getForces=True)
+            if file[i] is not None:
+                with open('systems/alanine-dipeptide-implicit-forces/'+file[i]+'.xml') as infile:
+                    state2 = XmlSerializer.deserialize(infile.read())
+                for f1, f2, in zip(state1.getForces().value_in_unit(kilojoules_per_mole/nanometer), state2.getForces().value_in_unit(kilojoules_per_mole/nanometer)):
+                    diff = norm(f1-f2)
+                    self.assertTrue(diff < 0.1 or diff/norm(f1) < 1e-4)
+
     def test_ProgrammaticForceField(self):
         """Test building a ForceField programmatically."""