Port #896 to biopolymer-topology-refactor branch (#1198)

* Port PR #896 to biopolymer-topology-refactor branch

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

* Fix GBSA tests

* Fix logic of calling setUseSwitchingFunction

* Condense switching function logic, update tests

* Update release history

* Remove an annotation

* Add comments about switching function with LJ-PME

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

docs/releasehistory.md

@@ -73,6 +73,8 @@ 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.
 
 ### Behaviors changed and bugfixes
 

openff/toolkit/tests/test_forcefield.py

@@ -3647,6 +3647,12 @@ def test_freesolv_gbsa_energies(
         ff = ForceField(
             "test_forcefields/test_forcefield.offxml", off_gbsas[gbsa_model]
         )
+
+        # OpenMM 7.7 and older don't properly handle parsing prmtop files with a GBSA model and
+        # a switching function; this bug may be fixed but for testing purposes, simply turn off
+        # the switching function
+        ff["vdW"].switch_width = 0.0 * unit.angstrom
+
         off_top = molecule.to_topology()
         if is_periodic:
             off_top.box_vectors = (
@@ -3695,6 +3701,13 @@ def test_freesolv_gbsa_energies(
             amber_nb_method = openmm.app.forcefield.NoCutoff
             amber_cutoff = None
 
+        # OpenMM will process False to mean "don't use one" and a quantity as "use with this distance"
+        switching_distance: Union[
+            bool, openmm_unit.Quantity
+        ] = off_nonbonded_force.getUseSwitchingFunction()
+        if switching_distance:
+            switching_distance = off_nonbonded_force.getSwitchingDistance()
+
         (
             amber_omm_system,
             amber_omm_topology,
@@ -3707,6 +3720,7 @@ def test_freesolv_gbsa_energies(
             nonbondedCutoff=amber_cutoff,
             gbsaModel="ACE",
             implicitSolventKappa=0.0,
+            switchDistance=switching_distance,
         )
 
         # Retrieve the GBSAForce from both the AMBER and OpenForceField systems
@@ -3782,8 +3796,6 @@ def test_freesolv_gbsa_energies(
 
         off_nonbonded_force.setReactionFieldDielectric(1.0)
 
-        # TODO: look into if switching distance is relevant #882
-
         # Create Contexts
         integrator = openmm.VerletIntegrator(1.0 * openmm_unit.femtoseconds)
         platform = Platform.getPlatformByName("Reference")

openff/toolkit/tests/test_parameters.py

@@ -1743,6 +1743,70 @@ 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),
+            (15.0 * unit.angstrom, 5.0 * unit.angstrom, True, 10.0 * unit.angstrom),
+            (9.0 * unit.angstrom, 0.0 * unit.angstrom, False, 0.0 * unit.angstrom),
+        ],
+    )
+    @pytest.mark.parametrize(
+        "method",
+        # It's possible that this test should not cover PME (LJ-PME), see comments in parameters.py
+        ["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()
+
+        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() == expected_use
+
+        if expected_use:
+
+            assert numpy.isclose(
+                nonbonded_force.getSwitchingDistance().value_in_unit(
+                    openmm_unit.angstrom
+                ),
+                expected_switching_distance.m_as(unit.angstrom),
+            )
+
     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

@@ -3647,7 +3647,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,
@@ -3684,8 +3684,15 @@ def create_force(self, system, topology, **kwargs):
             else:
                 force.setNonbondedMethod(openmm.NonbondedForce.PME)
                 force.setUseDispersionCorrection(True)
+
                 force.setCutoffDistance(to_openmm(self.cutoff))
 
+        # This applies a switching function whether the vdW method is LJ-PME or cut-off. It's not clear if this is a
+        # valid combination of settings, if this is something that all engines would support, or if it even has an
+        # effect. In OpenMM, it can be applied but it might not have any effect. The SMIRNOFF spec offers no clear
+        # guidance on this combination, but it is possible that is may be revised in the future to do so.
+        self._apply_switching_function(force)
+
         # Iterate over all defined Lennard-Jones types, allowing later matches to override earlier ones.
         atom_matches = self.find_matches(topology)
 
@@ -3706,6 +3713,15 @@ def create_force(self, system, topology, **kwargs):
             atom_matches, topology, [(atom,) for atom in topology.atoms]
         )
 
+    def _apply_switching_function(self, force):
+        """Apply a switching function to a NonbondedForce if self.switch_width is nonzero."""
+        from openff.units.openmm import to_openmm
+
+        if self.switch_width.m > 0:
+            switching_distance = to_openmm(self.cutoff - self.switch_width)
+            force.setSwitchingDistance(switching_distance)
+            force.setUseSwitchingFunction(True)
+
 
 class ElectrostaticsHandler(_NonbondedHandler):
     """Handles SMIRNOFF ``<Electrostatics>`` tags.