Merge pull request #1805 from swails/feature/gbsa-parameter

Make the nonpolar part of GBSA available via an argument to createSystem for both Amber and CHARMM files.

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

@@ -161,7 +161,7 @@ class AmberPrmtopFile(object):
                      implicitSolventKappa=None, temperature=298.15*u.kelvin,
                      soluteDielectric=1.0, solventDielectric=78.5,
                      removeCMMotion=True, hydrogenMass=None, ewaldErrorTolerance=0.0005,
-                     switchDistance=0.0*u.nanometer):
+                     switchDistance=0.0*u.nanometer, gbsaModel='ACE'):
         """Construct an OpenMM System representing the topology described by this
         prmtop file.
 
@@ -208,6 +208,11 @@ class AmberPrmtopFile(object):
             turned on for Lennard-Jones interactions. If the switchDistance is 0
             or evaluates to boolean False, no switching function will be used.
             Values greater than nonbondedCutoff or less than 0 raise ValueError
+        gbsaModel : str='ACE'
+            The SA model used to model the nonpolar solvation component of GB
+            implicit solvent models. If GB is active, this must be 'ACE' or None
+            (the latter indicates no SA model will be used). Other values will
+            result in a ValueError
 
         Returns
         -------
@@ -273,7 +278,7 @@ class AmberPrmtopFile(object):
                         nonbondedCutoff=nonbondedCutoff, nonbondedMethod=methodMap[nonbondedMethod],
                         flexibleConstraints=False, gbmodel=implicitString, soluteDielectric=soluteDielectric,
                         solventDielectric=solventDielectric, implicitSolventKappa=implicitSolventKappa,
-                        rigidWater=rigidWater, elements=self.elements)
+                        rigidWater=rigidWater, elements=self.elements, gbsaModel=gbsaModel)
 
         if hydrogenMass is not None:
             for atom1, atom2 in self.topology.bonds():

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

@@ -678,7 +678,8 @@ class CharmmPsfFile(object):
                      hydrogenMass=None,
                      ewaldErrorTolerance=0.0005,
                      flexibleConstraints=True,
-                     verbose=False):
+                     verbose=False,
+                     gbsaModel=None):
         """Construct an OpenMM System representing the topology described by the
         prmtop file. You MUST have loaded a parameter set into this PSF before
         calling createSystem. If not, AttributeError will be raised. ValueError
@@ -733,10 +734,16 @@ class CharmmPsfFile(object):
             Are our constraints flexible or not?
         verbose : bool=False
             Optionally prints out a running progress report
+        gbsaModel : str=None
+            Can be ACE (to use the ACE solvation model) or None. Other values
+            raise a ValueError
         """
         # Load the parameter set
         self.loadParameters(params.condense())
         hasbox = self.topology.getUnitCellDimensions() is not None
+        # Check GB input parameters
+        if implicitSolvent is not None and gbsaModel not in ('ACE', None):
+            raise ValueError('gbsaModel must be ACE or None')
         # Set the cutoff distance in nanometers
         cutoff = None
         if nonbondedMethod is not ff.NoCutoff:
@@ -1191,19 +1198,19 @@ class CharmmPsfFile(object):
                 implicitSolventKappa = implicitSolventKappa.value_in_unit(
                                             (1.0/u.nanometer).unit)
             if implicitSolvent is HCT:
-                gb = GBSAHCTForce(solventDielectric, soluteDielectric, None,
+                gb = GBSAHCTForce(solventDielectric, soluteDielectric, gbsaModel,
                                   cutoff, kappa=implicitSolventKappa)
             elif implicitSolvent is OBC1:
-                gb = GBSAOBC1Force(solventDielectric, soluteDielectric, None,
+                gb = GBSAOBC1Force(solventDielectric, soluteDielectric, gbsaModel,
                                    cutoff, kappa=implicitSolventKappa)
             elif implicitSolvent is OBC2:
-                gb = GBSAOBC2Force(solventDielectric, soluteDielectric, None,
+                gb = GBSAOBC2Force(solventDielectric, soluteDielectric, gbsaModel,
                                    cutoff, kappa=implicitSolventKappa)
             elif implicitSolvent is GBn:
-                gb = GBSAGBnForce(solventDielectric, soluteDielectric, None,
+                gb = GBSAGBnForce(solventDielectric, soluteDielectric, gbsaModel,
                                   cutoff, kappa=implicitSolventKappa)
             elif implicitSolvent is GBn2:
-                gb = GBSAGBn2Force(solventDielectric, soluteDielectric, None,
+                gb = GBSAGBn2Force(solventDielectric, soluteDielectric, gbsaModel,
                                    cutoff, kappa=implicitSolventKappa)
             gb_parms = gb.getStandardParameters(self.topology)
             for atom, gb_parm in zip(self.atom_list, gb_parms):

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

@@ -579,7 +579,8 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
           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):
+          EwaldErrorTolerance=None, flexibleConstraints=True, rigidWater=True, elements=None,
+          gbsaModel='ACE'):
     """
     Create an OpenMM System from an Amber prmtop file.
 
@@ -603,6 +604,7 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
       verbose (boolean) - if True, print out information on progress (default: False)
       flexibleConstraints (boolean) - if True, flexible bonds will be added in addition ot constrained bonds
       rigidWater (boolean=True) If true, water molecules will be fully rigid regardless of the value passed for the shake argument
+      gbsaModel (str='ACE') The string representing the SA model to use for GB calculations. Must be 'ACE' or None
 
     NOTES
 
@@ -648,6 +650,9 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
         warnings.warn("1-4 scaling parameters in topology file are being ignored. "
             "This is not recommended unless you know what you are doing.")
 
+    if gbmodel is not None and gbsaModel not in ('ACE', None):
+        raise ValueError('gbsaModel must be ACE or None')
+
     has_1264 = 'LENNARD_JONES_CCOEF' in prmtop._raw_data.keys()
     if has_1264:
         parm_ccoef = [float(x) for x in prmtop._raw_data['LENNARD_JONES_CCOEF']]
@@ -974,20 +979,20 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
             if units.is_quantity(cutoff):
                 cutoff = cutoff.value_in_unit(units.nanometers)
         if gbmodel == 'HCT':
-            gb = customgb.GBSAHCTForce(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
+            gb = customgb.GBSAHCTForce(solventDielectric, soluteDielectric, gbsaModel, cutoff, implicitSolventKappa)
         elif gbmodel == 'OBC1':
-            gb = customgb.GBSAOBC1Force(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
+            gb = customgb.GBSAOBC1Force(solventDielectric, soluteDielectric, gbsaModel, cutoff, implicitSolventKappa)
         elif gbmodel == 'OBC2':
             if implicitSolventKappa > 0:
-                gb = customgb.GBSAOBC2Force(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
+                gb = customgb.GBSAOBC2Force(solventDielectric, soluteDielectric, gbsaModel, cutoff, implicitSolventKappa)
             else:
                 gb = mm.GBSAOBCForce()
                 gb.setSoluteDielectric(soluteDielectric)
                 gb.setSolventDielectric(solventDielectric)
         elif gbmodel == 'GBn':
-            gb = customgb.GBSAGBnForce(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
+            gb = customgb.GBSAGBnForce(solventDielectric, soluteDielectric, gbsaModel, cutoff, implicitSolventKappa)
         elif gbmodel == 'GBn2':
-            gb = customgb.GBSAGBn2Force(solventDielectric, soluteDielectric, 'ACE', cutoff, implicitSolventKappa)
+            gb = customgb.GBSAGBn2Force(solventDielectric, soluteDielectric, gbsaModel, cutoff, implicitSolventKappa)
         else:
             raise ValueError("Illegal value specified for implicit solvent model")
         if isinstance(gb, mm.GBSAOBCForce):

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -88,8 +88,8 @@ class TestAmberPrmtopFile(unittest.TestCase):
         parameter.
 
         """
-        for implicitSolvent_value in [HCT, OBC1, OBC2, GBn]:
-            system = prmtop2.createSystem(implicitSolvent=implicitSolvent_value)
+        for implicitSolvent_value, gbsa in zip([HCT, OBC1, OBC2, GBn], ['ACE', None, 'ACE', None]):
+            system = prmtop2.createSystem(implicitSolvent=implicitSolvent_value, gbsaModel=gbsa)
             forces = system.getForces()
             if implicitSolvent_value in set([HCT, OBC1, GBn]):
                 force_type = CustomGBForce

wrappers/python/tests/TestCharmmFiles.py

@@ -61,7 +61,7 @@ class TestCharmmFiles(unittest.TestCase):
         """Test implicit solvent using the implicitSolvent parameter.
 
         """
-        system = self.psf_v.createSystem(self.params, implicitSolvent=OBC2)
+        system = self.psf_v.createSystem(self.params, implicitSolvent=OBC2, gbsaModel='ACE')
         self.assertTrue(any(isinstance(f, CustomGBForce) for f in system.getForces()))
 
     def test_ImplicitSolventParameters(self):