CREST-based TS conformer search adapter - #921
Conversation
8e09a0a to
45a03ee
Compare
45a03ee to
39d7c01
Compare
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #921 +/- ##
==========================================
+ Coverage 65.67% 66.09% +0.41%
==========================================
Files 121 125 +4
Lines 41015 41886 +871
Branches 10548 10738 +190
==========================================
+ Hits 26938 27685 +747
- Misses 11064 11141 +77
- Partials 3013 3060 +47
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:
|
e6a47d8 to
174199b
Compare
d2b0575 to
f4ae037
Compare
b31f635 to
b2e2076
Compare
b2e2076 to
89ba170
Compare
| from arc.species.species import ARCSpecies, TSGuess, SpeciesError, colliding_atoms | ||
| from arc.species.zmat import get_parameter_from_atom_indices, remove_zmat_atom_0, up_param, xyz_to_zmat | ||
| from arc.species.vectors import calculate_angle | ||
| from arc.job.adapters.ts.seed_hub import get_ts_seeds |
31231d0 to
d70196a
Compare
c988137 to
bba78a3
Compare
23beab3 to
03321d8
Compare
LeenFahoum
left a comment
There was a problem hiding this comment.
Looks good overall! I left a few comments.
| CREST_REACTIVE_CORE_SIZES = {'H_Abstraction': 3, | ||
| 'carbonyl_based_hydrolysis': 4, | ||
| 'ether_hydrolysis': 4, | ||
| 'nitrile_hydrolysis': 4, |
There was a problem hiding this comment.
Should we also add XY_Addition_MultipleBond: 4 here? The CREST constraints for this family use all four *1-*4 atoms, so it looks like its reactive core size should also be 4. If so, the [2, 2] case in test_tiny_system_gate_skips_crest_only_for_whole_molecule_core should also expect True
There was a problem hiding this comment.
@LeenFahoum Yup! it is my intention but right now #978 needs to go in before this PR if we want to have XY_Addition_MultipleBond because right now ARC loads RMG-DB as default so it does not load in that family. The #978 PR is to change default to all and then I can adjust this PR to include it too
There was a problem hiding this comment.
I have added XY_Addition_MultipleBond: 4 and flipped [2,2] case to assertTrue. Also added [3, 2] (5 atoms, one spectator, gated) and [6,2] (C2H4 + HCl at 8 atoms, not gated) so the boundary is pinned from both sides.
Also, with #978, t has been merged but the RMG-DB load is still default and not all so XY_Addition_MultipleBond won't be identified unless the user changes in their settings to load all
| if os.path.isfile(stale_best_path): | ||
| try: | ||
| os.remove(stale_best_path) | ||
| except OSError as e: |
There was a problem hiding this comment.
If removing the previous crest_best.xyz fails, could we return None here instead of continuing?
Otherwise process_completed_jobs() may later pick up the old geometry and treat it as the result of the new run, since it only checks whether the file exists
There was a problem hiding this comment.
True, fixed it. process_completed_jobs() accepts a geometry on os.path.exists(crest_best.xyz) alone, so a surviving stale file would be reported as the new run's result. crest_ts_conformer_search() now returns None when the removal fails, before writing any input.
I did change one idee though, the caller was doing break on None return, which would have dropped the remaining seed for the reaction. An undeletable file is a property of that directory, so now it's continue. Also made two tests to cover one being aborted and the other seeds still being prepared and then submitted.
| if job_info["status"] not in TERMINAL_JOB_STATUSES: | ||
| try: | ||
| job_info["status"] = check_job_status(job_id) | ||
| except Exception as e: |
There was a problem hiding this comment.
A temporary scheduler error could mark the job as errored even if it’s still running. I think its better to just log the exception here and keep the current status unchanged
There was a problem hiding this comment.
Good catch, the exception is now logged and the status is left as it was, so the job is polled again on the next cycle.
| selected_hydrogen = None | ||
| selected_heavy_atoms = None | ||
| for hydrogen_index in hydrogen_indices: | ||
| heavy_atoms = sorted( |
There was a problem hiding this comment.
This filter removes all H atoms, even though a free H can be valid in H-abstraction. I think it would be better to use the same free-vs-bound H check as _is_valid_h_abs_atom_assignment() so bonded H atoms are still excluded
There was a problem hiding this comment.
Nice find! _get_h_abs_atoms_from_xyz now admits the free H atoms as candidate partners, using the same _is_free_h_atom check as _is_valid_h_abs_atom_assignment, so R1H + R2 <> R1 + R2H recolves whilst a hydrogen bound to something else is still excluded as a spectator. CH4 + H test for this.
b7c237f to
afd7821
Compare
afd7821 to
82c91d3
Compare
| Returns: | ||
| List[str]: The families of both hydrolysis parameter sets. | ||
| """ | ||
| from arc.job.adapters.ts.heuristics import FAMILY_SETS |
82c91d3 to
ca6f8c0
Compare
ca6f8c0 to
8c18690
Compare
The absolute H-abstraction reactive-core distances BREAKING_BOND_LENGTH = 1.29 A and FORMING_BOND_LENGTH = 1.36 A are the medians of the breaking (A-H) and forming (H-B) bond lengths over a dataset of 1670 DFT-optimized H-abstraction transition states. That dataset is dominated by C, N and O partner atoms, which is why the absolute distances are applied only when the partner is one of those elements (ABSOLUTE_BOND_LENGTH_PARTNERS); any other partner keeps the multiplicative 1.2 x equilibrium-bond stretch of the well.
8c18690 to
f2edae5
Compare
Adds CREST as a TS-search adapter. CREST does not search for a saddle itself — it takes the TS
seeds ARC already generates heuristically and runs constrained GFN2-xTB metadynamics around them,
so the reactive core is pinned while the spectator degrees of freedom are conformationally sampled.
The result is a better-converged starting geometry for the subsequent DFT TS optimisation.
Commits
CREST: settings, constants, and installermake install-crest,angstrom_to_bohrCREST: TS-search adapter (H-abstraction) + seed hubCREST: XY-addition four-center seed builderCREST: document the TS-search adapterdocs/source/TS_search.rstEvery file is touched by exactly one commit.
Family support and gating
CREST is opt-in. An adapter runs only if it is both registered for the family in
ts_adapters_by_rmg_familyand present in the globalts_adapterslist;crestis deliberatelyabsent from the default list, alongside
goflowandrits, because it needs its own environmentand an external binary. Enable it per run:
Constraints are implemented for:
H_Abstractioncarbonyl_based_hydrolysis,ether_hydrolysis,nitrile_hydrolysis— ARC-native families, alwaysclassified, so these are what make CREST reachable in a normal run
XY_Addition_MultipleBond— registered forcrestonly, so the four-centre seed builder isreachable through CREST and through no default adapter. This family's group definition also
requires a halogen at
*4([F1s,Cl1s,Br1s]), so RMG files it under thehalogensfamily setrather than
default; in a default runARCReaction.familyisNonefor these reactions and theseed builder is not reached at all. It needs
rmg_family_set: 'halogens'or'all'.Chemistry decisions
H-abstraction reactive core is element-aware. ARC applies
r1_stretch = r2_stretch = 1.2tothe reference single-bond length, which makes the core element-dependent. For a carbon, nitrogen or
oxygen partner this PR uses absolute targets instead —
BREAKING_BOND_LENGTH = 1.29andFORMING_BOND_LENGTH = 1.36Å, the medians over 1670 DFT-optimised H-abstraction transition states,a set whose donor and acceptor are heavy atoms in all but 62 records. Every other partner element
keeps the multiplicative path, because those medians describe heavy-partner abstraction and are
wrong elsewhere: for an H-atom partner 1.36 Å sits about 0.46 Å out along the reaction coordinate
from a saddle near 0.90 Å, and for S, Si, P, Br and I the absolute target is shorter than the
equilibrium bond. The two bonds are decided independently, so one seed can take the absolute target
on the side facing a carbon and the multiplicative stretch on the side facing a bromine.
ABSOLUTE_BOND_LENGTH_PARTNERS = ('C', 'N', 'O')andDEFAULT_BOND_STRETCH = 1.2name the rule.Passing
r1_stretch/r2_stretchexplicitly still multiplies, for every element.Breaking / forming bond, in Å:
mainBehaviour is unchanged from
mainfor every partner outside C/N/O. For C/N/O the forming bond is0.05–0.21 Å longer than
mainproduces; a GFN2-xTB Hessian on the CH₄ + OH seed moves from−462.8 cm⁻¹ to −775.1 cm⁻¹, one imaginary mode in both cases. This is a change to the default
heuristicsadapter and so affects runs that never enable CREST.Hydrolysis uses all six pairwise distances among the electrophilic centre, leaving group, water
oxygen and transferring water hydrogen. Four atoms have 6 internal degrees of freedom and there are
C(4,2) = 6 pairwise distances, so this is exactly determining; a smaller set leaves the reactive core
free to relax away from the saddle during sampling.
Changes outside the CREST paths
These affect runs that never use CREST:
arc/constants.py— removed a duplicatebohr_to_angstrom = 0.529177that shadoweda0 * 1e10;added
angstrom_to_bohras its reciprocal.arc/scheduler.py—successful_methodsnow credits every entry in a guess'smethod_sources,deduplicated case-insensitively, so a geometry found by several adapters credits all of them.
arc/species/converter.py— newreorder_xyz_string; empty input now raisesConverterErrorrather than
IndexError, and blank interior lines are skipped per the documented contract.arc/job/adapters/ts/heuristics.py— the element-aware reactive-core rule above.Testing
106 new tests against
origin/main:arc/job/adapters/ts/crest_test.pyarc/job/adapters/ts/xy_addition_test.pyarc/settings/crest_test.pyarc/job/adapters/ts/heuristics_test.pyarc/species/converter_test.pyarc/scheduler_test.pyVerified: 884 passed, 32 skipped, 0 failed, 90 subtests across
arc/job/adapters/ts/andarc/settings/crest_test.py.The suite was validated by mutation testing rather than by coverage percentage — each new test was
checked to fail under a plausible corruption of the code it covers (flipped comparisons, swapped
atom indices, ±10% on physical constants, disabled guards, removed unit conversions). Notable
results: the
coords.refÅ→Bohr conversion and the TURBOMOLE column order are now pinned (both werepreviously mutable with the suite green), as are
--chrg/--uhf, the*3/*4orientation theasymmetric XY factors depend on, and the properness of the antiparallel rotation branch (returning
−Ithere is a mirror image, det = −1, i.e. a chirality inversion).Upgrade note
arc/constants.pxdgainsangstrom_to_bohr. A stale compiledarc.constantsshadows the source andlacks the symbol, so re-run
make compileafter pulling — otherwise the first CREST job failswith an
AttributeErrorat runtime rather than at import.Known limitations
XY_Addition_MultipleBondis only classified withrmg_family_set: 'halogens'or'all', per above.at
MAX_CREST_SEEDSmost-distinct seeds, with the number dropped logged. Distinctness is measuredon interatomic distance matrices, so a rigidly translated or rotated duplicate is recognised as one.
CREST jobs do. This is how every TS adapter behaves —
heuristics,autotst,gcn,kinbot,xtb_gsmandorca_neball run to completion insideexecute_incore— and it follows from thestage's contract, since the scheduler asks for guesses and needs them back before it proceeds.
CREST differs from the others in submitting its compute to the queue rather than running it in the
driver's own allocation, and in bounding it: a job is given up on after
MAX_CREST_WALL_TIME_SECONDSof running orMAX_CREST_QUEUE_WAIT_SECONDSwithout starting.written and no angle is checked, because restraining A–H, H–B and A–B fixes all three sides of the
reactive triangle and so determines the angle. Measured on three H-abstraction cases run on a
cluster (CH₄ + OH, C₂H₆ + HO₂, CH₃OH + CH₃; 21 geometries across
crest_bestand the ensembles),the restrained distances drift at most 0.006 Å from the seed and the A–H–B angle moves at most
9.2°, usually under 3°. The 0.15 Å tolerance is therefore about 25× the observed drift.
submit_job()inarc/job/local.pysignals failure as either(None, None)or('errored', '').The adapter handles both.