Re-combine electrostatics terms in some plugins

devtools/conda-envs/test_env.yaml

@@ -19,5 +19,6 @@ dependencies:
   - pytest
   - pytest-cov
   - pytest-xdist
+  - pytest-randomly
 
   - mypy =1.2

smirnoff_plugins/_tests/handlers/test_nonbonded.py

@@ -547,7 +547,9 @@ def test_axilrodteller_energies():
 
     distances = [3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 7.0]
 
-    energies = [0.1 * 11 / 8 * (r / 10) ** (-9) for r in distances] * unit.kilojoule_per_mole
+    energies = [
+        0.1 * 11 / 8 * (r / 10) ** (-9) for r in distances
+    ] * unit.kilojoule_per_mole
 
     for energy, distance in zip(energies, distances):
         omm_context.setPositions(

smirnoff_plugins/collections/nonbonded.py

@@ -87,6 +87,58 @@ def create(
 
         return handler
 
+    def _recombine_electrostatics_1_4(
+        self,
+        system: openmm.System,
+    ):
+        """
+        Re-combine the CustomBondForce holding 1-4 electrostatics interactions into the NonbondedForce holding the
+        pairwise intermolecule electrostatics interactions. Leave intact the CustomNonbondedForce and CustomBondForce
+        holding the intermolecular and intramolecular (1-4) vdW interactions.
+
+        Parameters
+        ----------
+        system : openmm.System
+            The system to modify.
+        See for context: https://github.com/openforcefield/openff-interchange/issues/863
+        """
+        import re
+
+        for force_index in range(system.getNumForces()):
+            force = system.getForce(force_index)
+
+            if isinstance(force, openmm.CustomBondForce):
+                if bool(
+                    re.match(
+                        "138.*qq/r",
+                        force.getEnergyFunction(),
+                    )
+                ):
+                    electrostatics_14_force_index = force_index
+                    continue
+
+            elif isinstance(force, openmm.NonbondedForce):
+                electrostatics_force = force
+
+        for exception_index in range(electrostatics_force.getNumExceptions()):
+            particle1, particle2, *_ = electrostatics_force.getExceptionParameters(
+                exception_index
+            )
+
+            charge1 = electrostatics_force.getParticleParameters(particle1)[0]
+            charge2 = electrostatics_force.getParticleParameters(particle2)[0]
+
+            electrostatics_force.setExceptionParameters(
+                exception_index,
+                particle1,
+                particle2,
+                charge1 * charge2 * self.scale_14,
+                0.0,
+                0.0,
+            )
+
+        system.removeForce(electrostatics_14_force_index)
+
 
 class SMIRNOFFDampedBuckingham68Collection(_NonbondedPlugin):
     """Collection storing damped Buckingham potentials."""
@@ -170,6 +222,16 @@ def modify_parameters(
             for name in self.potential_parameters()
         }
 
+    def modify_openmm_forces(
+        self,
+        interchange: Interchange,
+        system: openmm.System,
+        add_constrained_forces: bool,
+        constrained_pairs: Set[Tuple[int, ...]],
+        particle_map: Dict[Union[int, "VirtualSiteKey"], int],
+    ):
+        self._recombine_electrostatics_1_4(system)
+
 
 class SMIRNOFFDoubleExponentialCollection(_NonbondedPlugin):
     """Handler storing vdW potentials as produced by a SMIRNOFF force field."""
@@ -238,6 +300,16 @@ def modify_parameters(
             ),
         }
 
+    def modify_openmm_forces(
+        self,
+        interchange: Interchange,
+        system: openmm.System,
+        add_constrained_forces: bool,
+        constrained_pairs: Set[Tuple[int, ...]],
+        particle_map: Dict[Union[int, "VirtualSiteKey"], int],
+    ):
+        self._recombine_electrostatics_1_4(system)
+
 
 class SMIRNOFFDampedExp6810Collection(_NonbondedPlugin):
     """