Fixed bugs setting up nonbonded and GB forces

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

@@ -135,8 +135,6 @@ class AmberPrmtopFile(object):
             raise ValueError('Illegal value for nonbonded method')
         if not self._prmtop.getIfBox() and nonbondedMethod in (ff.CutoffPeriodic, ff.Ewald, ff.PME):
             raise ValueError('Illegal nonbonded method for a non-periodic system')
-        if nonbondedMethod == ff.NoCutoff:
-            nonbondedCutoff = None
         constraintMap = {None:None,
                          ff.HBonds:'h-bonds',
                          ff.AllBonds:'all-bonds',

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

@@ -676,6 +676,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
             min_box_width = min([boxX / units.nanometers, boxY / units.nanometers, boxZ / units.nanometers])
             CLEARANCE_FACTOR = 0.97 # reduce the cutoff to be a bit smaller than 1/2 smallest box length            
             nonbondedCutoff = units.Quantity((min_box_width * CLEARANCE_FACTOR) / 2.0, units.nanometers)
+        if nonbondedMethod != 'NoCutoff':
             force.setCutoffDistance(nonbondedCutoff)
         
         # Set nonbonded method.
@@ -789,15 +790,15 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
             symbls = prmtop.getAtomTypes()
         gb_parms = prmtop.getGBParms(symbls)
         if gbmodel == 'HCT':
-            gb = customgb.GBSAHCTForce(solventDielectric, soluteDielectric, True)
+            gb = customgb.GBSAHCTForce(solventDielectric, soluteDielectric, 'ACE')
         elif gbmodel == 'OBC1':
-            gb = customgb.GBSAOBC1Force(solventDielectric, soluteDielectric, True)
+            gb = customgb.GBSAOBC1Force(solventDielectric, soluteDielectric, 'ACE')
         elif gbmodel == 'OBC2':
             gb = mm.GBSAOBCForce()
             gb.setSoluteDielectric(soluteDielectric)
             gb.setSolventDielectric(solventDielectric)
         elif gbmodel == 'GBn':
-            gb = customgb.GBSAGBnForce(solventDielectric, soluteDielectric, True)
+            gb = customgb.GBSAGBnForce(solventDielectric, soluteDielectric, 'ACE')
         else:
             raise Exception("Illegal value specified for implicit solvent model")
         for iAtom in range(prmtop.getNumAtoms()):

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

@@ -217,6 +217,8 @@ def GBSAHCTForce(solventDielectric=78.5, soluteDielectric=1, SA=None):
     custom.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce.SingleParticle)
     if SA=='ACE':
         custom.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6", CustomGBForce.SingleParticle)
+    elif SA is not None:
+        raise ValueError('Unknown surface area method: '+SA)
     custom.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
                            "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
 
@@ -248,6 +250,8 @@ def GBSAOBC1Force(solventDielectric=78.5, soluteDielectric=1, SA=None):
     custom.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce.SingleParticle)
     if SA=='ACE':
         custom.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6", CustomGBForce.SingleParticle)
+    elif SA is not None:
+        raise ValueError('Unknown surface area method: '+SA)
     custom.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
                            "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
 
@@ -279,6 +283,8 @@ def GBSAOBC2Force(solventDielectric=78.5, soluteDielectric=1, SA=None):
     custom.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce.SingleParticle)
     if SA=='ACE':
         custom.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6", CustomGBForce.SingleParticle)
+    elif SA is not None:
+        raise ValueError('Unknown surface area method: '+SA)
     custom.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
                            "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)
 
@@ -329,6 +335,8 @@ def GBSAGBnForce(solventDielectric=78.5, soluteDielectric=1, SA=None):
     custom.addEnergyTerm("-0.5*138.935485*(1/soluteDielectric-1/solventDielectric)*q^2/B", CustomGBForce.SingleParticle)
     if SA=='ACE':
         custom.addEnergyTerm("28.3919551*(radius+0.14)^2*(radius/B)^6", CustomGBForce.SingleParticle)
+    elif SA is not None:
+        raise ValueError('Unknown surface area method: '+SA)
     custom.addEnergyTerm("-138.935485*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;"
                            "f=sqrt(r^2+B1*B2*exp(-r^2/(4*B1*B2)))", CustomGBForce.ParticlePairNoExclusions)