diff --git a/steamroll/steamroll.py b/steamroll/steamroll.py index 7c1c478..64a0177 100644 --- a/steamroll/steamroll.py +++ b/steamroll/steamroll.py @@ -12,6 +12,7 @@ from rdkit.Chem import rdDetermineBonds from rdkit.Geometry import Point3D +from .utils import strip_to_connectivity from .xyz2mol.xyz2mol import xyz2mol from .xyz2mol_tmc.xyz2mol_local import xyz2AC_obabel as xyz2ac_obabel from .xyz2mol_tmc.xyz2mol_tmc import TRANSITION_METALS_NUM, get_tmc_mol @@ -138,16 +139,7 @@ def _from_smiles_and_coords( raw.AddConformer(raw_conf, assignId=True) rdDetermineBonds.DetermineConnectivity(raw) - # Strip charges, isotopes, and bond orders to match raw's bare atoms. - query = Chem.RWMol(template) - for atom in query.GetAtoms(): - atom.SetFormalCharge(0) - atom.SetIsotope(0) - for bond in query.GetBonds(): - bond.SetBondType(Chem.BondType.SINGLE) - bond.SetIsAromatic(False) - - match = raw.GetSubstructMatch(query) + match = raw.GetSubstructMatch(strip_to_connectivity(template)) if not match or len(match) != n: raise ValueError("Could not find a valid atom mapping between SMILES and XYZ") @@ -254,7 +246,7 @@ def to_rdkit( os.unlink(xyz_file) if rdkm is None: raise SteamrollConversionError("xyz2mol_tm returned no molecule") - return Chem.AddHs(rdkm) + return remove_hydrogens(rdkm) if remove_Hs else Chem.AddHs(rdkm, addCoords=True) def _topology_ok(mol: Chem.rdchem.Mol) -> bool: return smiles is None or _smiles_matches(mol, smiles) diff --git a/steamroll/utils.py b/steamroll/utils.py new file mode 100644 index 0000000..f5cad95 --- /dev/null +++ b/steamroll/utils.py @@ -0,0 +1,25 @@ +"""Shared RDKit utilities for steamroll.""" + +from rdkit import Chem + + +def strip_to_connectivity(mol: Chem.rdchem.Mol) -> Chem.RWMol: + """Return a copy of mol with all bond orders set to SINGLE and charges/isotopes cleared. + + Useful for connectivity-only substructure matching where bond order or charge + differences would otherwise prevent a match. + + Args: + mol: molecule to strip + + Returns: + RWMol with all bonds set to SINGLE, charges and isotopes cleared + """ + rwmol = Chem.RWMol(mol) + for atom in rwmol.GetAtoms(): + atom.SetFormalCharge(0) + atom.SetIsotope(0) + for bond in rwmol.GetBonds(): + bond.SetBondType(Chem.BondType.SINGLE) + bond.SetIsAromatic(False) + return rwmol diff --git a/steamroll/xyz2mol_tmc/xyz2mol_tmc.py b/steamroll/xyz2mol_tmc/xyz2mol_tmc.py index ad30926..0264ae1 100644 --- a/steamroll/xyz2mol_tmc/xyz2mol_tmc.py +++ b/steamroll/xyz2mol_tmc/xyz2mol_tmc.py @@ -19,6 +19,7 @@ from rdkit.Chem import GetPeriodicTable, rdchem, rdEHTTools, rdmolops from rdkit.Chem.MolStandardize import rdMolStandardize +from ..utils import strip_to_connectivity from .xyz2mol_local import ( AC2mol, chiral_stereo_check, @@ -667,6 +668,18 @@ def get_tmc_mol(xyz_file, overall_charge, with_stereo=False): Chem.SanitizeMol(tmc_mol) if with_stereo: chiral_stereo_check(tmc_mol) + + # Re-attach 3D coords lost in the SMILES roundtrip; strip bonds on both sides to match. + emol_mol = emol.GetMol() + if emol_mol.GetNumConformers() > 0: + match = strip_to_connectivity(emol_mol).GetSubstructMatch(strip_to_connectivity(tmc_mol)) + if match: + emol_conf = emol_mol.GetConformer() + conf = Chem.Conformer(tmc_mol.GetNumAtoms()) + for i, emol_idx in enumerate(match): + conf.SetAtomPosition(i, emol_conf.GetAtomPosition(emol_idx)) + tmc_mol.AddConformer(conf, assignId=True) + return tmc_mol diff --git a/tests/data/bromobenzene_distorted.xyz b/tests/data/bromobenzene_distorted.xyz new file mode 100644 index 0000000..b603dff --- /dev/null +++ b/tests/data/bromobenzene_distorted.xyz @@ -0,0 +1,14 @@ +12 +charge: 0; smiles: Brc1ccccc1 +Br 0.000 2.700 0.000 +C 0.000 1.400 0.000 +C 1.212 0.700 0.000 +C 1.212 -0.700 0.000 +C 0.000 -1.400 0.000 +C -1.212 -0.700 0.000 +C -1.212 0.700 0.000 +H 2.147 1.240 0.000 +H 2.147 -1.240 0.000 +H 0.000 -2.480 0.000 +H -2.147 -1.240 0.000 +H -2.147 1.240 0.000 diff --git a/tests/data/fe_pyridone_complex.xyz b/tests/data/fe_pyridone_complex.xyz new file mode 100644 index 0000000..c35d352 --- /dev/null +++ b/tests/data/fe_pyridone_complex.xyz @@ -0,0 +1,24 @@ +22 +name: Initial Structure; charge: 0; multiplicity: 1; generated_by: Rowan; timestamp: 2026-03-24 20:30:48; +N 1.24092555 1.77685010 -0.92464799 +C -0.05854858 1.68824053 -0.52946717 +C -0.43100920 0.55965042 0.17266804 +Fe -2.25838900 0.37456831 0.76066321 +C -3.21333218 1.50765657 -0.35333976 +C -2.68190145 -1.35837996 0.41204101 +C 0.43446308 -0.46198651 0.48630235 +C 1.74682081 -0.33479285 0.06707576 +C 2.74360943 -1.41847837 0.37866789 +C 2.10741663 0.80095631 -0.63570112 +O 3.45530796 0.86325085 -1.02980232 +H -0.77100676 2.48022151 -0.76147252 +H -2.74460077 2.10026026 -1.13416147 +H -4.28243208 1.58213091 -0.22053052 +H -3.62964702 -1.52671921 -0.14915533 +H -2.72617674 -1.99033189 1.32893932 +H -1.88709903 -1.84501624 -0.22044833 +H 0.12072585 -1.34531713 1.04180932 +H 2.71787429 -2.18230247 -0.43716633 +H 2.37594128 -1.92637289 1.29662430 +H 3.71285415 -0.93651640 0.49768966 +H 4.02820349 1.59242821 -0.65070683 diff --git a/tests/data/naphthalene.xyz b/tests/data/naphthalene.xyz new file mode 100644 index 0000000..229d998 --- /dev/null +++ b/tests/data/naphthalene.xyz @@ -0,0 +1,20 @@ +18 +charge: 0; smiles: c1ccc2ccccc2c1 +C 0.000 0.718 0.000 +C 0.000 -0.718 0.000 +C 1.230 1.399 0.000 +C 2.459 0.718 0.000 +C 2.459 -0.718 0.000 +C 1.230 -1.399 0.000 +C -1.230 -1.399 0.000 +C -2.459 -0.718 0.000 +C -2.459 0.718 0.000 +C -1.230 1.399 0.000 +H 1.230 2.480 0.000 +H 3.393 1.247 0.000 +H 3.393 -1.247 0.000 +H 1.230 -2.480 0.000 +H -1.230 -2.480 0.000 +H -3.393 -1.247 0.000 +H -3.393 1.247 0.000 +H -1.230 2.480 0.000 diff --git a/tests/test_steamroll.py b/tests/test_steamroll.py index edad50a..eeccf75 100644 --- a/tests/test_steamroll.py +++ b/tests/test_steamroll.py @@ -11,48 +11,7 @@ from steamroll.steamroll import ATOMIC_NUMBERS, SteamrollTopologyMismatchError, fragment, to_rdkit _BROMOBENZENE_SMILES = "Brc1ccccc1" - -# Naphthalene: all-carbon fused ring; a pathological case for centroid-based atom mapping. _NAPHTHALENE_SMILES = "c1ccc2ccccc2c1" -_NAPHTHALENE_ATOMIC_NUMBERS = [6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 1, 1, 1, 1, 1, 1, 1, 1] -_NAPHTHALENE_COORDS = [ - [0.000, 0.718, 0.000], # C4a (junction) - [0.000, -0.718, 0.000], # C8a (junction) - [1.230, 1.399, 0.000], # C1 - [2.459, 0.718, 0.000], # C2 - [2.459, -0.718, 0.000], # C3 - [1.230, -1.399, 0.000], # C4 - [-1.230, -1.399, 0.000], # C5 - [-2.459, -0.718, 0.000], # C6 - [-2.459, 0.718, 0.000], # C7 - [-1.230, 1.399, 0.000], # C8 - [1.230, 2.480, 0.000], # H on C1 - [3.393, 1.247, 0.000], # H on C2 - [3.393, -1.247, 0.000], # H on C3 - [1.230, -2.480, 0.000], # H on C4 - [-1.230, -2.480, 0.000], # H on C5 - [-3.393, -1.247, 0.000], # H on C6 - [-3.393, 1.247, 0.000], # H on C7 - [-1.230, 2.480, 0.000], # H on C8 -] - -# Bromobenzene with C-Br shrunk to 1.3 Å (normal ~1.9 Å); pulls Br into -# bonding range of the ortho carbons, causing DetermineConnectivity to mis-bond. -_BROMOBENZENE_DISTORTED_ATOMIC_NUMBERS = [35, 6, 6, 6, 6, 6, 6, 1, 1, 1, 1, 1] -_BROMOBENZENE_DISTORTED_COORDS = [ - [0.000, 2.700, 0.000], # Br — 1.3 Å from C1 (distorted) - [0.000, 1.400, 0.000], # C1 - [1.212, 0.700, 0.000], # C2 - [1.212, -0.700, 0.000], # C3 - [0.000, -1.400, 0.000], # C4 - [-1.212, -0.700, 0.000], # C5 - [-1.212, 0.700, 0.000], # C6 - [2.147, 1.240, 0.000], # H on C2 - [2.147, -1.240, 0.000], # H on C3 - [0.000, -2.480, 0.000], # H on C4 - [-2.147, -1.240, 0.000], # H on C5 - [-2.147, 1.240, 0.000], # H on C6 -] HERE = Path(__file__).parent @@ -124,8 +83,7 @@ def test_smiles_distorted_halogen() -> None: With SMILES, the correct topology is recovered: Br has exactly 1 bond and all heavy-atom coordinates are preserved. """ - atomic_numbers = _BROMOBENZENE_DISTORTED_ATOMIC_NUMBERS - coordinates = _BROMOBENZENE_DISTORTED_COORDS + atomic_numbers, coordinates, _ = read_xyz(DATA_DIR / "bromobenzene_distorted.xyz") ref_smiles = Chem.MolToSmiles(Chem.MolFromSmiles(_BROMOBENZENE_SMILES), isomericSmiles=False) # Without SMILES: wrong topology @@ -157,13 +115,9 @@ def test_smiles_distorted_halogen() -> None: def test_smiles_fused_ring() -> None: """SMILES-based conversion correctly maps atoms in fused ring systems.""" + atomic_numbers, coordinates, _ = read_xyz(DATA_DIR / "naphthalene.xyz") ref_smiles = Chem.MolToSmiles(Chem.MolFromSmiles(_NAPHTHALENE_SMILES), isomericSmiles=False) - rdkm = to_rdkit( - _NAPHTHALENE_ATOMIC_NUMBERS, - _NAPHTHALENE_COORDS, - smiles=_NAPHTHALENE_SMILES, - remove_Hs=False, - ) + rdkm = to_rdkit(atomic_numbers, coordinates, smiles=_NAPHTHALENE_SMILES, remove_Hs=False) assert Chem.MolToSmiles(Chem.RemoveHs(rdkm), isomericSmiles=False) == ref_smiles # Verify bond lengths match naphthalene geometry; wrong atom assignments @@ -218,3 +172,17 @@ def test_smiles_high_symmetry(smiles: str) -> None: assert Chem.MolToSmiles(Chem.RemoveHs(rdkm), isomericSmiles=False) == Chem.MolToSmiles( ref, isomericSmiles=False ) + + +def test_tmc_conformer_preserved() -> None: + """to_rdkit preserves 3D coordinates for transition metal complexes. + + Previously get_tmc_mol discarded coordinates by roundtripping through SMILES, + returning a mol with no conformer. + """ + atomic_numbers, coordinates, charge = read_xyz(DATA_DIR / "fe_pyridone_complex.xyz") + rdkm = to_rdkit(atomic_numbers, coordinates, charge=charge, remove_Hs=True) + assert rdkm.GetNumConformers() == 1 + conf = rdkm.GetConformer() + positions = [conf.GetAtomPosition(i) for i in range(rdkm.GetNumAtoms())] + assert not all(p.x == 0.0 and p.y == 0.0 and p.z == 0.0 for p in positions)