Properly set the cutoff distance in HTF nonbonded forces

openfe/protocols/openmm_rfe/_rfe_utils/relative.py

@@ -900,25 +900,21 @@ def _add_nonbonded_force_terms(self):
         # Select functional form based on nonbonded method.
         # TODO: check _nonbonded_custom_ewald and _nonbonded_custom_cutoff
         # since they take arguments that are never used...
+        r_cutoff = self._old_system_forces['NonbondedForce'].getCutoffDistance()
+        sterics_energy_expression = self._nonbonded_custom(self._softcore_LJ_v2)
         if self._nonbonded_method in [openmm.NonbondedForce.NoCutoff]:
             sterics_energy_expression = self._nonbonded_custom(
                 self._softcore_LJ_v2)
         elif self._nonbonded_method in [openmm.NonbondedForce.CutoffPeriodic,
                                         openmm.NonbondedForce.CutoffNonPeriodic]:
             epsilon_solvent = self._old_system_forces['NonbondedForce'].getReactionFieldDielectric()
-            r_cutoff = self._old_system_forces['NonbondedForce'].getCutoffDistance()
-            sterics_energy_expression = self._nonbonded_custom(
-                self._softcore_LJ_v2)
             standard_nonbonded_force.setReactionFieldDielectric(
                 epsilon_solvent)
             standard_nonbonded_force.setCutoffDistance(r_cutoff)
         elif self._nonbonded_method in [openmm.NonbondedForce.PME,
                                         openmm.NonbondedForce.Ewald]:
             [alpha_ewald, nx, ny, nz] = self._old_system_forces['NonbondedForce'].getPMEParameters()
             delta = self._old_system_forces['NonbondedForce'].getEwaldErrorTolerance()
-            r_cutoff = self._old_system_forces['NonbondedForce'].getCutoffDistance()
-            sterics_energy_expression = self._nonbonded_custom(
-                self._softcore_LJ_v2)
             standard_nonbonded_force.setPMEParameters(alpha_ewald, nx, ny, nz)
             standard_nonbonded_force.setEwaldErrorTolerance(delta)
             standard_nonbonded_force.setCutoffDistance(r_cutoff)
@@ -939,6 +935,8 @@ def _add_nonbonded_force_terms(self):
 
         sterics_custom_nonbonded_force = openmm.CustomNonbondedForce(
             total_sterics_energy)
+        # Match cutoff from non-custom NB forces
+        sterics_custom_nonbonded_force.setCutoffDistance(r_cutoff)
 
         if self._softcore_LJ_v2:
             sterics_custom_nonbonded_force.addGlobalParameter(

openfe/tests/protocols/test_openmm_equil_rfe_protocols.py

@@ -547,6 +547,45 @@
         assert e1 == e2 == 0
 
 
+@pytest.mark.slow
+@pytest.mark.parametrize('cutoff',
+    [1.0 * unit.nanometer,
+     12.0 * unit.angstrom,
+     0.9 * unit.nanometer]
+)
+def test_dry_run_ligand_system_cutoff(
+    cutoff, benzene_system, toluene_system, benzene_to_toluene_mapping, tmpdir
+):
+    """
+    Test that the right nonbonded cutoff is propagated to the hybrid system.
+    """
+    settings = openmm_rfe.RelativeHybridTopologyProtocol.default_settings()
+    settings.solvation_settings.solvent_padding = 1.5 * unit.nanometer
+    settings.system_settings.nonbonded_cutoff = cutoff
+
+    protocol = openmm_rfe.RelativeHybridTopologyProtocol(
+            settings=settings,
+    )
+    dag = protocol.create(
+        stateA=benzene_system,
+        stateB=toluene_system,
+        mapping={'ligand': benzene_to_toluene_mapping},
+    )
+    dag_unit = list(dag.protocol_units)[0]
+
+    with tmpdir.as_cwd():
+        sampler = dag_unit.run(dry=True)['debug']['sampler']
+        hs = sampler._factory.hybrid_system
+
+        nbfs = [f for f in hs.getForces() if
+                isinstance(f, CustomNonbondedForce) or
+                isinstance(f, NonbondedForce)]
+
+        for f in nbfs:
+            f_cutoff = from_openmm(f.getCutoffDistance())
+            assert f_cutoff == cutoff
+
+
 @pytest.mark.flaky(reruns=3)  # bad minimisation can happen
 def test_dry_run_user_charges(benzene_modifications, tmpdir):
     """