Merge pull request #179 from swails/master

Implement GBn2 (same as igb==8 in Amber), a significantly improved GB model

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

@@ -61,6 +61,11 @@ class GBn(object):
         return 'GBn'
 GBn = GBn()
 
+class GBn2(object):
+    def __repr__(self):
+        return 'GBn2'
+GBn2 = GBn2()
+
 class AmberPrmtopFile(object):
     """AmberPrmtopFile parses an AMBER prmtop file and constructs a Topology and (optionally) an OpenMM System from it."""
 
@@ -69,6 +74,7 @@ class AmberPrmtopFile(object):
         top = Topology()
         ## The Topology read from the prmtop file
         self.topology = top
+        self.elements = []
 
         # Load the prmtop file
 
@@ -96,7 +102,14 @@ class AmberPrmtopFile(object):
             if atomName in atomReplacements:
                 atomName = atomReplacements[atomName]
 
-            # Try to guess the element.
+            # Get the element from the prmtop file if available
+            if prmtop.has_atomic_number:
+                try:
+                    element = elem.Element.getByAtomicNumber(int(prmtop._raw_data['ATOMIC_NUMBER'][index]))
+                except KeyError:
+                    element = None
+            else:
+                # Try to guess the element from the atom name.
 
                 upper = atomName.upper()
                 if upper.startswith('CL'):
@@ -110,7 +123,9 @@ class AmberPrmtopFile(object):
                         element = elem.get_by_symbol(atomName[0])
                     except KeyError:
                         element = None
+
             top.addAtom(atomName, element, r)
+            self.elements.append(element)
 
         # Add bonds to the topology
 
@@ -137,7 +152,7 @@ class AmberPrmtopFile(object):
          - constraints (object=None) Specifies which bonds angles should be implemented with constraints.
            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, or GBn.
+         - implicitSolvent (object=None) If not None, the implicit solvent model to use.  Allowed values are HCT, OBC1, OBC2, GBn, or GBn2.
          - 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
@@ -167,19 +182,22 @@ class AmberPrmtopFile(object):
             raise ValueError('Illegal value for constraints')
         if implicitSolvent is None:
             implicitString = None
-        elif implicitSolvent == HCT:
+        elif implicitSolvent is HCT:
             implicitString = 'HCT'
-        elif implicitSolvent == OBC1:
+        elif implicitSolvent is OBC1:
             implicitString = 'OBC1'
-        elif implicitSolvent == OBC2:
+        elif implicitSolvent is OBC2:
             implicitString = 'OBC2'
-        elif implicitSolvent == GBn:
+        elif implicitSolvent is GBn:
             implicitString = 'GBn'
+        elif implicitSolvent is GBn2:
+            implicitString = 'GBn2'
         else:
             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)
+                                                 soluteDielectric=soluteDielectric, solventDielectric=solventDielectric, 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

@@ -52,6 +52,7 @@ except:
 
 import simtk.unit as units
 import simtk.openmm
+from simtk.openmm.app import element as elem
 import customgbforces as customgb
 
 #=============================================================================================
@@ -472,34 +473,74 @@ class PrmtopLoader(object):
             self._excludedAtoms.append(atomList)
         return self._excludedAtoms
 
-    def getGBParms(self, symbls=None):
+    def getGBParms(self, gbmodel, elements):
         """Return list giving GB params, Radius and screening factor"""
-        try:
-            return self._gb_List
-        except AttributeError:
-            pass
-        self._gb_List=[]
-        radii=self._raw_data["RADII"]
-        screen=self._raw_data["SCREEN"]
+        gb_List=[]
+        radii=[float(r) for r in self._raw_data["RADII"]]
+        screen=[float(s) for s in self._raw_data["SCREEN"]]
+        if gbmodel == 'GBn2':
+            alpha = [0 for i in self._raw_data['RADII']]
+            beta = [0 for i in self._raw_data['RADII']]
+            gamma = [0 for i in self._raw_data['RADII']]
         # Update screening parameters for GBn if specified
-        if symbls:
-            for (i, symbl) in enumerate(symbls):
-                if symbl[0] in ('c', 'C'):
+        if gbmodel == 'GBn':
+            if elements is None:
+                raise Exception('GBn model requires element information')
+            for i, element in enumerate(elements):
+                if element is elem.carbon:
                     screen[i] = 0.48435382330
-                elif symbl[0] in ('h', 'H'):
+                elif element is elem.hydrogen:
                     screen[i] = 1.09085413633
-                elif symbl[0] in ('n', 'N'):
+                elif element is elem.nitrogen:
                     screen[i] = 0.700147318409
-                elif symbl[0] in ('o', 'O'):
+                elif element is elem.oxygen:
                     screen[i] = 1.06557401132
-                elif symbl[0] in ('s', 'S'):
+                elif element is elem.sulfur:
                     screen[i] = 0.602256336067
                 else:
                     screen[i] = 0.5
+        if gbmodel == 'GBn2':
+            if elements is None:
+                raise Exception('GBn2 model requires element information')
+            for i, element in enumerate(elements):
+                if element is elem.carbon:
+                    screen[i] = 1.058554
+                    alpha[i] = 0.733756
+                    beta[i] = 0.506378
+                    gamma[i] = 0.205844
+                elif element is elem.hydrogen:
+                    screen[i] = 1.425952
+                    alpha[i] = 0.788440
+                    beta[i] = 0.798699
+                    gamma[i] = 0.437334
+                elif element is elem.nitrogen:
+                    screen[i] = 0.733599
+                    alpha[i] = 0.503364
+                    beta[i] = 0.316828
+                    gamma[i] = 0.192915
+                elif element is elem.oxygen:
+                    screen[i] = 1.061039
+                    alpha[i] = 0.867814
+                    beta[i] = 0.876635
+                    gamma[i] = 0.387882
+                elif element is elem.sulfur:
+                    screen[i] = -0.703469
+                    alpha[i] = 0.867814
+                    beta[i] = 0.876635
+                    gamma[i] = 0.387882
+                else: # not optimized
+                    screen[i] = 0.5
+                    alpha[i] = 1.0
+                    beta[i] = 0.8
+                    gamma[i] = 4.85
         lengthConversionFactor = units.angstrom.conversion_factor_to(units.nanometer)
-        for iAtom in range(len(radii)):
-            self._gb_List.append((float(radii[iAtom])*lengthConversionFactor, float(screen[iAtom])))
-        return self._gb_List
+        if gbmodel == 'GBn2':
+            for rad, scr, alp, bet, gam in zip(radii, screen, alpha, beta, gamma):
+                gb_List.append((rad*lengthConversionFactor, scr, alp, bet, gam))
+        else:
+            for rad, scr in zip(radii, screen):
+                gb_List.append((rad*lengthConversionFactor, scr))
+        return gb_List
 
     def getBoxBetaAndDimensions(self):
         """Return periodic boundary box beta angle and dimensions"""
@@ -516,11 +557,17 @@ class PrmtopLoader(object):
     def has_scee_scnb(self):
         return ("SCEE_SCALE_FACTOR" in self._raw_data and "SCNB_SCALE_FACTOR" in self._raw_data)
 
+    @property
+    def has_atomic_number(self):
+        return 'ATOMIC_NUMBER' in self._raw_data
+
 #=============================================================================================
 # 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, EwaldErrorTolerance=None, flexibleConstraints=True, rigidWater=True):
+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,
+                    EwaldErrorTolerance=None, flexibleConstraints=True, rigidWater=True, elements=None):
     """
     Create an OpenMM System from an Amber prmtop file.
 
@@ -813,15 +860,12 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
     if gbmodel is not None:
         if verbose: print "Adding GB parameters..."
         charges = prmtop.getCharges()
-        symbls = None
         cutoff = None
         if nonbondedMethod != 'NoCutoff':
             cutoff = nonbondedCutoff
             if units.is_quantity(cutoff):
                 cutoff = cutoff.value_in_unit(units.nanometers)
-        if gbmodel == 'GBn':
-            symbls = prmtop.getAtomTypes()
-        gb_parms = prmtop.getGBParms(symbls)
+        gb_parms = prmtop.getGBParms(gbmodel, elements)
         if gbmodel == 'HCT':
             gb = customgb.GBSAHCTForce(solventDielectric, soluteDielectric, 'ACE', cutoff)
         elif gbmodel == 'OBC1':
@@ -832,11 +876,16 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
             gb.setSolventDielectric(solventDielectric)
         elif gbmodel == 'GBn':
             gb = customgb.GBSAGBnForce(solventDielectric, soluteDielectric, 'ACE', cutoff)
+        elif gbmodel == 'GBn2':
+            gb = customgb.GBSAGBn2Force(solventDielectric, soluteDielectric, 'ACE', cutoff)
         else:
             raise Exception("Illegal value specified for implicit solvent model")
         for iAtom in range(prmtop.getNumAtoms()):
             if gbmodel == 'OBC2':
                 gb.addParticle(charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1])
+            elif gbmodel == 'GBn2':
+                gb.addParticle([charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1],
+                                gb_parms[iAtom][2], gb_parms[iAtom][3], gb_parms[iAtom][4]])
             else:
                 gb.addParticle([charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1]])
         system.addForce(gb)

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

@@ -8,7 +8,7 @@ Medical Research, grant U54 GM072970. See https://simtk.org.
 
 Portions copyright (c) 2012 University of Virginia and the Authors.
 Authors: Christoph Klein, Michael R. Shirts
-Contributors:
+Contributors: Jason M. Swails
 
 Permission is hereby granted, free of charge, to any person obtaining a
 copy of this software and associated documentation files (the "Software"),
@@ -29,8 +29,9 @@ OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
 USE OR OTHER DEALINGS IN THE SOFTWARE.
 """
 
+from __future__ import division
+
 from simtk.openmm import CustomGBForce
-import sys, pdb, pickle
 
 d0=[2.26685,2.32548,2.38397,2.44235,2.50057,2.55867,2.61663,2.67444,
     2.73212,2.78965,2.84705,2.9043,2.96141,3.0184,3.07524,3.13196,
@@ -212,11 +213,11 @@ def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=Non
 
     custom = CustomGBForce()
 
-    custom.addPerParticleParameter("q");
-    custom.addPerParticleParameter("radius");
-    custom.addPerParticleParameter("scale");
-    custom.addGlobalParameter("solventDielectric", solventDielectric);
-    custom.addGlobalParameter("soluteDielectric", soluteDielectric);
+    custom.addPerParticleParameter("q")
+    custom.addPerParticleParameter("radius")
+    custom.addPerParticleParameter("scale")
+    custom.addGlobalParameter("solventDielectric", solventDielectric)
+    custom.addGlobalParameter("soluteDielectric", soluteDielectric)
     custom.addGlobalParameter("offset", 0.009)
     custom.addComputedValue("I", "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
                                   "U=r+sr2;"
@@ -237,11 +238,11 @@ def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=No
 
     custom = CustomGBForce()
 
-    custom.addPerParticleParameter("q");
-    custom.addPerParticleParameter("radius");
-    custom.addPerParticleParameter("scale");
-    custom.addGlobalParameter("solventDielectric", solventDielectric);
-    custom.addGlobalParameter("soluteDielectric", soluteDielectric);
+    custom.addPerParticleParameter("q")
+    custom.addPerParticleParameter("radius")
+    custom.addPerParticleParameter("scale")
+    custom.addGlobalParameter("solventDielectric", solventDielectric)
+    custom.addGlobalParameter("soluteDielectric", soluteDielectric)
     custom.addGlobalParameter("offset", 0.009)
     custom.addComputedValue("I",  "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
                                   "U=r+sr2;"
@@ -262,11 +263,11 @@ def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=No
 
     custom = CustomGBForce()
 
-    custom.addPerParticleParameter("q");
-    custom.addPerParticleParameter("radius");
-    custom.addPerParticleParameter("scale");
-    custom.addGlobalParameter("solventDielectric", solventDielectric);
-    custom.addGlobalParameter("soluteDielectric", soluteDielectric);
+    custom.addPerParticleParameter("q")
+    custom.addPerParticleParameter("radius")
+    custom.addPerParticleParameter("scale")
+    custom.addGlobalParameter("solventDielectric", solventDielectric)
+    custom.addGlobalParameter("soluteDielectric", soluteDielectric)
     custom.addGlobalParameter("offset", 0.009)
     custom.addComputedValue("I",  "step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
                                   "U=r+sr2;"
@@ -296,12 +297,12 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=Non
 
     custom = CustomGBForce()
 
-    custom.addPerParticleParameter("q");
-    custom.addPerParticleParameter("radius");
-    custom.addPerParticleParameter("scale");
+    custom.addPerParticleParameter("q")
+    custom.addPerParticleParameter("radius")
+    custom.addPerParticleParameter("scale")
 
-    custom.addGlobalParameter("solventDielectric", solventDielectric);
-    custom.addGlobalParameter("soluteDielectric", soluteDielectric);
+    custom.addGlobalParameter("solventDielectric", solventDielectric)
+    custom.addGlobalParameter("soluteDielectric", soluteDielectric)
     custom.addGlobalParameter("offset", 0.009)
     custom.addGlobalParameter("neckScale", 0.361825)
     custom.addGlobalParameter("neckCut", 0.68)
@@ -323,3 +324,50 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=Non
                               "psi=I*or; or=radius-offset", CustomGBForce.SingleParticle)
     _createEnergyTerms(custom, SA, cutoff)
     return custom
+
+"""
+Amber Equivalents: igb = 8
+"""
+def GBSAGBn2Force(solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=None):
+
+
+    """
+    Indexing for tables:
+        input: radius1, radius2
+        index = (radius2*200-20)*21 + (radius1*200-20)
+        output: index of desired value in row-by-row, 1D version of Tables 3 & 4
+    """
+
+
+    custom = CustomGBForce()
+
+    custom.addPerParticleParameter("q")
+    custom.addPerParticleParameter("radius")
+    custom.addPerParticleParameter("scale")
+    custom.addPerParticleParameter("alpha")
+    custom.addPerParticleParameter("beta")
+    custom.addPerParticleParameter("gamma")
+
+    custom.addGlobalParameter("solventDielectric", solventDielectric)
+    custom.addGlobalParameter("soluteDielectric", soluteDielectric)
+    custom.addGlobalParameter("offset", 0.0195141)
+    custom.addGlobalParameter("neckScale", 0.826836)
+    custom.addGlobalParameter("neckCut", 0.68)
+
+    custom.addFunction("getd0", d0, 0, 440)
+    custom.addFunction("getm0", m0, 0, 440)
+
+    custom.addComputedValue("I",  "Ivdw+neckScale*Ineck;"
+                                  "Ineck=step(radius1+radius2+neckCut-r)*getm0(index)/(1+100*(r-getd0(index))^2+0.3*1000000*(r-getd0(index))^6);"
+                                  "index = (radius2*200-20)*21 + (radius1*200-20);"
+                                  "Ivdw=step(r+sr2-or1)*0.5*(1/L-1/U+0.25*(r-sr2^2/r)*(1/(U^2)-1/(L^2))+0.5*log(L/U)/r);"
+                                  "U=r+sr2;"
+                                  "L=max(or1, D);"
+                                  "D=abs(r-sr2);"
+                                  "sr2 = scale2*or2;"
+                                  "or1 = radius1-offset; or2 = radius2-offset", CustomGBForce.ParticlePairNoExclusions)
+
+    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)
+    return custom