Fix conformer search correctness bugs and speed up hot pure-Python paths - #927
Open
alongd wants to merge 7 commits into
Open
Fix conformer search correctness bugs and speed up hot pure-Python paths#927alongd wants to merge 7 commits into
alongd wants to merge 7 commits into
Conversation
![github-advanced-security[bot]](https://avatars.githubusercontent.com/in/57789?s=60&v=4){kind=small}
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #927 +/- ##
==========================================
+ Coverage 63.24% 63.42% +0.17%
==========================================
Files 114 114
Lines 38318 38347 +29
Branches 10026 10033 +7
==========================================
+ Hits 24234 24321 +87
+ Misses 11169 11114 -55
+ Partials 2915 2912 -3
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
rdkit_conf_from_mol() ran a full AllChem.EmbedMolecule() distance-geometry embedding and then immediately overwrote every atom position from the given xyz, so the embedding result was always discarded. Build the conformer directly instead, the way conformers.embed_rdkit() already does. This function sits on the conformer search hot path: it is called once per torsion per sampled dihedral inside change_dihedrals_and_force_field_it(), which itself runs up to MAX_COMBINATION_ITERATIONS times. Measured 6.5x-70x faster on real species (ethane 12x, a 57-atom species 70x). Behavior change: previously, when the embedding failed (returns -1 for some strained structures), no conformer existed and the function returned None, which silently skipped the requested dihedral change at conformers.py:397 and conformers.py:761. A conformer is now always returned, so those dihedral changes are always applied.
Contributor
There was a problem hiding this comment.
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.
cluster_confs_by_rmsd() called compare_confs(), which rebuilt BOTH conformers' distance matrices from scratch on every pairwise comparison, so each conformer's matrix was recomputed once per representative rather than once. Precompute each conformer's upper-triangular distance matrix up front and feed it through compare_confs()' existing skip_conversion/dmat1/dmat2 path. Measured 3.7x-8x faster (10-150 conformers). xyz_to_dmat() built the coordinate array twice per call, and each build ran check_xyz_dict() again; build it once. 1.1x-1.4x faster. distance_matrix() materialized both an (N, M, 3) difference array and an (N, M, 3) squared copy; square in place to drop the second temporary. Output is bitwise identical (np.einsum was tried first and rejected: it differs from the original by 1-2 ULP).
Replace the per-pair Python double loop with a single numpy pairwise computation. The z-sorted early-break is preserved exactly as a mask (|dz| > 4.0 excludes the pair), the inclusive/exclusive distance bounds are unchanged, and bonds are still added in z-sorted (i, j > i) order so atom edge insertion order is untouched. Bond sets are identical on every species tested. The win is modest on real geometries (~3-6%): real molecules have genuine z-spread, so the original's early-break already pruned most of the inner loop. It only looks dramatic on artificially planar coordinates, where that break never fires. Adds the first direct tests for connect_the_dots(), covering 1 atom, H2, methane, a heteroatom ring, separated fragments and an empty molecule.
The MMFF optimization retry loop only incremented its counter in the `else` clause, so a call that raised on every attempt spun forever: `j` never advanced and `v` never changed. Break out of the loop on failure instead, and narrow the bare `except`, which also swallowed KeyboardInterrupt and SystemExit. Guard `ff is not None` before evaluating its energy, and import AllChem explicitly. The module called `Chem.AllChem.*` in five places while only importing `Chem`, which happened to work solely because converter.py imports the submodule first and thereby sets the attribute on the parent package. Without the bare `except` swallowing it, that latent AttributeError would otherwise escape into the conformer pipeline.
… the dedup scan Two correctness bugs in the iterative conformer search. Both change ARC's scientific output. 1. change_dihedrals_and_force_field_it() called get_force_field_energies() with optimize=True unconditionally, so the energy it returned always belonged to the FF-optimized geometry. When its own `optimize` argument was False it nevertheless stored the *unoptimized* torsion geometry alongside that energy, so every conformer it produced carried an energy describing a different structure. conformers_combinations_by_lowest_conformer() is exactly that caller, and get_lowest_confs() then ranked geometries by energies that were not theirs. Measured on ethylamine, the reported energy was off by up to 0.49 kcal/mol from the true energy of the reported geometry. Always store the matching optimized pair. The `optimize` argument is removed rather than left as a no-op: the FF optimization was performed either way, so it never saved any work, and its only effect was to decide whether the stored geometry was a lie. Conformers are now local minima labelled with their own energies. 2. The duplicate scan in conformers_combinations_by_lowest_conformer() treated compare_confs_fl(), a cheap first/last-atom-distance prefilter, as if a negative result ended the search. It only means the one candidate under test is not a duplicate, and the list is in generation order with no ordering invariant to exploit, so `break` skipped every conformer past the first dissimilar one and let duplicates through. Use `continue`, and carry the candidate's distance matrix across iterations, since compare_confs_fl() returns it as None precisely when the prefilter rejects. The conformer count for the deduce_new_conformers test drops from 9 to 6: those three were duplicates that the truncated scan had been admitting.
conformers_combinations_by_lowest_conformer() is an iterative descent: each round samples dihedrals from the lowest conformer found so far. It updated base_energy but never reassigned base_xyz, so all MAX_COMBINATION_ITERATIONS rounds re-sampled dihedrals from the original geometry and regenerated the same conformers, which were then discarded as duplicates. The loop-bottom exit could not stop this: it required the round's lowest conformer to compare equal to the frozen base geometry, which it essentially never does. So the loop always ran its full budget doing redundant force field work. Reassign base_xyz whenever the round improves on the base energy, and stop as soon as a round fails to improve. Also round base_energy at initialisation: per-conformer energies are stored via round(energy, 3), so comparing them against a full-precision base mixed precisions and made the convergence test unreliable. The old `if not newest_conformer_list: newest_conformer_list = [lowest_conf_i]` fallback is replaced by a plain break. It reused the previous round's lowest conformer, which is None on the first iteration. Measured via ARCSpecies.generate_conformers(n_confs=10) against official/main, the conformers produced are bitwise 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).
…ol to kcal/mol RDKit's MMFF/UFF and OpenBabel's MMFF94/UFF/GAFF implementations all report energies in kcal/mol, but DE_THRESHOLD and the de_threshold/e_confs docstrings described them as kJ/mol. The thresholds are compared directly against those energies, so the documented units were off by a factor of ~4.2.
8 tasks
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
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,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 calledget_force_field_energies()withoptimize=True, so the energy it returned described the FF-optimized geometry. But whenits own
optimizeargument wasFalseit stored the unoptimized torsion geometrynext to that energy.
conformers_combinations_by_lowest_conformer()is exactly thatcaller, 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
optimizeargument is kept as a deprecated, ignored keyword rather than removed:dropping it would have silently rebound a legacy positional
Falseontoforce_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 meansonly that that one conformer is not a duplicate, but the loop treated it as
break, soevery 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 sincecompare_confs_fl()returns it asNoneprecisely when the prefilter rejects.Because
continueremoves the short-circuit, this scan became a full pass per candidate, sothe 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 everycomparison. On hexanediol that cuts
xyz_to_dmat()builds by 37% (76 -> 48) with identicalconformers returned.
3. The iterative conformer descent never advanced its base geometry
conformers_combinations_by_lowest_conformer()updatedbase_energybut never reassignedbase_xyz, so allMAX_COMBINATION_ITERATIONS(25) rounds re-sampled dihedrals from theoriginal 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_xyznow advances whenever a round improves on the base energy, and the loop stopswhen 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 foreverThe MMFF retry loop incremented its counter only in the
elseclause, so a call thatraised on every attempt spun forever with no progress. The bare
exceptalso swallowedKeyboardInterrupt/SystemExit, and a failed optimization was still reported as anoptimized 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]againstxyzs == [X0, X1, X2]. Both call sites zip them, so X1 was handed E2 and X2 was droppedentirely. Non-converged conformers are now excluded from both lists together.
This also required importing
AllChemexplicitly: the module calledChem.AllChem.*infive places while importing only
Chem, which worked solely becauseconverter.pyimportsthe submodule first and thereby sets the attribute on the parent package. With the bare
exceptgone, that latentAttributeErrorwould have escaped into the conformer pipeline.Blast radius of 1-3
generate_conformers(n_confs=10)), the delivered conformers areidentical 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).
(
combination_threshold=10,n_confs=50,e_confs=50), fewer conformers aredelivered (e.g. hexanediol 21 -> 15,
OCC(O)CC(O)CO50 -> 15) because the duplicates thetruncated 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 generationis
EmbedMultipleConfs(..., randomSeed=1), a fixed seed, and it is used only by the earlierrandom-conformer stage, not here - this path is handed an explicit
xyz. So if a round failsto improve,
base_xyzis unchanged, and the next round re-samples the same dihedrals from thesame 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 thanthe 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 winIt ran a full ETKDG
AllChem.EmbedMolecule()and then overwrote every atom position fromthe 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
IndexErrorand a long one was silently truncated).Molecule.connect_the_dots()- vectorizedReplaces 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):
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, whichremove_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 overcoords.shape[1]rather thanhardcoding x/y/z so it stays equivalent to the original
sum((c1 - c2) ** 2)for anydimensionality. Non-finite coordinates now raise
ValueError: previously every comparisonagainst
NaNevaluatedFalse, control fell through to theelsebranch, and a bond wasadded between every such pair, silently producing a nonsense fully-bonded molecule.
cluster_confs_by_rmsd()- compute each fingerprint onceEach conformer's
np.triu(dmat)is computed once and passed throughcompare_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 coordinatearray 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.einsumwas tried first and rejected - it differsby 1-2 ULP, and this function is on a production path via
perceive.py:78). It shows nomeasurable speed change; it is included for the reduced allocation only, not as a win.
Units
DE_THRESHOLD,e_confsand thede_thresholddocstrings described force field energiesas 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_conformersconformer count, because ithad encoded the truncated dedup scan.
Notes for reviewers
connect_the_dots()is reachable only vias_bonds_mol_from_xyz()(test-only callers)and
Molecule.from_xyz()(no callers), andcluster_confs_by_rmsd()has no productioncall 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.MMFFOptimizeMoleculeConfs(), reusing one RDKitMolacross the dihedral loop, theO(N^2)
.index()into_rdkit_mol(), andcolliding_atoms().