Merge pull request #1484 from swails/gbneck_errchk

Fix error check for GBneck parameters

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

@@ -335,7 +335,7 @@ def _mbondi3_radii(topology):
 
 def _createEnergyTerms(force, solventDielectric, soluteDielectric, SA, cutoff, kappa, offset):
     # 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
@@ -535,7 +535,7 @@ class GBSAGBnForce(CustomAmberGBForce):
 
     def addParticle(self, parameters):
         parameters = CustomAmberGBForce.addParticle(self, parameters)
-        if parameters[1] < 0.1 or parameters[1] > 0.2:
+        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):

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -361,5 +361,17 @@ class TestAmberPrmtopFile(unittest.TestCase):
             # Make sure it says something about chamber
             self.assertTrue('chamber' in str(e).lower())
 
+    def testGBneckRadii(self):
+        """ Tests that GBneck radii limits are correctly enforced """
+        from simtk.openmm.app.internal.customgbforces import GBSAGBnForce
+        f = GBSAGBnForce()
+        # Make sure legal parameters do not raise
+        f.addParticle([0, 0.1, 0.5])
+        f.addParticle([0, 0.2, 0.5])
+        f.addParticle([0, 0.15, 0.5])
+        # Now make sure that out-of-range parameters *do* raise
+        self.assertRaises(ValueError, lambda: f.addParticle([0, 0.9, 0.5]))
+        self.assertRaises(ValueError, lambda: f.addParticle([0, 0.21, 0.5]))
+
 if __name__ == '__main__':
     unittest.main()