Skip to content

Commit ad7c1d5

Browse files
authored
Merge pull request #772 from Steinbeck-Lab/development
Add multi-frame XYZ conversion features and refactor types
2 parents 7126421 + 9b85d18 commit ad7c1d5

12 files changed

Lines changed: 1150 additions & 268 deletions

File tree

‎app/modules/toolkits/helpers.py‎

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -70,3 +70,84 @@ def parse_SMILES(smiles: str, framework: str = "rdkit", standardize: bool = Fals
7070
return Chem.MolFromMolBlock(mol_str)
7171
except Exception:
7272
raise InvalidInputException(name="smiles", value=smiles)
73+
74+
75+
def _is_plausible_atom_row(line: str) -> bool:
76+
"""Return True if *line* looks like a valid XYZ atom row.
77+
78+
A plausible row, after stripping, has at least 4 whitespace-separated
79+
tokens where the first token starts with an ASCII letter (element symbol)
80+
and tokens 2-4 are parseable as floats (the x, y, z coordinates).
81+
"""
82+
parts = line.strip().split()
83+
if len(parts) < 4:
84+
return False
85+
if not parts[0][0].isascii() or not parts[0][0].isalpha():
86+
return False
87+
try:
88+
float(parts[1])
89+
float(parts[2])
90+
float(parts[3])
91+
except ValueError:
92+
return False
93+
return True
94+
95+
96+
def split_xyz_frames(xyz_data: str) -> list[str]:
97+
"""Split a multi-frame XYZ string into a list of single-frame XYZ strings.
98+
99+
XYZ format: each record is ``<atom_count>\\n<comment>\\n<atom_count atom rows>``.
100+
A multi-record file simply concatenates these blocks.
101+
102+
The splitter is permissive:
103+
104+
* Accepts CRLF or LF line endings (output is always LF).
105+
* Recovers the next valid frame after a truncated or garbage block by
106+
advancing one line at a time rather than aborting the whole parse.
107+
* Validates each candidate atom row before accepting a block, so a fake
108+
``atom_count`` cannot greedily consume real lines from the next frame.
109+
* Drops zero-atom frames (degenerate).
110+
* Normalizes the count line (strips surrounding whitespace).
111+
* Each emitted frame is itself a valid XYZ string ending in ``\\n``.
112+
113+
Args:
114+
xyz_data: Raw text, possibly containing multiple concatenated XYZ records.
115+
116+
Returns:
117+
A list of self-contained single-frame XYZ strings. Empty list if the
118+
input contains no recognizable frames.
119+
"""
120+
if not xyz_data or not xyz_data.strip():
121+
return []
122+
123+
lines = xyz_data.replace("\r\n", "\n").replace("\r", "\n").split("\n")
124+
frames: list[str] = []
125+
i = 0
126+
n = len(lines)
127+
while i < n:
128+
token = lines[i].strip()
129+
if not token.isdecimal():
130+
i += 1
131+
continue
132+
atom_count = int(token)
133+
if atom_count == 0:
134+
# Skip the count line + the comment line; emit nothing.
135+
i += 2
136+
continue
137+
# Need atom_count + 2 (count + comment + atoms) lines starting at i.
138+
end = i + 2 + atom_count
139+
if end > n:
140+
# Truncated frame: advance by one line and keep scanning.
141+
i += 1
142+
continue
143+
# Validate every candidate atom row before accepting the block.
144+
atom_rows = lines[i + 2 : end]
145+
if not all(_is_plausible_atom_row(row) for row in atom_rows):
146+
i += 1
147+
continue
148+
# Normalize the count line; keep all other lines verbatim.
149+
block_lines = [lines[i].strip()] + lines[i + 1 : end]
150+
# Each emitted frame ends in exactly one newline.
151+
frames.append("\n".join(block_lines) + "\n")
152+
i = end
153+
return frames

‎app/modules/toolkits/rdkit_wrapper.py‎

Lines changed: 30 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -1038,50 +1038,42 @@ def convert_xyz_to_mol(
10381038
embed_chiral: bool = True,
10391039
use_huckel: bool = False,
10401040
cov_factor: float = 1.3,
1041-
) -> Chem.Mol:
1042-
"""Parse an XYZ block and assign connectivity + bond orders via xyz2mol.
1041+
) -> tuple[Chem.Mol, str]:
1042+
"""Parse an XYZ block and assign bonds via a two-tier RDKit pipeline.
10431043
1044-
Reads coordinates with ``Chem.MolFromXYZBlock`` (which yields a Mol with
1045-
only atoms and a 3D conformer, no bonds), copies it via ``Chem.Mol(raw)``
1046-
so the original is preserved, then runs ``rdDetermineBonds.DetermineBonds``
1047-
to perceive connectivity *and* bond orders in a single charge-aware pass.
1044+
Tier 1: ``rdDetermineBonds.DetermineBonds`` (xyz2mol) -- full bond-order
1045+
perception, charge-aware. Works for organic/main-group molecules.
1046+
Tier 2 (only if Tier 1 raises): ``rdDetermineBonds.DetermineConnectivity``
1047+
-- VdW connect-the-dots, all bonds order 1, but works for any element
1048+
including transition metals. Information-poor but accurate-by-construction.
10481049
10491050
Args:
1050-
xyz_data (str): Plain-text XYZ block (atom count, comment line, then
1051-
``element x y z`` rows). Must include the header lines.
1052-
charge (int): Net molecular charge. Required for non-neutral species.
1053-
Defaults to 0.
1054-
allow_charged_fragments (bool): If True, place formal charges on
1055-
atoms; if False, assign radical electrons instead. Defaults to
1056-
True (matches xyz2mol's default behavior).
1057-
embed_chiral (bool): Embed chirality information from 3D coords.
1058-
Defaults to True.
1059-
use_huckel (bool): Use extended Hückel theory for bond perception
1060-
instead of the van der Waals method. More accurate for
1061-
unusual valences but slower. Defaults to False.
1062-
cov_factor (float): Multiplier on covalent radii for VdW
1063-
connect-the-dots. Defaults to 1.3.
1051+
xyz_data: Plain-text XYZ block including the count + comment header.
1052+
charge: Net molecular charge (Tier 1 only). Defaults to 0.
1053+
allow_charged_fragments: xyz2mol option (Tier 1 only). Defaults to True.
1054+
embed_chiral: xyz2mol option (Tier 1 only). Defaults to True.
1055+
use_huckel: Use extended Huckel theory in Tier 1. Defaults to False.
1056+
cov_factor: Covalent-radius multiplier (both tiers). Defaults to 1.3.
10641057
10651058
Returns:
1066-
Chem.Mol: RDKit molecule with 3D conformer and perceived bonds.
1059+
A tuple ``(mol, method)`` where ``method`` is one of
1060+
``"bond_orders"`` (Tier 1) or ``"connectivity_only"`` (Tier 2).
10671061
10681062
Raises:
1069-
ValueError: If the XYZ block cannot be parsed or bond perception
1070-
fails (malformed input, unknown elements, geometry that cannot
1071-
be reconciled with the requested charge).
1063+
ValueError: If the XYZ block cannot be parsed, or if both tiers fail.
10721064
"""
10731065
if not xyz_data or not xyz_data.strip():
10741066
raise ValueError("Empty XYZ data.")
10751067

10761068
raw_mol = Chem.MolFromXYZBlock(xyz_data)
10771069
if raw_mol is None:
10781070
raise ValueError("Failed to parse XYZ block.")
1079-
10801071
if raw_mol.GetNumAtoms() == 0:
10811072
raise ValueError("XYZ block contains no atoms.")
10821073

1083-
mol = Chem.Mol(raw_mol)
1074+
# Tier 1: full bond-order perception via xyz2mol.
10841075
try:
1076+
mol = Chem.Mol(raw_mol)
10851077
rdDetermineBonds.DetermineBonds(
10861078
mol,
10871079
charge=charge,
@@ -1090,7 +1082,16 @@ def convert_xyz_to_mol(
10901082
useHueckel=use_huckel,
10911083
covFactor=cov_factor,
10921084
)
1093-
except (ValueError, RuntimeError) as exc:
1094-
raise ValueError(f"Bond perception failed: {exc}") from exc
1085+
return mol, "bond_orders"
1086+
except (ValueError, RuntimeError) as tier1_exc:
1087+
tier1_err = str(tier1_exc)
10951088

1096-
return mol
1089+
# Tier 2: connectivity-only (VdW connect-the-dots). Works for any element.
1090+
try:
1091+
mol = Chem.Mol(raw_mol)
1092+
rdDetermineBonds.DetermineConnectivity(mol, covFactor=cov_factor)
1093+
return mol, "connectivity_only"
1094+
except (ValueError, RuntimeError) as tier2_exc:
1095+
raise ValueError(
1096+
f"Bond perception failed (tier 1: {tier1_err}; tier 2: {tier2_exc})"
1097+
) from tier2_exc

0 commit comments

Comments
 (0)