Speedup for molecule isomorphism checking (#2114)

* potentially quick speedup for molecule isomorphism

* silence yappy test and handle nx graphs

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* fix atom counting with negative numbered Hs

* change test to expect new (equivalent) ordering for atom mapping

* add new error, test, and update releasenotes

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* fix greater-than/less-than mismatch

* privatize Molecule.are_isomorphic.to_networks, unconditionally deepcopy input, polish docstrings

* Update releasenotes for 0.18.0 release

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Author: j-wags

docs/releasehistory.md

@@ -7,20 +7,19 @@ Releases follow the `major.minor.micro` scheme recommended by [PEP440](https://w
 * `micro` increments represent bugfix releases or improvements in documentation
 
 
-## Current development
+## 0.18.0
 
 ### API-breaking changes
+- [PR #2114](https://github.com/openforcefield/openff-toolkit/pull/2114): Removes `Molecule.are_isomorphic.to_networkx` from the public API (by renaming to `_to_networkx`). We are nearly certain that nobody is using this, and it is DIFFERENT from the much more useful `Molecule.to_networkx`. This change is motivated by the observation that further behavior changes to this method can relieve critical performance bottlenecks, and that we shouldn't repeatedly change the behavior of public methods.  
 
 ### Behavior changes
+- [PR #2114](https://github.com/openforcefield/openff-toolkit/pull/2114): Fixes [Issue #2035](https://github.com/openforcefield/openff-toolkit/issues/2035) (and touches on several others) by improving the runtime of Molecule.are_isomorphic/is_isomorphic_with, which can greatly improve interchange creation and export runtime.  
+- [PR #2114](https://github.com/openforcefield/openff-toolkit/pull/2114): May change the specific mapping of symmetric atoms returned by `Molecule.are_isomorphic` (to a different but chemically equivalent mapping) 
 
 ### Bugfixes
 
 - [PR #2115](https://github.com/openforcefield/openff-toolkit/pull/2115): Avoid unnecessary conformer coordinate lookups which were causing performance issues with large molecules. Fixes [Issue #1855](https://github.com/openforcefield/openff-toolkit/issues/1844).
 
-### New features
-
-### Improved documentation and warnings
-
 
 ## 0.17.1
 

openff/toolkit/_tests/test_molecule.py

@@ -59,6 +59,7 @@
 from openff.toolkit.utils import get_data_file_path
 from openff.toolkit.utils.exceptions import (
     AtomMappingWarning,
+    BadMoleculeAssumptionError,
     HierarchyIteratorNameConflictError,
     IncompatibleShapeError,
     IncompatibleTypeError,
@@ -1521,7 +1522,7 @@ def test_isomorphic_general(self):
         ethanol_reverse = create_reversed_ethanol()
         assert ethanol.is_isomorphic_with(ethanol_reverse) is True
         # check a reference mapping between ethanol and ethanol_reverse matches that calculated
-        ref_mapping = {0: 8, 1: 7, 2: 6, 3: 3, 4: 4, 5: 5, 6: 1, 7: 2, 8: 0}
+        ref_mapping = {0: 8, 1: 7, 2: 6, 3: 5, 4: 4, 5: 3, 6: 2, 7: 1, 8: 0}
         assert Molecule.are_isomorphic(ethanol, ethanol_reverse, return_atom_map=True)[1] == ref_mapping
         # check matching with nx.Graph atomic numbers and connectivity only
         assert (
@@ -1669,6 +1670,16 @@ def test_isomorphic_perumtations(self, inputs):
             is inputs["result"]
         )
 
+    def test_isomorphic_speedup_assumptions(self):
+        """
+        This test ensures that the assumptions needed by the isomorphism checking speedup in
+        https://github.com/openforcefield/openff-toolkit/pull/2114 are policed
+        """
+        mol1 = Molecule.from_mapped_smiles("[F:1][H-:2][H-:3][Cl:4]")
+        mol2 = Molecule.from_mapped_smiles("[Cl:1][H-:2][H-:3][F:4]")
+        with pytest.raises(BadMoleculeAssumptionError):
+            assert mol1.is_isomorphic_with(mol2)
+
     @requires_openeye
     def test_strip_atom_stereochemistry(self):
         """Test the basic behavior of strip_atom_stereochemistry"""

openff/toolkit/topology/molecule.py

@@ -55,6 +55,7 @@
 from openff.toolkit.utils.constants import DEFAULT_AROMATICITY_MODEL
 from openff.toolkit.utils.exceptions import (
     AtomMappingWarning,
+    BadMoleculeAssumptionError,
     BondExistsError,
     HierarchyIteratorNameConflictError,
     HierarchySchemeNotFoundException,
@@ -2131,20 +2132,44 @@ def edge_match_func(x, y):
             edge_match_func = None  # type: ignore
 
         # Here we should work out what data type we have, also deal with lists?
-        def to_networkx(data: FrozenMolecule | nx.Graph) -> nx.Graph:
-            """For the given data type, return the networkx graph"""
+        def _to_networkx(data: FrozenMolecule | nx.Graph) -> nx.Graph:
+            """
+            For the given data type, return the networkx graph. Note that the graphs returned from this method
+            will have negative-numbered hydrogens (eg. the first hydrogen on a heavy atom will have atomic number
+            -1, the second -2, etc). This is a performance hack to reduce the number of potential symmetric matches
+            that must be evaluated during isomorphism checks.
+            """
             if strip_pyrimidal_n_atom_stereo:
                 SMARTS = "[N+0X3:1](-[*])(-[*])(-[*])"
 
+            data = deepcopy(data)
             if isinstance(data, FrozenMolecule):
                 # Molecule class instance
                 if strip_pyrimidal_n_atom_stereo:
                     # Make a copy of the molecule so we don't modify the original
-                    data = deepcopy(data)
                     data.strip_atom_stereochemistry(SMARTS, toolkit_registry=toolkit_registry)
+                for atom in data.atoms:
+                    h_counter = -1
+                    for neighbor in atom.bonded_atoms:
+                        if neighbor.atomic_number < 0:
+                            raise BadMoleculeAssumptionError(
+                                f"Molecule {data} appears to violate an assumption of the OpenFF Toolkit "
+                                f"(likely that an H can't have two bonds). Please check that your molecule is "
+                                f"what you think it is, and if it is indeed what you intend, open an issue at "
+                                f"https://github.com/openforcefield/openff-toolkit/issues"
+                            )
+                        if neighbor.atomic_number == 1:
+                            neighbor._atomic_number = h_counter
+                            h_counter -= 1
                 return data.to_networkx()
 
             elif isinstance(data, nx.Graph):
+                for node in data:
+                    h_counter = -1
+                    for neighbor in data.neighbors(node):
+                        if data.nodes[neighbor]["atomic_number"] == 1:
+                            data.nodes[neighbor]["atomic_number"] = h_counter
+                            h_counter -= 1
                 return data
 
             else:
@@ -2154,8 +2179,8 @@ def to_networkx(data: FrozenMolecule | nx.Graph) -> nx.Graph:
                     f"or networkx.Graph representation of the molecule."
                 )
 
-        mol1_netx = to_networkx(mol1)
-        mol2_netx = to_networkx(mol2)
+        mol1_netx = _to_networkx(mol1)
+        mol2_netx = _to_networkx(mol2)
 
         from networkx.algorithms.isomorphism import GraphMatcher
 
@@ -5665,6 +5690,11 @@ def _atom_nums_to_hill_formula(atom_nums: list[int]) -> str:
     SYMBOLS_ = deepcopy(SYMBOLS)
     SYMBOLS_[0] = "X"
 
+    atom_nums = deepcopy(atom_nums)
+    for idx in range(len(atom_nums)):
+        if atom_nums[idx] < 0:
+            atom_nums[idx] = 1
+
     atom_symbol_counts = Counter(SYMBOLS_[atom_num] for atom_num in atom_nums)
 
     formula = []

openff/toolkit/utils/__init__.py

@@ -15,6 +15,7 @@
     RDKIT_AVAILABLE,
     AmberToolsToolkitWrapper,
     AntechamberNotFoundError,
+    BadMoleculeAssumptionError,
     BuiltInToolkitWrapper,
     ChargeCalculationError,
     ChargeMethodUnavailableError,

openff/toolkit/utils/exceptions.py

@@ -760,3 +760,7 @@ def __init__(self, res_name, mol_bond, query_bond, reason):
         msg += f"\t{reason}"
         super().__init__(msg)
         self.msg = msg
+
+
+class BadMoleculeAssumptionError(OpenFFToolkitException):
+    """Exception raised when a molecule violates assumptions made by the OpenFF Toolkit's implementation"""

openff/toolkit/utils/toolkits.py

@@ -34,6 +34,7 @@
     "RDKIT_AVAILABLE",
     "AmberToolsToolkitWrapper",
     "AntechamberNotFoundError",
+    "BadMoleculeAssumptionError",
     "BuiltInToolkitWrapper",
     "ChargeCalculationError",
     "ChargeMethodUnavailableError",
@@ -72,6 +73,7 @@
 )
 from openff.toolkit.utils.exceptions import (
     AntechamberNotFoundError,
+    BadMoleculeAssumptionError,
     ChargeCalculationError,
     ChargeMethodUnavailableError,
     GAFFAtomTypeWarning,