Add code supporting input with and without units for CharmmPsfFile.

Add relevant tests (that used to fail). This also fixes a subtle bug where
switching functions were not applied when NBFIXes were specified.

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

@@ -102,6 +102,15 @@ class _ZeroDict(dict):
                 return [0, 0], []
             return 0, []
 
+def _strip_optunit(thing, unit):
+    """
+    Strips optional units, converting to specified unit type. If no unit
+    present, it just returns the number
+    """
+    if u.is_quantity(thing):
+        return thing.value_in_unit(unit)
+    return thing
+
 # ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 _resre = re.compile(r'(\d+)([a-zA-Z]*)')
@@ -1026,10 +1035,11 @@ class CharmmPsfFile(object):
             # See if we need to use a switching function
             if switchDistance and nonbondedMethod is not ff.NoCutoff:
                 # make sure it's legal
-                if switchDistance >= nonbondedCutoff:
+                if (_strip_optunit(switchDistance, u.nanometer) >=
+                        _strip_optunit(nonbondedCutoff, u.nanometer)):
                     raise ValueError('switchDistance is too large compared '
                                      'to the cutoff!')
-                if abs(switchDistance) != switchDistance:
+                if _strip_optunit(switchDistance, u.nanometer) < 0:
                     # Detects negatives for both Quantity and float
                     raise ValueError('switchDistance must be non-negative!')
                 force.setUseSwitchingFunction(True)
@@ -1070,10 +1080,11 @@ class CharmmPsfFile(object):
             # See if we need to use a switching function
             if switchDistance and nonbondedMethod is not ff.NoCutoff:
                 # make sure it's legal
-                if switchDistance >= nonbondedCutoff:
+                if (_strip_optunit(switchDistance, u.nanometer) >=
+                        _strip_optunit(nonbondedCutoff, u.nanometer)):
                     raise ValueError('switchDistance is too large compared '
                                      'to the cutoff!')
-                if abs(switchDistance) != switchDistance:
+                if _strip_optunit(switchDistance, u.nanometer) < 0:
                     # Detects negatives for both Quantity and float
                     raise ValueError('switchDistance must be non-negative!')
                 force.setUseSwitchingFunction(True)
@@ -1156,9 +1167,13 @@ class CharmmPsfFile(object):
                 raise ValueError('Unrecognized nonbonded method')
             if switchDistance and nonbondedMethod is not ff.NoCutoff:
                 # make sure it's legal
-                if switchDistance >= nonbondedCutoff:
+                if (_strip_optunit(switchDistance, u.nanometer) >=
+                        _strip_optunit(nonbondedCutoff, u.nanometer)):
                     raise ValueError('switchDistance is too large compared '
                                      'to the cutoff!')
+                if _strip_optunit(switchDistance, u.nanometer) < 0:
+                    # Detects negatives for both Quantity and float
+                    raise ValueError('switchDistance must be non-negative!')
                 cforce.setUseSwitchingFunction(True)
                 cforce.setSwitchingDistance(switchDistance)
             for i in lj_idx_list:

wrappers/python/tests/TestCharmmFiles.py

@@ -4,6 +4,7 @@ from simtk.openmm.app import *
 from simtk.openmm import *
 from simtk.unit import *
 import simtk.openmm.app.element as elem
+import warnings
 
 class TestCharmmFiles(unittest.TestCase):
 
@@ -90,7 +91,6 @@ class TestCharmmFiles(unittest.TestCase):
 
     def test_NBFIX(self):
         """Tests CHARMM systems with NBFIX Lennard-Jones modifications"""
-        import warnings
         warnings.filterwarnings('ignore', category=CharmmPSFWarning)
         psf = CharmmPsfFile('systems/ala3_solv.psf')
         crd = CharmmCrdFile('systems/ala3_solv.crd')
@@ -122,6 +122,34 @@ class TestCharmmFiles(unittest.TestCase):
         self.assertEqual(len(list(psf.topology.residues())), 20169)
         self.assertEqual(len(list(psf.topology.bonds())), 46634)
 
+    def testSystemOptions(self):
+        """ Test various options in CharmmPsfFile.createSystem """
+        warnings.filterwarnings('ignore', category=CharmmPSFWarning)
+        psf = CharmmPsfFile('systems/ala3_solv.psf')
+        crd = CharmmCrdFile('systems/ala3_solv.crd')
+        params = CharmmParameterSet('systems/par_all36_prot.prm',
+                                    'systems/toppar_water_ions.str')
+        # Box dimensions (found from bounding box)
+        psf.setBox(32.7119500*angstroms, 32.9959600*angstroms, 33.0071500*angstroms)
+
+        # Check some illegal options
+        self.assertRaises(ValueError, lambda:
+                    psf.createSystem(params, nonbondedMethod=5))
+        self.assertRaises(TypeError, lambda:
+                    psf.createSystem(params, nonbondedMethod=PME,
+                                     nonbondedCutoff=1*radian)
+        )
+        self.assertRaises(TypeError, lambda:
+                    psf.createSystem(params, nonbondedMethod=PME,
+                                     switchDistance=1*radian)
+        )
+
+        # Check what should be some legal options
+        psf.createSystem(params, nonbondedMethod=PME, switchDistance=0.8,
+                         nonbondedCutoff=1.2)
+        psf.createSystem(params, nonbondedMethod=PME, switchDistance=0.8,
+                         nonbondedCutoff=1.2*nanometer)
+
     def test_ImplicitSolventForces(self):
         """Compute forces for different implicit solvent types, and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""
         solventType = [HCT, OBC1, OBC2, GBn, GBn2]