This is version 1 of F0002, and its results are outdated. Version 2 (2026-10-22) corrected them: Rerun with corrected code on the whole F0001 benchmark set (31,914 Reaction-QM reactions). Version 1's counts came from searches with the two errors listed here and should not be used. the path search passed coordinates to geodesic interpolation in the wrong unit (bohr for ångström) (fixed in git commit 5b18bda). reference transition states and IRC ends were not checked strictly (convergence, frequencies, minimized ends, stereochemistry) (fixed in git commit 0594999). Read the current version.
Status: verified · Version 1 · Jan Estrada Pabon (ORCID) · 2026-10-05 · no DOI yet
Finding. Searching with mepd channels on the g-xTB surface found a verified route whose highest transition state lies more than 1 kcal/mol below the dataset's own transition state (re-converged at g-xTB) for 103 of 500 Reaction-QM reactions (median 6.02 kcal/mol lower) and for 12 of 66 Geodesic-set reactions (median 15.72 kcal/mol lower).
mepd is a tool for mapping reaction mechanisms. Given a reactant and a product, its search, mepd channels,
looks for the transition state between them: the top of the
energy hill the reaction has to climb. A lower hill means a faster reaction, so the lowest transition state
connecting the two is the one that matters.
We test the search on reactions whose transition state is already known: one set from Reaction-QM, and one from the Geodesic MLP dataset, most of whose reactions come from a benchmark by Jónsson and coworkers. The search runs on g-xTB, a fast approximate quantum-chemistry method, so each known transition state was first re-converged on g-xTB (Finding F0001). The question here: does the search find a route that is lower than the known one?
For each reaction, mepd channels built several starting arrangements of the reactant and product and searched for paths between them on g-xTB.2 Every transition state it reports was checked: optimized, then followed downhill to confirm it joins the requested reactant and product, either directly or through intermediates.3
We compared the energy of the best route's highest transition state with that of the known transition state, both on g-xTB. A route counts as lower when it is more than 1 kcal/mol below.4 The sets were 500 Reaction-QM reactions and 66 Geodesic-set reactions.5
Version 2 is in progress. We are running the same search on every reaction of the F0001 benchmark set whose re-converged transition state leads back to both of the dataset's endpoints: 31,914 Reaction-QM reactions. The numbers and routes below are preliminary and are updated as the search runs (last update 2026-10-20, 1,346 reactions searched).6
| Reactions searched so far | 1,346 of 31,914 |
|---|---|
| Search found a verified route | 1,165 |
| Lower than the known transition state | 393 (29.2% of those searched) |
| … as a single step | 228 (median 2.85 kcal/mol lower) |
| … through intermediates | 165 (median 10.03 kcal/mol lower) |
| About the same (within 1 kcal/mol) | 670 |
| Only higher routes | 102 |
| No route found (no transition state connecting the two ends) | 180 (13.4%) |
| Stopped at the time limit | 1 |
| Largest amount lower | 87.84 kcal/mol |
The ten single-step and ten multistep routes furthest below the known transition state so far, beside the dataset's own paths:
Version 1's results, on 500 Reaction-QM and 66 Geodesic-set reactions, are in Details.1
| Reaction-QM | Geodesic set | |
|---|---|---|
| Reactions | 500 | 66 |
| Search found a verified route | 278 | 60 |
| Lower than the known transition state | 103 | 12 |
| About the same (within 1 kcal/mol) | 123 | 47 |
| Only higher routes | 52 | 1 |
| No route found | 222 | 6 |
| Median amount lower | 6.02 kcal/mol | 15.72 kcal/mol |
| Largest amount lower | 47.78 kcal/mol (RXN_0000151161) | 34.62 kcal/mol (EW_system18) |
So for one in five reactions in each set, the known transition state is not the lowest one on the g-xTB surface: 20.6% of the Reaction-QM set and 18.2% of the Geodesic set. Some lower routes are single steps (58 and 9), others go through intermediates (45 and 3), including both of the largest.3
These are counts of what one search budget found, so they are a floor, not a ceiling.7 Whether the lower routes stay lower at a higher level of theory was not tested here.
mepd channels with five RDKit conformers per side, direct paths only, five path searches per mechanism and a
30-minute search budget (45-minute hard limit) per reaction, one core each, GSM with 20 nodes as the path method.
Reaction-QM endpoints are docked multi-fragment complexes (xtb aISS) built from the dataset's separate fragments,
and atoms are mapped between the two ends by mepd's SLAPMapper with the endpoint-RMSD metric. The Geodesic set was run
twice, once mapped the same way and once without atom mapping, since its start, end and TS already share one atom
order; the better of the two counts. 24 Reaction-QM and 2 Geodesic searches hit the time limit.
Only routes mepd itself verified count. A single-step channel is one TS whose IRC reconnects the requested reactant and product species (bond graph and stereochemistry). A multistep alternate channel is a chain of such TSs through intermediates from reactant to product; its height is its highest TS. Routes that leave to other species are excluded. Geodesic routes with an IRC end that has a bridging hydrogen (H within 1.45 Å of two heavy atoms) were also excluded, since bond perception misreads them; this affected 3 reactions.
The gap is the g-xTB energy of the reference TS minus that of the best route's highest TS; both structures have
the same atoms, so no reactant reference is needed. Outcomes: lower (gap > 1 kcal/mol), matched (within ±1),
higher (only routes more than 1 kcal/mol above). All 115 recorded gaps reproduce from the stored geometries
within 0.001 kcal/mol (reproduce/quick_check.py).
Reaction-QM: 500 reactions from F0001's survivors, in two batches of 250. Batch 1 was drawn from survivors whose g-xTB barrier lay within 5 kcal/mol of B3LYP; batch 2 uniformly from the remaining survivors whose two ends differ in bond graph. Lower routes turned up at a similar rate in both (49 and 54). Geodesic: 66 of the dataset's 97 reactions, those whose B3LYP-D3(BJ)/def2-SVP TS survives at g-xTB and whose two ends differ in bond graph. All are neutral closed-shell singlets except two Geodesic anions.
The search settings are those of version 1, with four differences. The reactions are all 31,914 of F0001's benchmark set whose re-converged transition state's IRC, with both ends minimized, reproduces both dataset endpoints (bond graph and stereochemistry); those minimized IRC ends are the reactant and product the search starts from. Each search may run for 2 hours before it is stopped. The searches run on mepd git commits 0594999 and 9aee085; the second only uses less memory, and gives the same routes and energies. The reactant reference (the floor) is the lowest of the reactant and its search conformers; an IRC end of a route found can lie slightly below it, so the barriers shown are measured from the floor, while the amounts lower do not depend on it (both transition states share it). "No route found" means that no transition state connecting the requested reactant and product was found within the budget; it is not evidence that none exists. The per-reaction outcomes so far are in data/v2_search_outcomes.csv.
Level of theory. g-xTB 2.0.1 throughout: path search (GSM), TS optimization and IRC (Sella through ASE and mepd's g-xTB engine); gas phase. Reference TSs are the datasets' DFT TSs re-converged at g-xTB (Reaction-QM: F0001; Geodesic: same procedure).
Software.
Chemical domain.
Settings the result depends on.
Tested on 566. Reaction-QM: 500 of F0001's 29,833 survivors. Batch 1 (250) was drawn from survivors whose g-xTB barrier was within 5 kcal/mol of B3LYP; batch 2 (250) was drawn uniformly from all other survivors whose two ends differ in bond graph. Lower routes were found at a similar rate in both (49 and 54). Geodesic: all 66 of the dataset's 97 reactions whose B3LYP TS survives at g-xTB and whose two ends differ in bond graph.
Not tested.
Known exceptions.
Reproducibility: statistical. The recorded energy gaps reproduce from the stored geometries within 0.05 kcal/mol (reproduce/quick_check.py). Rerunning the searches is not expected to find the same routes: they depend on random conformers, timing and the code version.
| File | Size | sha256 | What it is |
|---|---|---|---|
data/scores.csv |
58 kB | 41eeadc272ff… |
One row per reaction in both sets: outcome, energies of the reference TS and of the best route found, the gap. |
data/lower_ts.xyz.gz |
98 kB | 01c4879ab338… |
Multi-frame XYZ (Å) for the 115 reactions with a lower route: the reference g-xTB TS, then each TS of the best route. Comment line: reaction id, role (reference | found | found_stepKofN), energy for the reference. |
data/inputs.tar.gz |
0.8 MB | 41b2292ae6a7… |
The two reaction sets: set tables plus start, end and reference transition-state structures of all 566 reactions; unpacks into reproduce/rerun/ |
data/scores.csv columns:
check.py recomputes all 46 numbers in this report from these files and compares them with expected_numbers.json.
Software. Python 3.12, mepd, g-xTB 2.0.1, patched molecularGSM, RDKit, Sella; one core per search.
mepd version. Git commit 672b509, the start of the commit hash that identifies this exact version of the mepd code (git checkout 672b509).
Cost. Reaction-QM sweep about 500 reactions × up to 45 min on one core each; Geodesic 2 × 66. Recomputing the stored energies: a few seconds.
Download. F0002-v1.zip holds everything below: the data, the scripts and the input files.
| File | What it is |
|---|---|
check.py |
Recomputes every number in this report from data/. Needs only Python 3. |
reproduce/quick_check.py |
Recalculates the g-xTB energy of every stored transition state and compares the amounts lower with the table. Needs g-xTB 2.0.1. |
reproduce/figures.py |
Redraws the figure. Needs matplotlib. |
reproduce/build_tables.py |
Builds data/scores.csv and data/lower_ts.xyz.gz from scored searches. |
data/inputs.tar.gz |
The two reaction sets: start, end and reference transition-state structures for all 566 reactions. |
reproduce/rerun/ |
The scripts that run and score the searches (src/benchmark/), the step that combines the two Geodesic searches, and the mepd input file mepd_profile.toml. |
Check the numbers (a few seconds):
python check.py
python reproduce/quick_check.py --jobs 3 # set GXTB=/path/to/xtb if g-xTB is not on PATH as gxtb-xtb
Rerun the searches. This needs mepd at git commit 672b509 (mepd is not public yet) and the mepd command on PATH. It also needs g-xTB 2.0.1 on PATH as gxtb-xtb, and the patched GSM from jandestrada/molecularGSM (branch mepd), whose path goes in mepd_profile.toml. Each search takes up to 45 minutes on one core: 500 for Reaction-QM and 66 twice for Geodesic. --jobs sets how many run at once.
cd reproduce/rerun
tar xzf ../../data/inputs.tar.gz
# From doi:10.5281/zenodo.18551029, put B3LYPD3_TZVP.h5 and B3LYPD3_TZVP_reaction_info.csv in data/reaction_qm/
export PYTHONPATH=$PWD/src OMP_NUM_THREADS=1
SEARCH="--n-conformers 5 --workers 1 --jobs 8 --timeout 2700 -- --direct-only --no-parallel
--pairs-per-mechanism 5 --rmsd-cutoff 0.1 --search-budget 1800
--hessian-minimum-frequency-cutoff -20 --hessian-minima-rescue-displacement 0.3"
python -m benchmark.reaction_qm_channels --inputs mepd_profile.toml --dataset sets/reaction_qm_500.csv \
--output reaction_qm_runs --atom-mapping-metric endpoint-rmsd $SEARCH
python -m benchmark.reaction_qm_channels --inputs mepd_profile.toml --dataset sets/geodesic_66.csv \
--output geodesic_runs_mapped --atom-mapping-metric endpoint-rmsd $SEARCH
python -m benchmark.reaction_qm_channels --inputs mepd_profile.toml --dataset sets/geodesic_66.csv \
--output geodesic_runs_unmapped --no-atom-mapping --max-pairs 5 --atom-mapping-metric rmsd-geodesic $SEARCH
Score the searches and rebuild the tables:
python -m benchmark.score_channels --runs reaction_qm_runs --dataset sets/reaction_qm_500.csv \
--inputs mepd_profile.toml --out reaction_qm_scores.json
python -m benchmark.score_channels --runs geodesic_runs_mapped --floor-roots geodesic_runs_unmapped \
--dataset sets/geodesic_66.csv --inputs mepd_profile.toml --out geodesic_mapped_scores.json
python -m benchmark.score_channels --runs geodesic_runs_unmapped --floor-roots geodesic_runs_mapped \
--dataset sets/geodesic_66.csv --inputs mepd_profile.toml --out geodesic_unmapped_scores.json
python combine_geodesic.py
cd ../..
python reproduce/build_tables.py --out rerun_data
Then compare rerun_data/ with data/. These are the scripts the original searches used, tidied for reuse: files and folders renamed, and set-building code that is not needed (the sets are in data/inputs.tar.gz) removed. A dry run produces exactly the same mepd channels commands as the original searches.
Reproductions.
Reconstructed: written up after the run; some of it is inferred. These searches ran before mepd recorded its own provenance, so the details were pieced together afterward from the run folders' logs. The code is mepd git commit 672b509. Whether the Geodesic searches ran before its one change to atom mapping was applied (that is, with its parent, git commit 28b963f) was not recorded. Some run folders used in scoring have since been deleted; the amounts lower do not depend on them, because each is the difference between two transition-state energies.
Data this builds on.
Cited.
Software used.
Other Findings.