chematic is a cheminformatics library I've been writing in pure Rust, with bindings for Python and WASM — molecule parsing, property calculation, similarity search, 3D generation. Its design goal is compatibility with RDKit, the C++ library that's the de facto standard in the field.
The piece I'd been stuck on was aromaticity detection: deciding whether a given ring in a molecule counts as aromatic. The default logic, built by inferring local rules from cases where chematic disagreed with RDKit, had plateaued at 99.44%. I set that approach aside, built a separate experimental engine by porting RDKit's actual source, and got that one to 100.0000%. It also turned up an unrelated bug I hadn't been looking for.
What makes this hard
A ring like benzene gets extra stability because its electrons delocalize across the whole ring instead of staying in fixed double bonds. Chemists call this aromatic. Once a SMILES parser or SSSR (smallest set of smallest rings — the standard ring-perception algorithm) has found the rings in a molecule, something still has to decide whether each one satisfies the conditions for aromaticity, including Hückel's rule (a ring is aromatic if it has 4n+2 π electrons for some integer n).
The trouble is that this decision isn't unique for fused ring systems (rings sharing an edge or atom) or rings containing heteroatoms. RDKit itself ships four different aromaticity models — AROMATICITY_DEFAULT (identical to AROMATICITY_RDKIT), AROMATICITY_SIMPLE, AROMATICITY_MDL, and AROMATICITY_MMFF94 — and which one you pick changes the answer. "Correct" here is less a physical fact than an agreement about which tool's, which model's judgment you're matching. chematic targets RDKit's default model.
Round one: reverse-engineering rules from failures
The original implementation worked by collecting molecules where chematic's output diverged from RDKit's, then working backward to local rules — evaluating whether a given atom contributes to aromaticity on a per-candidate-ring, per-component basis.
That approach topped out at 99.44% atom agreement and 98.82% bond agreement against the production baseline (tracked in docs/rdkit_compat.md). The failure pattern was understood, too: fused ring systems involving a bridgehead nitrogen (a nitrogen shared between two rings) would sometimes borrow that nitrogen's lone pair for one ring even when the other ring the nitrogen belongs to was invalid. Rings with a carbonyl group, like tropone or 2-pyridone, classified correctly, so the bug wasn't in carbonyl handling — it was specifically in how bridgehead nitrogens got evaluated.
A concrete failing case: C1=Cc2ccccc2C2=NCCCN12, a benzene ring fused through a bridgehead nitrogen to a non-aromatic ring made of three sp3 carbons.
In the ring that shouldn't be aromatic, four ring carbons each contribute one π electron from a ring double bond, and the bridgehead nitrogen contributes two more from its lone pair — six total, which satisfies 4n+2 with n=1, so the whole ring gets marked aromatic. But the nitrogen's other ring has nothing legitimate to lend a lone pair to. The same lone-pair rule works correctly for indolizine, where both fused rings really are aromatic, so the rule itself wasn't wrong — the code just never checked whether the ring receiving the electrons was actually valid before lending them.
Reading RDKit's source instead of guessing more rules
crates/chematic-perception/src/rdkit_parity.rs ports getAtomDonorTypeArom, countAtomElec, isAtomCandForArom, applyHuckel, and applyHuckelToFused from RDKit's Code/GraphMol/Aromaticity.cpp, line by line, pinned to one commit — the payoff of giving up on inferring further special cases from failures and just checking what RDKit's own code does:
RDKit release: Release_2026_03_4
RDKit commit: e89c9f656a694fab4105139844cba88d2e013354
The BSD 3-Clause license text and copyright notice are preserved verbatim in THIRD_PARTY_NOTICES.md.
The difference the source revealed: RDKit computes each atom's ElectronDonorType (which kind of electron, if any, it can donate toward aromaticity) exactly once for the whole molecule. Whether a multiple bond counts toward electron donation depends only on whether that atom sits in some SSSR ring anywhere in the molecule — not on which candidate ring is currently being evaluated. chematic's existing implementation recomputed this per candidate ring and per component instead. The bridgehead-nitrogen false positive traced directly back to that per-ring recomputation.
Result: an exact match on 4,999 of 5,000 molecules
I validated the ported engine first against a 55-molecule diagnostic corpus:
false_positive corpus fixed: 33/33
false_negative corpus fixed: 5/5
negative_control maintained: 17/17
RDKit atom-flag agreement: 100.00% (1214/1214)
Then against the full 5,000-molecule benchmark corpus, comparing atom by atom and bond by bond with RDKit's actual output. One molecule got excluded for an unrelated, pre-existing kekulize() limitation; the remaining 4,999 gave:
Comparable molecules: 4,999 / 5,000
Atoms: 138,635 / 138,635
Bonds: 150,004 / 150,004
Unexplained aromaticity differences: 0
100.0000% on both atoms and bonds — down from 99.44% / 98.82% under the old rule-inference approach, with zero differences left unaccounted for.
The new engine ships as a new Rust API, kept separate from the existing assign_aromaticity_ex / apply_aromaticity_ex functions because it can fail in ways those can't:
pub enum AromaticityError {
KekulizationFailed { reason: String },
InternalInvariantViolation { reason: String },
}
pub fn assign_aromaticity_rdkit_parity_experimental(
mol: &Molecule,
) -> Result<AromaticityModel, AromaticityError>;
The existing API guarantees it won't fail regardless of whether you pick the Hückel or RdkitLike model. The new engine assumes its input is already kekulized (aromatic bonds converted to explicit alternating single/double bonds) and can hit known kekulize() limitations, so forcing it into the old no-fail contract wasn't an option — it gets its own function instead. It's Rust-only and marked experimental for now; it hasn't reached the Python bindings yet.
A nastier bug turned up along the way
Checking round-trip stability of aromaticity — whether canonicalizing a SMILES string twice in a row gives the same result — across a 94-case corpus surfaced something unrelated. chematic_core::apply_kekule() silently drops stereo_neighbor_order, internal bookkeeping about a stereocenter's neighbor ordering, whenever it rebuilds a molecule.
No error, no panic. The chirality flag itself (the @ / @@ designation) survives, but the neighbor-ordering data used to interpret it is gone. Code that later reinterprets @ / @@ against a different ordering — canonical atom rank, or CIP (Cahn-Ingold-Prelog, the naming system behind R/S descriptors) digraph position — can end up reporting the same stereocenter with the opposite configuration.
Not all 94 cases traced back to this bug. Reclassifying after the fix, only 9 were actually caused by it. 19 were a separate, already-known issue where a single bond at a ring fusion got incorrectly promoted to an aromatic bond. The remaining 66 were canonical-SMILES-writer formatting variance in complex fused ring systems, unrelated to aromaticity and out of scope for this pass.
I grepped the rest of the codebase for the same shape of bug mechanically:
| Function | Symptom |
|---|---|
apply_kekule |
Drops stereo configuration on rebuild (where this started) |
enumerate_stereoisomers |
Drops existing stereocenter info; newly assigned stereocenters sometimes flip CIP label between S and R |
transfer_hydrogen_aromatic / clone_mol
|
Same loss. clone_mol got deleted in favor of the Molecule: Clone impl — callers now just use .clone()
|
invert_stereocenter |
Inversion itself didn't work at all (below), plus the same data loss |
transfer_hydrogen |
Same pattern, on the keto-enol tautomerization path |
invert_stereocenter — which is supposed to flip the configuration at a given atom — was worse. It only had handling for 2D wedge bonds. A molecule whose stereochemistry was specified purely through SMILES @ / @@ fell into a branch meant for "no annotation present," which did nothing. No error raised. This is reproducible against the old release still sitting on PyPI:
>>> import chematic # chematic==0.4.29
>>> mol = chematic.from_smiles("C[C@@H](C(=O)O)N") # L-alanine
>>> mol.cip_stereo()
[{'atom_idx': 1, 'descriptor': 'S'}]
>>> inverted = mol.invert_stereocenter(1)
>>> inverted.cip_stereo()
[{'atom_idx': 1, 'descriptor': 'S'}] # unchanged
The same code against the fixed 0.4.30:
>>> import chematic # chematic==0.4.30
>>> mol.invert_stereocenter(1).cip_stereo()
[{'atom_idx': 1, 'descriptor': 'R'}] # S -> R
Every one of these fixes went through a revert-and-confirm-red step before merging: temporarily undo the fix and check that a test actually fails, then reapply it. A bug that raises neither an error nor a panic is easy for a test suite to paper over, so the process was to reproduce the broken state first and only then fix it, rather than trust a test written after the fact.
Where this stands
Aromaticity detection has landed on the Rust side first, as assign_aromaticity_rdkit_parity_experimental — it isn't in the Python bindings yet. The stereochemistry fixes, including invert_stereocenter and enumerate_stereoisomers, already shipped in 0.4.30.
Stacking rules on top of failing test cases has a ceiling. What actually got the error rate to zero was switching to reading the reference implementation and reproducing it line by line. Finding a bug that hands back a wrong stereocenter without ever raising an error is the kind of thing that only turns up when a verification pass gets carried all the way through instead of stopping once the numbers look good.
Top comments (0)