Merge pull request #938 from openforcefield/test-nonbonded-combinations

Error out when combining mismatched non-bonded settings

.github/workflows/ci.yaml

@@ -103,7 +103,7 @@ jobs:
       if: always()
       run: |
         python -m pytest $COV openff/interchange/ \
-          -r fExs -nauto --durations=10 \
+          -r fExs -n logical --durations=10 \
           -m "slow or not slow" \
           --ignore=openff/interchange/_tests/energy_tests/test_energies.py
 

.github/workflows/examples.yaml

@@ -70,7 +70,7 @@ jobs:
     - name: Run example notebooks
       if: always()
       run: |
-        python -m pytest --nbval-lax --dist loadscope -nauto \
+        python -m pytest --nbval-lax --dist loadscope -n logical \
           examples/ \
           --ignore=examples/deprecated/ \
           --ignore=examples/lammps/ \

devtools/conda-envs/examples_env.yaml

@@ -5,9 +5,9 @@ channels:
 dependencies:
   # Core
   - python
-  - numpy =1
-  - pydantic =1
-  - openmm >=7.6
+  - numpy
+  - pydantic =2
+  - openmm
   # OpenFF stack
   - openff-toolkit =0.15.2
   - openff-models
@@ -16,8 +16,9 @@ dependencies:
   - dgl =1
   # Optional features
   - unyt
-  - mbuild
+  - mbuild =0.17
   - foyer >=0.12.1
+  - gmso =0.12
   # Examples
   - jax
   - mdtraj

docs/releasehistory.md

@@ -16,6 +16,7 @@ Please note that all releases prior to a version 1.0.0 are considered pre-releas
 * #933 Fixes #934 in which atom order was sometimes mangled in `Interchange.from_openmm`.
 * #929 A warning is raised when positions are not passed to `Interchange.from_openmm`.
 * #930 Adds `additional_forces` argument to `create_openmm_simulation`.
+* #938 An error is raised when non-bonded settings do not match when using `Interchange.combine`.
 
 ## 0.3.23 - 2024-03-07
 

openff/interchange/_tests/unit_tests/components/test_interchange.py

@@ -1,13 +1,10 @@
-from copy import deepcopy
-
 import numpy
 import pytest
-from openff.toolkit.topology import Molecule, Topology
+from openff.toolkit import Molecule, Topology, unit
 from openff.toolkit.typing.engines.smirnoff.parameters import (
     ElectrostaticsHandler,
     ParameterHandler,
 )
-from openff.units import unit
 from openff.utilities.testing import skip_if_missing
 
 from openff.interchange import Interchange
@@ -100,85 +97,6 @@ def test_box_setter(self):
         with pytest.raises(ValidationError):
             tmp.box = [2, 2, 3, 90, 90, 90]
 
-    @skip_if_missing("openmm")
-    def test_basic_combination(self, monkeypatch, sage_unconstrained):
-        """Test basic use of Interchange.__add__() based on the README example"""
-        monkeypatch.setenv("INTERCHANGE_EXPERIMENTAL", "1")
-
-        mol = Molecule.from_smiles("C")
-        mol.generate_conformers(n_conformers=1)
-        top = Topology.from_molecules([mol])
-
-        interchange = Interchange.from_smirnoff(sage_unconstrained, top)
-
-        interchange.box = [4, 4, 4] * numpy.eye(3)
-        interchange.positions = mol.conformers[0]
-
-        # Copy and translate atoms by [1, 1, 1]
-        other = Interchange()
-        other = deepcopy(interchange)
-        other.positions += 1.0 * unit.nanometer
-
-        combined = interchange.combine(other)
-
-        # Just see if it can be converted into OpenMM and run
-        get_openmm_energies(combined)
-
-    def test_parameters_do_not_clash(self, monkeypatch, sage_unconstrained):
-        monkeypatch.setenv("INTERCHANGE_EXPERIMENTAL", "1")
-
-        thf = Molecule.from_smiles("C1CCOC1")
-        ace = Molecule.from_smiles("CC(=O)O")
-
-        thf.generate_conformers(n_conformers=1)
-        ace.generate_conformers(n_conformers=1)
-
-        def make_interchange(molecule: Molecule) -> Interchange:
-            molecule.generate_conformers(n_conformers=1)
-            interchange = Interchange.from_smirnoff(
-                force_field=sage_unconstrained,
-                topology=[molecule],
-            )
-            interchange.positions = molecule.conformers[0]
-
-            return interchange
-
-        thf_interchange = make_interchange(thf)
-        ace_interchange = make_interchange(ace)
-        complex_interchange = thf_interchange.combine(ace_interchange)
-
-        thf_vdw = thf_interchange["vdW"].get_system_parameters()
-        ace_vdw = ace_interchange["vdW"].get_system_parameters()
-        add_vdw = complex_interchange["vdW"].get_system_parameters()
-
-        numpy.testing.assert_equal(numpy.vstack([thf_vdw, ace_vdw]), add_vdw)
-
-        # TODO: Ensure the de-duplication is maintained after exports
-
-    def test_positions_setting(self, monkeypatch, sage):
-        """Test that positions exist on the result if and only if
-        both input objects have positions."""
-
-        monkeypatch.setenv("INTERCHANGE_EXPERIMENTAL", "1")
-
-        ethane = Molecule.from_smiles("CC")
-        methane = Molecule.from_smiles("C")
-
-        ethane_interchange = Interchange.from_smirnoff(
-            sage,
-            [ethane],
-        )
-        methane_interchange = Interchange.from_smirnoff(sage, [methane])
-
-        ethane.generate_conformers(n_conformers=1)
-        methane.generate_conformers(n_conformers=1)
-
-        assert (methane_interchange.combine(ethane_interchange)).positions is None
-        methane_interchange.positions = methane.conformers[0]
-        assert (methane_interchange.combine(ethane_interchange)).positions is None
-        ethane_interchange.positions = ethane.conformers[0]
-        assert (methane_interchange.combine(ethane_interchange)).positions is not None
-
     def test_input_topology_not_modified(self, sage):
         molecule = Molecule.from_smiles("CCO")
         molecule.generate_conformers(n_conformers=1)

openff/interchange/_tests/unit_tests/operations/test_combine.py

@@ -0,0 +1,137 @@
+import copy
+
+import numpy
+import pytest
+from openff.toolkit import ForceField, Molecule, Topology, unit
+from openff.utilities.testing import skip_if_missing
+
+from openff.interchange import Interchange
+from openff.interchange.drivers import get_openmm_energies
+from openff.interchange.exceptions import (
+    CutoffMismatchError,
+    SwitchingFunctionMismatchError,
+)
+
+
+class TestCombine:
+    @skip_if_missing("openmm")
+    def test_basic_combination(self, monkeypatch, sage_unconstrained):
+        """Test basic use of Interchange.__add__() based on the README example"""
+        monkeypatch.setenv("INTERCHANGE_EXPERIMENTAL", "1")
+
+        mol = Molecule.from_smiles("C")
+        mol.generate_conformers(n_conformers=1)
+        top = Topology.from_molecules([mol])
+
+        interchange = Interchange.from_smirnoff(sage_unconstrained, top)
+
+        interchange.box = [4, 4, 4] * numpy.eye(3)
+        interchange.positions = mol.conformers[0]
+
+        # Copy and translate atoms by [1, 1, 1]
+        other = Interchange()
+        other = copy.deepcopy(interchange)
+        other.positions += 1.0 * unit.nanometer
+
+        combined = interchange.combine(other)
+
+        # Just see if it can be converted into OpenMM and run
+        get_openmm_energies(combined)
+
+    def test_parameters_do_not_clash(self, monkeypatch, sage_unconstrained):
+        monkeypatch.setenv("INTERCHANGE_EXPERIMENTAL", "1")
+
+        thf = Molecule.from_smiles("C1CCOC1")
+        ace = Molecule.from_smiles("CC(=O)O")
+
+        thf.generate_conformers(n_conformers=1)
+        ace.generate_conformers(n_conformers=1)
+
+        def make_interchange(molecule: Molecule) -> Interchange:
+            molecule.generate_conformers(n_conformers=1)
+            interchange = Interchange.from_smirnoff(
+                force_field=sage_unconstrained,
+                topology=[molecule],
+            )
+            interchange.positions = molecule.conformers[0]
+
+            return interchange
+
+        thf_interchange = make_interchange(thf)
+        ace_interchange = make_interchange(ace)
+        complex_interchange = thf_interchange.combine(ace_interchange)
+
+        thf_vdw = thf_interchange["vdW"].get_system_parameters()
+        ace_vdw = ace_interchange["vdW"].get_system_parameters()
+        add_vdw = complex_interchange["vdW"].get_system_parameters()
+
+        numpy.testing.assert_equal(numpy.vstack([thf_vdw, ace_vdw]), add_vdw)
+
+        # TODO: Ensure the de-duplication is maintained after exports
+
+    def test_positions_setting(self, monkeypatch, sage):
+        """Test that positions exist on the result if and only if
+        both input objects have positions."""
+
+        monkeypatch.setenv("INTERCHANGE_EXPERIMENTAL", "1")
+
+        ethane = Molecule.from_smiles("CC")
+        methane = Molecule.from_smiles("C")
+
+        ethane_interchange = Interchange.from_smirnoff(
+            sage,
+            [ethane],
+        )
+        methane_interchange = Interchange.from_smirnoff(sage, [methane])
+
+        ethane.generate_conformers(n_conformers=1)
+        methane.generate_conformers(n_conformers=1)
+
+        assert (methane_interchange.combine(ethane_interchange)).positions is None
+        methane_interchange.positions = methane.conformers[0]
+        assert (methane_interchange.combine(ethane_interchange)).positions is None
+        ethane_interchange.positions = ethane.conformers[0]
+        assert (methane_interchange.combine(ethane_interchange)).positions is not None
+
+    @pytest.mark.parametrize("handler", ["Electrostatics", "vdW"])
+    def test_error_mismatched_cutoffs(
+        self,
+        monkeypatch,
+        sage,
+        basic_top,
+        handler,
+    ):
+
+        monkeypatch.setenv("INTERCHANGE_EXPERIMENTAL", "1")
+
+        sage_modified = ForceField("openff-2.1.0.offxml")
+
+        sage_modified[handler].cutoff *= 1.5
+
+        with pytest.raises(
+            CutoffMismatchError,
+            match=f"{handler} cutoffs do not match",
+        ):
+            sage.create_interchange(basic_top).combine(
+                sage_modified.create_interchange(basic_top),
+            )
+
+    def test_error_mismatched_switching_function(
+        self,
+        monkeypatch,
+        sage,
+        basic_top,
+    ):
+
+        monkeypatch.setenv("INTERCHANGE_EXPERIMENTAL", "1")
+
+        sage_modified = ForceField("openff-2.1.0.offxml")
+
+        sage_modified["vdW"].switch_width *= 0.0
+
+        with pytest.raises(
+            SwitchingFunctionMismatchError,
+        ):
+            sage.create_interchange(basic_top).combine(
+                sage_modified.create_interchange(basic_top),
+            )

openff/interchange/components/interchange.py

@@ -29,7 +29,6 @@
     InvalidTopologyError,
     MissingParameterHandlerError,
     MissingPositionsError,
-    UnsupportedCombinationError,
     UnsupportedExportError,
 )
 from openff.interchange.operations.minimize import (
@@ -946,80 +945,9 @@ def __add__(self, other: "Interchange") -> "Interchange":
     @experimental
     def combine(self, other: "Interchange") -> "Interchange":
         """Combine two Interchange objects. This method is unstable and not yet unsafe."""
-        from openff.interchange.components.toolkit import _combine_topologies
+        from openff.interchange.operations._combine import _combine
 
-        warnings.warn(
-            "Interchange object combination is experimental and likely to produce "
-            "strange results. Any workflow using this method is not guaranteed to "
-            "be suitable for production. Use with extreme caution and thoroughly "
-            "validate results!",
-            stacklevel=2,
-        )
-
-        self_copy = copy.deepcopy(self)
-
-        self_copy.topology = _combine_topologies(self.topology, other.topology)
-        atom_offset = self.topology.n_atoms
-
-        if "Electrostatics" in self_copy.collections:
-            self_copy["Electrostatics"]._charges = dict()
-            self_copy["Electrostatics"]._charges_cached = False
-
-        if "Electrostatics" in other.collections:
-            other["Electrostatics"]._charges = dict()
-            other["Electrostatics"]._charges_cached = False
-
-        for handler_name, handler in other.collections.items():
-            # TODO: Actually specify behavior in this case
-            try:
-                self_handler = self_copy.collections[handler_name]
-            except KeyError:
-                self_copy.collections[handler_name] = handler
-                warnings.warn(
-                    f"'other' Interchange object has handler with name {handler_name} not "
-                    f"found in 'self,' but it has now been added.",
-                    stacklevel=2,
-                )
-                continue
-
-            for top_key, pot_key in handler.key_map.items():
-                new_atom_indices = tuple(
-                    idx + atom_offset for idx in top_key.atom_indices
-                )
-                new_top_key = top_key.__class__(**top_key.dict())
-                try:
-                    new_top_key.atom_indices = new_atom_indices
-                except ValueError:
-                    assert len(new_atom_indices) == 1
-                    new_top_key.this_atom_index = new_atom_indices[0]
-
-                self_handler.key_map.update({new_top_key: pot_key})
-                if handler_name == "Constraints":
-                    self_handler.potentials.update(
-                        {pot_key: handler.potentials[pot_key]},
-                    )
-                else:
-                    self_handler.potentials.update(
-                        {pot_key: handler.potentials[pot_key]},
-                    )
-
-            self_copy.collections[handler_name] = self_handler
-
-        if self_copy.positions is not None and other.positions is not None:
-            new_positions = np.vstack([self_copy.positions, other.positions])
-            self_copy.positions = new_positions
-        else:
-            warnings.warn(
-                "Setting positions to None because one or both objects added together were missing positions.",
-            )
-            self_copy.positions = None
-
-        if not np.all(self_copy.box == other.box):
-            raise UnsupportedCombinationError(
-                "Combination with unequal box vectors is not curretnly supported",
-            )
-
-        return self_copy
+        return _combine(self, other)
 
     def __repr__(self) -> str:
         periodic = self.box is not None