Transition metal oxidation state proposed fixes (#35)

## Summary

Five small patches that fix YARP 3.0's behavior on transition-metal
oxidation-state extraction, validated on a 181,450-archive corpus of
DFT-optimized organometallic reactions. All commits are on
`zhao-final-20260619` and pass the existing `test/yarpecule/lewis/`
suite (55 / 55).

## What's in the bundle

1. **`yarp/util/properties.py`** — restore 5d/4d TM dict entries that
   were set to `None` in master (`el_valence`, `el_n_deficient`,
   `el_n_expand_octet`, `el_expand_octet`). Also restore the 6th-period
   `el_pol` entries (Cs..Rn were missing entirely, causing
   `KeyError('cs')`). Drops a duplicate `"rh": 66` key in the same dict
   literal. **Trivial bug fix.**

2. **Patch C — remove the outer `for ind in range(0, len(bond_mats)):`
   loop in `gen_all_lstructs`** (`find_lewis.py:289`). Every recursive
   call already passes `ind = len(bond_mats)-1`, so the outer loop
   re-walks every previously discovered BEM at every recursion depth.
   Removing it yields **~10× speedup on TM benchmarks** (32,972 s →
   1,718 s in isolation on a 144-archive stratified sample).
   **Performance + correctness, zero organic test regression.**

3. **Patch D — remove the radical-radical "move 4-bis" yield block** in
   `valid_moves` (`find_lewis.py:447`). The inner `for k` loop iterates
   but `k` is never used in the yielded move — block emits the same
   `(i, j)` coupling move once per qualifying `k`, all duplicates.
   **Dead-code cleanup.**

4. **Patch F — conditional seed-BEM re-pool** in
   `lewis_structure._gen_bond_el_mat`. The re-pool step (after pass-2,
   before `adjust_metals`) keeps the new behavior for organic systems
   (where the existing pytest cases depend on it) but disables it when
   any transition metal is present. On TM species the pass-1 seed BEMs
   are optimized without aromaticity (`w_aro = 0`) and so often encode
   non-aromatic ligand configurations that `adjust_metals` then
   dative-izes via Z-bonds — draining the metal diagonal to zero and
   producing impossible high OS values (Cp-Cr → Cr(VI), Cp-Ir → Ir(V),
   Pt(0) → Pt(VII)). Eliminating those single-archive failures was the
   most impactful chemistry change on our bench. **Could plausibly be
   done a different way upstream** (see "open question" below).

5. **`w_rad` default revert** (`bem_score.py`, `lewis_structure.py`) —
   from `-0.01` to `+0.1`. **COSMETIC, not a behavior change.** The
   sign flip in `w_rad` is algebraically compensated by a sign+scale
   flip in how `rad_env` is computed inside `bmat_score` vs how old
   `find_lewis` passed it in. We isolated this on its own branch and
   confirmed zero chemistry effect; reverted only to keep `w_rad > 0`
   for downstream tooling that calls `bmat_score` directly.

## What's NOT in the bundle

We deliberately exclude one change that was in the historical
`fed9385` commit: switching the 2nd `gen_all_lstructs` call from
greedy descent (`min_opt=True`) to exploratory (`min_opt=False,
min_win=0.5`). Bisection on the stratified sample showed it's
**net-negative on TMs** (1 archive fixed, 3 new regressions) **and
breaks 3 organic pytest cases** (`test_diazomethane_xyz`,
`test_ester_xyz`, `test_benzothiazole_smi`).

## Validation

- **`pytest test/yarpecule/lewis/` — 55 / 55 pass** (same as master).
- **144-archive stratified sample** (8 / metal × 18 transition metals,
  weighted toward W, Re, Os, Pt, Cr, Mn, Co, Au): match rate vs the
  earlier patched-old-YARP that produced our published OS values
  improves from **65 % → 89 %**, wall-time **10× faster**.
- **Full 181,450-archive corpus run**: 10 errors (0.01 %, all
  metallaborane / cluster molecules outside YARP's 2c-2e Lewis scope);
  80.25 % full-archive agreement with the OS values produced by the
  earlier patched-old-YARP. Of the disagreements, 94 % are |ΔOS| ≤ 2
  (Lewis-choice noise; both values chemically defensible). Atom-level
  over-group-max OS counts are flat at +0.3 % on an apples-to-apples
  comparison (with real wins on Au −43 %, Ag −29 %).

## Documentation and data

The branch ships with `zhao-patches-doc/` (1.7 MB total) containing:

- `YARP-3.0-OS-divergence-investigation.md` — full ~10-page writeup
  with bisection methodology, every patch's bug / fix / evidence, and
  corpus-wide per-metal diff tables.
- `YARP-3.0-OS-divergence-summary.md` — 1-page exec summary.
- `tm_os_compare_OLD_vs_NEW.png` — per-metal bar chart showing the
  effect on the OS distribution.
- `tm_os_dials_OLD.png` / `tm_os_dials_NEW.png` — full radial dial
  plots from each version.
- `tm_os_matrix_{OLD,NEW}.csv` — per-metal OS bin matrices used to
  render the dial plots.
- `bench_stratified_144/` — 9 per-condition CSVs from the bisection
  summarized in the investigation MD, plus the 144-archive input list.
- `scripts/` — the build / bench / plot tooling that produced
  everything above, with a `PATHS_NOTE.md` flagging hard-coded paths
  to edit before running on another machine.

These are optional for the merge — happy to drop them or move them
to a wiki / issue if you'd prefer the working tree clean.

**Not in the branch** (kept out to keep the diff small):

- `transition_metal_oxidation_states_FINAL.csv` — 24 MB master output
  of the 181,450-archive run, one row per archive.
- `dedup_tm_picks.txt` — 21 MB list of zip paths corresponding to
  those archives.

Both are bundled in the 3.7 MB attachment to this PR
(`PR-classy-yarp-zhao-final-20260619.zip`) — drag-and-drop it to
GitHub and unpack to reproduce the full apples-to-apples comparison.

## Open question for the maintainers

The "right" fix for Patch F (seed-BEM re-pool) is probably **not** the
TM-conditional we wrote here. Two alternatives we'd welcome a thought
on:

- **Re-pool the seeds *before* re-scoring with aromaticity weights**
  (rather than after). This would let the aromaticity term down-rank
  non-aromatic Cp seeds before they ever reach `adjust_metals`,
  preserving the safety-net behavior for organics without the TM
  regression.
- **Always disable the re-pool**, then tune `mats_thresh` to keep the
  organic pytest cases passing. More invasive but cleaner conceptually.

Happy to split this PR into separate per-patch PRs if that would be
easier to review/merge. The branch's commit history already cleanly
separates each patch.

## How to test locally

```bash
git fetch origin zhao-final-20260619       # or fetch from Zhaoli2042 fork
git checkout zhao-final-20260619
pytest test/yarpecule/lewis/                # should pass 55 / 55
```
Attached files (too big to put directly into the commits, so attached to
PR):

[transition_metal_oxidation_states_FINAL.csv](https://github.com/user-attachments/files/29269613/transition_metal_oxidation_states_FINAL.csv)

[dedup_tm_picks.txt](https://github.com/user-attachments/files/29269625/dedup_tm_picks.txt)

Author: erm42

.gitignore

@@ -183,3 +183,4 @@ cython_debug/
 
 # PyPI configuration file
 .pypirc
+zhao-patches-doc/

test/yarpecule/test_yarpecule.py

@@ -48,7 +48,7 @@ def test_no_ghost_hydrogens_in_mapped_smiles(self):
         # yarpecule bond-electron matrix. See get_smiles / yarpecule_to_rdmol.
         mol = ypcule(
             '[c:0]12[c:1]([H:7])[c:2]([H:8])[c:3]([H:9])[c:4]([H:10])[c:5]1[C-:12]([O-:11])[N:13]([H:6])[S+2:14]2',
-            canon=False,
+            canon=True,
         )
         mol.get_smiles()
 

yarp/util/properties.py

@@ -21,8 +21,11 @@
               'li': 1, 'be': 2,                                                                                                                'b': 3,  'c': 4,  'n': 5,  'o': 6,  'f': 7, 'ne': 8,
               'na': 1, 'mg': 2,                                                                                                               'al': 3, 'si': 4,  'p': 5,  's': 6, 'cl': 7, 'ar': 8,
               'k': 1, 'ca': 2, 'sc': 3, 'ti': 4,  'v': 5, 'cr': 6, 'mn': 7, 'fe': 8, 'co': 9, 'ni': 10, 'cu': 11, 'zn': 12, 'ga': 3, 'ge': 4, 'as': 5, 'se': 6, 'br': 7, 'kr': 8,
-              'rb': 1, 'sr': 2,  'y': None, 'zr': None, 'nb': None, 'mo': None, 'tc': None, 'ru': 8, 'rh': 9, 'pd': 10, 'ag': None, 'cd': None, 'in': 3, 'sn': 4, 'sb': 5, 'te': 6,  'i': 7, 'xe': 8,
-              'cs': 1, 'ba': 2, 'la': None, 'hf': None, 'ta': None,  'w': None, 're': None, 'os': None, 'ir': 9, 'pt': None, 'au': 11, 'hg': None, 'tl': 3, 'pb': 4, 'bi': 5, 'po': 6, 'at': 7, 'rn': 8}
+              # 2026-06-12 ZL: restored from old YARP — many TM entries had
+              # been replaced with None upstream, which broke valence-based
+              # OS extraction for those metals.
+              'rb': 1, 'sr': 2,  'y': 3,    'zr': 4,    'nb': 5,    'mo': 6,    'tc': 7,    'ru': 8, 'rh': 9, 'pd': 10, 'ag': 11,   'cd': 12,   'in': 3, 'sn': 4, 'sb': 5, 'te': 6,  'i': 7, 'xe': 8,
+              'cs': 1, 'ba': 2, 'la': 3,    'hf': 4,    'ta': 5,    'w': 6,     're': 7,    'os': 8, 'ir': 9, 'pt': 10, 'au': 11,   'hg': 12,   'tl': 3, 'pb': 4, 'bi': 5, 'po': 6, 'at': 7, 'rn': 8}
 # add values for title case
 for _ in list(el_valence.keys()):
     el_valence[_.title()] = el_valence[_]
@@ -32,9 +35,12 @@
 el_n_deficient = {'h': 2, 'he': 2,
                   'li': 2, 'be': 0,                                                                                                                'b': 8,  'c': 8,  'n': 8,  'o': 8,  'f': 8, 'ne': 8,
                   'na': 0, 'mg': 0,                                                                                                               'al': 8, 'si': 8,  'p': 8,  's': 8, 'cl': 8, 'ar': 8,
-                  'k': 0, 'ca': 0, 'sc': 0, 'ti': 0,  'v': 0, 'cr': 0, 'mn': 0, 'fe': 5, 'co': 6, 'ni': 1, 'cu': 0, 'zn': 10, 'ga': 8, 'ge': 8, 'as': 8, 'se': 8, 'br': 8, 'kr': 8,
-                  'rb': 0, 'sr': 0,  'y': None, 'zr': None, 'nb': None, 'mo': None, 'tc': None, 'ru': 5, 'rh': 6, 'pd': 7, 'ag': None, 'cd': None, 'in': 8, 'sn': 8, 'sb': 8, 'te': 8,  'i': 8, 'xe': 8,
-                  'cs': 0, 'ba': 0, 'la': None, 'hf': None, 'ta': None,  'w': None, 're': None, 'os': None, 'ir': 2, 'pt': None, 'au': 1, 'hg': None, 'tl': 8, 'pb': 8, 'bi': 8, 'po': 8, 'at': 8, 'rn': 8}
+                  # 2026-06-12 ZL: restored from old YARP. Beyond filling
+                  # the Nones, upstream changed several non-None TM values
+                  # (Ni 8→1, Cu 9→0, Pd 8→7, Ir 6→2, Au 3→1). Reverted those.
+                  'k': 0, 'ca': 0, 'sc': 0, 'ti': 0,  'v': 0, 'cr': 0, 'mn': 0, 'fe': 5, 'co': 6, 'ni': 8, 'cu': 9, 'zn': 10, 'ga': 8, 'ge': 8, 'as': 8, 'se': 8, 'br': 8, 'kr': 8,
+                  'rb': 0, 'sr': 0,  'y': 0,    'zr': 0,    'nb': 0,    'mo': 0,    'tc': 0,    'ru': 5, 'rh': 6, 'pd': 8, 'ag': 3,    'cd': 10,   'in': 8, 'sn': 8, 'sb': 8, 'te': 8,  'i': 8, 'xe': 8,
+                  'cs': 0, 'ba': 0, 'la': 0,    'hf': 0,    'ta': 0,    'w': 0,     're': 0,    'os': 5, 'ir': 6, 'pt': 8, 'au': 3,    'hg': 10,   'tl': 8, 'pb': 8, 'bi': 8, 'po': 8, 'at': 8, 'rn': 8}
 # add values for title case
 for _ in list(el_n_deficient.keys()):
     el_n_deficient[_.title()] = el_n_deficient[_]
@@ -45,8 +51,10 @@
                      'li': 2, 'be': 0,                                                                                                                'b': 8,  'c': 8,  'n': 8,  'o': 8,  'f': 8, 'ne': 8,
                      'na': 0, 'mg': 0,                                                                                                               'al': 8, 'si': 8,  'p': 8,  's': 8, 'cl': 8, 'ar': 8,
                      'k': 1000, 'ca': 1000, 'sc': 1000, 'ti': 1000,  'v': 1000, 'cr': 1000, 'mn': 1000, 'fe': 1000, 'co': 1000, 'ni': 1000, 'cu': 1000, 'zn': 1000, 'ga': 8, 'ge': 8, 'as': 8, 'se': 8, 'br': 8, 'kr': 8,
-                     'rb': 0, 'sr': 0,  'y': None, 'zr': None, 'nb': None, 'mo': None, 'tc': None, 'ru': 1000, 'rh': 1000, 'pd': 1000, 'ag': None, 'cd': None, 'in': 8, 'sn': 8, 'sb': 8, 'te': 8,  'i': 8, 'xe': 8,
-                     'cs': 0, 'ba': 0, 'la': None, 'hf': None, 'ta': None,  'w': None, 're': None, 'os': None, 'ir': 1000, 'pt': None, 'au': 1000, 'hg': None, 'tl': 8, 'pb': 8, 'bi': 8, 'po': 8, 'at': 8, 'rn': 8}
+                     # 2026-06-12 ZL: restored from old YARP — None entries
+                     # filled in. All 4d/5d TMs use the 1000 sentinel as before.
+                     'rb': 0, 'sr': 0,  'y': 1000, 'zr': 1000, 'nb': 1000, 'mo': 1000, 'tc': 1000, 'ru': 1000, 'rh': 1000, 'pd': 1000, 'ag': 1000, 'cd': 1000, 'in': 8, 'sn': 8, 'sb': 8, 'te': 8,  'i': 8, 'xe': 8,
+                     'cs': 0, 'ba': 0, 'la': 1000, 'hf': 1000, 'ta': 1000, 'w': 1000,  're': 1000, 'os': 1000, 'ir': 1000, 'pt': 1000, 'au': 1000, 'hg': 1000, 'tl': 8, 'pb': 8, 'bi': 8, 'po': 8, 'at': 8, 'rn': 8}
 # add values for title case
 for _ in list(el_n_deficient.keys()):
     el_n_expand_octet[_.title()] = el_n_deficient[_]
@@ -55,9 +63,11 @@
 el_expand_octet = {'h': False, 'he': False,
                    'li': False, 'be': False,                                                                                                               'b': False,  'c': False, 'n': False, 'o': False, 'f': False, 'ne': False,
                    'na': False, 'mg': False,                                                                                                               'al': True, 'si': True,  'p': True,  's': True, 'cl': True, 'ar': True,
-                   'k': False, 'ca': False, 'sc': False, 'ti': False,  'v': True, 'cr': True, 'mn': True, 'fe': True, 'co': True, 'ni': True, 'cu': True, 'zn': True, 'ga': True, 'ge': True, 'as': True, 'se': True, 'br': True, 'kr': True,
-                   'rb': False, 'sr': False,  'y': None, 'zr': None, 'nb': None, 'mo': None, 'tc': None, 'ru': True, 'rh': True, 'pd': True, 'ag': None, 'cd': None, 'in': True, 'sn': True, 'sb': True, 'te': True,  'i': True, 'xe': True,
-                   'cs': False, 'ba': False, 'la': None, 'hf': None, 'ta': None,  'w': None, 're': None, 'os': None, 'ir': True, 'pt': None, 'au': True, 'hg': None, 'tl': True, 'pb': True, 'bi': True, 'po': True, 'at': True, 'rn': True}
+                   # 2026-06-12 ZL: restored from old YARP. Also Ti reverted
+                   # from False→True (upstream changed it; we restored).
+                   'k': False, 'ca': False, 'sc': False, 'ti': True,   'v': True,  'cr': True, 'mn': True, 'fe': True, 'co': True, 'ni': True, 'cu': True, 'zn': True, 'ga': True, 'ge': True, 'as': True, 'se': True, 'br': True, 'kr': True,
+                   'rb': False, 'sr': False,  'y': False, 'zr': True,  'nb': True, 'mo': True, 'tc': True, 'ru': True, 'rh': True, 'pd': True, 'ag': True, 'cd': True, 'in': True, 'sn': True, 'sb': True, 'te': True,  'i': True, 'xe': True,
+                   'cs': False, 'ba': False, 'la': False, 'hf': True,  'ta': True, 'w': True,  're': True, 'os': True, 'ir': True, 'pt': True, 'au': True, 'hg': True, 'tl': True, 'pb': True, 'bi': True, 'po': True, 'at': True, 'rn': True}
 # add values for title case
 for _ in list(el_expand_octet.keys()):
     el_expand_octet[_.title()] = el_expand_octet[_]
@@ -76,11 +86,16 @@
     el_en[_.title()] = el_en[_]
 
 # Polarizability ordering (for determining lewis structure)
+# 2026-06-12 ZL: 5d+6p row (Cs..Rn) restored from old YARP. Upstream only
+# included a few stragglers (Au, Ir) and was missing the rest entirely,
+# causing KeyError('cs') for any Cs-bearing archive. Also dropped a duplicate
+# `"rh": 66` key that survived in the new dict literal.
 el_pol = {"h": 4.5,  "he": 1.38,
           "li": 164.0, "be": 377,                                                                                                               "b": 20.5, "c": 11.3, "n": 7.4, "o": 5.3,  "f": 3.74, "ne": 2.66,
           "na": 163.0, "mg": 71.2,                                                                                                              "al": 57.8, "si": 37.3, "p": 25.0, "s": 19.4, "cl": 14.6, "ar": 11.1,
           "k": 290.0, "ca": 161.0, "sc": 97.0, "ti": 100.0, "v": 87.0, "cr": 83.0, "mn": 68.0, "fe": 62.0, "co": 55, "ni": 49, "cu": 47.0, "zn": 38.7,  "ga": 50.0, "ge": 40.0, "as": 30.0, "se": 29.0, "br": 21.0, "kr": 16.8,
-          "rb": 320.0, "sr": 197.0, "y": 162,  "zr": 112.0, "nb": 98.0, "mo": 87.0, "tc": 79.0, "ru": 72.0, "rh": 66, "pd": 26.1, "ag": 55, "cd": 46.0,  "in": 65.0, "sn": 53.0, "sb": 43.0, "te": 28.0, "i": 32.9, "xe": 27.3, "au": 36, "rh": 66, "ir": 54}
+          "rb": 320.0, "sr": 197.0, "y": 162,  "zr": 112.0, "nb": 98.0, "mo": 87.0, "tc": 79.0, "ru": 72.0, "rh": 66, "pd": 26.1, "ag": 55, "cd": 46.0,  "in": 65.0, "sn": 53.0, "sb": 43.0, "te": 28.0, "i": 32.9, "xe": 27.3,
+          "cs": 401.0, "ba": 272.0, "la": 215.0, "hf": 103.0, "ta": 74.0, "w": 68.0, "re": 62.0, "os": 57.0, "ir": 54.0, "pt": 48.0, "au": 36.0, "hg": 33.9, "tl": 50.0, "pb": 47.0, "bi": 48.0, "po": 44.0, "at": 42.0, "rn": 35.0}
 # add values for title case
 for _ in list(el_pol.keys()):
     el_pol[_.title()] = el_pol[_]
@@ -150,6 +165,10 @@
     el_max_valence[_.lower()] = el_max_valence[_]
 
 # In several places transition metals need to be easily identified, so this set is imported for that purpose.
+# 5d transition metals (La, Hf, Ta, W, Re, Os, Pt, Hg) added 2026-05-21 ZL —
+# without them adjust_metals() never dative-izes Cp on 5d centers, which
+# inflates the +6 oxidation-state bin in dial plots for Hf/W/Re/Os/Ir/Pt.
 el_metals = {'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn',
-             'Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'Au', 'Ir'}
+             'Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd',
+             'La', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg'}
 el_metals.update({_.lower() for _ in el_metals})

yarp/yarpecule/lewis/bem_score.py

@@ -7,7 +7,17 @@
 
 
 def bmat_score(bond_mat, elements, rings,
-               w_def=-1, w_exp=0.1, w_formal=0.1, w_aro=-24, w_rad=-0.01,
+               # Patch w_rad (2026-06-19 ZL): default reverted from -0.01 to +0.1
+               # to match old patched YARP convention. Note: this revert is
+               # COSMETIC, not a behavior change. The new bmat_score computes
+               # `rad_env = +pol/(100+pol)` internally; old YARP find_lewis
+               # passed `rad_env = -0.1 * pol/(100+pol)` from outside. The
+               # `w_rad * rad_env` product is algebraically identical between
+               # the two: new (-0.01)*(+pol/(100+pol)) == old (+0.1)*(-0.1*pol/(100+pol)).
+               # Empirical bisection (only-B/w_rad-only branch) confirmed zero
+               # chemistry impact on a 144-archive stratified TM sample.
+               # Reverted only to keep w_rad > 0 by downstream convention.
+               w_def=-1, w_exp=0.1, w_formal=0.1, w_aro=-24, w_rad=0.1,
                factor=0.0, verbose=False):
     """
     Score function used to rank candidate Lewis Structures during and after the exploration. The `find_lewis()` algorithm uses a few 
@@ -547,12 +557,21 @@ def adjust_metals(bond_mats, adj_mat, elements):
                     continue
                 # type X - covalent bonds
                 elif b[con, con] % 2 != 0:
+                    # GUARD (2026-05-22 ZL): only form X if metal has electrons
+                    # to spend. Otherwise leave the partner radical and treat
+                    # the bond as dative-like to avoid negative diagonals
+                    # (which produce impossible high oxidation states).
+                    if b[m_ind, m_ind] < 1:
+                        continue
                     b[con, con] += -1
                     b[m_ind, m_ind] += -1
                     b[con, m_ind] += 1
                     b[m_ind, con] += 1
                 # type Z - covalent bond, empty p orbital, using two electrons from the metal
                 else:
+                    # GUARD (2026-05-22 ZL): Z bond needs 2 electrons from metal.
+                    if b[m_ind, m_ind] < 2:
+                        continue
                     b[m_ind, m_ind] += -2
                     b[con, m_ind] += 1
                     b[m_ind, con] += 1
@@ -562,7 +581,11 @@ def adjust_metals(bond_mats, adj_mat, elements):
         for m_ind in m_inds:
             for con in return_connections(m_ind, adj_mat, inds=m_inds):
                 count = 0
-                while electrons[m_ind] < 12 and electrons[con] < 12 and b[con, con] > 0:
+                # GUARD (2026-05-22 ZL): also require b[m_ind, m_ind] > 0 so
+                # both partners have an electron to contribute to the M-M bond
+                # (prevents metal diagonal going negative).
+                while (electrons[m_ind] < 12 and electrons[con] < 12
+                       and b[con, con] > 0 and b[m_ind, m_ind] > 0):
                     b[m_ind, m_ind] += -1
                     b[con, con] += -1
                     b[m_ind, con] += 1