"""CTfile multiple groups: explicit copies, parent atoms and multiplier.
MDL CTfile Formats (December 1999), pp. 2-15/16 requires every copy in the
atom list. KET draws the parent once. Nested repeats expand inside out;
atom and bond stereo are retained.
Crossing groups without a nesting order are not guessed.
"""
from __future__ import annotations
import itertools
from typing import Any
[docs]
class MultipleGroupRefused(ValueError):
"""The bracket cannot be faithfully expressed as a CTfile multiple group."""
def _groups(mol: Any) -> list[dict[str, Any]]:
from rdkit import Chem
return [
{
"atoms": list(g.GetAtoms()),
"parent": list(g.GetParentAtoms()),
"props": {
k: g.GetProp(k)
for k in g.GetPropNames()
if k not in ("index", "DATAFIELDS", "PARENT")
},
}
for g in Chem.GetMolSubstanceGroups(mol)
]
def _write_groups(mol: Any, groups: list[dict[str, Any]]) -> None:
from rdkit import Chem
Chem.ClearMolSubstanceGroups(mol)
for info in groups:
group = Chem.CreateMolSubstanceGroup(mol, info["props"]["TYPE"])
for key, value in info["props"].items():
group.SetProp(key, value)
group.SetAtoms(info["atoms"])
group.SetParentAtoms(info["parent"])
atoms = set(info["atoms"])
group.SetBonds(
[
b.GetIdx()
for b in mol.GetBonds()
if (b.GetBeginAtomIdx() in atoms) != (b.GetEndAtomIdx() in atoms)
]
)
def _orient(mol: Any, desired: dict[int, list[int]]) -> None:
"""Preserve tetrahedral orientation after graph edits change neighbor order."""
from rdkit import Chem
for index, neighbors in desired.items():
atom = mol.GetAtomWithIdx(index)
if atom.GetChiralTag() == Chem.ChiralType.CHI_UNSPECIFIED:
continue
actual = [a.GetIdx() for a in atom.GetNeighbors()]
if set(actual) != set(neighbors):
raise MultipleGroupRefused("repeat connectivity does not preserve a stereo center")
order = [actual.index(a) for a in neighbors]
if sum(a > b for i, a in enumerate(order) for b in order[i + 1 :]) % 2:
atom.InvertChirality()
[docs]
def expand_multiple_groups(mol: Any) -> Any:
"""Expand inner repeats before outer ones; retain atom and bond stereo."""
from rdkit import Chem
from rdkit.Chem import rdDepictor
groups = _groups(mol)
if not any(g["props"]["TYPE"] == "MUL" and not g["parent"] for g in groups):
return mol
result = Chem.RWMol(mol)
Chem.ClearMolSubstanceGroups(result)
for i, group in enumerate(groups):
first = set(group["atoms"])
for other in groups[i + 1 :]:
second = set(other["atoms"])
if first & second and not (first < second or second < first):
raise MultipleGroupRefused("overlapping groups without a strict nesting order")
while True:
pending = [g for g in groups if g["props"]["TYPE"] == "MUL" and not g["parent"]]
if not pending:
break
group = min(pending, key=lambda g: len(g["atoms"]))
parent = list(group["atoms"])
inside = set(parent)
count = int(group["props"].get("MULT", group["props"].get("LABEL", "1")))
if count < 1 or result.GetNumAtoms() + (count - 1) * len(parent) > 10_000:
raise MultipleGroupRefused(
"expanded multiple group exceeds 10,000 atoms or has invalid count"
)
source = result.GetMol()
internal = [
b
for b in source.GetBonds()
if b.GetBeginAtomIdx() in inside and b.GetEndAtomIdx() in inside
]
crossings = [
b
for b in source.GetBonds()
if (b.GetBeginAtomIdx() in inside) != (b.GetEndAtomIdx() in inside)
]
if len(crossings) not in (0, 2):
raise MultipleGroupRefused("a multiple group needs zero or two crossing bonds")
if crossings and crossings[0].GetBondType() != crossings[1].GetBondType():
raise MultipleGroupRefused("multiple-group crossing bond orders differ")
reachable = {parent[0]}
while True:
more = reachable | {
a
for b in internal
if b.GetBeginAtomIdx() in reachable or b.GetEndAtomIdx() in reachable
for a in (b.GetBeginAtomIdx(), b.GetEndAtomIdx())
}
if more == reachable:
break
reachable = more
if reachable != inside:
raise MultipleGroupRefused("the multiple-group parent must be connected")
copies = [{a: a for a in parent}]
desired = {a.GetIdx(): [n.GetIdx() for n in a.GetNeighbors()] for a in source.GetAtoms()}
stereo = [
(b.GetBeginAtomIdx(), b.GetEndAtomIdx(), b.GetStereo(), list(b.GetStereoAtoms()))
for b in source.GetBonds()
if b.GetStereoAtoms()
]
for _ in range(1, count):
mapping = {a: result.AddAtom(Chem.Atom(source.GetAtomWithIdx(a))) for a in parent}
for a in parent:
desired[mapping[a]] = [
mapping.get(n.GetIdx(), n.GetIdx())
for n in source.GetAtomWithIdx(a).GetNeighbors()
]
for bond in internal:
left, right = mapping[bond.GetBeginAtomIdx()], mapping[bond.GetEndAtomIdx()]
result.AddBond(left, right, bond.GetBondType())
result.GetBondBetweenAtoms(left, right).SetBondDir(bond.GetBondDir())
if bond.GetStereoAtoms():
stereo.append(
(
left,
right,
bond.GetStereo(),
[mapping.get(a, a) for a in bond.GetStereoAtoms()],
)
)
copies.append(mapping)
replacements = {}
if crossings and count > 1:
ends = [
(b.GetBeginAtomIdx(), b.GetEndAtomIdx())
if b.GetBeginAtomIdx() in inside
else (b.GetEndAtomIdx(), b.GetBeginAtomIdx())
for b in crossings
]
(head, before), (tail, after) = ends
result.RemoveBond(tail, after)
for i, mapping in enumerate(copies):
for atom, old, new in [
(mapping[head], before, copies[i - 1][tail] if i else before),
(mapping[tail], after, copies[i + 1][head] if i + 1 < count else after),
]:
desired[atom] = [new if n == old else n for n in desired[atom]]
replacements[(atom, old)] = new
desired[after] = [copies[-1][tail] if n == tail else n for n in desired[after]]
replacements[(after, tail)] = copies[-1][tail]
for previous, following in zip(copies, copies[1:], strict=False):
result.AddBond(previous[tail], following[head], crossings[0].GetBondType())
result.AddBond(copies[-1][tail], after, crossings[1].GetBondType())
_orient(result, desired)
for left, right, kind, refs in stereo:
bond = result.GetBondBetweenAtoms(left, right)
if bond is None:
raise MultipleGroupRefused(
"a stereogenic crossing bond requires explicit linkage geometry"
)
bond.SetStereoAtoms(
replacements.get((left, refs[0]), refs[0]),
replacements.get((right, refs[1]), refs[1]),
)
bond.SetStereo(kind)
added = [mapping[a] for mapping in copies[1:] for a in parent]
children = [g for g in groups if set(g["atoms"]) < inside]
for other in groups:
if inside < set(other["atoms"]):
other["atoms"].extend(added)
for mapping in copies[1:]:
groups.extend(
{
"atoms": [mapping[a] for a in child["atoms"]],
"parent": [mapping[a] for a in child["parent"]],
"props": dict(child["props"]),
}
for child in children
)
group["parent"] = parent
group["atoms"].extend(added)
value = result.GetMol()
Chem.SanitizeMol(value)
for bond in value.GetBonds():
if bond.GetBondDir() in (Chem.BondDir.ENDUPRIGHT, Chem.BondDir.ENDDOWNRIGHT):
bond.SetBondDir(Chem.BondDir.NONE)
Chem.AssignStereochemistry(value, cleanIt=False, force=True)
Chem.SetDoubleBondNeighborDirections(value)
_write_groups(value, groups)
value.RemoveAllConformers()
rdDepictor.Compute2DCoords(value)
return value
[docs]
def collapse_multiple_groups(mol: Any) -> Any:
"""Verify each outer parent by stereo-aware, mapped-graph re-expansion."""
from rdkit import Chem
mol = Chem.Mol(mol)
Chem.AssignStereochemistry(mol, cleanIt=True, force=True)
groups = _groups(mol)
pending = [g for g in groups if g["props"]["TYPE"] == "MUL" and g["parent"]]
if not pending:
return mol
group = max(pending, key=lambda g: len(g["atoms"]))
parent, inside = group["parent"], set(group["atoms"])
try:
count = int(group["props"].get("MULT", group["props"].get("LABEL", "1")))
except ValueError:
return mol
if not set(parent) <= inside or len(inside) != len(parent) * count or count < 1:
return mol
remove = inside - set(parent)
edges = [
(b.GetBeginAtomIdx(), b.GetEndAtomIdx(), b.GetBondType())
for b in mol.GetBonds()
if (b.GetBeginAtomIdx() in inside) != (b.GetEndAtomIdx() in inside)
]
if len(edges) not in (0, 2):
return mol
kept = [i for i in range(mol.GetNumAtoms()) if i not in remove]
indices = {old: new for new, old in enumerate(kept)}
original = Chem.Mol(mol)
for atom in original.GetAtoms():
atom.SetAtomMapNum(atom.GetIdx() + 1 if atom.GetIdx() in indices else 0)
target = Chem.MolToSmiles(original)
external = [(b if a in inside else a, a if a in inside else b, order) for a, b, order in edges]
changed = [edge for edge in external if edge[1] in remove]
if len(parent) ** len(changed) > 10_000:
return mol # foreign record verification is bounded, never guessed
for replacements in itertools.product(parent, repeat=len(changed)):
candidate = Chem.RWMol(mol)
Chem.ClearMolSubstanceGroups(candidate)
desired = {a: [n.GetIdx() for n in mol.GetAtomWithIdx(a).GetNeighbors()] for a in kept}
relink = {}
valid = True
for (outside, old, order), replacement in zip(changed, replacements, strict=True):
if candidate.GetBondBetweenAtoms(outside, replacement) is not None:
valid = False
break
candidate.RemoveBond(outside, old)
candidate.AddBond(outside, replacement, order)
desired[outside] = [replacement if n == old else n for n in desired[outside]]
relink[(outside, old)] = replacement
# The deleted copy chain identifies which parent edge this restores.
reachable = {old}
while True:
more = reachable | {
n.GetIdx()
for a in reachable
for n in mol.GetAtomWithIdx(a).GetNeighbors()
if n.GetIdx() in remove
}
if more == reachable:
break
reachable = more
lost = [n for n in desired[replacement] if n in reachable]
if len(lost) != 1:
valid = False
break
desired[replacement] = [outside if n == lost[0] else n for n in desired[replacement]]
relink[(replacement, lost[0])] = outside
if not valid:
continue
stereo = []
for bond in mol.GetBonds():
left, right = bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()
refs = list(bond.GetStereoAtoms())
if refs and left in indices and right in indices:
stereo.append(
(
left,
right,
bond.GetStereo(),
[
relink.get((left, refs[0]), refs[0]),
relink.get((right, refs[1]), refs[1]),
],
)
)
for atom in sorted(remove, reverse=True):
candidate.RemoveAtom(atom)
try:
_orient(
candidate,
{indices[a]: [indices[n] for n in neighbors] for a, neighbors in desired.items()},
)
for left, right, kind, refs in stereo:
bond = candidate.GetBondBetweenAtoms(indices[left], indices[right])
bond.SetStereoAtoms(indices[refs[0]], indices[refs[1]])
bond.SetStereo(kind)
small = candidate.GetMol()
Chem.SanitizeMol(small)
chosen: dict[str, Any] = {
"atoms": [indices[a] for a in parent],
"parent": [],
"props": dict(group["props"]),
}
chosen["props"]["LABEL"] = str(count)
_write_groups(small, [chosen])
for old, new in indices.items():
small.GetAtomWithIdx(new).SetAtomMapNum(old + 1)
expanded = expand_multiple_groups(small)
for atom in expanded.GetAtoms():
if atom.GetIdx() >= small.GetNumAtoms():
atom.SetAtomMapNum(0)
Chem.AssignStereochemistry(expanded, cleanIt=True, force=True)
if Chem.MolToSmiles(expanded) != target:
continue
for atom in small.GetAtoms():
atom.SetAtomMapNum(mol.GetAtomWithIdx(kept[atom.GetIdx()]).GetAtomMapNum())
retained = [chosen]
for other in groups:
if other is group:
continue
atoms = set(other["atoms"])
if atoms <= remove:
continue # copied child, restored by outer expansion
if atoms & remove and not inside < atoms:
valid = False
break
retained.append(
{
"atoms": [indices[a] for a in other["atoms"] if a in indices],
"parent": [indices[a] for a in other["parent"] if a in indices],
"props": dict(other["props"]),
}
)
if not valid:
continue
_write_groups(small, retained)
small.RemoveAllConformers()
return collapse_multiple_groups(small)
except (ValueError, RuntimeError, KeyError):
continue
return mol