preserve 3D conformer through TMC conversion

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)

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

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
 
 

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

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

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

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)