Merge pull request #249 from swails/master

Add a salt concentration term to the custom GB forces

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

@@ -141,8 +141,10 @@ class AmberPrmtopFile(object):
             top.setUnitCellDimensions(tuple(x.value_in_unit(unit.nanometer) for x in prmtop.getBoxBetaAndDimensions()[1:4])*unit.nanometer)
 
     def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
-                     constraints=None, rigidWater=True, implicitSolvent=None, soluteDielectric=1.0, solventDielectric=78.5, removeCMMotion=True,
-                     hydrogenMass=None, ewaldErrorTolerance=0.0005):
+                     constraints=None, rigidWater=True, implicitSolvent=None,
+                     implicitSolventKappa=0.0*(1/unit.nanometer),
+                     soluteDielectric=1.0, solventDielectric=78.5,
+                     removeCMMotion=True, hydrogenMass=None, ewaldErrorTolerance=0.0005):
         """Construct an OpenMM System representing the topology described by this prmtop file.
 
         Parameters:
@@ -153,6 +155,7 @@ class AmberPrmtopFile(object):
            Allowed values are None, HBonds, AllBonds, or HAngles.
          - rigidWater (boolean=True) If true, water molecules will be fully rigid regardless of the value passed for the constraints argument
          - implicitSolvent (object=None) If not None, the implicit solvent model to use.  Allowed values are HCT, OBC1, OBC2, GBn, or GBn2.
+         - implicitSolventKappa (float=0.0*1/unit.nanometer) The Debye-screening parameter corresponding to ionic strength used for implicit solvent
          - soluteDielectric (float=1.0) The solute dielectric constant to use in the implicit solvent model.
          - solventDielectric (float=78.5) The solvent dielectric constant to use in the implicit solvent model.
          - removeCMMotion (boolean=True) If true, a CMMotionRemover will be added to the System
@@ -196,8 +199,8 @@ class AmberPrmtopFile(object):
             raise ValueError('Illegal value for implicit solvent model')
         sys = amber_file_parser.readAmberSystem(prmtop_loader=self._prmtop, shake=constraintString, nonbondedCutoff=nonbondedCutoff,
                                                  nonbondedMethod=methodMap[nonbondedMethod], flexibleConstraints=False, gbmodel=implicitString,
-                                                 soluteDielectric=soluteDielectric, solventDielectric=solventDielectric, rigidWater=rigidWater,
-                                                 elements=self.elements)
+                                                 soluteDielectric=soluteDielectric, solventDielectric=solventDielectric, implicitSolventKappa=implicitSolventKappa,
+                                                 rigidWater=rigidWater, elements=self.elements)
         if hydrogenMass is not None:
             for atom1, atom2 in self.topology.bonds():
                 if atom1.element == elem.hydrogen:

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

@@ -565,8 +565,10 @@ class PrmtopLoader(object):
 # AMBER System builder (based on, but not identical to, systemManager from 'zander')
 #=============================================================================================
 
-def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmodel=None, soluteDielectric=1.0, solventDielectric=78.5,
-                    nonbondedCutoff=None, nonbondedMethod='NoCutoff', scee=None, scnb=None, mm=None, verbose=False,
+def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmodel=None,
+          soluteDielectric=1.0, solventDielectric=78.5,
+          implicitSolventKappa=0.0*(1/units.nanometer), nonbondedCutoff=None,
+          nonbondedMethod='NoCutoff', scee=None, scnb=None, mm=None, verbose=False,
           EwaldErrorTolerance=None, flexibleConstraints=True, rigidWater=True, elements=None):
     """
     Create an OpenMM System from an Amber prmtop file.
@@ -580,6 +582,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
       gbmodel (String) - if 'OBC', OBC GBSA will be used; if 'GBVI', GB/VI will be used (default: None)
       soluteDielectric (float) - The solute dielectric constant to use in the implicit solvent model (default: 1.0)
       solventDielectric (float) - The solvent dielectric constant to use in the implicit solvent model (default: 78.5)
+      implicitSolventKappa (float) - The Debye screening parameter corresponding to implicit solvent ionic strength
       nonbondedCutoff (float) - if specified, will set nonbondedCutoff (default: None)
       scnb (float) - 1-4 Lennard-Jones scaling factor (default: taken from prmtop or 1.2 if not present there)
       scee (float) - 1-4 electrostatics scaling factor (default: taken from prmtop or 2.0 if not present there)
@@ -858,6 +861,9 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
 
     # Add GBSA model.
     if gbmodel is not None:
+        # Convert implicitSolventKappa to nanometers if it is a unit.
+        if units.is_quantity(implicitSolventKappa):
+            implicitSolventKappa = implicitSolventKappa.value_in_unit(1/units.nanometers)
         if verbose: print "Adding GB parameters..."
         charges = prmtop.getCharges()
         cutoff = None
@@ -867,17 +873,20 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
                 cutoff = cutoff.value_in_unit(units.nanometers)
         gb_parms = prmtop.getGBParms(gbmodel, elements)
         if gbmodel == 'HCT':
-            gb = customgb.GBSAHCTForce(solventDielectric, soluteDielectric, 'ACE', cutoff)
+            gb = customgb.GBSAHCTForce(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
         elif gbmodel == 'OBC1':
-            gb = customgb.GBSAOBC1Force(solventDielectric, soluteDielectric, 'ACE', cutoff)
+            gb = customgb.GBSAOBC1Force(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
         elif gbmodel == 'OBC2':
+            if implicitSolventKappa > 0:
+                gb = customgb.GBSAOBC2Force(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
+            else:
                 gb = mm.GBSAOBCForce()
                 gb.setSoluteDielectric(soluteDielectric)
                 gb.setSolventDielectric(solventDielectric)
         elif gbmodel == 'GBn':
-            gb = customgb.GBSAGBnForce(solventDielectric, soluteDielectric, 'ACE', cutoff)
+            gb = customgb.GBSAGBnForce(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
         elif gbmodel == 'GBn2':
-            gb = customgb.GBSAGBn2Force(solventDielectric, soluteDielectric, 'ACE', cutoff)
+            gb = customgb.GBSAGBn2Force(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
         else:
             raise Exception("Illegal value specified for implicit solvent model")
         for iAtom in range(prmtop.getNumAtoms()):

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

@@ -189,16 +189,24 @@ for i in range (len(d0)):
     m0[i]=m0[i]*10
 
 
-def _createEnergyTerms(force, SA, cutoff):
+def _createEnergyTerms(force, SA, cutoff, kappa):
     # Add the energy terms to the CustomGBForce.  These are identical for all the GB models.
     
-    force.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce.SingleParticle)
+    if kappa > 0:
+        force.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-exp(-kappa*B)/solventDielectric)*q^2/B",
+                CustomGBForce.SingleParticle)
+    elif kappa < 0:
+        # Do kappa check here to avoid repeating code everywhere
+        raise ValueError('kappa/ionic strength must be >= 0')
+    else:
+        force.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B",
+                CustomGBForce.SingleParticle)
     if SA=='ACE':
         force.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6", CustomGBForce.SingleParticle)
     elif SA is not None:
         raise ValueError('Unknown surface area method: '+SA)
     if cutoff is None:
-        force.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
+        force.addEnergyTerm("-138.935485*(1/soluteDielectric-exp(-kappa*f)/solventDielectric)*q1*q2/f;"
                             "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
     else:
         force.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2*(1/f-"+str(1/cutoff)+");"
@@ -209,13 +217,15 @@ Amber Equivalent: igb = 1
 """
 
 
-def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=None):
+def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None,
+                 cutoff=None, kappa=0.0):
 
     custom = CustomGBForce()
 
     custom.addPerParticleParameter("q")
     custom.addPerParticleParameter("radius")
     custom.addPerParticleParameter("scale")
+    if kappa > 0: custom.addGlobalParameter('kappa', kappa)
     custom.addGlobalParameter("solventDielectric", solventDielectric)
     custom.addGlobalParameter("soluteDielectric", soluteDielectric)
     custom.addGlobalParameter("offset", 0.009)
@@ -228,19 +238,21 @@ def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=Non
 
     custom.addComputedValue("B", "1/(1/or-I);"
                                   "or=radius-offset", CustomGBForce.SingleParticle)
-    _createEnergyTerms(custom, SA, cutoff)
+    _createEnergyTerms(custom, SA, cutoff, kappa)
     return custom
 
 """
 Amber Equivalents: igb = 2
 """
-def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=None):
+def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None,
+                  cutoff=None, kappa=0.0):
 
     custom = CustomGBForce()
 
     custom.addPerParticleParameter("q")
     custom.addPerParticleParameter("radius")
     custom.addPerParticleParameter("scale")
+    if kappa > 0: custom.addGlobalParameter('kappa', kappa)
     custom.addGlobalParameter("solventDielectric", solventDielectric)
     custom.addGlobalParameter("soluteDielectric", soluteDielectric)
     custom.addGlobalParameter("offset", 0.009)
@@ -253,19 +265,21 @@ def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=No
 
     custom.addComputedValue("B", "1/(1/or-tanh(0.8*psi+2.909125*psi^3)/radius);"
                                   "psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
-    _createEnergyTerms(custom, SA, cutoff)
+    _createEnergyTerms(custom, SA, cutoff, kappa)
     return custom
 
 """
 Amber Equivalents: igb = 5
 """
-def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=None):
+def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None,
+                  cutoff=None, kappa=0.0):
 
     custom = CustomGBForce()
 
     custom.addPerParticleParameter("q")
     custom.addPerParticleParameter("radius")
     custom.addPerParticleParameter("scale")
+    if kappa > 0: custom.addGlobalParameter('kappa', kappa)
     custom.addGlobalParameter("solventDielectric", solventDielectric)
     custom.addGlobalParameter("soluteDielectric", soluteDielectric)
     custom.addGlobalParameter("offset", 0.009)
@@ -278,13 +292,14 @@ def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=No
 
     custom.addComputedValue("B", "1/(1/or-tanh(psi-0.8*psi^2+4.85*psi^3)/radius);"
                                   "psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
-    _createEnergyTerms(custom, SA, cutoff)
+    _createEnergyTerms(custom, SA, cutoff, kappa)
     return custom
 
 """
 Amber Equivalents: igb = 7
 """
-def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=None):
+def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None,
+                  cutoff=None, kappa=0.0):
 
 
     """
@@ -301,6 +316,7 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=Non
     custom.addPerParticleParameter("radius")
     custom.addPerParticleParameter("scale")
 
+    if kappa > 0: custom.addGlobalParameter('kappa', kappa)
     custom.addGlobalParameter("solventDielectric", solventDielectric)
     custom.addGlobalParameter("soluteDielectric", soluteDielectric)
     custom.addGlobalParameter("offset", 0.009)
@@ -322,13 +338,14 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=Non
 
     custom.addComputedValue("B", "1/(1/or-tanh(1.09511284*psi-1.907992938*psi^2+2.50798245*psi^3)/radius);"
                               "psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
-    _createEnergyTerms(custom, SA, cutoff)
+    _createEnergyTerms(custom, SA, cutoff, kappa)
     return custom
 
 """
 Amber Equivalents: igb = 8
 """
-def GBSAGBn2Force(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=None):
+def GBSAGBn2Force(solventDielectric=78.5, soluteDielectric=1, SA=None,
+                  cutoff=None, kappa=0.0):
 
 
     """
@@ -348,6 +365,7 @@ def GBSAGBn2Force(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=No
     custom.addPerParticleParameter("beta")
     custom.addPerParticleParameter("gamma")
 
+    if kappa > 0: custom.addGlobalParameter('kappa', kappa)
     custom.addGlobalParameter("solventDielectric", solventDielectric)
     custom.addGlobalParameter("soluteDielectric", soluteDielectric)
     custom.addGlobalParameter("offset", 0.0195141)
@@ -369,5 +387,5 @@ def GBSAGBn2Force(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=No
 
     custom.addComputedValue("B", "1/(1/or-tanh(alpha*psi-beta*psi^2+gamma*psi^3)/radius);"
                               "psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
-    _createEnergyTerms(custom, SA, cutoff)
+    _createEnergyTerms(custom, SA, cutoff, kappa)
     return custom