Port PR #896 to biopolymer-topology-refactor branch

Co-authored-by: SimonBoothroyd <[email protected]>

docs/releasehistory.md

@@ -73,6 +73,11 @@ print(value_roundtrip)
 - [PR #1192](https://github.com/openforcefield/openforcefield/pull/1192): Add re-exports for core classes to the new
   `openff.toolkit.app` module and re-exports for parameter types to the new `openff.toolkit.topology.parametertypes` module.
   This does not affect existing paths and gives some new, easier to remember paths to core objects.
+- [PR #1198](https://github.com/openforcefield/openforcefield/pull/1198): Ensure the vdW switch
+  width is correctly set and enabled. If the switch width is greater than the cutoff distance, it is
+  considered invalid and an exception
+  ([`InvalidSwitchingDistanceError`](openff.toolkit.utils.exceptions.InvalidSwitchingDistanceError))
+  is raised.
 
 ### Behaviors changed and bugfixes
 

openff/toolkit/tests/test_parameters.py

@@ -40,6 +40,7 @@
     DuplicateParameterError,
     IncompatibleParameterError,
     IncompatibleUnitError,
+    InvalidSwitchingDistanceError,
     MissingIndexedAttributeError,
     NotEnoughPointsForInterpolationError,
     ParameterLookupError,
@@ -1743,6 +1744,74 @@ class TestvdWType:
     Test the behavior of vdWType
     """
 
+    @pytest.mark.parametrize(
+        "cutoff,switch_width,expected_use,expected_switching_distance",
+        [
+            (9.0 * unit.angstrom, 1.0 * unit.angstrom, True, 8.0 * unit.angstrom),
+            (11.0 * unit.angstrom, 2.0 * unit.angstrom, True, 9.0 * unit.angstrom),
+            (9.0 * unit.angstrom, 9.0 * unit.angstrom, False, 0.0 * unit.angstrom),
+            (9.0 * unit.angstrom, 10.0 * unit.angstrom, False, 0.0 * unit.angstrom),
+        ],
+    )
+    @pytest.mark.parametrize(
+        "method",
+        ["PME", "cutoff"],
+    )
+    def test_switch_width(
+        self, cutoff, switch_width, expected_use, expected_switching_distance, method
+    ):
+        """Test that create_force works on a vdWHandler which has a switch width
+        specified.
+        """
+
+        import openmm
+        from openmm import unit as openmm_unit
+
+        # Create a dummy topology containing only argon and give it a set of
+        # box vectors.
+        topology = Molecule.from_smiles("[Ar]").to_topology()
+        topology.box_vectors = unit.Quantity(numpy.eye(3) * 20 * unit.angstrom)
+
+        # create a VdW handler with only parameters for argon.
+        vdw_handler = vdWHandler(
+            version=0.3,
+            cutoff=cutoff,
+            switch_width=switch_width,
+            method=method,
+        )
+        vdw_handler.add_parameter(
+            {
+                "smirks": "[#18:1]",
+                "epsilon": 1.0 * unit.kilojoules_per_mole,
+                "sigma": 1.0 * unit.angstrom,
+            }
+        )
+
+        omm_sys = openmm.System()
+
+        if expected_use:
+
+            vdw_handler.create_force(omm_sys, topology)
+
+            nonbonded_force = [
+                force
+                for force in omm_sys.getForces()
+                if isinstance(force, openmm.NonbondedForce)
+            ][0]
+
+            assert nonbonded_force.getUseSwitchingFunction()
+
+            assert numpy.isclose(
+                nonbonded_force.getSwitchingDistance().value_in_unit(
+                    openmm_unit.angstrom
+                ),
+                expected_switching_distance.m_as(unit.angstrom),
+            )
+
+        else:
+            with pytest.raises(InvalidSwitchingDistanceError, match=str(cutoff.m)):
+                vdw_handler.create_force(omm_sys, topology)
+
     def test_sigma_rmin_half(self):
         """Test the setter/getter behavior or sigma and rmin_half"""
         from openff.toolkit.typing.engines.smirnoff.parameters import vdWHandler

openff/toolkit/typing/engines/smirnoff/parameters.py

@@ -92,6 +92,7 @@
     FractionalBondOrderInterpolationMethodUnsupportedError,
     IncompatibleParameterError,
     IncompatibleUnitError,
+    InvalidSwitchingDistanceError,
     MissingIndexedAttributeError,
     NonintegralMoleculeChargeException,
     NotEnoughPointsForInterpolationError,
@@ -3647,7 +3648,7 @@ def check_handler_compatibility(self, other_handler):
         """
         float_attrs_to_compare = ["scale12", "scale13", "scale14", "scale15"]
         string_attrs_to_compare = ["potential", "combining_rules", "method"]
-        unit_attrs_to_compare = ["cutoff"]
+        unit_attrs_to_compare = ["cutoff", "switch_width"]
 
         self._check_attributes_are_equal(
             other_handler,
@@ -3676,6 +3677,17 @@ def create_force(self, system, topology, **kwargs):
                 force.setCutoffDistance(to_openmm(self.cutoff))
                 force.setEwaldErrorTolerance(1.0e-4)
 
+                if self.switch_width >= self.cutoff:
+                    raise InvalidSwitchingDistanceError(
+                        "Attempted to generate a NonbondedForce with a switching width greater than the or equal to the cutoff "
+                        f"distance. Found switch_width={self.switch_width} and cutoff={self.cutoff}."
+                    )
+                else:
+                    force.setSwitchingDistance(
+                        to_openmm(self.cutoff - self.switch_width)
+                    )
+                    force.setUseSwitchingFunction(self.switch_width < self.cutoff)
+
         # If method is cutoff, then we currently support openMM's PME for periodic system and NoCutoff for nonperiodic
         elif self.method == "cutoff":
             # If we're given a nonperiodic box, we always set NoCutoff. Later we'll add support for CutoffNonPeriodic
@@ -3684,8 +3696,20 @@ def create_force(self, system, topology, **kwargs):
             else:
                 force.setNonbondedMethod(openmm.NonbondedForce.PME)
                 force.setUseDispersionCorrection(True)
+
                 force.setCutoffDistance(to_openmm(self.cutoff))
 
+                if self.switch_width >= self.cutoff:
+                    raise InvalidSwitchingDistanceError(
+                        "Attempted to generate a NonbondedForce with a switching width greater than the or equal to the cutoff "
+                        f"distance. Found switch_width={self.switch_width} and cutoff={self.cutoff}."
+                    )
+                else:
+                    force.setSwitchingDistance(
+                        to_openmm(self.cutoff - self.switch_width)
+                    )
+                    force.setUseSwitchingFunction(self.switch_width < self.cutoff)
+
         # Iterate over all defined Lennard-Jones types, allowing later matches to override earlier ones.
         atom_matches = self.find_matches(topology)
 

openff/toolkit/utils/exceptions.py

@@ -243,6 +243,12 @@ class SMIRNOFFSpecUnimplementedError(OpenFFToolkitException):
     """
 
 
+class InvalidSwitchingDistanceError(OpenFFToolkitException):
+    """
+    Exception for when a force field specifies a switch distance greater than the cutoff.
+    """
+
+
 class FractionalBondOrderInterpolationMethodUnsupportedError(OpenFFToolkitException):
     """
     Exception for when an unsupported fractional bond order interpolation assignment method is called.