Repository navigation
Conversation
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #1044 +/- ##
==========================================
+ Coverage 65.66% 67.09% +1.42%
==========================================
Files 121 123 +2
Lines 41073 43741 +2668
Branches 10566 11170 +604
==========================================
+ Hits 26972 29347 +2375
- Misses 11085 11236 +151
- Partials 3016 3158 +142
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
10e4352 to
d545be9
Compare
| if uni_is_product: | ||
| try: | ||
| atom_map = cached_map_rxn(rxn=rxn, product_dict_index=i, map_cache=map_cache) | ||
| atom_map = cached_map_rxn(rxn=rxn, product_dict_index=i, map_cache=map_cache, |
There was a problem hiding this comment.
the requested atom map can still be ignored in some cases. when a molecule splits into fragments, parts of the code still choose their own atom correspondence. so you ask for map a, and you may instead receive another map, reported as successful
There was a problem hiding this comment.
You're right, and it was worse than partial — fixed in 9a8d039 / d8879b8.
Where it was ignored: the dissociation route. interpolate_addition only consulted the map inside if uni_is_product:; in the else branch the recipe bonds are already in uni ordering, so the map was read nowhere. A path's changing bonds come from the RMG recipe, not from the map, so every requested map built every path. For CCOO• → C2H4 + HO2, two maps with genuinely different reaction centers — H4 eliminated (broken (0,4),(1,2), formed (3,4)) vs H5 ((0,5),(1,2), (3,5)) — returned byte-identical 7-guess sets, every one reported successful.
The fix makes the requested center select the path: get_requested_reactive_bonds() reads the center from the map (reactant ordering, converted through the map when the uni species is the product) and paths changing bonds the map does not are skipped. Subset rather than equality, so an enumerated map carrying extra spectator-permutation bonds still selects the path whose center it contains; both index-space readings are accepted since reverse-discovered paths report recipe bonds in uni ordering already. The two maps above now give different guesses (3 each, each from its own path). Inert when no map is requested.
Where I could not reproduce it: the addition route. My first reading of the data agreed with you here too, and it was wrong — the map pairs that returned identical geometries share a reaction center in the reactant well, so identical output is correct. Comparing centers rather than geometries is what separated the two cases.
On fragment matching specifically: map_and_verify_fragments does re-derive a correspondence by subgraph isomorphism, independent of the requested map. But its species assignment only feeds product_groups, which is then sorted by fragment size, so I could not make it change a guess geometry in either direction. The place the map genuinely cannot be honored is the combinatorial-fragmentation fallback, which has no template paths to select between — that now logs a warning that its guesses may describe a different correspondence, rather than letting them pass as the caller's map.
| atom_map = cached_map_rxn(rxn=rxn, product_dict_index=i, map_cache=map_cache, | ||
| forced_atom_map=forced_atom_map) |
There was a problem hiding this comment.
The isomerization route still builds every path under the requested map. forced_atom_map replaces the per-path map here, but bb/fb still come from each path's recipe, and nothing compares them with the requested centre. That is the same gap 9a8d039 closed for additions.
For the test reaction sBu → nBu (6 product_dicts), I requested the cluster representative whose centre is the 1,3-shift of H4 (formed (2,4), broken (0,4)). The adapter returned 15 successful guesses. 8 of them move a C3 methyl hydrogen (H10/H11/H12) to C2, which is the 1,2-shift centre; only 5 move H4.
The bespoke strategies that run before this call (ring-closure bespoke, bespoke_retro_da_bicyclic), _try_bespoke_family_fallbacks and the direct-contraction supplement also read recipe bonds of all paths and never see the map.
Suggested fix: compute the requested centre once with get_formed_and_broken_bonds() under the map. For an isomerization no index conversion is needed. Then skip any path whose recipe bb ∪ fb does not match it, both here and in the bespoke and fallback branches. If no path matches, warn and produce nothing for that map.
| if forced_atom_map is not None and cut_lists: | ||
| # This strategy has no template paths to select between: the cuts come from combinatorial bond | ||
| # breaking and the fragment correspondence from subgraph isomorphism, neither of which the | ||
| # requested map constrains. Guesses from here may describe a correspondence other than the one | ||
| # asked for, so say so rather than letting them pass silently as the caller's map. | ||
| logger.warning(f'The linear TS search adapter is building addition/dissociation TS guesses for ' | ||
| f'{rxn.label} by combinatorial bond cutting, which does not consult the requested ' | ||
| f'atom map; the resulting guesses may describe a different correspondence.') |
There was a problem hiding this comment.
The fragmentation supplement still returns a guess for the wrong hydrogen under the requested map, and it is recorded as success=True. With MAP_H5 from the new test, the w=0.65 guess comes from frag_fallback (bb [(1,2)]): the template path for H5 was rejected at that weight for colliding atoms. In that guess d(O3,H4)=2.59 Å and d(O3,H5)=4.03 Å, so it is the H4 elimination, identical to MAP_H4's w=0.65 guess. So "3 guesses each, each from its own path" is not quite what happens.
The warning also fires once per weight per map (6 times in a two-map run on this reaction), whenever cut_lists is non-empty, whether or not Strategy 2 contributes anything. It also says the map does not constrain these cuts, but for a dissociation the requested centre's broken bonds are the cut.
Suggestion: when a map is requested, drop Strategy-2 records that are inconsistent with requested_reactive (or, if they can't be checked, don't record them under that map), and warn once per job, only when such a record is actually kept. Outside this hunk, the XY_elimination builder (early return [ts_family]), the Retroene builder (product_dicts[0]) and the Baeyer-Villiger_step2 builder also ignore the map.
| if requested_reactive is None: | ||
| return True | ||
| for bonds in (reactive_bonds_uni, reactive_bonds_r): | ||
| if bonds and {canonical_bond(a, b) for a, b in bonds}.issubset(requested_reactive): |
There was a problem hiding this comment.
I think a subset test is too loose here. Relabelling equivalent atoms on the same centre (e.g. swapping the H's of one methyl) changes no bond, so it adds nothing to formed ∪ broken. Extra bonds only appear when the map moves atoms between different neighbours, and such a map describes additional bond breaking and forming: a different, non-elementary correspondence. With a subset test that map still selects the recipe path and the guess gets attributed to it, which is the silent substitution this PR is trying to stop.
Unless there is a case I'm missing, an equality check of the recipe's bb ∪ fb against the map's formed ∪ broken is the right test, and a map that matches no path should be reported and skipped. That would also be the guard against element-preserving but chemically impossible maps, which atom_map_fits_rxn cannot catch.
| """ | ||
| if requested_reactive is None: | ||
| return True | ||
| for bonds in (reactive_bonds_uni, reactive_bonds_r): |
There was a problem hiding this comment.
Accepting either reading means comparing bonds from two different index spaces. For an addition, reactive_bonds_r uses multi-species ordering while requested_reactive uses uni ordering, so a forward path can pass the check by a coincidence of indices. product_dict['discovered_in_reverse'] already records which reading applies (interpolate_isomerization passes it into _build_path_context). Could you pass it in and compare only the one correct set?
| def test_requested_center_selects_the_path(self): | ||
| """Test that two maps eliminating a different H give different guesses rather than identical ones.""" | ||
| xyzs_h4 = self._guess_xyzs(self.MAP_H4) | ||
| xyzs_h5 = self._guess_xyzs(self.MAP_H5) | ||
| self.assertGreater(len(xyzs_h4), 0) | ||
| self.assertGreater(len(xyzs_h5), 0) | ||
| self.assertNotEqual(xyzs_h4, xyzs_h5) |
There was a problem hiding this comment.
This only checks that the two maps give different outputs, so it can't tell the right path from the wrong one. I inverted the filter at line 2129 of linear.py, so it keeps exactly the paths the map does not describe, and every new test still passed.
Something like this would pin it: for every MAP_H4 guess, d(O3,H4) < d(O3,H5), and the reverse for MAP_H5. Right now that assertion fails on the w=0.65 guess (see the comment on the fragmentation supplement), which is exactly the case to catch.
| except Exception as e: | ||
| logger.debug(f'Linear (rxn={rxn.label}): could not determine well symbols to validate an atom map: {e}.') | ||
| return False | ||
| if len(atom_map) != len(r_symbols) or len(r_symbols) != len(p_symbols): | ||
| return False |
There was a problem hiding this comment.
When this returns False because of an exception, the caller still reports "wrong number of atoms or mapping atoms onto different elements", which won't be true. Consider returning a reason, or warning here with the exception text.
| if not path_matches_requested_center(list(bb_conc) + list(fb_conc), | ||
| list(sb_conc) + list(cb_conc), | ||
| requested_reactive): | ||
| logger.debug(f'Linear addition (rxn={rxn.label}): the concerted path changes bonds that the ' | ||
| f'requested atom map does not; skipping path.') | ||
| continue |
There was a problem hiding this comment.
Test gaps: no test notices if this concerted-path filter is removed. Deleting the rxn._atom_map = list(forced_atom_map) line at 1051 is also invisible to the tests, so the claim that the fallback paths honour the requested map is unverified. Hard-coding the species count to 1 in get_well_symbols is invisible too, since there is no A+A case; C2H6 <=> CH3 + CH3 works and would make a cheap test. The addition direction (uni_is_product=True) of get_requested_reactive_bonds has no test either.
| # A map that cannot index the center bonds disables filtering rather than raising. | ||
| self.assertIsNone(get_requested_reactive_bonds(rxn, [0, 1], uni_is_product=True)) |
There was a problem hiding this comment.
[0, 1] doesn't reach the IndexError branch this comment describes; it fails earlier, inside get_formed_and_broken_bonds, with a ValueError from list.index. Since the adapter validates the map's length before getting here, the IndexError branch at linear.py 1858-1861 looks unreachable and could be dropped.
| This is the reaction center the caller asked for. The changing bonds of an RMG template path come | ||
| from the recipe and not from the atom map, so without comparing against this center a requested | ||
| map has almost no influence on which path (which leaving or migrating atom, say) gets built: the | ||
| map is otherwise only consulted to convert bond indices, and only when the unimolecular species is | ||
| the product. ``get_formed_and_broken_bonds`` reports in reactant ordering, which is already the | ||
| unimolecular ordering for a dissociation and is mapped through the requested map for an addition. |
There was a problem hiding this comment.
Minor style point (project convention): this docstring, and the one for path_matches_requested_center at 1871-1875, explain why the function exists rather than what it does. The same goes for the new inline # comments at 1044-1045, 1049-1050, 1053-1055, 2092-2094 and 2538-2541, and in the tests at 6340, 6415-6418, 6466, 6471 and 6477-6482. Could the docstrings describe behaviour only, and the rationale move to the commit message?
| return all(p_symbols[p_index] == r_symbol for r_symbol, p_index in zip(r_symbols, atom_map)) | ||
|
|
||
|
|
||
|
|
There was a problem hiding this comment.
Nit: three blank lines before class LinearAdapter (PEP 8 E303). Same before class TestLinearAdapterAtomMapInAddition and before if __name__ in the test file (6409-6411 and 6484-6486).
Summary
LinearAdapterderives an atom map per reaction path viamap_rxnand never reads a map set on the reaction, so pre-setting one changed nothing and was not reported. This addsLinearAdapter(atom_map=...), taking one map or several: each requested map replaces the derived one for a full pass over the weight grid, so a TS guess can be requested for a specific reactants-to-products correspondence. Requested maps are validated (permutation, atom count, element-preserving) and unusable ones are reported rather than silently dropped. A pre-setrxn.atom_mapstill does not take part in the derivation - other adapters populate it as a side effect of reading it - but the adapter now warns that it is unused and points at the new argument, and the docstrings say guesses are generated per product dictionary.Closes #1042