Fix conformer search correctness bugs and speed up hot pure-Python paths

[alongd]

What this is

Seven commits. Four of them change ARC's scientific output and are described first;
the behaviour-preserving speedups follow. Everything was measured on one machine against
main (b0dd288c7) using a pristine baseline worktree, interleaved A/B where possible,

= 15 repetitions, reporting medians.

Behaviour changes (please review these first)

1. Conformer geometries were paired with energies belonging to a different structure

change_dihedrals_and_force_field_it() always called get_force_field_energies() with
optimize=True, so the energy it returned described the FF-optimized geometry. But when
its own optimize argument was False it stored the unoptimized torsion geometry
next to that energy. conformers_combinations_by_lowest_conformer() is exactly that
caller, so every conformer it produced carried an energy that was not its own, and
get_lowest_confs() ranked geometries by foreign energies.

Measured on ethylamine, the reported energy was off by up to 0.49 kcal/mol from the
true energy of the reported geometry. It is now 0.000000 on every sample.

The optimize argument is kept as a deprecated, ignored keyword rather than removed:
dropping it would have silently rebound a legacy positional False onto force_field.

2. The duplicate scan stopped at the first non-duplicate

compare_confs_fl() is a cheap first/last-atom-distance prefilter. A negative result means
only that that one conformer is not a duplicate, but the loop treated it as break, so
every conformer past the first dissimilar one went unchecked and duplicates were admitted.
The list is in generation order with no ordering invariant that would justify the early
exit. Now continue, carrying the candidate's distance matrix across iterations since
compare_confs_fl() returns it as None precisely when the prefilter rejects.

Because continue removes the short-circuit, this scan became a full pass per candidate, so
the first/last-atom distance and the candidate's distance matrix are now computed once per
candidate and threaded into compare_confs_fl() rather than being rebuilt on every
comparison. On hexanediol that cuts xyz_to_dmat() builds by 37% (76 -> 48) with identical
conformers returned.

3. The iterative conformer descent never advanced its base geometry

conformers_combinations_by_lowest_conformer() updated base_energy but never reassigned
base_xyz, so all MAX_COMBINATION_ITERATIONS (25) rounds re-sampled dihedrals from the
original geometry and regenerated the same conformers, which were then discarded as
duplicates. The loop-bottom exit required the round's lowest conformer to compare equal to
the frozen base, which essentially never happens, so the full budget was always spent on
redundant force-field work.

base_xyz now advances whenever a round improves on the base energy, and the loop stops
when a round fails to improve. Energies are no longer rounded to 3 decimals before the
internal comparison, which had made the convergence test degenerate against a 1e-3
tolerance.

4. rdkit_force_field() could hang forever

The MMFF retry loop incremented its counter only in the else clause, so a call that
raised on every attempt spun forever with no progress. The bare except also swallowed
KeyboardInterrupt/SystemExit, and a failed optimization was still reported as an
optimized conformer. Fixed, narrowed, and the conformer is now skipped so the existing
UFF/OpenBabel fallback can engage.

The same function's UFF fallback had the identical geometry/energy pairing defect as #1: it
appended the xyz unconditionally but the energy only on convergence, so with three conformers
where 0 and 2 converge and 1 does not you got energies == [E0, E2] against
xyzs == [X0, X1, X2]. Both call sites zip them, so X1 was handed E2 and X2 was dropped
entirely. Non-converged conformers are now excluded from both lists together.

This also required importing AllChem explicitly: the module called Chem.AllChem.* in
five places while importing only Chem, which worked solely because converter.py imports
the submodule first and thereby sets the attribute on the parent package. With the bare
except gone, that latent AttributeError would have escaped into the conformer pipeline.

Blast radius of 1-3

• With default settings (generate_conformers(n_confs=10)), the delivered conformers are
  identical for ethanol, n-butanol, n-hexanol, iso-octane, glycerol and hexyl radical,
  while wall time drops 1.6x-2.3x on the species with several rotors (n-hexanol
  11.61s -> 5.11s, glycerol 3.03s -> 1.37s, hexyl radical 2.68s -> 1.63s).
• Under deliberately aggressive settings that force this code path
  (combination_threshold=10, n_confs=50, e_confs=50), fewer conformers are
  delivered (e.g. hexanediol 21 -> 15, OCC(O)CC(O)CO 50 -> 15) because the duplicates the
  truncated scan had been admitting are now caught. In every case tested the lowest-energy
  conformer is preserved to within 3e-4 kcal/mol.

Does the greedy stop cost us the global minimum?

The goal of this search is to find the global minimum, so the only acceptance criterion that
matters is whether the new early stop can return a higher minimum than before. It cannot,
and this is provable rather than merely measured.

The round is deterministic in base_xyz. The only randomness anywhere in conformer generation
is EmbedMultipleConfs(..., randomSeed=1), a fixed seed, and it is used only by the earlier
random-conformer stage, not here - this path is handed an explicit xyz. So if a round fails
to improve, base_xyz is unchanged, and the next round re-samples the same dihedrals from the
same geometry and produces bit-identical conformers, which then dedup away. Continuing cannot
reach a new basin; it is exactly the redundancy described above. The conformers generated by
the final, non-improving round are still kept - they are appended before the loop exits - so
nothing is discarded either.

Verified empirically across 18 species (alcohols, diols, polyols, glycols, amines,
diamines, an acid, an ester, a thiol, an aldehyde, a branched alkane and a radical), under
both default settings and forced settings (combination_threshold=10, n_confs=50,
e_confs=50): the baseline ran its full 25 rounds and never found a lower minimum than
the early-terminating branch. In all 36 comparisons the largest deviation is 4e-4 kcal/mol,
which is the de-rounding described above, not a search difference.

Performance (behaviour-preserving)

rdkit_conf_from_mol() - the largest win

It ran a full ETKDG AllChem.EmbedMolecule() and then overwrote every atom position from
the supplied xyz, so the entire embedding was computed and thrown away - once per torsion
per sampled dihedral, inside the 25-iteration loop above. It now builds the conformer
directly. Also Set3D(True), AddConformer(assignId=True), and a coordinate-count check
(previously a short xyz raised a raw IndexError and a long one was silently truncated).

species	atoms	before	after	speedup
ethane	8	0.609 ms	0.069 ms	8.9x
isobutane	14	0.772 ms	0.088 ms	8.7x
1-hexanol	21	1.543 ms	0.125 ms	12.4x
methylcyclohexane	21	1.032 ms	0.131 ms	7.9x
anthracene	24	2.413 ms	0.176 ms	13.7x
branched dodecanol	57	18.137 ms	0.343 ms	52.9x

Molecule.connect_the_dots() - vectorized

Replaces the O(N^2) Python double loop with a numpy pairwise formulation. Measured
compiled against compiled (this module is cythonized, so timing a Python copy of the
old function would measure the Cython compiler rather than the change):

species	atoms	before	after	speedup
naphthalene	18	0.281 ms	0.164 ms	1.7-2.1x
hexanol	21	0.391 ms	0.188 ms	1.7-2.2x
hexyl radical	22	0.433 ms	0.179 ms	2.4x
octane	26	0.462 ms	0.179 ms	2.5-3.1x
iso-octane	26	0.529 ms	0.182 ms	2.8-2.9x
decalin	28	0.506 ms	0.186 ms	2.6-2.7x
pentacontane	152	10.354 ms	0.545 ms	18.6-20.1x

The function's "delete all bonds and set them again" preamble never worked: get_bonds()
returns a {atom: bond} dict, so iterating it yielded atoms, which remove_edge() rejects.
Any call on an already-bonded molecule raised. It is fixed here (collect the bonds once,
deduplicating, then remove them) since this PR is the one rewriting and documenting the
function, and connect_the_dots() is now covered by an idempotency test.

Two deliberate details: the squared distance accumulates one axis at a time rather than
materializing an (N, N, 3) cube, and it iterates over coords.shape[1] rather than
hardcoding x/y/z so it stays equivalent to the original sum((c1 - c2) ** 2) for any
dimensionality. Non-finite coordinates now raise ValueError: previously every comparison
against NaN evaluated False, control fell through to the else branch, and a bond was
added between every such pair, silently producing a nonsense fully-bonded molecule.

cluster_confs_by_rmsd() - compute each fingerprint once

Each conformer's np.triu(dmat) is computed once and passed through compare_confs()
instead of both distance matrices being rebuilt on every pairwise comparison.
3.7x (10 conformers) rising to 11.1x (150). xyz_to_dmat() builds its coordinate
array once instead of twice: 1.11x-1.46x.

common.distance_matrix() squares in place rather than allocating a second (N, N, 3)
temporary. This is bitwise identical (np.einsum was tried first and rejected - it differs
by 1-2 ULP, and this function is on a production path via perceive.py:78). It shows no
measurable speed change; it is included for the reduced allocation only, not as a win.

Units

DE_THRESHOLD, e_confs and the de_threshold docstrings described force field energies
as kJ/mol. RDKit's MMFF/UFF and OpenBabel's MMFF94/UFF/GAFF all report kcal/mol, and these
thresholds are compared directly against those energies, so the documented units were off
by a factor of ~4.2. Comments only; no behaviour change.

Testing

Full unit suite: 2493 passed, 3 skipped (the 3 skips are pre-existing, UMA env
unavailable), up from 2477 on main. 16 new tests.

One existing expectation changed: test_deduce_new_conformers conformer count, because it
had encoded the truncated dedup scan.

Notes for reviewers

• Two of the three functions this started from are not on ARC's production path:
  connect_the_dots() is reachable only via s_bonds_mol_from_xyz() (test-only callers)
  and Molecule.from_xyz() (no callers), and cluster_confs_by_rmsd() has no production
  call sites. Their speedups are real but do not move ARC's wall-clock. The wins that do are
  rdkit_conf_from_mol() and the conformer descent fix.
• Deliberately left out, as separate concerns: batching MMFF via
  MMFFOptimizeMoleculeConfs(), reusing one RDKit Mol across the dihedral loop, the
  O(N^2) .index() in to_rdkit_mol(), and colliding_atoms().

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 64.34%. Comparing base (06c6ce6) to head (3c2fa61).
⚠️ Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #927      +/-   ##
==========================================
+ Coverage   64.20%   64.34%   +0.13%     
==========================================
  Files         119      119              
  Lines       39578    39601      +23     
  Branches    10266    10269       +3     
==========================================
+ Hits        25412    25482      +70     
+ Misses      11196    11137      -59     
- Partials     2970     2982      +12     

Flag	Coverage Δ
functionaltests	64.34% <ø> (+0.13%)	⬆️
unittests	64.34% <ø> (+0.13%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[Copilot]

Pull request overview

This PR fixes multiple correctness issues in ARC’s conformer generation/deduplication pipeline (several of which affected scientific output), and then applies targeted performance improvements to hot pure-Python/numpy paths (RDKit conformer construction, distance-matrix usage, bond perception).

Changes:

• Fixes conformer energy/geometry mispairing and dedup scan logic, and makes the iterative “lowest conformer” descent actually advance and terminate on non-improvement.
• Hardens RDKit force-field optimization against infinite retry/hang scenarios and keeps xyz/energy lists index-aligned across fallbacks.
• Speeds up key paths by avoiding redundant coordinate/distance-matrix construction and vectorizing Molecule.connect_the_dots().

Reviewed changes

Copilot reviewed 8 out of 8 changed files in this pull request and generated 2 comments.

Show a summary per file

File	Description
arc/species/converter.py	Avoid redundant coordinate array builds; construct RDKit conformers directly from provided xyz; add precompute hooks to first/last-atom prefilter; reduce repeated dmat work in RMSD clustering.
arc/species/converter_test.py	Adds regression tests for xyz_to_dmat, pinned-coordinate conformer construction, and compare_confs_fl precompute equivalence.
arc/species/conformers.py	Correctness fixes in conformer descent/dedup and energy handling; deprecates/ignores legacy optimize arg in dihedral+FF helper; imports and uses AllChem explicitly; hardens rdkit_force_field.
arc/species/conformers_test.py	Updates/expands tests to pin the corrected conformer search behavior and prevent RDKit FF hang/misalignment regressions.
arc/molecule/molecule.py	Vectorizes connect_the_dots(), fixes bond-removal preamble, and adds explicit validation for bad coordinates.
arc/molecule/molecule_test.py	Adds comprehensive connect_the_dots() tests including idempotency, NaN/inf handling, and reference-algorithm equivalence.
arc/common.py	Reduces allocation in distance_matrix() by squaring in-place.
arc/common_test.py	Adds unit tests for distance_matrix() including trivial, normal, and error cases.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[calvinp0]

It seems there are two cases of output is None but this currently only handles one of them.

1. optimize=False, never run UFF and no energies to report
2. optimize=True, MMFF gives us nothing and then UFF was called and raise and try_ob=False. So output is still None but we return it with geoms and zero energies but are treated as if it were optimized

[alongd]

Good catch — fixed.

You're right that output is None was covering two unrelated situations. Case 1 (optimize=False) is the intended one: no optimization was asked for, so geometries are kept and there are no energies to report. Case 2 is the bug: optimization was requested, UFF raised, and with try_ob=False the exception was swallowed, leaving output at None and dropping straight into the same branch — so unoptimized geometries were emitted as though they had been optimized, with an empty energies list.

The failure path now logs a warning and returns empty lists instead, so by the time control reaches the loop below, output is None can only mean "optimization was not requested" — which is what the comment there claims.

Two tests pin the pair apart so they can't be conflated again:

• test_rdkit_force_field_uff_failure_without_openbabel_reports_nothing — asserts both lists come back empty (fails before the fix; it returned all 3 geometries).
• test_rdkit_force_field_unoptimized_keeps_geometries_in_uff_fallback — pins that optimize=False still yields every geometry and no energies, so the fix doesn't over-correct into dropping the legitimate case.

[calvinp0]

So since we are not optimizing, we are then getting the FF-relaced structures instead which I suppose means many sampling points fall into the same energy well (and thus duplicates are thrown away) which in turn gets us ferew, different conformers? This was intended?

[calvinp0]

Why are we dropping the optimization being True for the conformer generation?

[alongd]

Answering both of your comments in this thread together, since they have the same root.

optimize never controlled whether optimization happened. The inner call was hardcoded:

xyz_, energy = get_force_field_energies(..., optimize=True, ...)   # always optimizes
if energy and xyz_:
    energies.append(energy[0])            # ALWAYS the relaxed energy
    if optimize:
        xyzs.append(xyz_[0])              # relaxed geometry
    else:
        xyzs.append(xyz_dihedrals)        # rigid geometry

The flag only chose which geometry was returned alongside an always-relaxed energy. This call site passed optimize=False, so it stored the rigid geometry paired with the relaxed structure's energy — the two described different molecules. That mispairing was worth 0.49 kcal/mol on the case I measured. So the argument is gone rather than pinned to True because there was never a path here that skipped optimization; keeping the flag would only have preserved the ability to request an inconsistent pair.

On getting fewer conformers — yes, intended. The sampling points are seeds for the minimization, not measurement points. Several seeds relaxing into one well means those seeds located the same minimum, and dropping the duplicates is the correct outcome rather than a loss. For the NCC torsion in your other comment, the six seeds collapse to three because that torsion genuinely has three staggered minima.

The acceptance criterion is whether the search still reaches the global minimum, not how many conformers survive. I checked that across 18 species × 2 modes: there is no species where this branch's lowest conformer is higher in energy than main's.

[calvinp0]

So the dihedral we asked for will not be the dihedral of the xyz since the xyz is a relaxed structure. We may then need to re-measure the dihedral?
On NCC, torsion (1,2,3,8):

requested	rigid (old optimize=False xyz)	relaxed (new xyz)	drift	energy
0	360.00	59.91	59.9	-6.150
60	60.00	59.91	0.1	-6.150
120	120.00	179.53	59.5	-6.150
180	180.00	179.53	0.5	-6.150
240	240.00	299.16	59.2	-6.150
300	300.00	299.16	0.8	-6.150

Or we could add a torsion constraintt to the MMFF optimisation for the torsions being scanned?

[alongd]

Reproduced your table exactly, to the decimal (0→59.91, 120→179.53, 240→299.16). Two things worth adding, then the fix.

The flat profile predates this PR. All six points report the same energy (−6.150) — and they do on main too, because main also took its energy from the relaxed structure while returning the rigid geometry. So main had geometry at the requested angle paired with an energy belonging to a different structure; this branch has geometry and energy agreeing, with the dihedral field now being the odd one out. Neither version ever produced a real torsion profile, so what you've found is a labelling defect rather than a regression in the scan.

Re-measuring the dihedral looks like the obvious fix but breaks the consumer. plotter.py:1088 matches the stored dihedral against the sampling point within 0.1° and then calls min(y); storing measured values drops most of those matches and can empty y outright. So in this PR I've left dihedral as the seed label and documented it as such in change_dihedrals_and_force_field_it — a seed for the minimization, not a measurement of the returned geometry, whose actual dihedral may differ because a well's minimum doesn't generally sit on the sampling grid.

Your constraint suggestion is the right answer, and it's what I'm implementing — in #932. That branch already does exactly this for ring puckers, in optimize_conformer_with_frozen_ring, whose docstring names the identical failure mode: the seed's pucker, "which an unconstrained optimization would otherwise erase", is held through stage 1 while everything else relaxes, then stage 2 releases the constraint and polishes into that basin's true local minimum. Torsion seeds have the same problem and should get the same treatment (ff.MMFFAddTorsionConstraint for stage 1).

One refinement on the form it should take: a constrained scan — reporting the constrained energies — would compare badly against base_energy in the descent loop, which is unconstrained. Constrained energies are always ≥ their unconstrained counterparts, so lowest_conf_i['FF energy'] < base_energy - BASE_CONFORMER_ENERGY_TOL would be biased against advancing and could stall the search. The two-stage form avoids that, because stage 2's energy is a genuine unconstrained minimum and stays directly comparable.

Worth flagging what this actually buys, since the ethylamine case doesn't show it: the real cost of releasing the constraint immediately is that a rigid rotation into a sterically strained arrangement can throw the minimizer clean out of the basin you were aiming at, landing it on a conformer you already have — where it is silently deduplicated as a success. That failure grows with substituent bulk, i.e. exactly the molecules where conformer search earns its keep. #932 has the global-minimum-recovery metric to measure whether the two-stage version recovers basins that are currently lost.

[alongd]

Correcting my reply above: I said I'd implement the constrained pre-relax in #932. I measured it first, and the measurement doesn't support doing it — so we're not, and I'd rather post the numbers than quietly drop it.

One-stage (current) vs two-stage (constrained pre-relax, then release and polish), over every rotor ARC detects, seeded every 30°:

species	rotors	distinct geometries 1→2 stage	global min
n-butane	3	7 → 7	−5.076 = −5.076
2,3-dimethylbutane	5	11 → 11	6.401 = 6.401
neopentyl alcohol	5	11 → 11	12.045 = 12.045
2,2,3,3-tetramethylbutane	7	43 → 38	29.762 = 29.762
glycerol	5	8 → 10	34.319 = 34.319
1-hexanol	6	13 → 13	−2.417 = −2.417
1,2-dichloroethane	1	3 → 3	5.014 = 5.014
ethylene glycol	3	6 → 7	16.065 = 16.065
iso-octane	7	20 → 25	18.432 = 18.432
n-heptane	6	13 → 13	−5.680 = −5.680

The global minimum is identical in all ten — not lower in a single case. The constraint changes which basins get visited (more in three species, fewer in one), but not the answer the search exists to produce.

It also specifically contradicts the mechanism I proposed to you. I argued that steric bulk makes the released minimizer overshoot its intended basin, so crowded molecules should lose minima. The most crowded species here, 2,2,3,3-tetramethylbutane, finds fewer distinct geometries with the constraint, not more.

Against that, the cost is concrete: two FF minimizations per seed instead of one, in the hottest loop of the conformer search — which would hand back a large share of the 1.6–2.3× end-to-end speedup this PR is for.

So the torsion path stays as it is, with dihedral documented as a seed label. The ring case in #932 remains different on purpose: a pucker seed is a distinct chemical state that an unconstrained relax collapses, so holding it changes which minima exist at all; a torsion seed is only a starting point on the way to minima that are reachable regardless.

Probe and raw output are kept outside the repo (torsion_two_stage_probe.py / torsion_two_stage_results.txt) — happy to paste the script or run it over any species set you think would break the conclusion. Thanks for pushing on this one; it was worth measuring.

[calvinp0]

Can dmat1_ be None? I may have missed it but it appears compare_confs_fl returns the dmat1 it was passed unchanged on the dissimilar path, and a computed one on a similar path.

[alongd]

Correct — that guard was dead code, and it's removed.

compare_confs_fl() returns the dmat1 it was handed on the dissimilar path and a freshly computed one on the similar path, so dmat1_ is None only when dmat1 was already None. if dmat1_ is not None: dmat1 = dmat1_ was therefore equivalent to a plain assignment in every case, and it now unpacks straight into dmat1, with a comment recording why that can't clobber an already-computed matrix.

[MaroGol]

conf can no longer be None on any path - the function either raises or returns a conformer - so the if conf is not None: guards at conformers.py:401 and :784 are now dead. Same situation as the dmat1_ guard you dropped earlier in this PR.

[alongd]

Correct on both counts, and there was a third one — fixed.

rdkit_conf_from_mol() used to reach a None through the embed:

try:
    AllChem.EmbedMolecule(rd_mol)
except:
    pass
conf = None
if rd_mol.GetNumConformers():        # 0 when the embed failed -> conf stays None
    conf = rd_mol.GetConformer(id=0)

The bare except: pass meant a failed embed left zero conformers and the function returned (None, rd_mol). That is what the guards were for. Now the conformer is constructed directly and always added, so the function either raises ConverterError or returns a real Conformer — rd_mol.GetConformer(id=conf_id) raises ValueError on a bad id rather than returning None, so there is no surviving path to a None.

Verified rather than assumed, since "no path returns None" is the kind of claim that is cheap to check and expensive to get wrong:

NCC / C / O / c1ccccc1 / CC(C)(C)C(C)(C)C / [OH] / O=C=O   ->  conf is None? False   (all)
GetConformer(valid id) is None?  False
GetConformer(bogus id)        ->  raises ValueError
coord/atom count mismatch     ->  raises ConverterError

Removed: the guards at conformers.py:401 and :784 — neither had an else, so dropping them is behaviour-preserving.

The third: plotter.py:118, if conf is None: return False in show_sticks(), same shape and equally dead. Worth noting it was the only way show_sticks() returned False for a failed embed; with the embed gone, the failure it guarded cannot occur.

The invariant is pinned by an existing test rather than left implicit — test_rdkit_conf_from_mol_coordinates_pinned asserts assertIsNotNone(conf) across four species including C12(CC3)CC(CC1)CC3C2, a cage chosen precisely because ETKDG fails to converge on it. That is the exact input that produced conf = None on main, so if anyone reintroduces a None path the test catches it.

Folded into Drop the discarded ETKDG embed in rdkit_conf_from_mol() rather than added as a follow-up commit, since that is the commit that made the guards dead — so it is now self-contained: removes the None-producing path, its dead guards, and pins the invariant.

One thing I deliberately did not change, to flag rather than silently fix: show_sticks() catches (ValueError, AttributeError), and ConverterError subclasses Exception, not ValueError — so a malformed xyz propagates out of show_sticks() instead of returning False. That is pre-existing (main raises ConverterError for a wrong-type xyz too) and orthogonal to this PR, so it seemed wrong to widen the scope here. Happy to fix it in this PR if you'd prefer.

[MaroGol]

Apart from the one comment I had, the optimizations look good.