Transition metal oxidation state patch

yarp/util/properties.py

@@ -150,6 +150,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

@@ -547,12 +547,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 +571,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

yarp/yarpecule/lewis/find_lewis.py

@@ -200,9 +200,13 @@ def gen_init(obj_fun, adj_mat, elements, rings, q):
         yield obj_fun(bond_mat), bond_mat, reactive
 
 
+# Defaults rolled back 2026-06-12 ZL: upstream bumped N_score=100 → 1000 and
+# counter=0 → 100 (which by itself would break immediately if N_score is also
+# 100). The production OS recalculation pipeline ran with N_score=100, so the
+# old behavior is restored here.
 def gen_all_lstructs(obj_fun, bond_mats, scores, hashes, elements,
                      reactive, rings, ring_atoms, bridgeheads, seps, min_score,
-                     ind=0, counter=100, N_score=1000, N_max=10000, min_opt=False, min_win=False):
+                     ind=0, counter=0, N_score=100, N_max=10000, min_opt=False, min_win=False):
     """ 
     A generator for Lewis search algorithm that recursively applies a set of valid bond-electron moves to find all relevant resonance structures. 
 

yarp/yarpecule/lewis/lewis_structure.py

@@ -255,12 +255,14 @@ def obj_fun(x): return bmat_score(x, elements, self._rings,
                 seed_scores += [score]
                 seed_bond_mats += [bond_mat]
                 seed_hashes.add(bmat_hash(bond_mat))
+                # N_score=100 rolled back 2026-06-12 ZL — matches production
+                # OS-recalculation run that fed the published dial plot.
                 seed_bond_mats, seed_scores, _, _, _ = gen_all_lstructs(obj_fun, seed_bond_mats, seed_scores, seed_hashes,
                                                                         elements, reactive, self._rings, ring_atoms, bridgeheads,
                                                                         # allow all charge transfers in first pass
                                                                         seps=np.zeros([len(elements), len(elements)]),
                                                                         min_score=seed_scores[0], ind=len(seed_bond_mats)-1,
-                                                                        N_score=1000, N_max=10000, min_opt=True)
+                                                                        N_score=100, N_max=10000, min_opt=True)
 
         # Update objective function to include (anti)aromaticity considerations
         def obj_fun(x): return bmat_score(x, elements, self._rings,
@@ -279,13 +281,14 @@ def obj_fun(x): return bmat_score(x, elements, self._rings,
         hashes = set([bmat_hash(seed_bond_mats[0])])
 
         # Next round of BEM searching
+        # N_score=100 rolled back 2026-06-12 ZL (see above).
         bond_mats, scores, hashes, _, _ = gen_all_lstructs(obj_fun, bond_mats, scores,
                                                            hashes, elements, reactive,
                                                            self._rings, ring_atoms, bridgeheads,
                                                            # set according to local_opt flag
                                                            seps,
                                                            min_score=min(scores), ind=len(bond_mats)-1,
-                                                           N_score=1000, N_max=10000, min_opt=True)
+                                                           N_score=100, N_max=10000, min_opt=True)
 
         # Collect all discovered BEMs
         for i, bem in enumerate(seed_bond_mats):