Implement GBn2 (equivalent to igb==8 in Amber) via CustomGBForce. [1] The energy

computed with OpenMM matches that computed with Amber for the new model as well
as the existing models do (see table below). GBn2 is a much better GB model than
either the OBC models or the original GB neck by all indications...

This commit also does a few other things:
   -  Made the necessary changes to allow users to select GBn2 as the
      implicitSolvent model in AmberPrmtop.createSystem()

   -  Since AmberTools 12, Amber topology files have been written out with a
      ATOMIC_NUMBER section. When available, use this data to initialize element
      information in the Topology object generated by AmberPrmtopFile (when not
      available, try to determine elements from the atom name).

   -  GB parameters for GBn (and the new GBn2) should be determined based on the
      type of element. While this was approximated using the atom name, the
      better approach (and the one used in Amber since the addition of
      ATOMIC_NUMBER to the topology file) is to set the screening parameters
      from the element information.

[1] H. Nguyen, et al., J. Chem. Theory Comput. 9 (4) pp. 2020 - 2034

I also checked all of the GB model energies computed via OpenMM with those
computed via Amber, keeping all settings the same (I disabled the surface area
component in OpenMM and set rgbmax=cut=1000 in Amber, used full double
precision/CUDA in OpenMM, and used no constraints).

(all energies are kcal/mol)

             |  OpenMM    |    Amber
--------------------------------------
HCT (igb=1)  | -2284.0383 | -2283.9012
OBC1 (igb=2) | -2334.0054 | -2333.8666
OBC2 (igb=5) | -2220.8282 | -2220.6934
GBn (igb=7)* | -2155.8072 | -2155.6449
GBn2 (igb=8) | -2228.1681 | -2227.9147  <-- new

*Energies do not match OpenMM 5.2 due to bug fix in commit f7b98d01

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,21 +102,30 @@ class AmberPrmtopFile(object):
             if atomName in atomReplacements:
                 atomName = atomReplacements[atomName]
 
-            # Try to guess the element.
-
-            upper = atomName.upper()
-            if upper.startswith('CL'):
-                element = elem.chlorine
-            elif upper.startswith('NA'):
-                element = elem.sodium
-            elif upper.startswith('MG'):
-                element = elem.magnesium
-            else:
+            # Get the element from the prmtop file if available
+            if prmtop.has_atomic_number:
                 try:
-                    element = elem.get_by_symbol(atomName[0])
+                    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'):
+                    element = elem.chlorine
+                elif upper.startswith('NA'):
+                    element = elem.sodium
+                elif upper.startswith('MG'):
+                    element = elem.magnesium
+                else:
+                    try:
+                        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
 
@@ -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