Status: verified · Version 1 · Jan Estrada Pabon (ORCID) · 2026-10-15 · no DOI yet
Finding. On 30 benchmark reactions (3–56 atoms; 60 calls at 12 and 30 images), mepd's geodesic interpolation with the numba kernel produced bit-identical paths in all 60 calls and ran 2.1 times faster in total (26.0 → 12.3 s); the per-call median gain was 2.4 times, falling to 1.5 times for the one 56-atom reaction.
To find how a reaction happens, path-search methods such as the nudged elastic band and the growing string method start from a guess of the path between reactant and product. mepd makes that guess with geodesic interpolation (Zhu, Thompson and Martínez, 2019): it finds the shortest path between the two structures in a space where bonds stretch and bend smoothly, so atoms don't collide halfway.
Each geodesic path comes from a least-squares fit, solved with SciPy. A search may call it many times; GSM with geodesic node placement calls it at every step that adds or moves a node. The question here: can that solve be made faster without changing a single path?
The least-squares solver SciPy uses for this problem was rewritten step by step, with its inner loop compiled by Numba and doing the same floating-point operations in the same order. The matrices it needs are now assembled with array operations instead of Python loops.1
Both versions of mepd interpolated the same 30 reactions, 15 from a geodesic benchmark dataset and 15 from Reaction-QM, at 12 and at 30 images: 60 calls each.2 The two versions ran alternately, three rounds each, on one thread, and every path was compared coordinate by coordinate.3
All 60 paths were identical, to the last bit, between the two versions and across all three rounds.
| Before | With numba | |
|---|---|---|
| All 60 calls (median of 3 rounds) | 26.0 s | 12.3 s |
| Median call | 0.38 s | 0.15 s |
Overall the interpolation ran 2.1 times faster. Per call the gain had a median of 2.4 times, ranging from 1.3 to 3.4 times. It shrinks with molecule size: the median was 2.6 times for the 18 reactions with up to 15 atoms and 2.1 times for the larger ones. For the single 56-atom reaction it was 1.5 times at 12 images and 1.3 times at 30. Longer paths gained less too: 2.5 times at 12 images, 2.0 times at 30.4
mepd/geodesic_interpolation2/_fast_lsq.py re-implements the path SciPy's least_squares takes here
(method trf, no bounds, lsmr trust-region solver, callable sparse Jacobian), with the LSMR iterations in
a Numba-compiled kernel. Other options, or a missing Numba, fall back to scipy.optimize.least_squares.
The Jacobians are now assembled vectorized, directly in compressed sparse row form, and one Jacobian
evaluation used only for a sparsity pattern SciPy ignores with a callable Jacobian was dropped. The kernel
relies on private SciPy helpers (scipy.optimize._lsq.common); the mepd test test_geodesic_fast_lsq
compares it with SciPy.
From each source set (66 Geodesic-set reactions whose reference TS survives at g-xTB; a 500-reaction
Reaction-QM subset), reactions whose start and end differ by 1–4 bonds were sorted by atom count and 15 were
taken at evenly spaced ranks (reproduce/original_run/select_set30.py). The set covers sizes evenly in rank,
not in atoms: 18 of 30 have at most 15 atoms and only one has more than 26. All are neutral singlets with H, B,
C, N, O, F, Si, P, S and Cl. The set was built to benchmark GSM variants; here it only supplies endpoints.
reproduce/time_geodesic.py calls mepd.chainhelpers.run_geodesic with only the number of images, so the
interpolation's own defaults apply (friction 0.01; mepd's path methods usually pass 0.001).5 Calls ran
one after another in one process with OMP_NUM_THREADS=1, after one untimed warm-up call (Numba compiles on
first use; that compile time is excluded). Rounds alternated before, after, before, … in separate processes;
totals varied by at most 6.2% between rounds of the same version. The machine was a shared laptop (Intel Core
Ultra 9 290HX Plus, 24 logical CPUs) running other jobs, at a load average of about 8–10; absolute times
depend on that, the ratio less so. Paths are compared as the exact decimal text of every float64 coordinate.
The gain falls with atom count and with images, so the totals here are weighted toward small molecules. Why it falls was not investigated; molecules above 26 atoms are represented by one reaction.
Level of theory. None; geodesic interpolation is a geometric construction with no energies.
Software.
Chemical domain.
Settings the result depends on.
Tested on 30. Per source set, reactions with 1–4 bond changes between start and end, sorted by atom count, 15 taken at evenly spaced ranks. This follows the size distribution of the source sets, which is weighted to small molecules: one reaction above 26 atoms.
Not tested.
Reproducibility: statistical. Paths: identical to data/paths.xyz.gz within 1e-8 Bohr per coordinate (bit-identical on the original machine). Timing: total speed-up at least 1.5x on a rerun (reproduce/quick_check.py); absolute times depend on the machine and its load.
| File | Size | sha256 | What it is |
|---|---|---|---|
data/reactions.csv |
1 kB | 6ed8bff6b065… |
The 30 reactions. |
data/endpoints.xyz.gz |
20 kB | 5f4dfa31ec72… |
Multi-frame XYZ (Angstrom): start and end of each reaction. Comment line: idx, reaction_id, role (start | end), charge=, mult=. |
data/paths.xyz.gz |
0.5 MB | e64801224054… |
Every image of the 60 interpolated paths (Bohr, as mepd holds them; exact float64 decimal text), from the first "after" round. Comment line: idx, nimages=, image=. |
data/timings.csv |
32 kB | 3e6e8c385a0f… |
One row per call, version and round. |
data/reactions.csv columns:
data/timings.csv columns:
check.py recomputes all 23 numbers in this report from these files and compares them with expected_numbers.json.
Software. Python 3.12.14, numpy 2.5.3, scipy 1.18.1, numba 0.68.0, llvmlite 0.50.0 (mepd's uv environment).
mepd version. Git commit b366826, the start of the commit hash that identifies this exact version of the mepd code (git checkout b366826).
mepd version. Git commit 1b8d63f, the start of the commit hash that identifies this exact version of the mepd code (git checkout 1b8d63f).
Cost. About 4 minutes for all six rounds; the quick check about a minute.
Download. F0003-v1.zip holds everything below: the data and the scripts.
| File | What it is |
|---|---|
check.py |
Recomputes every number in this report from data/. Needs only Python 3. |
reproduce/quick_check.py |
Reruns both versions once on all 60 calls, checks every path against data/paths.xyz.gz and reports the speed-up. About a minute. |
reproduce/time_geodesic.py |
Times one version on all 60 calls and writes its paths. |
reproduce/build_tables.py |
Builds data/timings.csv and data/paths.xyz.gz from six time_geodesic.py runs. |
reproduce/figures.py |
Redraws the figure. Needs matplotlib. |
data/endpoints.xyz.gz |
Start and end structures of the 30 reactions, the only input. |
reproduce/original_run/select_set30.py |
How the 30 reactions were picked from the two source datasets. |
Check the numbers (a few seconds):
python check.py
Rerun. This needs two mepd checkouts, at commit 1b8d63f (before) and b366826 (with numba). mepd is not
public yet. Use one Python environment with mepd's dependencies, including numba.
MEPD_BEFORE=/path/to/mepd-1b8d63f MEPD_AFTER=/path/to/mepd-b366826 python reproduce/quick_check.py
# the full timing, three alternating rounds (about 4 minutes):
for r in 1 2 3; do for v in before after; do
PYTHONPATH=/path/to/mepd-$v OMP_NUM_THREADS=1 python reproduce/time_geodesic.py $v$r rerun
done; done
python reproduce/build_tables.py rerun
The paths should match data/paths.xyz.gz exactly on the same software versions. The times depend on the machine
and what else it is running; compare the ratios.
Reproductions.
Reconstructed: written up after the run; some of it is inferred. This is a timing script, not an mepd command, so mepd's provenance.json was not written. Code commits (both clean), package versions, CPU and run order were recorded at run time in provenance/run_environment.json. The machine load is approximate (load average read before and after the runs). The reaction selection script was recovered from the session that ran it on 2026-10-13.
Cited.
Data this builds on.
Software used.