|
1 | 1 | """SMILES input and graph component helpers.""" |
2 | 2 |
|
3 | 3 | from rdkit import Chem |
| 4 | +from rdkit.Chem import rdCIPLabeler |
4 | 5 |
|
5 | 6 | from .molecule import Molecule |
6 | 7 |
|
@@ -29,31 +30,53 @@ def read_smiles(smiles: str) -> Molecule: |
29 | 30 | atom_metadata = _atom_metadata(rdmol) |
30 | 31 |
|
31 | 32 | Chem.AssignStereochemistry(rdmol, force=True, cleanIt=True) |
32 | | - chiral_centers = dict(Chem.FindMolChiralCenters(rdmol, includeUnassigned=False)) |
| 33 | + legacy_centers = dict(Chem.FindMolChiralCenters(rdmol, includeUnassigned=False)) |
| 34 | + # The legacy labeler decides *which* centers carry an absolute descriptor: |
| 35 | + # centers it leaves unassigned are ring-symmetry-dependent ones that the |
| 36 | + # scoped small-ring / cis-trans fallbacks must handle (OPSIN rejects |
| 37 | + # absolute R/S there). But its CIP *values* mis-rank ligands in |
| 38 | + # deep-sphere comparisons (e.g. ring-closure duplicate atoms) and invert |
| 39 | + # 3-coordinate sulfur centers, so the label itself comes from |
| 40 | + # rdCIPLabeler, which overwrites the _CIPCode properties. |
| 41 | + try: |
| 42 | + rdCIPLabeler.AssignCIPLabels(rdmol) |
| 43 | + except Exception: |
| 44 | + pass # unsanitized input; keep the legacy labels |
| 45 | + chiral_centers = {} |
| 46 | + for idx, legacy_code in legacy_centers.items(): |
| 47 | + atom = rdmol.GetAtomWithIdx(idx) |
| 48 | + code = atom.GetProp("_CIPCode") if atom.HasProp("_CIPCode") else legacy_code |
| 49 | + if code in ("R", "S"): |
| 50 | + chiral_centers[idx] = code |
33 | 51 |
|
34 | 52 | for atom in rdmol.GetAtoms(): |
35 | 53 | stereo = chiral_centers.get(atom.GetIdx()) |
36 | | - if stereo and atom.GetSymbol() == "S" and atom.GetTotalDegree() == 3: |
37 | | - stereo = "R" if stereo == "S" else "S" |
38 | 54 | raw_stereo = _raw_tetrahedral_stereo(atom) if not stereo else None |
39 | 55 | mol.add_atom( |
40 | 56 | symbol=atom.GetSymbol(), |
41 | 57 | idx=atom.GetIdx(), |
42 | 58 | charge=atom.GetFormalCharge(), |
43 | 59 | stereo=stereo, |
44 | 60 | raw_stereo=raw_stereo, |
| 61 | + cip_code=atom.GetProp("_CIPCode") if atom.HasProp("_CIPCode") else None, |
45 | 62 | is_aromatic=atom_metadata[atom.GetIdx()]["is_aromatic"], |
46 | 63 | explicit_h_count=atom_metadata[atom.GetIdx()]["explicit_h_count"], |
47 | 64 | total_h_count=atom_metadata[atom.GetIdx()]["total_h_count"], |
48 | 65 | ) |
49 | 66 |
|
50 | 67 | for bond in rdmol.GetBonds(): |
51 | | - stereo = None |
52 | | - st = bond.GetStereo() |
53 | | - if st == Chem.rdchem.BondStereo.STEREOE: |
54 | | - stereo = "E" |
55 | | - elif st == Chem.rdchem.BondStereo.STEREOZ: |
56 | | - stereo = "Z" |
| 68 | + # rdCIPLabeler rewrites bond enums to STEREOCIS/STEREOTRANS but stamps |
| 69 | + # the authoritative E/Z label as _CIPCode; fall back to the enums when |
| 70 | + # the labeler did not run. |
| 71 | + stereo = bond.GetPropsAsDict().get("_CIPCode") |
| 72 | + if stereo not in ("E", "Z"): |
| 73 | + st = bond.GetStereo() |
| 74 | + if st == Chem.rdchem.BondStereo.STEREOE: |
| 75 | + stereo = "E" |
| 76 | + elif st == Chem.rdchem.BondStereo.STEREOZ: |
| 77 | + stereo = "Z" |
| 78 | + else: |
| 79 | + stereo = None |
57 | 80 |
|
58 | 81 | in_small_ring = any(bond.IsInRingSize(i) for i in range(3, 8)) |
59 | 82 |
|
|
0 commit comments