Precompute GBNeck terms

The previous version of this code would compute
the values of d0 and m0 for every iteration. This
was unnecessary as these values are fixed at the
start of the calculation by the radii of the particles
involved. This resulted in unnecessary computation
and memory access that dramatically slowed simulations
using the GBn or GBn2 solvation models.

The new version pre-computes the values based on the
radii that are present in the system. Linear interpolation
is used, which is consistent with Amber. The previous
version of OpenMM used cubic, which gives slightly
different results.

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

@@ -1217,6 +1217,7 @@ class CharmmPsfFile(object):
             else:
                 raise ValueError('Illegal nonbonded method for use with GBSA')
             gb.setForceGroup(self.GB_FORCE_GROUP)
+            gb.finalize()
             system.addForce(gb)
             force.setReactionFieldDielectric(1.0) # applies to NonbondedForce
 

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

@@ -1006,8 +1006,8 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
         for i, (r, s) in enumerate(zip(radii, screen)):
             if abs(r - gb_parms[i][0]) > 1e-4 or abs(s - gb_parms[i][1]) > 1e-4:
                 if not warned:
-                    warnings.warn('Non-optimal GB parameters detected for GB '
-                                  'model %s' % gbmodel)
+                    warnings.warn(
+                        'Non-optimal GB parameters detected for GB model %s' % gbmodel)
                     warned = True
             gb_parms[i][0], gb_parms[i][1] = r, s
 
@@ -1019,6 +1019,7 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
                                 gb_parm[2], gb_parm[3], gb_parm[4]])
             else:
                 gb.addParticle([charge, gb_parm[0], gb_parm[1]])
+        gb.finalize()
         system.addForce(gb)
         if nonbondedMethod == 'NoCutoff':
             gb.setNonbondedMethod(mm.NonbondedForce.NoCutoff)

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-2015 University of Virginia and the Authors.
 Authors: Christoph Klein, Michael R. Shirts
-Contributors: Jason M. Swails, Peter Eastman
+Contributors: Jason M. Swails, Peter Eastman, Justin L. MacCallum
 
 Permission is hereby granted, free of charge, to any person obtaining a
 copy of this software and associated documentation files (the "Software"),
@@ -29,19 +29,26 @@ OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
 USE OR OTHER DEALINGS IN THE SOFTWARE.
 """
 
+
 from __future__ import division, absolute_import
 
 from collections import defaultdict
 import copy
+import numpy as np
+from scipy import interpolate
 from simtk.openmm.app import element as E
-from simtk.openmm import CustomGBForce, Continuous2DFunction
+from simtk.openmm import CustomGBForce, Discrete2DFunction
 import simtk.unit as u
 
+
 def strip_unit(value, unit):
+    """Strip off any units and return value in unit"""
     if not u.is_quantity(value):
         return value
     return value.value_in_unit(unit)
 
+
+# GBn and GBn2 lookup tables
 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,
     3.18854,3.24498,3.30132,3.35752,3.4136,
@@ -106,8 +113,6 @@ d0=[2.26685,2.32548,2.38397,2.44235,2.50057,2.55867,2.61663,2.67444,
     3.6629,3.71694,3.77095,3.82489,3.8788,3.93265,3.98646,4.04022,
     4.09395,4.14762,4.20126,4.25485,4.3084]
 
-
-
 m0=[0.0381511,0.0338587,0.0301776,0.027003,0.0242506,0.0218529,
     0.0197547,0.0179109,0.0162844,0.0148442,0.0135647,0.0124243,
     0.0114047,0.0104906,0.00966876,0.008928,0.0082587,0.00765255,
@@ -192,10 +197,10 @@ m0=[0.0381511,0.0338587,0.0301776,0.027003,0.0242506,0.0218529,
     0.0332033,0.030322,0.0277596,0.0254732,0.0234266,0.0215892,
     0.0199351,0.018442,0.0170909,0.0158654,0.0147514,0.0137365,
     0.0128101,0.0119627,0.0111863]
-# Rescale to nm
-for i in range (len(d0)):
-    d0[i]=d0[i]/10
-    m0[i]=m0[i]*10
+# Rescale to nanometers
+d0 = [d / 10. for d in d0]
+m0 = [m * 10. for m in m0]
+
 
 def _get_bonded_atom_list(topology):
     """ Returns a list of atoms bonded to each other atom in a dict """
@@ -205,6 +210,7 @@ def _get_bonded_atom_list(topology):
         bondeds[a2].append(a1)
     return bondeds
 
+
 def _is_carboxylateO(atom, all_bonds):
     if atom is not E.oxygen: return False
     bondeds = all_bonds[atom]
@@ -224,6 +230,7 @@ def _is_carboxylateO(atom, all_bonds):
     # If we got here, must be a carboxylate
     return True
 
+
 def _bondi_radii(topology):
     """ Sets the bondi radii """
     radii = [0.0 for atom in topology.atoms()]
@@ -250,6 +257,7 @@ def _bondi_radii(topology):
             radii[i] = 1.5
     return radii  # converted to nanometers above
 
+
 def _mbondi_radii(topology):
     """ Sets the mbondi radii """
     radii = [0.0 for atom in topology.atoms()]
@@ -286,6 +294,7 @@ def _mbondi_radii(topology):
             radii[i] = 1.5
     return radii  # converted to nanometers above
 
+
 def _mbondi2_radii(topology):
     """ Sets the mbondi2 radii """
     radii = [0.0 for atom in topology.atoms()]
@@ -320,6 +329,7 @@ def _mbondi2_radii(topology):
             radii[i] = 1.5
     return radii  # Converted to nanometers above
 
+
 def _mbondi3_radii(topology):
     """ Sets the mbondi3 radii """
     radii = _mbondi2_radii(topology)
@@ -333,9 +343,13 @@ def _mbondi3_radii(topology):
                 radii[i] = 1.17
     return radii  # Converted to nanometers above
 
+
 def _createEnergyTerms(force, solventDielectric, soluteDielectric, SA, cutoff, kappa, offset):
-    # Add the energy terms to the CustomGBForce.  These are identical for all the GB models.
+    """Add the energy terms to the CustomGBForce.
+
+    These are identical for all the GB models.
 
+    """
     params = "; solventDielectric=%.16g; soluteDielectric=%.16g; kappa=%.16g; offset=%.16g" % (solventDielectric, soluteDielectric, kappa, offset)
     if cutoff is not None:
         params += "; cutoff=%.16g" % cutoff
@@ -367,6 +381,7 @@ def _createEnergyTerms(force, solventDielectric, soluteDielectric, SA, cutoff, k
             force.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2*(1/f-"+str(1/cutoff)+");"
                                 "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))"+params, CustomGBForce.ParticlePairNoExclusions)
 
+
 _SCREEN_PARAMETERS = { # normal, GBn, GBn2
         E.hydrogen : (0.85, 1.09085413633, 1.425952),
         E.carbon : (0.72, 0.48435382330, 1.058554),
@@ -377,33 +392,52 @@ _SCREEN_PARAMETERS = { # normal, GBn, GBn2
         E.sulfur : (0.96, 0.602256336067, -0.703469),
         None : (0.8, 0.5, 0.5) # default
 }
+
 _SCREEN_PARAMETERS[E.deuterium] = _SCREEN_PARAMETERS[E.hydrogen]
 def _screen_parameter(atom):
     if atom.element in _SCREEN_PARAMETERS:
         return _SCREEN_PARAMETERS[atom.element]
     return _SCREEN_PARAMETERS[None]
 
-class CustomAmberGBForce(CustomGBForce):
 
+class CustomAmberGBForceBase(CustomGBForce):
+    """Base class for all of the Amber custom GB forces.
+
+    Should not be instantiated directly, use one of its
+    derived classes instead.
+
+    """
     OFFSET = 0.009
     RADIUS_ARG_POSITION = 1
     SCREEN_POSITION = 2
 
-    def __init__(self, *args, **kwargs):
-        raise NotImplementedError('Cannot instantiate ABC')
+    def __init__(self):
+        CustomGBForce.__init__(self)
+        self.parameters = []
 
     def addParticle(self, params):
-        params = copy.deepcopy(params)
-        params[self.RADIUS_ARG_POSITION] = strip_unit(params[self.RADIUS_ARG_POSITION], u.nanometer) - self.OFFSET
-        params[self.SCREEN_POSITION] *= params[self.RADIUS_ARG_POSITION]
-        CustomGBForce.addParticle(self, params)
-        return params
+        """Add a particle to the force
+
+        Particles are added in order. The number of particles
+        added must match the number of particles in the system.
+
+        Parameters
+        ----------
+        params : list
+            A list of parameters to add to the force. The meaning
+            parameters depends on the model.
 
-    def setParticleParameters(self, idx, params):
+        Returns
+        -------
+        list
+            The list of parameters after stripping off units and
+            modifying SCREEN.
+
+        """
         params = copy.deepcopy(params)
         params[self.RADIUS_ARG_POSITION] = strip_unit(params[self.RADIUS_ARG_POSITION], u.nanometer) - self.OFFSET
         params[self.SCREEN_POSITION] *= params[self.RADIUS_ARG_POSITION]
-        CustomGBForce.addParticle(self, params)
+        self.parameters.append(params)
         return params
 
     @staticmethod
@@ -421,17 +455,50 @@ class CustomAmberGBForce(CustomGBForce):
             List of all parameters needed for this GB model. These can be passed
             to addParticle or setParticleParameters after the charge is inserted
             at the beginning of the list
+
         """
-        raise NotImplementedError('Must override getStandardParameters in subclasses')
+        raise NotImplementedError(
+            'getStandardParameters must be defined in derived classes')
 
-"""
-Amber Equivalent: igb = 1
-"""
-class GBSAHCTForce(CustomAmberGBForce):
+    def finalize(self):
+        """Finalize this force so it can be added to a system.
 
+        This method must be called before the force is added
+        to the system.
+
+        """
+        self._addParticles()
+
+    def _addParticles(self):
+        for params in self.parameters:
+            CustomGBForce.addParticle(self, params)
+
+
+class GBSAHCTForce(CustomAmberGBForceBase):
+    """This class is equivalent to Amber ``igb=1``
+
+    The list of parameters to ``addParticle`` is: ``[q, or, sr]``.
+
+    Parameters
+    ----------
+    solventDielectric: float
+        Dielectric constant for the solvent
+    soluteDielectric: float
+        Dielectric constant for the solute
+    SA: string or None
+        Surface area model to use
+    cutoff: float or Quantity or None
+        Cutoff distance to use. If float, value is in nm. If ``None``,
+        then no cutoffs are used.
+    kappa: float or Quantity
+        Debye kappa parameter related to modelling salt in GB. It has
+        units of 1 / length with 1 / nanometer assumed if a float
+        is given. A value of zero corresponds to zero salt concentration.
+
+    """
     def __init__(self, solventDielectric=78.5, soluteDielectric=1, SA=None,
                  cutoff=None, kappa=0.0):
-        CustomGBForce.__init__(self)
+        CustomAmberGBForceBase.__init__(self)
 
         self.addPerParticleParameter("q")
         self.addPerParticleParameter("or") # Offset radius
@@ -447,20 +514,53 @@ class GBSAHCTForce(CustomAmberGBForce):
 
     @staticmethod
     def getStandardParameters(topology):
+        """ Gets list of standard parameters for this GB model based on an input Topology
+
+        Parameters
+        ----------
+        topology : simtk.openmm.app.Topology
+            Topology of the system to get parameters for
+
+        Returns
+        -------
+        list of float
+            List of all parameters needed for this GB model. These can be passed
+            to addParticle or setParticleParameters after the charge is inserted
+            at the beginning of the list
+
+        """
         radii = [[x/10] for x in _mbondi_radii(topology)]
         for i, atom in enumerate(topology.atoms()):
             radii[i].append(_screen_parameter(atom)[0])
         return radii
 
-"""
-Amber Equivalents: igb = 2
-"""
-class GBSAOBC1Force(CustomAmberGBForce):
 
+class GBSAOBC1Force(CustomAmberGBForceBase):
+    """This class is equivalent to Amber ``igb=2``
+
+    The list of parameters to ``addParticle`` is: ``[q, or, sr]``.
+
+    Parameters
+    ----------
+    solventDielectric: float
+        Dielectric constant for the solvent
+    soluteDielectric: float
+        Dielectric constant for the solute
+    SA: string or None
+        Surface area model to use
+    cutoff: float or Quantity or None
+        Cutoff distance to use. If float, value is in nm. If ``None``,
+        then no cutoffs are used.
+    kappa: float or Quantity
+        Debye kappa parameter related to modelling salt in GB. It has
+        units of 1 / length with 1 / nanometer assumed if a float
+        is given. A value of zero corresponds to zero salt concentration.
+
+    """
     def __init__(self, solventDielectric=78.5, soluteDielectric=1, SA=None,
                  cutoff=None, kappa=0.0):
 
-        CustomGBForce.__init__(self)
+        CustomAmberGBForceBase.__init__(self)
 
         self.addPerParticleParameter("q")
         self.addPerParticleParameter("or") # Offset radius
@@ -476,20 +576,55 @@ class GBSAOBC1Force(CustomAmberGBForce):
 
     @staticmethod
     def getStandardParameters(topology):
+        """ Gets list of standard parameters for this GB model based on an input Topology
+
+        Parameters
+        ----------
+        topology : simtk.openmm.app.Topology
+            Topology of the system to get parameters for
+
+        Returns
+        -------
+        list of float
+            List of all parameters needed for this GB model. These can be passed
+            to addParticle or setParticleParameters after the charge is inserted
+            at the beginning of the list
+
+        """
         radii = [[x/10] for x in _mbondi2_radii(topology)]
         for i, atom in enumerate(topology.atoms()):
             radii[i].append(_screen_parameter(atom)[0])
         return radii
 
-"""
-Amber Equivalents: igb = 5
-"""
-class GBSAOBC2Force(GBSAOBC1Force):
 
+class GBSAOBC2Force(GBSAOBC1Force):
+    """This class is equivalent to Amber ``igb=5``
+
+    The list of parameters to ``addParticle`` is: ``[q, or, sr]``.
+
+    Parameters
+    ----------
+    solventDielectric: float
+        Dielectric constant for the solvent
+    soluteDielectric: float
+        Dielectric constant for the solute
+    SA: string or None
+        Surface area model to use
+    cutoff: float or Quantity or None
+        Cutoff distance to use. If float, value is in nm. If ``None``,
+        then no cutoffs are used.
+    kappa: float or Quantity
+        Debye kappa parameter related to modelling salt in GB. It has
+        units of 1 / length with 1 / nanometer assumed if a float
+        is given. A value of zero corresponds to zero salt concentration.
+
+    """
     def __init__(self, solventDielectric=78.5, soluteDielectric=1, SA=None,
                  cutoff=None, kappa=0.0):
 
-        CustomGBForce.__init__(self)
+        # Note that this is not GBSAOBC1Force, as our initialization
+        # is different. We inherit for getStandardParameters.
+        CustomAmberGBForceBase.__init__(self)
 
         self.addPerParticleParameter("q")
         self.addPerParticleParameter("or") # Offset radius
@@ -503,90 +638,168 @@ class GBSAOBC2Force(GBSAOBC1Force):
                                      "psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
         _createEnergyTerms(self, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
 
-"""
-Amber Equivalents: igb = 7
-"""
-class GBSAGBnForce(CustomAmberGBForce):
-
-    def __init__(self, solventDielectric=78.5, soluteDielectric=1, SA=None,
-                 cutoff=None, kappa=0.0):
 
-        CustomGBForce.__init__(self)
-
-        self.addPerParticleParameter("q")
-        self.addPerParticleParameter("or") # Offset radius
-        self.addPerParticleParameter("sr") # Scaled offset radius
-
-        self.addTabulatedFunction("getd0", Continuous2DFunction(21, 21, d0, 0.1, 0.2, 0.1, 0.2))
-        self.addTabulatedFunction("getm0", Continuous2DFunction(21, 21, m0, 0.1, 0.2, 0.1, 0.2))
+class GBSAGBnForce(CustomAmberGBForceBase):
+    """This class is equivalent to Amber ``igb=7``
+
+    The list of parameters to ``addParticle`` is: ``[q, or, sr]``.
+
+    Parameters
+    ----------
+    solventDielectric: float
+        Dielectric constant for the solvent
+    soluteDielectric: float
+        Dielectric constant for the solute
+    SA: string or None
+        Surface area model to use
+    cutoff: float or Quantity or None
+        Cutoff distance to use. If float, value is in nm. If ``None``,
+        then no cutoffs are used.
+    kappa: float or Quantity
+        Debye kappa parameter related to modelling salt in GB. It has
+        units of 1 / length with 1 / nanometer assumed if a float
+        is given. A value of zero corresponds to zero salt concentration.
+
+    """
+    OFFSET = 0.009
 
-        self.addComputedValue("I",  "Ivdw+neckScale*Ineck;"
-                                    "Ineck=step(radius1+radius2+neckCut-r)*getm0(radius1,radius2)/(1+100*(r-getd0(radius1,radius2))^2+0.3*1000000*(r-getd0(radius1,radius2))^6);"
-                                    "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);"
-                                    "radius1=or1+offset; radius2=or2+offset;"
-                                    "neckScale=0.361825; neckCut=0.68; offset=0.009", CustomGBForce.ParticlePairNoExclusions)
+    def __init__(self, solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=None, kappa=0.0):
 
-        self.addComputedValue("B", "1/(1/or-tanh(1.09511284*psi-1.907992938*psi^2+2.50798245*psi^3)/radius);"
-                                  "psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
-        _createEnergyTerms(self, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.009)
+        CustomAmberGBForceBase.__init__(self)
+        self.solventDielectric = solventDielectric
+        self.soluteDielectric = soluteDielectric
+        self.SA = SA
+        self.cutoff = cutoff
+        self.kappa = kappa
+        self._uniqueRadii = None
+        self._radiusToIndex = None
 
     def addParticle(self, parameters):
-        parameters = CustomAmberGBForce.addParticle(self, parameters)
+        parameters = CustomAmberGBForceBase.addParticle(self, parameters)
         if parameters[1] + self.OFFSET < 0.1 or parameters[1] + self.OFFSET > 0.2:
             raise ValueError('Radii must be between 1 and 2 Angstroms for neck lookup')
 
-    def setParticleParameters(self, idx, parameters):
-        parameters = CustomAmberGBForce.setParticleParameters(self, idx, parameters)
-        if parameters[1] < 0.1 or parameters[1] > 0.2:
-            raise ValueError('Radii must be between 1 and 2 Angstroms for neck lookup')
+    def finalize(self):
+        self._findUniqueRadii()
+        self._createRadiusToIndexMap()
+        self._addEnergyTerms()
+        self._addParticles()
 
     @staticmethod
     def getStandardParameters(topology):
+        """ Gets list of standard parameters for this GB model based on an input Topology
+
+        Parameters
+        ----------
+        topology : simtk.openmm.app.Topology
+            Topology of the system to get parameters for
+
+        Returns
+        -------
+        list of float
+            List of all parameters needed for this GB model. These can be passed
+            to addParticle or setParticleParameters after the charge is inserted
+            at the beginning of the list
+
+        """
         radii = [[x/10] for x in _bondi_radii(topology)]
         for i, atom in enumerate(topology.atoms()):
             radii[i].append(_screen_parameter(atom)[1])
         return radii
 
-"""
-Amber Equivalents: igb = 8
-"""
-class GBSAGBn2Force(GBSAGBnForce):
-
-    OFFSET = 0.0195141
-
-    def __init__(self, solventDielectric=78.5, soluteDielectric=1, SA=None,
-                 cutoff=None, kappa=0.0):
-
-        CustomGBForce.__init__(self)
-
+    def _findUniqueRadii(self):
+        radii = [p[self.RADIUS_ARG_POSITION] for p in self.parameters]
+        self._uniqueRadii = list(sorted(set(radii)))
+
+    def _createRadiusToIndexMap(self):
+        self._radiusToIndex = {r: i for i, r in enumerate(self._uniqueRadii)}
+
+    def _createUniqueTable(self, fullTable):
+        fullTable = np.array(fullTable).reshape((21, 21))
+        xs = np.linspace(0.1, 0.2, 21)
+        interp = interpolate.interp2d(xs, xs, fullTable.T, kind='linear')
+        table = []
+        for r1 in self._uniqueRadii:
+            for r2 in self._uniqueRadii:
+                table.append(interp(r1 + self.OFFSET, r2 + self.OFFSET)[0])
+        return table
+
+    def _addParticles(self):
+        for p in self.parameters:
+            radIndex = self._radiusToIndex[p[self.RADIUS_ARG_POSITION]]
+            CustomGBForce.addParticle(self, p + [radIndex])
+
+    def _addEnergyTerms(self):
         self.addPerParticleParameter("q")
         self.addPerParticleParameter("or") # Offset radius
         self.addPerParticleParameter("sr") # Scaled offset radius
-        self.addPerParticleParameter("alpha")
-        self.addPerParticleParameter("beta")
-        self.addPerParticleParameter("gamma")
+        self.addPerParticleParameter("radindex")
+
+        n = len(self._uniqueRadii)
+        m0Table = self._createUniqueTable(m0)
+        d0Table = self._createUniqueTable(d0)
+        self.addTabulatedFunction("getd0", Discrete2DFunction(n, n, d0Table))
+        self.addTabulatedFunction("getm0", Discrete2DFunction(n, n, m0Table))
+
+        self.addComputedValue("I", "Ivdw+neckScale*Ineck;"
+                                   "Ineck=step(radius1+radius2+neckCut-r)*getm0(radindex1,radindex2)/(1+100*(r-getd0(radindex1,radindex2))^2+"
+                                   "0.3*1000000*(r-getd0(radindex1,radindex2))^6);"
+                                   "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);"
+                                   "radius1=or1+offset; radius2=or2+offset;"
+                                   "neckScale=0.361825; neckCut=0.68; offset=0.009", CustomGBForce.ParticlePairNoExclusions)
 
-        self.addTabulatedFunction("getd0", Continuous2DFunction(21, 21, d0, 0.1, 0.2, 0.1, 0.2))
-        self.addTabulatedFunction("getm0", Continuous2DFunction(21, 21, m0, 0.1, 0.2, 0.1, 0.2))
+        self.addComputedValue("B", "1/(1/or-tanh(1.09511284*psi-1.907992938*psi^2+2.50798245*psi^3)/radius);"
+                                  "psi=I*or; radius=or+offset; offset=0.009", CustomGBForce.SingleParticle)
+        _createEnergyTerms(self, self.solventDielectric, self.soluteDielectric, self.SA, self.cutoff, self.kappa, self.OFFSET)
 
-        self.addComputedValue("I",  "Ivdw+neckScale*Ineck;"
-                                    "Ineck=step(radius1+radius2+neckCut-r)*getm0(radius1,radius2)/(1+100*(r-getd0(radius1,radius2))^2+0.3*1000000*(r-getd0(radius1,radius2))^6);"
-                                    "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);"
-                                    "radius1=or1+offset; radius2=or2+offset;"
-                                    "neckScale=0.826836; neckCut=0.68; offset=0.0195141", CustomGBForce.ParticlePairNoExclusions)
 
-        self.addComputedValue("B", "1/(1/or-tanh(alpha*psi-beta*psi^2+gamma*psi^3)/radius);"
-                                     "psi=I*or; radius=or+offset; offset=0.0195141", CustomGBForce.SingleParticle)
-        _createEnergyTerms(self, solventDielectric, soluteDielectric, SA, cutoff, kappa, 0.0195141)
+class GBSAGBn2Force(GBSAGBnForce):
+    """This class is equivalent to Amber ``igb=8``
+
+    The list of parameters to ``addParticle`` is: ``[q, or, sr, alpha, beta, gamma]``.
+
+    Parameters
+    ----------
+    solventDielectric: float
+        Dielectric constant for the solvent
+    soluteDielectric: float
+        Dielectric constant for the solute
+    SA: string or None
+        Surface area model to use
+    cutoff: float or Quantity or None
+        Cutoff distance to use. If float, value is in nm. If ``None``,
+        then no cutoffs are used.
+    kappa: float or Quantity
+        Debye kappa parameter related to modelling salt in GB. It has
+        units of 1 / length with 1 / nanometer assumed if a float
+        is given. A value of zero corresponds to zero salt concentration.
+
+    """
+    OFFSET = 0.0195141
+
+    def __init__(self, solventDielectric=78.5, soluteDielectric=1, SA=None, cutoff=None, kappa=0.0):
+        GBSAGBnForce.__init__(self, solventDielectric, soluteDielectric, SA, cutoff, kappa)
 
     @staticmethod
     def getStandardParameters(topology):
+        """ Gets list of standard parameters for this GB model based on an input Topology
+
+        Parameters
+        ----------
+        topology : simtk.openmm.app.Topology
+            Topology of the system to get parameters for
+
+        Returns
+        -------
+        list of float
+            List of all parameters needed for this GB model. These can be passed
+            to addParticle or setParticleParameters after the charge is inserted
+            at the beginning of the list
+
+        """
         radii = [[x/10] for x in _mbondi3_radii(topology)]
         for i, atom in enumerate(topology.atoms()):
             radii[i].append(_screen_parameter(atom)[2])
@@ -603,3 +816,32 @@ class GBSAGBn2Force(GBSAGBnForce):
             else:
                 radii[i].extend([0.8, 4.85, 0.5])
         return radii
+
+    def _addEnergyTerms(self):
+        self.addPerParticleParameter("q")
+        self.addPerParticleParameter("or") # Offset radius
+        self.addPerParticleParameter("sr") # Scaled offset radius
+        self.addPerParticleParameter("alpha")
+        self.addPerParticleParameter("beta")
+        self.addPerParticleParameter("gamma")
+        self.addPerParticleParameter("radindex")
+
+        n = len(self._uniqueRadii)
+        m0Table = self._createUniqueTable(m0)
+        d0Table = self._createUniqueTable(d0)
+        self.addTabulatedFunction("getd0", Discrete2DFunction(n, n, d0Table))
+        self.addTabulatedFunction("getm0", Discrete2DFunction(n, n, m0Table))
+
+        self.addComputedValue("I", "Ivdw+neckScale*Ineck;"
+                                   "Ineck=step(radius1+radius2+neckCut-r)*getm0(radindex1,radindex2)/(1+100*(r-getd0(radindex1,radindex2))^2+"
+                                   "0.3*1000000*(r-getd0(radindex1,radindex2))^6);"
+                                   "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);"
+                                   "radius1=or1+offset; radius2=or2+offset;"
+                                   "neckScale=0.826836; neckCut=0.68; offset=0.0195141", CustomGBForce.ParticlePairNoExclusions)
+
+        self.addComputedValue("B", "1/(1/or-tanh(alpha*psi-beta*psi^2+gamma*psi^3)/radius);"
+                                     "psi=I*or; radius=or+offset; offset=0.0195141", CustomGBForce.SingleParticle)
+        _createEnergyTerms(self, self.solventDielectric, self.soluteDielectric, self.SA, self.cutoff, self.kappa, self.OFFSET)