Merge pull request #1998 from peastman/switch

Added switchDistance option to ForceField.createSystem()

docs-source/usersguide/application.rst

@@ -792,6 +792,16 @@ The error tolerance is roughly equal to the fractional error in the forces due
 to truncating the Ewald summation.  If you do not specify it, a default value of
 0.0005 is used.
 
+Another optional parameter when using a cutoff is :code:`switchDistance`.  This
+causes Lennard-Jones interactions to smoothly go to zero over some finite range,
+rather than being sharply truncated at the cutoff distance.  This can improve
+energy conservation.  To use it, specify a distance at which the interactions
+should start being reduced.  For example:
+::
+
+    system = prmtop.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
+            switchDistance=0.9*nanometer)
+
 
 Nonbonded Forces for AMOEBA
 ---------------------------

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

@@ -1011,7 +1011,7 @@ class ForceField(object):
 
     def createSystem(self, topology, nonbondedMethod=NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
                      constraints=None, rigidWater=True, removeCMMotion=True, hydrogenMass=None, residueTemplates=dict(),
-                     ignoreExternalBonds=False, **args):
+                     ignoreExternalBonds=False, switchDistance=None, flexibleConstraints=False, **args):
         """Construct an OpenMM System representing a Topology with this force field.
 
         Parameters
@@ -1047,6 +1047,9 @@ class ForceField(object):
             not terminated properly.  This option can create ambiguities where multiple
             templates match the same residue.  If that happens, use the residueTemplates
             argument to specify which one to use.
+        switchDistance : float=None
+            The distance at which the potential energy switching function is turned on for
+            Lennard-Jones interactions. If this is None, no switching function will be used.
         flexibleConstraints : boolean=False
             If True, parameters for constrained degrees of freedom will be added to the System
         args
@@ -1059,6 +1062,8 @@ class ForceField(object):
         system
             the newly created System
         """
+        args['switchDistance'] = switchDistance
+        args['flexibleConstraints'] = flexibleConstraints
         data = ForceField._SystemData()
         data.atoms = list(topology.atoms())
         for atom in data.atoms:
@@ -2322,6 +2327,9 @@ class NonbondedGenerator(object):
             force.addParticle(values[0], values[1], values[2])
         force.setNonbondedMethod(methodMap[nonbondedMethod])
         force.setCutoffDistance(nonbondedCutoff)
+        if args['switchDistance'] is not None:
+            force.setUseSwitchingFunction(True)
+            force.setSwitchingDistance(args['switchDistance'])
         if 'ewaldErrorTolerance' in args:
             force.setEwaldErrorTolerance(args['ewaldErrorTolerance'])
         if 'useDispersionCorrection' in args:
@@ -2439,6 +2447,9 @@ class LennardJonesGenerator(object):
             self.force.setNonbondedMethod(mm.CustomNonbondedForce.CutoffNonPeriodic)
         else:
             raise AssertionError('Unrecognized nonbonded method [%s]' % nonbondedMethod)
+        if args['switchDistance'] is not None:
+            force.setUseSwitchingFunction(True)
+            force.setSwitchingDistance(args['switchDistance'])
 
         # Add the particles
 

wrappers/python/tests/TestForceField.py

@@ -76,6 +76,21 @@ class TestForceField(unittest.TestCase):
                     cutoff_distance = force.getCutoffDistance()
             self.assertEqual(cutoff_distance, cutoff_check)
 
+    def test_SwitchingDistance(self):
+        """Test that the switchDistance parameter is processed correctly."""
+
+        for switchDistance in [None, 0.9*nanometers]:
+            system = self.forcefield1.createSystem(self.pdb1.topology,
+                                                   nonbondedMethod=PME,
+                                                   switchDistance=switchDistance)
+            for force in system.getForces():
+                if isinstance(force, NonbondedForce):
+                    if switchDistance is None:
+                        self.assertFalse(force.getUseSwitchingFunction())
+                    else:
+                        self.assertTrue(force.getUseSwitchingFunction())
+                        self.assertEqual(switchDistance, force.getSwitchingDistance())
+
     def test_RemoveCMMotion(self):
         """Test both options (True and False) for the removeCMMotion parameter."""
         for b in [True, False]: