Merge pull request #1272 from swails/amber-switching

Amber switching

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

@@ -38,7 +38,7 @@ from simtk.openmm.app import PDBFile
 from simtk.openmm.app.internal import amber_file_parser
 from . import forcefield as ff
 from . import element as elem
-import simtk.unit as unit
+import simtk.unit as u
 import simtk.openmm as mm
 from simtk.openmm.app.internal.unitcell import computePeriodicBoxVectors
 
@@ -69,6 +69,15 @@ class GBn2(object):
         return 'GBn2'
 GBn2 = GBn2()
 
+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
+
 class AmberPrmtopFile(object):
     """AmberPrmtopFile parses an AMBER prmtop file and constructs a Topology and (optionally) an OpenMM System from it."""
 
@@ -145,12 +154,13 @@ class AmberPrmtopFile(object):
             box = prmtop.getBoxBetaAndDimensions()
             top.setPeriodicBoxVectors(computePeriodicBoxVectors(*(box[1:4] + box[0:1]*3)))
 
-    def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
+    def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*u.nanometer,
                      constraints=None, rigidWater=True, implicitSolvent=None,
-                     implicitSolventSaltConc=0.0*(unit.moles/unit.liter),
-                     implicitSolventKappa=None, temperature=298.15*unit.kelvin,
+                     implicitSolventSaltConc=0.0*(u.moles/u.liter),
+                     implicitSolventKappa=None, temperature=298.15*u.kelvin,
                      soluteDielectric=1.0, solventDielectric=78.5,
-                     removeCMMotion=True, hydrogenMass=None, ewaldErrorTolerance=0.0005):
+                     removeCMMotion=True, hydrogenMass=None, ewaldErrorTolerance=0.0005,
+                     switchDistance=0.0*u.nanometer):
         """Construct an OpenMM System representing the topology described by this
         prmtop file.
 
@@ -192,6 +202,11 @@ class AmberPrmtopFile(object):
             total mass the same.
         ewaldErrorTolerance : float=0.0005
             The error tolerance to use if nonbondedMethod is Ewald or PME.
+        switchDistance : float=0*nanometers
+            The distance at which the potential energy switching function is
+            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
 
         Returns
         -------
@@ -237,10 +252,10 @@ class AmberPrmtopFile(object):
             raise ValueError('Illegal value for implicit solvent model')
         # If implicitSolventKappa is None, compute it from the salt concentration
         if implicitSolvent is not None and implicitSolventKappa is None:
-            if unit.is_quantity(implicitSolventSaltConc):
-                implicitSolventSaltConc = implicitSolventSaltConc.value_in_unit(unit.moles/unit.liter)
-            if unit.is_quantity(temperature):
-                temperature = temperature.value_in_unit(unit.kelvin)
+            if u.is_quantity(implicitSolventSaltConc):
+                implicitSolventSaltConc = implicitSolventSaltConc.value_in_unit(u.moles/u.liter)
+            if u.is_quantity(temperature):
+                temperature = temperature.value_in_unit(u.kelvin)
             # The constant is 1 / sqrt( epsilon_0 * kB / (2 * NA * q^2 * 1000) )
             # where NA is avogadro's number, epsilon_0 is the permittivity of
             # free space, q is the elementary charge (this number matches
@@ -271,4 +286,17 @@ class AmberPrmtopFile(object):
                 force.setEwaldErrorTolerance(ewaldErrorTolerance)
         if removeCMMotion:
             sys.addForce(mm.CMMotionRemover())
+
+        if switchDistance and nonbondedMethod is not ff.NoCutoff:
+            # make sure it's legal
+            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!')
+            force.setUseSwitchingFunction(True)
+            force.setSwitchingDistance(switchDistance)
+
         return sys

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -296,6 +296,38 @@ class TestAmberPrmtopFile(unittest.TestCase):
                 diff = norm(f1-f2)
                 self.assertTrue(diff < 0.1 or diff/norm(f1) < 1e-4)
 
+    def testSwitchFunction(self):
+        """ Tests the switching function option in AmberPrmtopFile """
+        system = prmtop1.createSystem(nonbondedMethod=PME,
+                                      nonbondedCutoff=1*nanometer,
+                                      switchDistance=0.8*nanometer)
+        for force in system.getForces():
+            if isinstance(force, NonbondedForce):
+                self.assertTrue(force.getUseSwitchingFunction())
+                self.assertEqual(force.getSwitchingDistance(), 0.8*nanometer)
+                break
+        else:
+            assert False, 'Did not find expected nonbonded force!'
+
+        # Check error handling
+        system = prmtop1.createSystem(nonbondedMethod=PME,
+                                      nonbondedCutoff=1*nanometer)
+        for force in system.getForces():
+            if isinstance(force, NonbondedForce):
+                self.assertFalse(force.getUseSwitchingFunction())
+                break
+        else:
+            assert False, 'Did not find expected nonbonded force!'
+
+        self.assertRaises(ValueError, lambda:
+                prmtop1.createSystem(nonbondedMethod=PME,
+                    nonbondedCutoff=1*nanometer, switchDistance=-1)
+        )
+        self.assertRaises(ValueError, lambda:
+                prmtop1.createSystem(nonbondedMethod=PME,
+                    nonbondedCutoff=1*nanometer, switchDistance=1.2)
+        )
+
     def test_with_dcd_reporter(self):
         """Check that an amber simulation like the docs example works with a DCD reporter."""