From b049629e0780dd2a72fd2e6a2b95cce7901c679c Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 30 Jul 2026 19:22:30 +0100 Subject: [PATCH 1/8] Warn when a mapping stops short at a pairable attachment point. --- src/BioSimSpace/Align/_align.py | 157 +++++++++++++++++- .../Sandpit/Exscientia/Align/_align.py | 150 ++++++++++++++++- tests/Align/test_align.py | 13 ++ tests/Sandpit/Exscientia/Align/test_align.py | 13 ++ 4 files changed, 331 insertions(+), 2 deletions(-) diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 666cfe44..2706fcb8 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -896,7 +896,7 @@ def matchAtoms( """ if roi is None: - return _matchAtoms( + result = _matchAtoms( molecule0=molecule0, molecule1=molecule1, scoring_function=scoring_function, @@ -913,6 +913,58 @@ def matchAtoms( property_map1=property_map1, mcs_kwargs=mcs_kwargs, ) + + # Only check a mapping generated from the defaults. If the user has + # configured the MCS then both the baseline and our idea of a sensible + # mapping may not match their intent. This also stops the retry below + # from recursing, since it passes 'mcs_kwargs'. Skip when a prematch is + # given, since the retry could then fall back on Sire MCS, which + # ignores 'mcs_kwargs', making the comparison meaningless. + if not mcs_kwargs and not prematch: + best = result[0] if return_scores else result + if isinstance(best, list): + best = best[0] if best else {} + + # Attachment points where the MCS stopped on both sides. + flagged = ( + _flag_unmapped_attachments(molecule0, molecule1, best) if best else [] + ) + + if flagged: + # Retry with ring matching relaxed to see if it does better. + retry = matchAtoms( + molecule0, + molecule1, + scoring_function=scoring_function, + prematch=prematch, + timeout=timeout, + complete_rings_only=complete_rings_only, + max_scoring_matches=max_scoring_matches, + prune_perturbed_constraints=prune_perturbed_constraints, + prune_crossing_constraints=prune_crossing_constraints, + prune_atom_types=prune_atom_types, + property_map0=property_map0, + property_map1=property_map1, + mcs_kwargs={"ringMatchesRingOnly": False}, + ) + + # Only trust the retry if it extends the mapping, i.e. keeps + # every existing pair and adds sensible ones. + if ( + len(retry) > len(best) + and set(best.items()) <= set(retry.items()) + and _is_sensible_extension(molecule0, molecule1, best, retry) + ): + _warnings.warn( + f"Mapping leaves heavy atoms unmapped on both sides of " + f"atom(s) {flagged} " + f"in molecule0. Relaxing 'ringMatchesRingOnly' gives a " + f"common core of {len(retry)} rather than {len(best)}. " + f"Consider passing " + f"mcs_kwargs={{'ringMatchesRingOnly': False}}." + ) + + return result else: return _roiMatch( molecule0, @@ -927,6 +979,109 @@ def matchAtoms( ) +def _flag_unmapped_attachments(molecule0, molecule1, mapping): + """ + Internal function to find mapped atoms that have an unmapped heavy atom + neighbour of a common element in both molecules. These are attachment + points where the MCS stopped on both sides, which usually means that a + pairable atom was missed. + + Parameters + ---------- + + molecule0 : :class:`Molecule ` + The first molecule. + + molecule1 : :class:`Molecule ` + The second molecule. + + mapping : dict + The atom mapping between the two molecules. + + Returns + ------- + + flagged : [int] + The indices of the flagged atoms in molecule0. + """ + from sire.legacy import Mol as _SireMol + + mol0 = molecule0._getSireObject() + mol1 = molecule1._getSireObject() + conn0 = mol0.property("connectivity") + conn1 = mol1.property("connectivity") + + def _heavy_elements(mol, conn, idx, mapped): + elements = set() + for i in conn.connections_to(_SireMol.AtomIdx(idx)): + if i.value() not in mapped: + protons = mol.atom(i).property("element").num_protons() + if protons > 1: + elements.add(protons) + return elements + + mapped1 = set(mapping.values()) + flagged = [] + + for idx0, idx1 in mapping.items(): + elements0 = _heavy_elements(mol0, conn0, idx0, mapping) + if not elements0: + continue + elements1 = _heavy_elements(mol1, conn1, idx1, mapped1) + if elements0 & elements1: + flagged.append(idx0) + + return flagged + + +def _is_sensible_extension(molecule0, molecule1, mapping, extended): + """ + Internal function to test whether the atoms that 'extended' adds relative + to 'mapping' pair like with like. The MCS uses CompareAny, so it is free to + pair a heavy atom with a hydrogen, which grows the common core without + improving the mapping. + + Parameters + ---------- + + molecule0 : :class:`Molecule ` + The first molecule. + + molecule1 : :class:`Molecule ` + The second molecule. + + mapping : dict + The original atom mapping. + + extended : dict + The larger atom mapping to test. + + Returns + ------- + + is_sensible : bool + Whether the added atoms pair heavy with heavy and hydrogen with + hydrogen. + """ + from sire.legacy import Mol as _SireMol + + mol0 = molecule0._getSireObject() + mol1 = molecule1._getSireObject() + + for idx0, idx1 in extended.items(): + if idx0 not in mapping: + protons0 = ( + mol0.atom(_SireMol.AtomIdx(idx0)).property("element").num_protons() + ) + protons1 = ( + mol1.atom(_SireMol.AtomIdx(idx1)).property("element").num_protons() + ) + if (protons0 > 1) != (protons1 > 1): + return False + + return True + + def _matchAtoms( molecule0, molecule1, diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index 6ff60101..ce8a05ab 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py @@ -703,6 +703,109 @@ def generateNetwork( return edges, scores +def _flag_unmapped_attachments(molecule0, molecule1, mapping): + """ + Internal function to find mapped atoms that have an unmapped heavy atom + neighbour of a common element in both molecules. These are attachment + points where the MCS stopped on both sides, which usually means that a + pairable atom was missed. + + Parameters + ---------- + + molecule0 : :class:`Molecule ` + The first molecule. + + molecule1 : :class:`Molecule ` + The second molecule. + + mapping : dict + The atom mapping between the two molecules. + + Returns + ------- + + flagged : [int] + The indices of the flagged atoms in molecule0. + """ + from sire.legacy import Mol as _SireMol + + mol0 = molecule0._getSireObject() + mol1 = molecule1._getSireObject() + conn0 = mol0.property("connectivity") + conn1 = mol1.property("connectivity") + + def _heavy_elements(mol, conn, idx, mapped): + elements = set() + for i in conn.connections_to(_SireMol.AtomIdx(idx)): + if i.value() not in mapped: + protons = mol.atom(i).property("element").num_protons() + if protons > 1: + elements.add(protons) + return elements + + mapped1 = set(mapping.values()) + flagged = [] + + for idx0, idx1 in mapping.items(): + elements0 = _heavy_elements(mol0, conn0, idx0, mapping) + if not elements0: + continue + elements1 = _heavy_elements(mol1, conn1, idx1, mapped1) + if elements0 & elements1: + flagged.append(idx0) + + return flagged + + +def _is_sensible_extension(molecule0, molecule1, mapping, extended): + """ + Internal function to test whether the atoms that 'extended' adds relative + to 'mapping' pair like with like. The MCS uses CompareAny, so it is free to + pair a heavy atom with a hydrogen, which grows the common core without + improving the mapping. + + Parameters + ---------- + + molecule0 : :class:`Molecule ` + The first molecule. + + molecule1 : :class:`Molecule ` + The second molecule. + + mapping : dict + The original atom mapping. + + extended : dict + The larger atom mapping to test. + + Returns + ------- + + is_sensible : bool + Whether the added atoms pair heavy with heavy and hydrogen with + hydrogen. + """ + from sire.legacy import Mol as _SireMol + + mol0 = molecule0._getSireObject() + mol1 = molecule1._getSireObject() + + for idx0, idx1 in extended.items(): + if idx0 not in mapping: + protons0 = ( + mol0.atom(_SireMol.AtomIdx(idx0)).property("element").num_protons() + ) + protons1 = ( + mol1.atom(_SireMol.AtomIdx(idx1)).property("element").num_protons() + ) + if (protons0 > 1) != (protons1 > 1): + return False + + return True + + def defaultMCSOptions(): """ Return the default options used for the RDKit maximum common substructure @@ -948,7 +1051,9 @@ def matchAtoms( mol0 = molecule0._getSireObject() mol1 = molecule1._getSireObject() - # Convert the timeout to seconds and take the value as an integer. + # Convert the timeout to seconds and take the value as an integer, keeping + # the original for any onward call that expects a Time. + orig_timeout = timeout timeout = int(timeout.seconds().value()) # Use RDKkit to find the maximum common substructure. @@ -1076,6 +1181,49 @@ def matchAtoms( _prune_crossing_constraints(molecule0, molecule1, x) for x in mappings ] + # Warn if the mapping stopped short at an attachment point where a pairable + # atom exists. Only check a mapping generated from the defaults. If the + # user has configured the MCS then both the baseline and our idea of a + # sensible mapping may not match their intent. This also stops the retry + # below from recursing, since it passes 'mcs_kwargs'. Skip when a prematch + # is given, since the retry could then fall back on Sire MCS, which ignores + # 'mcs_kwargs', making the comparison meaningless. + if not mcs_kwargs and not prematch and mappings: + flagged = _flag_unmapped_attachments(molecule0, molecule1, mappings[0]) + + if flagged: + # Retry with ring matching relaxed to see if it does better. + retry = matchAtoms( + molecule0, + molecule1, + engine=engine, + scoring_function=scoring_function, + prematch=prematch, + timeout=orig_timeout, + complete_rings_only=complete_rings_only, + prune_perturbed_constraints=prune_perturbed_constraints, + prune_crossing_constraints=prune_crossing_constraints, + max_scoring_matches=max_scoring_matches, + property_map0=property_map0, + property_map1=property_map1, + mcs_kwargs={"ringMatchesRingOnly": False}, + ) + + # Only trust the retry if it extends the mapping, i.e. keeps every + # existing pair and adds sensible ones. + if ( + len(retry) > len(mappings[0]) + and set(mappings[0].items()) <= set(retry.items()) + and _is_sensible_extension(molecule0, molecule1, mappings[0], retry) + ): + _warnings.warn( + f"Mapping leaves heavy atoms unmapped on both sides of " + f"atom(s) {flagged} in molecule0. Relaxing " + f"'ringMatchesRingOnly' gives a common core of " + f"{len(retry)} rather than {len(mappings[0])}. Consider " + f"passing mcs_kwargs={{'ringMatchesRingOnly': False}}." + ) + if matches == 1: if return_scores: return (mappings[0], scores[0]) diff --git a/tests/Align/test_align.py b/tests/Align/test_align.py index ef7fc280..aa2ba2af 100644 --- a/tests/Align/test_align.py +++ b/tests/Align/test_align.py @@ -1,4 +1,5 @@ import sys +import warnings import pytest import sire as sr @@ -1397,3 +1398,15 @@ def test_mcs_kwargs_merge(ejm31, jmc28): # No ring is broken or made, so neither property is set. assert not sire_mol.has_property("ring_breaking_bonds") assert not sire_mol.has_property("ring_making_bonds") + + +def test_unmapped_attachment_warning(ejm31, jmc28): + # The default mapping stops at the carbonyl carbon (atom 17), leaving + # heavy atoms unmapped on both sides, so we should be told about it. + with pytest.warns(UserWarning, match="ringMatchesRingOnly"): + BSS.Align.matchAtoms(ejm31, jmc28) + + # No warning once the option has been set explicitly. + with warnings.catch_warnings(): + warnings.simplefilter("error") + BSS.Align.matchAtoms(ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False}) diff --git a/tests/Sandpit/Exscientia/Align/test_align.py b/tests/Sandpit/Exscientia/Align/test_align.py index 6c85a1ee..b6597eba 100644 --- a/tests/Sandpit/Exscientia/Align/test_align.py +++ b/tests/Sandpit/Exscientia/Align/test_align.py @@ -1,4 +1,5 @@ import sys +import warnings import pytest from sire.legacy.Maths import Vector @@ -968,3 +969,15 @@ def test_mcs_kwargs_merge(ejm31, jmc28): # No ring is broken or made, so neither property is set. assert not sire_mol.has_property("ring_breaking_bonds") assert not sire_mol.has_property("ring_making_bonds") + + +def test_unmapped_attachment_warning(ejm31, jmc28): + # The default mapping stops at the carbonyl carbon (atom 17), leaving + # heavy atoms unmapped on both sides, so we should be told about it. + with pytest.warns(UserWarning, match="ringMatchesRingOnly"): + BSS.Align.matchAtoms(ejm31, jmc28) + + # No warning once the option has been set explicitly. + with warnings.catch_warnings(): + warnings.simplefilter("error") + BSS.Align.matchAtoms(ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False}) From 8279dd577c35fe325c12daae21494e91db06ef4e Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 31 Jul 2026 12:34:01 +0100 Subject: [PATCH 2/8] Harden the mapping diagnostic against crashes and custom properties. --- src/BioSimSpace/Align/_align.py | 187 +++++++++++++----- .../Sandpit/Exscientia/Align/_align.py | 166 ++++++++++++---- tests/Align/test_align.py | 16 +- tests/Sandpit/Exscientia/Align/test_align.py | 16 +- 4 files changed, 291 insertions(+), 94 deletions(-) diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 2706fcb8..fc04962c 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -921,49 +921,76 @@ def matchAtoms( # given, since the retry could then fall back on Sire MCS, which # ignores 'mcs_kwargs', making the comparison meaningless. if not mcs_kwargs and not prematch: - best = result[0] if return_scores else result - if isinstance(best, list): - best = best[0] if best else {} - - # Attachment points where the MCS stopped on both sides. - flagged = ( - _flag_unmapped_attachments(molecule0, molecule1, best) if best else [] - ) - - if flagged: - # Retry with ring matching relaxed to see if it does better. - retry = matchAtoms( - molecule0, - molecule1, - scoring_function=scoring_function, - prematch=prematch, - timeout=timeout, - complete_rings_only=complete_rings_only, - max_scoring_matches=max_scoring_matches, - prune_perturbed_constraints=prune_perturbed_constraints, - prune_crossing_constraints=prune_crossing_constraints, - prune_atom_types=prune_atom_types, - property_map0=property_map0, - property_map1=property_map1, - mcs_kwargs={"ringMatchesRingOnly": False}, + # This is a diagnostic, so it must never be able to break a call + # that would otherwise have succeeded. + try: + best = result[0] if return_scores else result + if isinstance(best, list): + best = best[0] if best else {} + + # Attachment points where the MCS stopped on both sides. + flagged = ( + _flag_unmapped_attachments( + molecule0, molecule1, best, property_map0, property_map1 + ) + if best + else [] ) - # Only trust the retry if it extends the mapping, i.e. keeps - # every existing pair and adds sensible ones. - if ( - len(retry) > len(best) - and set(best.items()) <= set(retry.items()) - and _is_sensible_extension(molecule0, molecule1, best, retry) - ): - _warnings.warn( - f"Mapping leaves heavy atoms unmapped on both sides of " - f"atom(s) {flagged} " - f"in molecule0. Relaxing 'ringMatchesRingOnly' gives a " - f"common core of {len(retry)} rather than {len(best)}. " - f"Consider passing " - f"mcs_kwargs={{'ringMatchesRingOnly': False}}." + if flagged: + # Retry with ring matching relaxed to see if it does + # better. Note that the RDKit documentation implies that + # 'completeRingsOnly' forces 'ringMatchesRingOnly', which + # would make this a no-op. It doesn't: as of RDKit + # 2026.03.4 the relaxed search still returns a larger MCS + # with 'completeRingsOnly' enabled. If a future RDKit + # changes this, the feature will silently stop firing. + # + # 'matches' and 'return_scores' are passed explicitly so + # that 'retry' is always a plain dict. The comparison below + # relies on it. 'prematch' is omitted since the enclosing + # guard means it's always empty. + retry = matchAtoms( + molecule0, + molecule1, + scoring_function=scoring_function, + matches=1, + return_scores=False, + timeout=timeout, + complete_rings_only=complete_rings_only, + max_scoring_matches=max_scoring_matches, + prune_perturbed_constraints=prune_perturbed_constraints, + prune_crossing_constraints=prune_crossing_constraints, + prune_atom_types=prune_atom_types, + property_map0=property_map0, + property_map1=property_map1, + mcs_kwargs={"ringMatchesRingOnly": False}, ) + # Only trust the retry if it extends the mapping, i.e. + # keeps every existing pair and adds sensible ones. + if ( + len(retry) > len(best) + and set(best.items()) <= set(retry.items()) + and _is_sensible_extension( + molecule0, + molecule1, + best, + retry, + property_map0, + property_map1, + ) + ): + _warnings.warn( + f"Mapping leaves heavy atoms unmapped on both sides " + f"of atom(s) {_format_flagged(flagged)} in molecule0. " + f"Relaxing 'ringMatchesRingOnly' gives a common core " + f"of {len(retry)} rather than {len(best)}. Consider " + f"passing mcs_kwargs={{'ringMatchesRingOnly': False}}." + ) + except Exception as e: + _warnings.warn(f"Unable to check the quality of the mapping: {e}") + return result else: return _roiMatch( @@ -979,7 +1006,37 @@ def matchAtoms( ) -def _flag_unmapped_attachments(molecule0, molecule1, mapping): +def _format_flagged(flagged, max_show=5): + """ + Internal helper to format a list of atom indices for use in a warning + message, truncating so that the message stays readable for large + molecules. + + Parameters + ---------- + + flagged : [int] + The indices to format. + + max_show : int + The maximum number of indices to show. + + Returns + ------- + + string : str + The formatted indices. + """ + if len(flagged) <= max_show: + return str(flagged) + else: + shown = ", ".join(str(x) for x in flagged[:max_show]) + return f"[{shown}, ... ({len(flagged)} in total)]" + + +def _flag_unmapped_attachments( + molecule0, molecule1, mapping, property_map0={}, property_map1={} +): """ Internal function to find mapped atoms that have an unmapped heavy atom neighbour of a common element in both molecules. These are attachment @@ -998,6 +1055,16 @@ def _flag_unmapped_attachments(molecule0, molecule1, mapping): mapping : dict The atom mapping between the two molecules. + property_map0 : dict + A dictionary that maps "properties" in molecule0 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + property_map1 : dict + A dictionary that maps "properties" in molecule1 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + Returns ------- @@ -1008,33 +1075,42 @@ def _flag_unmapped_attachments(molecule0, molecule1, mapping): mol0 = molecule0._getSireObject() mol1 = molecule1._getSireObject() - conn0 = mol0.property("connectivity") - conn1 = mol1.property("connectivity") - def _heavy_elements(mol, conn, idx, mapped): + # Build the connectivity explicitly, since the molecules aren't guaranteed + # to have a stored "connectivity" property. + conn0 = _SireMol.Connectivity(mol0, _SireMol.CovalentBondHunter()) + conn1 = _SireMol.Connectivity(mol1, _SireMol.CovalentBondHunter()) + + element0 = property_map0.get("element", "element") + element1 = property_map1.get("element", "element") + + def _heavy_elements(mol, conn, idx, mapped, element): elements = set() for i in conn.connections_to(_SireMol.AtomIdx(idx)): if i.value() not in mapped: - protons = mol.atom(i).property("element").num_protons() + protons = mol.atom(i).property(element).num_protons() if protons > 1: elements.add(protons) return elements + mapped0 = set(mapping) mapped1 = set(mapping.values()) flagged = [] for idx0, idx1 in mapping.items(): - elements0 = _heavy_elements(mol0, conn0, idx0, mapping) + elements0 = _heavy_elements(mol0, conn0, idx0, mapped0, element0) if not elements0: continue - elements1 = _heavy_elements(mol1, conn1, idx1, mapped1) + elements1 = _heavy_elements(mol1, conn1, idx1, mapped1, element1) if elements0 & elements1: flagged.append(idx0) return flagged -def _is_sensible_extension(molecule0, molecule1, mapping, extended): +def _is_sensible_extension( + molecule0, molecule1, mapping, extended, property_map0={}, property_map1={} +): """ Internal function to test whether the atoms that 'extended' adds relative to 'mapping' pair like with like. The MCS uses CompareAny, so it is free to @@ -1056,6 +1132,16 @@ def _is_sensible_extension(molecule0, molecule1, mapping, extended): extended : dict The larger atom mapping to test. + property_map0 : dict + A dictionary that maps "properties" in molecule0 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + property_map1 : dict + A dictionary that maps "properties" in molecule1 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + Returns ------- @@ -1068,13 +1154,16 @@ def _is_sensible_extension(molecule0, molecule1, mapping, extended): mol0 = molecule0._getSireObject() mol1 = molecule1._getSireObject() + element0 = property_map0.get("element", "element") + element1 = property_map1.get("element", "element") + for idx0, idx1 in extended.items(): if idx0 not in mapping: protons0 = ( - mol0.atom(_SireMol.AtomIdx(idx0)).property("element").num_protons() + mol0.atom(_SireMol.AtomIdx(idx0)).property(element0).num_protons() ) protons1 = ( - mol1.atom(_SireMol.AtomIdx(idx1)).property("element").num_protons() + mol1.atom(_SireMol.AtomIdx(idx1)).property(element1).num_protons() ) if (protons0 > 1) != (protons1 > 1): return False diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index ce8a05ab..b6093733 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py @@ -703,7 +703,37 @@ def generateNetwork( return edges, scores -def _flag_unmapped_attachments(molecule0, molecule1, mapping): +def _format_flagged(flagged, max_show=5): + """ + Internal helper to format a list of atom indices for use in a warning + message, truncating so that the message stays readable for large + molecules. + + Parameters + ---------- + + flagged : [int] + The indices to format. + + max_show : int + The maximum number of indices to show. + + Returns + ------- + + string : str + The formatted indices. + """ + if len(flagged) <= max_show: + return str(flagged) + else: + shown = ", ".join(str(x) for x in flagged[:max_show]) + return f"[{shown}, ... ({len(flagged)} in total)]" + + +def _flag_unmapped_attachments( + molecule0, molecule1, mapping, property_map0={}, property_map1={} +): """ Internal function to find mapped atoms that have an unmapped heavy atom neighbour of a common element in both molecules. These are attachment @@ -722,6 +752,16 @@ def _flag_unmapped_attachments(molecule0, molecule1, mapping): mapping : dict The atom mapping between the two molecules. + property_map0 : dict + A dictionary that maps "properties" in molecule0 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + property_map1 : dict + A dictionary that maps "properties" in molecule1 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + Returns ------- @@ -732,33 +772,42 @@ def _flag_unmapped_attachments(molecule0, molecule1, mapping): mol0 = molecule0._getSireObject() mol1 = molecule1._getSireObject() - conn0 = mol0.property("connectivity") - conn1 = mol1.property("connectivity") - def _heavy_elements(mol, conn, idx, mapped): + # Build the connectivity explicitly, since the molecules aren't guaranteed + # to have a stored "connectivity" property. + conn0 = _SireMol.Connectivity(mol0, _SireMol.CovalentBondHunter()) + conn1 = _SireMol.Connectivity(mol1, _SireMol.CovalentBondHunter()) + + element0 = property_map0.get("element", "element") + element1 = property_map1.get("element", "element") + + def _heavy_elements(mol, conn, idx, mapped, element): elements = set() for i in conn.connections_to(_SireMol.AtomIdx(idx)): if i.value() not in mapped: - protons = mol.atom(i).property("element").num_protons() + protons = mol.atom(i).property(element).num_protons() if protons > 1: elements.add(protons) return elements + mapped0 = set(mapping) mapped1 = set(mapping.values()) flagged = [] for idx0, idx1 in mapping.items(): - elements0 = _heavy_elements(mol0, conn0, idx0, mapping) + elements0 = _heavy_elements(mol0, conn0, idx0, mapped0, element0) if not elements0: continue - elements1 = _heavy_elements(mol1, conn1, idx1, mapped1) + elements1 = _heavy_elements(mol1, conn1, idx1, mapped1, element1) if elements0 & elements1: flagged.append(idx0) return flagged -def _is_sensible_extension(molecule0, molecule1, mapping, extended): +def _is_sensible_extension( + molecule0, molecule1, mapping, extended, property_map0={}, property_map1={} +): """ Internal function to test whether the atoms that 'extended' adds relative to 'mapping' pair like with like. The MCS uses CompareAny, so it is free to @@ -780,6 +829,16 @@ def _is_sensible_extension(molecule0, molecule1, mapping, extended): extended : dict The larger atom mapping to test. + property_map0 : dict + A dictionary that maps "properties" in molecule0 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + property_map1 : dict + A dictionary that maps "properties" in molecule1 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + Returns ------- @@ -792,13 +851,16 @@ def _is_sensible_extension(molecule0, molecule1, mapping, extended): mol0 = molecule0._getSireObject() mol1 = molecule1._getSireObject() + element0 = property_map0.get("element", "element") + element1 = property_map1.get("element", "element") + for idx0, idx1 in extended.items(): if idx0 not in mapping: protons0 = ( - mol0.atom(_SireMol.AtomIdx(idx0)).property("element").num_protons() + mol0.atom(_SireMol.AtomIdx(idx0)).property(element0).num_protons() ) protons1 = ( - mol1.atom(_SireMol.AtomIdx(idx1)).property("element").num_protons() + mol1.atom(_SireMol.AtomIdx(idx1)).property(element1).num_protons() ) if (protons0 > 1) != (protons1 > 1): return False @@ -1189,41 +1251,63 @@ def matchAtoms( # is given, since the retry could then fall back on Sire MCS, which ignores # 'mcs_kwargs', making the comparison meaningless. if not mcs_kwargs and not prematch and mappings: - flagged = _flag_unmapped_attachments(molecule0, molecule1, mappings[0]) - - if flagged: - # Retry with ring matching relaxed to see if it does better. - retry = matchAtoms( - molecule0, - molecule1, - engine=engine, - scoring_function=scoring_function, - prematch=prematch, - timeout=orig_timeout, - complete_rings_only=complete_rings_only, - prune_perturbed_constraints=prune_perturbed_constraints, - prune_crossing_constraints=prune_crossing_constraints, - max_scoring_matches=max_scoring_matches, - property_map0=property_map0, - property_map1=property_map1, - mcs_kwargs={"ringMatchesRingOnly": False}, + # This is a diagnostic, so it must never be able to break a call that + # would otherwise have succeeded. + try: + best = mappings[0] + flagged = _flag_unmapped_attachments( + molecule0, molecule1, best, property_map0, property_map1 ) - # Only trust the retry if it extends the mapping, i.e. keeps every - # existing pair and adds sensible ones. - if ( - len(retry) > len(mappings[0]) - and set(mappings[0].items()) <= set(retry.items()) - and _is_sensible_extension(molecule0, molecule1, mappings[0], retry) - ): - _warnings.warn( - f"Mapping leaves heavy atoms unmapped on both sides of " - f"atom(s) {flagged} in molecule0. Relaxing " - f"'ringMatchesRingOnly' gives a common core of " - f"{len(retry)} rather than {len(mappings[0])}. Consider " - f"passing mcs_kwargs={{'ringMatchesRingOnly': False}}." + if flagged: + # Retry with ring matching relaxed to see if it does better. + # Note that the RDKit documentation implies that + # 'completeRingsOnly' forces 'ringMatchesRingOnly', which would + # make this a no-op. It doesn't: as of RDKit 2026.03.4 the + # relaxed search still returns a larger MCS with + # 'completeRingsOnly' enabled. If a future RDKit changes this, + # the feature will silently stop firing. + # + # 'matches' and 'return_scores' are passed explicitly so that + # 'retry' is always a plain dict. The comparison below relies + # on it. 'prematch' is omitted since the enclosing guard means + # it's always empty. + retry = matchAtoms( + molecule0, + molecule1, + engine=engine, + scoring_function=scoring_function, + matches=1, + return_scores=False, + timeout=orig_timeout, + complete_rings_only=complete_rings_only, + prune_perturbed_constraints=prune_perturbed_constraints, + prune_crossing_constraints=prune_crossing_constraints, + max_scoring_matches=max_scoring_matches, + property_map0=property_map0, + property_map1=property_map1, + mcs_kwargs={"ringMatchesRingOnly": False}, ) + # Only trust the retry if it extends the mapping, i.e. keeps + # every existing pair and adds sensible ones. + if ( + len(retry) > len(best) + and set(best.items()) <= set(retry.items()) + and _is_sensible_extension( + molecule0, molecule1, best, retry, property_map0, property_map1 + ) + ): + _warnings.warn( + f"Mapping leaves heavy atoms unmapped on both sides " + f"of atom(s) {_format_flagged(flagged)} in molecule0. " + f"Relaxing 'ringMatchesRingOnly' gives a common core " + f"of {len(retry)} rather than {len(best)}. Consider " + f"passing mcs_kwargs={{'ringMatchesRingOnly': False}}." + ) + except Exception as e: + _warnings.warn(f"Unable to check the quality of the mapping: {e}") + if matches == 1: if return_scores: return (mappings[0], scores[0]) diff --git a/tests/Align/test_align.py b/tests/Align/test_align.py index aa2ba2af..bcc62091 100644 --- a/tests/Align/test_align.py +++ b/tests/Align/test_align.py @@ -1406,7 +1406,19 @@ def test_unmapped_attachment_warning(ejm31, jmc28): with pytest.warns(UserWarning, match="ringMatchesRingOnly"): BSS.Align.matchAtoms(ejm31, jmc28) - # No warning once the option has been set explicitly. + # No warning once the option has been set explicitly. Only promote the + # warning we care about, so that unrelated warnings from RDKit or Sire + # don't fail the test. with warnings.catch_warnings(): - warnings.simplefilter("error") + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") BSS.Align.matchAtoms(ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False}) + + # The returned mapping must be unchanged by the check, since it is only + # meant to be an observation. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + checked = BSS.Align.matchAtoms(ejm31, jmc28) + unchecked = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() + ) + assert checked == unchecked diff --git a/tests/Sandpit/Exscientia/Align/test_align.py b/tests/Sandpit/Exscientia/Align/test_align.py index b6597eba..9b434931 100644 --- a/tests/Sandpit/Exscientia/Align/test_align.py +++ b/tests/Sandpit/Exscientia/Align/test_align.py @@ -977,7 +977,19 @@ def test_unmapped_attachment_warning(ejm31, jmc28): with pytest.warns(UserWarning, match="ringMatchesRingOnly"): BSS.Align.matchAtoms(ejm31, jmc28) - # No warning once the option has been set explicitly. + # No warning once the option has been set explicitly. Only promote the + # warning we care about, so that unrelated warnings from RDKit or Sire + # don't fail the test. with warnings.catch_warnings(): - warnings.simplefilter("error") + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") BSS.Align.matchAtoms(ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False}) + + # The returned mapping must be unchanged by the check, since it is only + # meant to be an observation. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + checked = BSS.Align.matchAtoms(ejm31, jmc28) + unchecked = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() + ) + assert checked == unchecked From 30895c10cacff6fde734751b29bb6ee52531fd1c Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 31 Jul 2026 12:53:14 +0100 Subject: [PATCH 3/8] Check the mapping before pruning, not after. --- src/BioSimSpace/Align/_align.py | 157 +++++++++--------- .../Sandpit/Exscientia/Align/_align.py | 50 +++--- 2 files changed, 104 insertions(+), 103 deletions(-) diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index fc04962c..40e20c86 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -896,7 +896,7 @@ def matchAtoms( """ if roi is None: - result = _matchAtoms( + return _matchAtoms( molecule0=molecule0, molecule1=molecule1, scoring_function=scoring_function, @@ -913,85 +913,6 @@ def matchAtoms( property_map1=property_map1, mcs_kwargs=mcs_kwargs, ) - - # Only check a mapping generated from the defaults. If the user has - # configured the MCS then both the baseline and our idea of a sensible - # mapping may not match their intent. This also stops the retry below - # from recursing, since it passes 'mcs_kwargs'. Skip when a prematch is - # given, since the retry could then fall back on Sire MCS, which - # ignores 'mcs_kwargs', making the comparison meaningless. - if not mcs_kwargs and not prematch: - # This is a diagnostic, so it must never be able to break a call - # that would otherwise have succeeded. - try: - best = result[0] if return_scores else result - if isinstance(best, list): - best = best[0] if best else {} - - # Attachment points where the MCS stopped on both sides. - flagged = ( - _flag_unmapped_attachments( - molecule0, molecule1, best, property_map0, property_map1 - ) - if best - else [] - ) - - if flagged: - # Retry with ring matching relaxed to see if it does - # better. Note that the RDKit documentation implies that - # 'completeRingsOnly' forces 'ringMatchesRingOnly', which - # would make this a no-op. It doesn't: as of RDKit - # 2026.03.4 the relaxed search still returns a larger MCS - # with 'completeRingsOnly' enabled. If a future RDKit - # changes this, the feature will silently stop firing. - # - # 'matches' and 'return_scores' are passed explicitly so - # that 'retry' is always a plain dict. The comparison below - # relies on it. 'prematch' is omitted since the enclosing - # guard means it's always empty. - retry = matchAtoms( - molecule0, - molecule1, - scoring_function=scoring_function, - matches=1, - return_scores=False, - timeout=timeout, - complete_rings_only=complete_rings_only, - max_scoring_matches=max_scoring_matches, - prune_perturbed_constraints=prune_perturbed_constraints, - prune_crossing_constraints=prune_crossing_constraints, - prune_atom_types=prune_atom_types, - property_map0=property_map0, - property_map1=property_map1, - mcs_kwargs={"ringMatchesRingOnly": False}, - ) - - # Only trust the retry if it extends the mapping, i.e. - # keeps every existing pair and adds sensible ones. - if ( - len(retry) > len(best) - and set(best.items()) <= set(retry.items()) - and _is_sensible_extension( - molecule0, - molecule1, - best, - retry, - property_map0, - property_map1, - ) - ): - _warnings.warn( - f"Mapping leaves heavy atoms unmapped on both sides " - f"of atom(s) {_format_flagged(flagged)} in molecule0. " - f"Relaxing 'ringMatchesRingOnly' gives a common core " - f"of {len(retry)} rather than {len(best)}. Consider " - f"passing mcs_kwargs={{'ringMatchesRingOnly': False}}." - ) - except Exception as e: - _warnings.warn(f"Unable to check the quality of the mapping: {e}") - - return result else: return _roiMatch( molecule0, @@ -1279,7 +1200,10 @@ def _matchAtoms( mol0 = molecule0._getSireObject() mol1 = molecule1._getSireObject() - # Convert the timeout to seconds and take the value as an integer. + # Convert the timeout to seconds and take the value as an integer. Keep + # the original, since the mapping check below re-enters this function, + # which expects a Time object. + orig_timeout = timeout timeout = int(timeout.seconds().value()) # Use RDKkit to find the maximum common substructure. @@ -1397,6 +1321,77 @@ def _matchAtoms( property_map1, ) + # Warn if the mapping stopped short at an attachment point where a pairable + # atom exists. This is done before the pruning below, since pruning deletes + # correctly mapped heavy atom pairs, which manufactures exactly the + # signature that the check looks for. Only check a mapping generated from + # the defaults. If the user has configured the MCS then both the baseline + # and our idea of a sensible mapping may not match their intent. This also + # stops the retry from recursing, since it passes 'mcs_kwargs'. Skip when a + # prematch is given, since the retry could then fall back on Sire MCS, + # which ignores 'mcs_kwargs', making the comparison meaningless. + if not mcs_kwargs and not prematch and mappings: + # This is a diagnostic, so it must never be able to break a call that + # would otherwise have succeeded. + try: + best = mappings[0] + + # Attachment points where the MCS stopped on both sides. + flagged = _flag_unmapped_attachments( + molecule0, molecule1, best, property_map0, property_map1 + ) + + if flagged: + # Retry with ring matching relaxed to see if it does better. + # Note that the RDKit documentation implies that + # 'completeRingsOnly' forces 'ringMatchesRingOnly', which would + # make this a no-op. It doesn't: as of RDKit 2026.03.4 the + # relaxed search still returns a larger MCS with + # 'completeRingsOnly' enabled. If a future RDKit changes this, + # the feature will silently stop firing. + # + # Pruning is disabled so that the retry is compared like for + # like against the unpruned mapping above. 'matches' and + # 'return_scores' are passed explicitly so that 'retry' is + # always a plain dict, which the comparison below relies on. + # 'prematch' is omitted since the enclosing guard means it's + # always empty. + retry = _matchAtoms( + molecule0=molecule0, + molecule1=molecule1, + scoring_function=scoring_function, + matches=1, + return_scores=False, + timeout=orig_timeout, + complete_rings_only=complete_rings_only, + max_scoring_matches=max_scoring_matches, + prune_perturbed_constraints=False, + prune_crossing_constraints=False, + prune_atom_types=False, + property_map0=property_map0, + property_map1=property_map1, + mcs_kwargs={"ringMatchesRingOnly": False}, + ) + + # Only trust the retry if it extends the mapping, i.e. keeps + # every existing pair and adds sensible ones. + if ( + len(retry) > len(best) + and set(best.items()) <= set(retry.items()) + and _is_sensible_extension( + molecule0, molecule1, best, retry, property_map0, property_map1 + ) + ): + _warnings.warn( + f"Mapping leaves heavy atoms unmapped on both sides " + f"of atom(s) {_format_flagged(flagged)} in molecule0. " + f"Relaxing 'ringMatchesRingOnly' gives a common core " + f"of {len(retry)} rather than {len(best)}. Consider " + f"passing mcs_kwargs={{'ringMatchesRingOnly': False}}." + ) + except Exception as e: + _warnings.warn(f"Unable to check the quality of the mapping: {e}") + # Optionally post-process the MCS for use with AMBER. if prune_perturbed_constraints: mappings = [ diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index b6093733..ec945082 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py @@ -1233,28 +1233,22 @@ def matchAtoms( property_map1, ) - # Optionally post-process the MCS. - if prune_perturbed_constraints: - mappings = [ - _prune_perturbed_constraints(molecule0, molecule1, x) for x in mappings - ] - if prune_crossing_constraints: - mappings = [ - _prune_crossing_constraints(molecule0, molecule1, x) for x in mappings - ] - # Warn if the mapping stopped short at an attachment point where a pairable - # atom exists. Only check a mapping generated from the defaults. If the - # user has configured the MCS then both the baseline and our idea of a - # sensible mapping may not match their intent. This also stops the retry - # below from recursing, since it passes 'mcs_kwargs'. Skip when a prematch - # is given, since the retry could then fall back on Sire MCS, which ignores - # 'mcs_kwargs', making the comparison meaningless. + # atom exists. This is done before the pruning below, since pruning deletes + # correctly mapped heavy atom pairs, which manufactures exactly the + # signature that the check looks for. Only check a mapping generated from + # the defaults. If the user has configured the MCS then both the baseline + # and our idea of a sensible mapping may not match their intent. This also + # stops the retry from recursing, since it passes 'mcs_kwargs'. Skip when a + # prematch is given, since the retry could then fall back on Sire MCS, + # which ignores 'mcs_kwargs', making the comparison meaningless. if not mcs_kwargs and not prematch and mappings: # This is a diagnostic, so it must never be able to break a call that # would otherwise have succeeded. try: best = mappings[0] + + # Attachment points where the MCS stopped on both sides. flagged = _flag_unmapped_attachments( molecule0, molecule1, best, property_map0, property_map1 ) @@ -1268,10 +1262,12 @@ def matchAtoms( # 'completeRingsOnly' enabled. If a future RDKit changes this, # the feature will silently stop firing. # - # 'matches' and 'return_scores' are passed explicitly so that - # 'retry' is always a plain dict. The comparison below relies - # on it. 'prematch' is omitted since the enclosing guard means - # it's always empty. + # Pruning is disabled so that the retry is compared like for + # like against the unpruned mapping above. 'matches' and + # 'return_scores' are passed explicitly so that 'retry' is + # always a plain dict, which the comparison below relies on. + # 'prematch' is omitted since the enclosing guard means it's + # always empty. retry = matchAtoms( molecule0, molecule1, @@ -1281,8 +1277,8 @@ def matchAtoms( return_scores=False, timeout=orig_timeout, complete_rings_only=complete_rings_only, - prune_perturbed_constraints=prune_perturbed_constraints, - prune_crossing_constraints=prune_crossing_constraints, + prune_perturbed_constraints=False, + prune_crossing_constraints=False, max_scoring_matches=max_scoring_matches, property_map0=property_map0, property_map1=property_map1, @@ -1308,6 +1304,16 @@ def matchAtoms( except Exception as e: _warnings.warn(f"Unable to check the quality of the mapping: {e}") + # Optionally post-process the MCS. + if prune_perturbed_constraints: + mappings = [ + _prune_perturbed_constraints(molecule0, molecule1, x) for x in mappings + ] + if prune_crossing_constraints: + mappings = [ + _prune_crossing_constraints(molecule0, molecule1, x) for x in mappings + ] + if matches == 1: if return_scores: return (mappings[0], scores[0]) From 0dc6b77a5eb544d6585e66779280092b34df2833 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 31 Jul 2026 13:09:53 +0100 Subject: [PATCH 4/8] Add negative, unit and branch coverage for the mapping check. --- tests/Align/test_align.py | 171 ++++++++++++++++++ tests/Sandpit/Exscientia/Align/test_align.py | 178 ++++++++++++++++++- 2 files changed, 348 insertions(+), 1 deletion(-) diff --git a/tests/Align/test_align.py b/tests/Align/test_align.py index bcc62091..6c51a44b 100644 --- a/tests/Align/test_align.py +++ b/tests/Align/test_align.py @@ -1422,3 +1422,174 @@ def test_unmapped_attachment_warning(ejm31, jmc28): ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() ) assert checked == unchecked + + +def test_unmapped_attachment_no_warning(ejm31): + # A molecule mapped to itself leaves nothing unmapped, so there is + # nothing to flag. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, ejm31) + + +@pytest.mark.skipif( + not has_antechamber or not has_tleap, + reason="Requires antechamber and tLEaP to be installed.", +) +def test_unmapped_attachment_no_warning_r_group(monkeypatch): + """ + Regression test for the check running on the pruned mapping. Pruning + deletes correctly mapped heavy atom pairs, which looks identical to an + MCS that stopped short, so these ordinary R-group edits used to be + flagged and pay for a second MCS search for nothing. + + No warning was ever emitted for them, since the gate rejected the retry, + so the flagged atoms are spied on directly rather than the warning. + """ + from BioSimSpace.Align import _align + + pairs = [ + ("Cc1ccccc1", "CCc1ccccc1"), # methyl -> ethyl + ("COc1ccccc1", "CCOc1ccccc1"), # methoxy -> ethoxy + ] + + flagged = [] + original = _align._flag_unmapped_attachments + + def _spy(*args, **kwargs): + result = original(*args, **kwargs) + flagged.append(result) + return result + + monkeypatch.setattr(_align, "_flag_unmapped_attachments", _spy) + + for smiles0, smiles1 in pairs: + molecule0 = BSS.Parameters.gaff2(smiles0).getMolecule() + molecule1 = BSS.Parameters.gaff2(smiles1).getMolecule() + + del flagged[:] + + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + warnings.filterwarnings("error", message=".*Unable to check.*") + BSS.Align.matchAtoms( + molecule0, + molecule1, + prune_perturbed_constraints=True, + prune_crossing_constraints=True, + ) + + # The check should run once, on the unpruned mapping, and find + # nothing. A second entry would mean the retry had run too. + assert flagged == [[]] + + +def test_unmapped_attachment_warning_other_branches(ejm31, jmc28): + # The check must cope with the alternative return shapes, and must not + # change what is returned in any of them. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + reference = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() + ) + + mappings = BSS.Align.matchAtoms(ejm31, jmc28, matches=5) + assert isinstance(mappings, list) + assert mappings[0] == reference + + mapping, score = BSS.Align.matchAtoms(ejm31, jmc28, return_scores=True) + assert mapping == reference + + mappings, scores = BSS.Align.matchAtoms( + ejm31, jmc28, matches=5, return_scores=True + ) + assert len(mappings) == len(scores) + assert mappings[0] == reference + + # A prematch skips the check, since the retry could fall back on the Sire + # MCS, which ignores 'mcs_kwargs'. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, jmc28, prematch={0: 0}) + + +def test_flag_unmapped_attachments(ejm31): + """ + Unit test for the attachment point check, using hand-built mappings so + that no MCS search is involved. + """ + from sire.legacy import Mol as _SireMol + + from BioSimSpace.Align._align import _flag_unmapped_attachments + + sire_mol = ejm31._sire_object + connectivity = _SireMol.Connectivity(sire_mol, _SireMol.CovalentBondHunter()) + + # Map the molecule onto itself. Nothing is unmapped, so nothing is flagged. + identity = {x: x for x in range(sire_mol.num_atoms())} + assert _flag_unmapped_attachments(ejm31, ejm31, identity) == [] + + # Find a heavy atom with a heavy atom neighbour. + for atom in sire_mol.atoms(): + idx = atom.index().value() + if atom.property("element").num_protons() == 1: + continue + neighbours = [ + i.value() + for i in connectivity.connections_to(_SireMol.AtomIdx(idx)) + if sire_mol.atom(i).property("element").num_protons() > 1 + ] + if neighbours: + break + + # Drop the neighbour from the mapping. It is now an unmapped heavy atom + # of the same element on both sides of 'idx', so 'idx' is flagged. + truncated = dict(identity) + del truncated[neighbours[0]] + assert idx in _flag_unmapped_attachments(ejm31, ejm31, truncated) + + # Hydrogens are ignored, so dropping one flags nothing. + hydrogen = next( + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() == 1 + ) + truncated = dict(identity) + del truncated[hydrogen] + assert _flag_unmapped_attachments(ejm31, ejm31, truncated) == [] + + +def test_is_sensible_extension(ejm31): + """ + Unit test for the heavy/hydrogen check on the atoms that the retry adds. + """ + from BioSimSpace.Align._align import _is_sensible_extension + + sire_mol = ejm31._sire_object + + heavy = [ + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() > 1 + ] + hydrogens = [ + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() == 1 + ] + + mapping = {heavy[0]: heavy[0]} + + # Heavy to heavy and hydrogen to hydrogen are both sensible. + extended = dict(mapping) + extended[heavy[1]] = heavy[1] + extended[hydrogens[0]] = hydrogens[0] + assert _is_sensible_extension(ejm31, ejm31, mapping, extended) + + # Pairing a heavy atom with a hydrogen is not. + extended = dict(mapping) + extended[heavy[1]] = hydrogens[0] + assert not _is_sensible_extension(ejm31, ejm31, mapping, extended) + + # Adding nothing is trivially sensible. + assert _is_sensible_extension(ejm31, ejm31, mapping, dict(mapping)) diff --git a/tests/Sandpit/Exscientia/Align/test_align.py b/tests/Sandpit/Exscientia/Align/test_align.py index 9b434931..ea3c78ae 100644 --- a/tests/Sandpit/Exscientia/Align/test_align.py +++ b/tests/Sandpit/Exscientia/Align/test_align.py @@ -7,7 +7,12 @@ from sire.legacy.Mol import AtomIdx, Element, PartialMolecule import BioSimSpace.Sandpit.Exscientia as BSS -from tests.Sandpit.Exscientia.conftest import has_antechamber, has_openff, url +from tests.Sandpit.Exscientia.conftest import ( + has_antechamber, + has_openff, + has_tleap, + url, +) @pytest.fixture(scope="session") @@ -993,3 +998,174 @@ def test_unmapped_attachment_warning(ejm31, jmc28): ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() ) assert checked == unchecked + + +def test_unmapped_attachment_no_warning(ejm31): + # A molecule mapped to itself leaves nothing unmapped, so there is + # nothing to flag. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, ejm31) + + +@pytest.mark.skipif( + not has_antechamber or not has_tleap, + reason="Requires antechamber and tLEaP to be installed.", +) +def test_unmapped_attachment_no_warning_r_group(monkeypatch): + """ + Regression test for the check running on the pruned mapping. Pruning + deletes correctly mapped heavy atom pairs, which looks identical to an + MCS that stopped short, so these ordinary R-group edits used to be + flagged and pay for a second MCS search for nothing. + + No warning was ever emitted for them, since the gate rejected the retry, + so the flagged atoms are spied on directly rather than the warning. + """ + from BioSimSpace.Sandpit.Exscientia.Align import _align + + pairs = [ + ("Cc1ccccc1", "CCc1ccccc1"), # methyl -> ethyl + ("COc1ccccc1", "CCOc1ccccc1"), # methoxy -> ethoxy + ] + + flagged = [] + original = _align._flag_unmapped_attachments + + def _spy(*args, **kwargs): + result = original(*args, **kwargs) + flagged.append(result) + return result + + monkeypatch.setattr(_align, "_flag_unmapped_attachments", _spy) + + for smiles0, smiles1 in pairs: + molecule0 = BSS.Parameters.gaff2(smiles0).getMolecule() + molecule1 = BSS.Parameters.gaff2(smiles1).getMolecule() + + del flagged[:] + + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + warnings.filterwarnings("error", message=".*Unable to check.*") + BSS.Align.matchAtoms( + molecule0, + molecule1, + prune_perturbed_constraints=True, + prune_crossing_constraints=True, + ) + + # The check should run once, on the unpruned mapping, and find + # nothing. A second entry would mean the retry had run too. + assert flagged == [[]] + + +def test_unmapped_attachment_warning_other_branches(ejm31, jmc28): + # The check must cope with the alternative return shapes, and must not + # change what is returned in any of them. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + reference = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() + ) + + mappings = BSS.Align.matchAtoms(ejm31, jmc28, matches=5) + assert isinstance(mappings, list) + assert mappings[0] == reference + + mapping, score = BSS.Align.matchAtoms(ejm31, jmc28, return_scores=True) + assert mapping == reference + + mappings, scores = BSS.Align.matchAtoms( + ejm31, jmc28, matches=5, return_scores=True + ) + assert len(mappings) == len(scores) + assert mappings[0] == reference + + # A prematch skips the check, since the retry could fall back on the Sire + # MCS, which ignores 'mcs_kwargs'. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, jmc28, prematch={0: 0}) + + +def test_flag_unmapped_attachments(ejm31): + """ + Unit test for the attachment point check, using hand-built mappings so + that no MCS search is involved. + """ + from sire.legacy import Mol as _SireMol + + from BioSimSpace.Sandpit.Exscientia.Align._align import _flag_unmapped_attachments + + sire_mol = ejm31._sire_object + connectivity = _SireMol.Connectivity(sire_mol, _SireMol.CovalentBondHunter()) + + # Map the molecule onto itself. Nothing is unmapped, so nothing is flagged. + identity = {x: x for x in range(sire_mol.num_atoms())} + assert _flag_unmapped_attachments(ejm31, ejm31, identity) == [] + + # Find a heavy atom with a heavy atom neighbour. + for atom in sire_mol.atoms(): + idx = atom.index().value() + if atom.property("element").num_protons() == 1: + continue + neighbours = [ + i.value() + for i in connectivity.connections_to(_SireMol.AtomIdx(idx)) + if sire_mol.atom(i).property("element").num_protons() > 1 + ] + if neighbours: + break + + # Drop the neighbour from the mapping. It is now an unmapped heavy atom + # of the same element on both sides of 'idx', so 'idx' is flagged. + truncated = dict(identity) + del truncated[neighbours[0]] + assert idx in _flag_unmapped_attachments(ejm31, ejm31, truncated) + + # Hydrogens are ignored, so dropping one flags nothing. + hydrogen = next( + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() == 1 + ) + truncated = dict(identity) + del truncated[hydrogen] + assert _flag_unmapped_attachments(ejm31, ejm31, truncated) == [] + + +def test_is_sensible_extension(ejm31): + """ + Unit test for the heavy/hydrogen check on the atoms that the retry adds. + """ + from BioSimSpace.Sandpit.Exscientia.Align._align import _is_sensible_extension + + sire_mol = ejm31._sire_object + + heavy = [ + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() > 1 + ] + hydrogens = [ + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() == 1 + ] + + mapping = {heavy[0]: heavy[0]} + + # Heavy to heavy and hydrogen to hydrogen are both sensible. + extended = dict(mapping) + extended[heavy[1]] = heavy[1] + extended[hydrogens[0]] = hydrogens[0] + assert _is_sensible_extension(ejm31, ejm31, mapping, extended) + + # Pairing a heavy atom with a hydrogen is not. + extended = dict(mapping) + extended[heavy[1]] = hydrogens[0] + assert not _is_sensible_extension(ejm31, ejm31, mapping, extended) + + # Adding nothing is trivially sensible. + assert _is_sensible_extension(ejm31, ejm31, mapping, dict(mapping)) From 63c5704f60117b766af9d3b3d3545453443a7657 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 31 Jul 2026 13:26:03 +0100 Subject: [PATCH 5/8] Only check the mapping where the advice can be acted on. --- src/BioSimSpace/Align/_align.py | 27 +++++++++-- .../Sandpit/Exscientia/Align/_align.py | 22 +++++++-- tests/Align/test_align.py | 46 +++++++++++++++++++ tests/Sandpit/Exscientia/Align/test_align.py | 28 +++++++++++ 4 files changed, 117 insertions(+), 6 deletions(-) diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 40e20c86..9f1c0580 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -744,6 +744,7 @@ def matchAtoms( property_map0={}, property_map1={}, mcs_kwargs={}, + _check_mapping=True, ): """ Find mappings between atom indices in molecule0 to those in molecule1. @@ -912,6 +913,7 @@ def matchAtoms( property_map0=property_map0, property_map1=property_map1, mcs_kwargs=mcs_kwargs, + _check_mapping=_check_mapping, ) else: return _roiMatch( @@ -1109,6 +1111,7 @@ def _matchAtoms( property_map0={}, property_map1={}, mcs_kwargs={}, + _check_mapping=True, ): import sys as _sys @@ -1330,9 +1333,12 @@ def _matchAtoms( # stops the retry from recursing, since it passes 'mcs_kwargs'. Skip when a # prematch is given, since the retry could then fall back on Sire MCS, # which ignores 'mcs_kwargs', making the comparison meaningless. - if not mcs_kwargs and not prematch and mappings: + if _check_mapping and not mcs_kwargs and not prematch and mappings: # This is a diagnostic, so it must never be able to break a call that - # would otherwise have succeeded. + # would otherwise have succeeded. The warning itself is emitted outside + # the guard, since it would otherwise be swallowed and re-reported as a + # failure whenever the user has promoted warnings to errors. + message = None try: best = mappings[0] @@ -1382,7 +1388,7 @@ def _matchAtoms( molecule0, molecule1, best, retry, property_map0, property_map1 ) ): - _warnings.warn( + message = ( f"Mapping leaves heavy atoms unmapped on both sides " f"of atom(s) {_format_flagged(flagged)} in molecule0. " f"Relaxing 'ringMatchesRingOnly' gives a common core " @@ -1392,6 +1398,9 @@ def _matchAtoms( except Exception as e: _warnings.warn(f"Unable to check the quality of the mapping: {e}") + if message is not None: + _warnings.warn(message) + # Optionally post-process the MCS for use with AMBER. if prune_perturbed_constraints: mappings = [ @@ -1732,6 +1741,9 @@ def _roiMatch( mapping = matchAtoms( res0_extracted, res1_extracted, + # The mapping check would report indices that are local to the + # extracted residue, not to molecule0 as its message claims. + _check_mapping=False, ) # Look up the absolute atom indices in the molecule if not using a custom ROI mapping. @@ -1955,6 +1967,9 @@ def _rmsdAlign(molecule0, molecule1, mapping=None, property_map0={}, property_ma molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) # Extract the Sire molecule from each BioSimSpace molecule. @@ -2157,6 +2172,9 @@ def _flexAlign( molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) # Convert the mapping to AtomIdx key:value pairs. @@ -2614,6 +2632,9 @@ def viewMapping( molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) molecule0 = rmsdAlign(molecule0, molecule1, mapping) diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index ec945082..554e002c 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py @@ -907,6 +907,7 @@ def matchAtoms( property_map0={}, property_map1={}, mcs_kwargs={}, + _check_mapping=True, ): """ Find mappings between atom indices in molecule0 to those in molecule1. @@ -1242,9 +1243,12 @@ def matchAtoms( # stops the retry from recursing, since it passes 'mcs_kwargs'. Skip when a # prematch is given, since the retry could then fall back on Sire MCS, # which ignores 'mcs_kwargs', making the comparison meaningless. - if not mcs_kwargs and not prematch and mappings: + if _check_mapping and not mcs_kwargs and not prematch and mappings: # This is a diagnostic, so it must never be able to break a call that - # would otherwise have succeeded. + # would otherwise have succeeded. The warning itself is emitted outside + # the guard, since it would otherwise be swallowed and re-reported as a + # failure whenever the user has promoted warnings to errors. + message = None try: best = mappings[0] @@ -1294,7 +1298,7 @@ def matchAtoms( molecule0, molecule1, best, retry, property_map0, property_map1 ) ): - _warnings.warn( + message = ( f"Mapping leaves heavy atoms unmapped on both sides " f"of atom(s) {_format_flagged(flagged)} in molecule0. " f"Relaxing 'ringMatchesRingOnly' gives a common core " @@ -1304,6 +1308,9 @@ def matchAtoms( except Exception as e: _warnings.warn(f"Unable to check the quality of the mapping: {e}") + if message is not None: + _warnings.warn(message) + # Optionally post-process the MCS. if prune_perturbed_constraints: mappings = [ @@ -1415,6 +1422,9 @@ def rmsdAlign(molecule0, molecule1, mapping=None, property_map0={}, property_map molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) # Extract the Sire molecule from each BioSimSpace molecule. @@ -1572,6 +1582,9 @@ def flexAlign( molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) # Convert the mapping to AtomIdx key:value pairs. @@ -1905,6 +1918,9 @@ def viewMapping( molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) molecule0 = rmsdAlign(molecule0, molecule1, mapping) diff --git a/tests/Align/test_align.py b/tests/Align/test_align.py index 6c51a44b..fa429504 100644 --- a/tests/Align/test_align.py +++ b/tests/Align/test_align.py @@ -1593,3 +1593,49 @@ def test_is_sensible_extension(ejm31): # Adding nothing is trivially sensible. assert _is_sensible_extension(ejm31, ejm31, mapping, dict(mapping)) + + +def test_unmapped_attachment_check_suppressed(ejm31, jmc28): + """ + The check should only fire where the user can act on its advice, i.e. + where 'mcs_kwargs' can be passed through. It should also never fire on + the ROI path, where the flagged indices would be local to the extracted + residue rather than to molecule0 as the message claims. + """ + # These functions don't take 'mcs_kwargs'. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.rmsdAlign(ejm31, jmc28) + BSS.Align.flexAlign(ejm31, jmc28) + + # 'merge' does, so the check stays on. + with pytest.warns(UserWarning, match="ringMatchesRingOnly"): + BSS.Align.merge(ejm31, jmc28, force=True) + + +def test_unmapped_attachment_check_suppressed_roi(protein_inputs): + # The ROI path maps each residue of interest separately, so any flagged + # indices would be local to that residue rather than to molecule0. + proteins, protein_mapping, roi = protein_inputs + p0 = BSS.IO.readMolecules( + BSS.IO.expand(BSS.tutorialUrl(), f"{proteins}_mut_peptide.pdb") + )[0] + p1 = BSS.IO.readMolecules( + BSS.IO.expand(BSS.tutorialUrl(), f"{proteins}_wt_peptide.pdb") + )[0] + + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + assert BSS.Align.matchAtoms(p0, p1, roi=roi) == protein_mapping + + +def test_unmapped_attachment_warning_not_swallowed(ejm31, jmc28): + # Promoting the warning to an error must surface the warning itself, not + # a report that the check failed. The warning is emitted outside the + # try/except that guards the check for exactly this reason. + # Anchored, since the wrapped "Unable to check ..." message quotes the + # original and would otherwise match too. + with pytest.raises(UserWarning, match=r"^Mapping leaves heavy atoms"): + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, jmc28) diff --git a/tests/Sandpit/Exscientia/Align/test_align.py b/tests/Sandpit/Exscientia/Align/test_align.py index ea3c78ae..0e413b82 100644 --- a/tests/Sandpit/Exscientia/Align/test_align.py +++ b/tests/Sandpit/Exscientia/Align/test_align.py @@ -1169,3 +1169,31 @@ def test_is_sensible_extension(ejm31): # Adding nothing is trivially sensible. assert _is_sensible_extension(ejm31, ejm31, mapping, dict(mapping)) + + +def test_unmapped_attachment_check_suppressed(ejm31, jmc28): + """ + The check should only fire where the user can act on its advice, i.e. + where 'mcs_kwargs' can be passed through. + """ + # These functions don't take 'mcs_kwargs'. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.rmsdAlign(ejm31, jmc28) + BSS.Align.flexAlign(ejm31, jmc28) + + # 'merge' does, so the check stays on. + with pytest.warns(UserWarning, match="ringMatchesRingOnly"): + BSS.Align.merge(ejm31, jmc28, force=True) + + +def test_unmapped_attachment_warning_not_swallowed(ejm31, jmc28): + # Promoting the warning to an error must surface the warning itself, not + # a report that the check failed. The warning is emitted outside the + # try/except that guards the check for exactly this reason. + # Anchored, since the wrapped "Unable to check ..." message quotes the + # original and would otherwise match too. + with pytest.raises(UserWarning, match=r"^Mapping leaves heavy atoms"): + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, jmc28) From ea14be897f42c2725d20a94b00d1a1872d0dacd3 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 31 Jul 2026 14:33:59 +0100 Subject: [PATCH 6/8] Restore the viewMapping check and tidy the private parameter. --- src/BioSimSpace/Align/_align.py | 19 +++++++++++-------- .../Sandpit/Exscientia/Align/_align.py | 18 ++++++++++-------- tests/Align/test_align.py | 8 ++++++-- tests/Sandpit/Exscientia/Align/test_align.py | 8 ++++++-- 4 files changed, 33 insertions(+), 20 deletions(-) diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 9f1c0580..7b1391d7 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -744,6 +744,7 @@ def matchAtoms( property_map0={}, property_map1={}, mcs_kwargs={}, + *, _check_mapping=True, ): """ @@ -1111,6 +1112,7 @@ def _matchAtoms( property_map0={}, property_map1={}, mcs_kwargs={}, + *, _check_mapping=True, ): import sys as _sys @@ -1334,10 +1336,11 @@ def _matchAtoms( # prematch is given, since the retry could then fall back on Sire MCS, # which ignores 'mcs_kwargs', making the comparison meaningless. if _check_mapping and not mcs_kwargs and not prematch and mappings: - # This is a diagnostic, so it must never be able to break a call that - # would otherwise have succeeded. The warning itself is emitted outside - # the guard, since it would otherwise be swallowed and re-reported as a - # failure whenever the user has promoted warnings to errors. + # This is a diagnostic, so a failure inside it is reported rather than + # raised. Note that this only holds while warnings are warnings: if the + # user has promoted them to errors then either notice below will raise + # out of here. The warning itself is emitted outside the guard, since + # it would otherwise be swallowed and re-reported as a failure. message = None try: best = mappings[0] @@ -1377,10 +1380,13 @@ def _matchAtoms( property_map0=property_map0, property_map1=property_map1, mcs_kwargs={"ringMatchesRingOnly": False}, + _check_mapping=False, ) # Only trust the retry if it extends the mapping, i.e. keeps - # every existing pair and adds sensible ones. + # every existing pair and adds sensible ones. The subset test is + # sensitive to symmetry: relabelling a symmetric ring the other + # way round discards the retry even though it is equivalent. if ( len(retry) > len(best) and set(best.items()) <= set(retry.items()) @@ -2632,9 +2638,6 @@ def viewMapping( molecule1, property_map0=property_map0, property_map1=property_map1, - # This function doesn't take 'mcs_kwargs', so the advice from the - # mapping check can't be acted on here. - _check_mapping=False, ) molecule0 = rmsdAlign(molecule0, molecule1, mapping) diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index 554e002c..d4be9b2e 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py @@ -907,6 +907,7 @@ def matchAtoms( property_map0={}, property_map1={}, mcs_kwargs={}, + *, _check_mapping=True, ): """ @@ -1244,10 +1245,11 @@ def matchAtoms( # prematch is given, since the retry could then fall back on Sire MCS, # which ignores 'mcs_kwargs', making the comparison meaningless. if _check_mapping and not mcs_kwargs and not prematch and mappings: - # This is a diagnostic, so it must never be able to break a call that - # would otherwise have succeeded. The warning itself is emitted outside - # the guard, since it would otherwise be swallowed and re-reported as a - # failure whenever the user has promoted warnings to errors. + # This is a diagnostic, so a failure inside it is reported rather than + # raised. Note that this only holds while warnings are warnings: if the + # user has promoted them to errors then either notice below will raise + # out of here. The warning itself is emitted outside the guard, since + # it would otherwise be swallowed and re-reported as a failure. message = None try: best = mappings[0] @@ -1287,10 +1289,13 @@ def matchAtoms( property_map0=property_map0, property_map1=property_map1, mcs_kwargs={"ringMatchesRingOnly": False}, + _check_mapping=False, ) # Only trust the retry if it extends the mapping, i.e. keeps - # every existing pair and adds sensible ones. + # every existing pair and adds sensible ones. The subset test is + # sensitive to symmetry: relabelling a symmetric ring the other + # way round discards the retry even though it is equivalent. if ( len(retry) > len(best) and set(best.items()) <= set(retry.items()) @@ -1918,9 +1923,6 @@ def viewMapping( molecule1, property_map0=property_map0, property_map1=property_map1, - # This function doesn't take 'mcs_kwargs', so the advice from the - # mapping check can't be acted on here. - _check_mapping=False, ) molecule0 = rmsdAlign(molecule0, molecule1, mapping) diff --git a/tests/Align/test_align.py b/tests/Align/test_align.py index fa429504..dca7b31b 100644 --- a/tests/Align/test_align.py +++ b/tests/Align/test_align.py @@ -1450,7 +1450,9 @@ def test_unmapped_attachment_no_warning_r_group(monkeypatch): pairs = [ ("Cc1ccccc1", "CCc1ccccc1"), # methyl -> ethyl + ("CCc1ccccc1", "CCCc1ccccc1"), # ethyl -> propyl ("COc1ccccc1", "CCOc1ccccc1"), # methoxy -> ethoxy + ("O=C(N)c1ccccc1", "O=C(N)c1ccccc1Cl"), # hydrogen -> chlorine ] flagged = [] @@ -1602,11 +1604,13 @@ def test_unmapped_attachment_check_suppressed(ejm31, jmc28): the ROI path, where the flagged indices would be local to the extracted residue rather than to molecule0 as the message claims. """ - # These functions don't take 'mcs_kwargs'. + # 'rmsdAlign' doesn't take 'mcs_kwargs'. 'flexAlign' is suppressed for the + # same reason, but isn't exercised here since it needs fkcombu. Nor is + # 'viewMapping', which keeps the check but returns early outside a + # notebook, before it ever reaches 'matchAtoms'. with warnings.catch_warnings(): warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") BSS.Align.rmsdAlign(ejm31, jmc28) - BSS.Align.flexAlign(ejm31, jmc28) # 'merge' does, so the check stays on. with pytest.warns(UserWarning, match="ringMatchesRingOnly"): diff --git a/tests/Sandpit/Exscientia/Align/test_align.py b/tests/Sandpit/Exscientia/Align/test_align.py index 0e413b82..6d265a16 100644 --- a/tests/Sandpit/Exscientia/Align/test_align.py +++ b/tests/Sandpit/Exscientia/Align/test_align.py @@ -1026,7 +1026,9 @@ def test_unmapped_attachment_no_warning_r_group(monkeypatch): pairs = [ ("Cc1ccccc1", "CCc1ccccc1"), # methyl -> ethyl + ("CCc1ccccc1", "CCCc1ccccc1"), # ethyl -> propyl ("COc1ccccc1", "CCOc1ccccc1"), # methoxy -> ethoxy + ("O=C(N)c1ccccc1", "O=C(N)c1ccccc1Cl"), # hydrogen -> chlorine ] flagged = [] @@ -1176,11 +1178,13 @@ def test_unmapped_attachment_check_suppressed(ejm31, jmc28): The check should only fire where the user can act on its advice, i.e. where 'mcs_kwargs' can be passed through. """ - # These functions don't take 'mcs_kwargs'. + # 'rmsdAlign' doesn't take 'mcs_kwargs'. 'flexAlign' is suppressed for the + # same reason, but isn't exercised here since it needs fkcombu. Nor is + # 'viewMapping', which keeps the check but returns early outside a + # notebook, before it ever reaches 'matchAtoms'. with warnings.catch_warnings(): warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") BSS.Align.rmsdAlign(ejm31, jmc28) - BSS.Align.flexAlign(ejm31, jmc28) # 'merge' does, so the check stays on. with pytest.warns(UserWarning, match="ringMatchesRingOnly"): From 155728e8303b5a9a8f21e8260a532843bd018ebd Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 31 Jul 2026 14:41:34 +0100 Subject: [PATCH 7/8] Note that the subset test is sensitive to equivalent relabellings. --- src/BioSimSpace/Align/_align.py | 4 ++-- src/BioSimSpace/Sandpit/Exscientia/Align/_align.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 7b1391d7..28060aef 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -1385,8 +1385,8 @@ def _matchAtoms( # Only trust the retry if it extends the mapping, i.e. keeps # every existing pair and adds sensible ones. The subset test is - # sensitive to symmetry: relabelling a symmetric ring the other - # way round discards the retry even though it is equivalent. + # sensitive to relabelling: an equivalent mapping that traverses + # a ring the other way, or permutes hydrogens, is discarded. if ( len(retry) > len(best) and set(best.items()) <= set(retry.items()) diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index d4be9b2e..db4675e3 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py @@ -1294,8 +1294,8 @@ def matchAtoms( # Only trust the retry if it extends the mapping, i.e. keeps # every existing pair and adds sensible ones. The subset test is - # sensitive to symmetry: relabelling a symmetric ring the other - # way round discards the retry even though it is equivalent. + # sensitive to relabelling: an equivalent mapping that traverses + # a ring the other way, or permutes hydrogens, is discarded. if ( len(retry) > len(best) and set(best.items()) <= set(retry.items()) From 934131b2e52cd921f24fff25b276944cd7fa5b8c Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 31 Jul 2026 17:47:55 +0100 Subject: [PATCH 8/8] Take the mapping check flag off the public matchAtoms signature. --- src/BioSimSpace/Align/_align.py | 10 +++------- src/BioSimSpace/FreeEnergy/_atm.py | 9 +++++++-- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 28060aef..72f0f50d 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -744,8 +744,6 @@ def matchAtoms( property_map0={}, property_map1={}, mcs_kwargs={}, - *, - _check_mapping=True, ): """ Find mappings between atom indices in molecule0 to those in molecule1. @@ -914,7 +912,6 @@ def matchAtoms( property_map0=property_map0, property_map1=property_map1, mcs_kwargs=mcs_kwargs, - _check_mapping=_check_mapping, ) else: return _roiMatch( @@ -1105,7 +1102,6 @@ def _matchAtoms( timeout=5 * _Units.Time.second, complete_rings_only=True, max_scoring_matches=1000, - roi=None, prune_perturbed_constraints=False, prune_crossing_constraints=False, prune_atom_types=False, @@ -1744,7 +1740,7 @@ def _roiMatch( ) mapping = None else: - mapping = matchAtoms( + mapping = _matchAtoms( res0_extracted, res1_extracted, # The mapping check would report indices that are local to the @@ -1968,7 +1964,7 @@ def _rmsdAlign(molecule0, molecule1, mapping=None, property_map0={}, property_ma # Get the best match atom mapping. else: - mapping = matchAtoms( + mapping = _matchAtoms( molecule0, molecule1, property_map0=property_map0, @@ -2173,7 +2169,7 @@ def _flexAlign( # Get the best match atom mapping. else: - mapping = matchAtoms( + mapping = _matchAtoms( molecule0, molecule1, property_map0=property_map0, diff --git a/src/BioSimSpace/FreeEnergy/_atm.py b/src/BioSimSpace/FreeEnergy/_atm.py index b1f126f2..e7d237e5 100644 --- a/src/BioSimSpace/FreeEnergy/_atm.py +++ b/src/BioSimSpace/FreeEnergy/_atm.py @@ -593,7 +593,10 @@ def _makeSystemFromThree(protein, ligand_bound, ligand_free, displacement): BioSimSpace._SireWrappers.System The system for the ATM simulation. """ - from ..Align import matchAtoms as _matchAtoms + # The private form is used so that the mapping check can be switched + # off below. It is otherwise identical to 'matchAtoms', which is a + # dispatcher on 'roi' and is never given one here. + from ..Align._align import _matchAtoms from ..Align import rmsdAlign as _rmsdAlign from ..Types import Vector as _Vector @@ -683,7 +686,9 @@ def _findTranslationVector(system, displacement, protein, ligand): out_of_protein = displacement.value() * initial_normal_vector return out_of_protein - mapping = _matchAtoms(ligand_free, ligand_bound) + # ATM doesn't expose 'mcs_kwargs', nor a way to pass in a mapping, so + # the advice from the mapping check can't be acted on here. + mapping = _matchAtoms(ligand_free, ligand_bound, _check_mapping=False) ligand_free_aligned = _rmsdAlign(ligand_free, ligand_bound, mapping) prot_lig1 = (protein + ligand_bound).toSystem()