A page of the mepd manual, from mepd git commit 7db52fb. mepd's settings are named as on the command line (--search-budget) and in Python (search_budget). mepd Findings

Reaction channels#

mepd channels searches for the ways a reactant can turn into a given product. You give it the two ends of a reaction (xyz files with the same atoms, SMILES, or one reaction SMILES). It returns every transition state (TS) it found and checked, sorted by what each one connects:

The same search runs from Python as mepd.channels(start, end, output=..., profile=..., **options). Its options are the fields of mepd.ChannelsOptions: each field is the command-line option of the same name, with underscores for dashes (n_conformers is --n-conformers), and has the same default.

A search works in its output folder and writes as it goes. A rerun into the same folder resumes it: finished path searches and TS optimizations are read back, not repeated.

What a channel search does#

  1. Ends. SMILES are turned into 3D structures and minimized.
  2. Atom mapping. The product's atoms are renumbered so that atom i is the same atom in the reactant and in the product.
  3. Conformers. Several conformers of each end are generated, minimized and deduplicated.
  4. Pairs. Every reactant conformer is paired with every product conformer.
  5. Mechanisms. Each pair is mapped again, once per mechanism (which bonds break and which form), and only the pairs that fit each mechanism best are kept.
  6. Path search. Each kept pair gets a path search (GSM, NEB or the profile's method). A path that turns out to go through an intermediate is split there, and each piece is searched again.
  7. TS guesses. The highest point of every final path piece is a TS guess.
  8. TS optimization and IRC. Each guess is optimized to a saddle point, checked to have exactly one imaginary frequency, and followed down both sides (IRC) to two minima.
  9. Sorting. The minima the IRCs reach are compared with the requested ends by their bonds and stereochemistry. Each step becomes a channel, part of a multistep route, or an off-target step.

The sections below follow these steps. Settings lists the options, with the values used in our published benchmarks.

The pipeline#

Ends#

An end given as SMILES is embedded in 3D (RDKit, with OpenBabel as a fallback). If either end was given as SMILES (or both as one reaction SMILES), both ends are minimized at the profile's level of theory before anything else, with minimize_ends (on by default). If that minimization fails, the search stops. Ends given as xyz files are used as they are at this stage.

Why: an embedded SMILES is a guess, not a minimum. The atom mapping below compares geometries, and strain from a poor embedding would distort that comparison.

Atom mapping#

A path search moves each atom from where it is in the reactant to where it is in the product, so both ends must list the same atoms in the same order. mepd checks the order with SLAPMapper, which finds the atom correspondences that break and form the fewest bonds (bond orders are not counted). With atom_mapping on (the default), the product's atoms are renumbered when a correspondence fits better than the given order. With it off (--no-atom-mapping), the given order is kept; use that when the two ends already share a correct atom order, e.g. a dataset whose ends came from one IRC.

This first mapping only fixes a consistent starting order. It does not decide the mechanism: every conformer pair is mapped again (see Mechanisms).

Conformers#

Each end is sampled separately (backend, RDKit by default):

Each side's final pool is written to conformers/start_pool.xyz and conformers/end_pool.xyz. A later run can reuse them with start_pool / end_pool.

Pairing conformers#

Every reactant conformer is paired with every product conformer (with pair_from = both, the default). max_pairs can cap the number of pairs (the first ones in that order are kept). By default nothing is capped here: the selection happens per mechanism, next.

Mechanisms#

For each conformer pair, SLAPMapper is run again on that pair's geometries. Its tied best correspondences are grouped into mechanisms: the set of bonds broken and formed, with each atom named by its element and symmetry class (so swapping a CH₂'s two hydrogens is the same mechanism). For a Claisen rearrangement, for example, the [3,3] shift and a [1,3] shift tie, because bond orders are not counted, and both are kept.

Within one mechanism there are usually many equivalent numberings (which of a methyl's hydrogens goes where). They are all scored against this pair's geometry with atom_mapping_metric, and the best one is kept. The pair's current numbering is scored as one more candidate. Scores choose the numbering within a mechanism. They never choose between mechanisms.

Why per pair: which hydrogen ends up where depends on the conformers. A poor assignment forces an unnecessary rotation, such as a CH₂ turning, onto the path.

Then, for each mechanism, only the pairs_per_mechanism pairs with the best scores get a path search (see pairs_per_mechanism). The list of searches, with each one's reactant conformer, product conformer, mechanism and score, is written to pair_mechanisms.json.

If the mapping fails for a pair, that pair is kept with its current numbering under the mechanism "unmapped". With atom_mapping off, this step is skipped entirely: every pair is searched with the given numbering, and pairs_per_mechanism does not apply.

Each kept (pair, mechanism) gets its own recursive path search (mepd's MSMEP) with the method of the profile (path_min_method: GSM, NEB and others; mepd's default is FSM-TS):

  1. A path is built between the two ends (geodesic interpolation) and optimized by the path method.
  2. The result is checked: is it one elementary step? - Every local minimum inside the path is minimized. If one relaxes to a structure that is neither end, the path is split there. - Otherwise, the two sides of the highest point are minimized (a quick approximate IRC). If they do not lead back to the two ends, the path is split at the minima found. - With validate_minima_with_hessian (on), a candidate intermediate must pass a frequency check before the path is split at it (see hessian_minimum_frequency_cutoff).
  3. Each piece of a split is searched the same way, recursively, until every piece is one elementary step, or is kept unsplit for one of these reasons: the same pair of ends keeps coming back without anything new (5 times), a piece returns to ends it is already inside (5 times), the search budget has run out, or the search was stopped.

Each pair's search is written to pairs/pair_<i>_<j>/tree/: one file per piece searched.

TS guesses#

Every final piece of every pair's tree gives one TS guess: the highest-energy image of its path, endpoints excluded. A piece kept unsplit (for example when the search budget ran out) also gives a guess. That guess is the highest point of a path that may hold more than one step.

TS optimization and IRC#

Each guess is optimized to a first-order saddle point (Sella for the g-xTB and other ASE-based engines; with g-xTB starting from its exact Hessian). By default the optimization converges when the largest force is below 0.001 eV/Å, within 1000 steps, and an optimization that runs out of steps is continued once from where it stopped. The converged structure must then have exactly one imaginary frequency below −10 cm⁻¹. Smaller negative values are taken as numerical noise along very soft modes.

The TS is then followed down both sides (IRC), and each end of the IRC is minimized to a true minimum: the minimization is frequency-checked, and a structure with an imaginary mode is pushed along it and minimized again, up to 3 times.

Why minimize the IRC ends: an IRC stops short of the minimum, typically about 2 kcal/mol above it, and sometimes on the way to a different species. Only a real minimum can be compared with the requested ends.

A guess whose optimization fails is recorded, so a rerun into the same folder does not try it again. A guess whose optimization or IRC fails, or whose IRC reaches the same species on both sides, contributes no step.

When more than one search runs at once (workers > 1), each pair's TS optimizations start as soon as that pair's search is finished, alongside the searches still running. Otherwise they all run after the last search.

Single-step and multistep routes#

Every IRC that ends in two different species is one elementary step between them. The steps form a graph whose nodes are species (compared as in What counts as a verified route), with the requested reactant and product as two of the nodes:

The same TS is often found from several pairs. Two TSs of the same step count as one when their energies agree within node_ene_thre (kcal/mol) and their geometries match: either all their interatomic distances agree within 0.1 bohr, or their snap-RMSD, also against the mirror image, is below node_rms_thre (bohr). Both thresholds come from the profile's [chain_inputs] (5 and 5 by default; 1 and 1 in our benchmarks). Each distinct TS of a step is kept, the lowest-energy member of each group written as ts.xyz with its irc.xyz, and members.txt lists the guesses that reached it. In a multistep route, the lowest TS of each step defines the route, and the step's other TSs are kept in step_<n>/alternate_ts_<m>/.

A route's barrier is set by its highest TS. ChannelsResult.lowest() returns the route (channel or multistep) whose highest TS is lowest.

The search warns when a step reaches a stereoisomer of a requested end instead of the end itself. With direct_only, or when nothing but off-target steps was found, it also warns for every off-target step: the result does not match the requested ends.

What counts as a verified route#

A step is verified when all of these hold:

  1. The TS optimization converged, and the structure has exactly one imaginary frequency (below −10 cm⁻¹).
  2. The IRC ran to both sides, and both ends were minimized to true minima (no imaginary frequency below −10 cm⁻¹, after up to 3 pushes off a saddle).
  3. The two minimized ends are different species.
  4. For a channel, one minimized end is the same species as the requested reactant and the other the same species as the requested product. For a multistep route, every step is verified, and the steps chain from the reactant through intermediates to the product.

"Same species" means the same molecule, whatever its conformation and atom numbering:

When the stereochemistry cannot be decided (bond orders could not be assigned, a SMILES failed to build, or 100 atoms or more), the structures count as not the same. RDKit decides first. When it cannot, OpenBabel's stereo SMILES may confirm "same", but never "different".

The comparison is made against the requested ends themselves (after SMILES minimization and renumbering), not against the conformer a search started from.

Settings#

Values used in our published benchmarks, and mepd's defaults:

Setting CLI option Benchmarks Default
n_conformers --n-conformers 5 (RDKit) 50
rmsd_cutoff --rmsd-cutoff 0.1 bohr 0.1 bohr
atom_mapping_metric --atom-mapping-metric endpoint-rmsd snap-gi-xtb
pairs_per_mechanism --pairs-per-mechanism 5 3
direct_only --direct-only / --allow-multistep on the profile's, else off
search_budget --search-budget 1800 s 0 (none)
Hard time limit (the benchmark driver's --timeout) 7200 s (none in mepd)
hessian_minimum_frequency_cutoff --hessian-minimum-frequency-cutoff −20 cm⁻¹ −20 cm⁻¹
hessian_minima_rescue_displacement --hessian-minima-rescue-displacement 0.3 bohr 0.3 bohr
workers --workers, -j 1 0 (automatic)
parallel --parallel / --no-parallel off on

The path search's own settings come from the profile (--inputs, -i); see GSM settings. The benchmark profile also set [chain_inputs] node_rms_thre = 1.0 and node_ene_thre = 1.0, [gi_inputs] nimages = 20, and g-xTB with one process per search.

n_conformers#

CLI --n-conformers, field n_conformers, default 50; 0 means no cap. Benchmarks: 5, with the RDKit backend.

The most conformers kept for each end, the input geometry included. The cap applies after RDKit's embeddings have been deduplicated, and picks the conformers by farthest-point selection rather than lowest energy. Minimizing can merge conformers afterwards, so a side can end with fewer. With 5 per side there are at most 5 × 5 = 25 conformer pairs before mechanisms are counted.

Why farthest-point: a pool cut by energy clusters around the lowest conformers. Spreading the picks covers more of the shapes the molecule can take, which is what pairing different conformers is for.

rmsd_cutoff#

CLI --rmsd-cutoff, field rmsd_cutoff, default 0.1 bohr. Benchmarks: 0.1.

Two conformers of the same end closer than this (snap-RMSD, in bohr) count as the same, and the later one is dropped. It is applied after sampling, again after minimization, and molecule by molecule for multi-molecule ends. Mirror images are merged with the same cutoff.

atom_mapping_metric#

CLI --atom-mapping-metric, field atom_mapping_metric, default snap-gi-xtb. Benchmarks: endpoint-rmsd.

How the numberings of one mechanism, and then the pairs of one mechanism, are compared (see Mechanisms). The mechanisms themselves always come from SLAPMapper.

Cost: endpoint RMSD is the cheapest metric. The path-based metrics build an interpolation for every pair × mechanism × numbering they score, which can use up a search budget before any path search starts. The price is that endpoint RMSD knows nothing about the path.

pairs_per_mechanism#

CLI --pairs-per-mechanism, field pairs_per_mechanism, default 3; 0 means every pair × every mechanism. Benchmarks: 5.

For each mechanism, how many conformer pairs get a path search. Every pair is scored once per mechanism (its best numbering under atom_mapping_metric). Within each mechanism the pairs are sorted by that score, lowest first, and the first pairs_per_mechanism are kept. With endpoint-rmsd, these are the pairs whose reactant and product conformers are closest in shape. Pairs whose score could not be computed go last.

Every mechanism gets its own searches, so a mechanism that suits only a few conformers is still tried. With 25 pairs and 2 mechanisms, pairs_per_mechanism = 5 gives 10 path searches.

It applies only with atom_mapping on and pair_from = both (the default). Pairs built with --pair-from are capped by --atom-mapping-candidates instead.

Why: most pairs of one mechanism repeat each other's search. A few pairs per mechanism cover the mechanisms at a fraction of the cost of all pairs. The closest pairs are not necessarily those that lead to the lowest TS, so a higher value (or 0) finds more when a missed channel matters.

direct_only#

CLI --direct-only / --allow-multistep, field direct_only, default: the profile's path_min_inputs.direct_only, else off. Benchmarks: on.

When a path search splits at intermediates, only the pieces with at least one end at a requested species are searched further. The requested species are the original reactant and product, in any conformer or stereoisomer. A piece between two other species is skipped. For example, A → B splitting into A → C, C → D, D → A′, A′ → B (A′ another conformer or stereoisomer of A) searches A → C, D → A′ and A′ → B, and skips C → D.

What it skips, and what it does not:

Why: to spend the compute on the requested pair. A piece between two other species cannot be a direct step between the requested ends.

search_budget#

CLI --search-budget, field search_budget, default 0 (no budget). Benchmarks: 1800 s.

A wall-clock time in seconds. It is counted from the start of the channel search proper: after a SMILES end has been minimized and the ends' atom mapping checked, and before conformers are sampled. So sampling, minimizing conformers and mapping pairs all count against it.

Once the time has passed:

stats.json records search_budget_reached.

Why: a slow reaction still returns the channels it found instead of being cut off with nothing.

The time taken after the budget is not bounded by mepd: the running optimizations finish, and every TS optimization and IRC still runs. A run can therefore take much longer than its budget.

Hard time limit#

Not a mepd setting. Our benchmark driver (benchmark.py search --timeout, 7200 s per reaction in the published runs; its default is 3600) runs each reaction's search in its own process and kills it when the time is up. A killed reaction is recorded as timeout, and no route is reported for it, even if some steps were already on disk. The search budget above is the graceful stop inside mepd; this limit is the hard stop around it.

hessian_minimum_frequency_cutoff#

CLI --hessian-minimum-frequency-cutoff, field hessian_minimum_frequency_cutoff, default −20 cm⁻¹. Benchmarks: −20. Used when validate_minima_with_hessian is on (CLI --validate-minima-with-hessian / -H, on by default for mepd channels).

When a path search is about to be split at an apparent intermediate, the intermediate's frequencies are computed. It counts as a minimum only if its lowest frequency is at or above this cutoff (imaginary frequencies are negative). One that fails is rescued if possible. Otherwise it is rejected: the path is not split there, and the elementary-step check goes on without it. A single atom always passes. A Hessian that cannot be computed counts as a failure.

This check is on the split points inside a path search only. TSs and IRC ends have their own frequency checks (see What counts as a verified route), with a cutoff of −10 cm⁻¹.

Why: a minimization can stop on a saddle point, where the forces also vanish. Splitting there would add an intermediate and a step that do not exist. The cutoff is below zero because very soft modes (rotations of groups, loose complexes) can come out slightly imaginary from numerical noise.

hessian_minima_rescue_displacement#

CLI --hessian-minima-rescue-displacement, field hessian_minima_rescue_displacement, default 0.3 bohr. Benchmarks: 0.3.

When an intermediate fails the frequency check, it is pushed along its lowest mode by this distance, first one way and then the other, and minimized again. The first result that passes the check replaces it. Larger pushes are not tried. If neither push passes, the intermediate is rejected.

Why only one size: with g-xTB, larger pushes mostly minimized back to the same soft saddle point.

workers#

CLI --workers, -j, field workers, default 0 (automatic). Benchmarks: 1.

How many pair searches run at once, each in its own process, and how many TS optimizations and IRCs run at once afterwards. With 0, it is the number of searches to run, capped at the number of CPU cores, and is worked out again for each stage. With more than one worker, each pair's TS optimizations start as soon as that pair's search is done.

With 1, the searches run one after another, and all TS optimizations wait for the last search. In the benchmarks, parallelism came from running several reactions at once, one core each.

parallel#

CLI --parallel / --no-parallel, field parallel, default on; with parallel_workers (--parallel-workers, default automatic: the CPU cores divided by workers). Benchmarks: off.

Within one pair's search, the pieces of a split are searched at the same time instead of one after another.

GSM settings#

The path method and its settings are not mepd channels options. They come from the profile (--inputs / -i, or profile= in Python), in path_min_method and [path_min_inputs]. Our benchmarks used the growing string method (GSM; mepd runs the compiled molecularGSM, with mepd's patches):

path_min_method = "GSM"

[path_min_inputs]
nnodes = 20
max_opt_iters = 80
early_stop_on_minima = false

[gi_inputs]
nimages = 20

nnodes#

Profile path_min_inputs.nnodes; mepd's default for GSM is 9. Benchmarks: 20.

The number of nodes on the string, endpoints included. By default (seed_with_geodesic_interpolation = true) GSM does not grow its string from the two ends. It starts from mepd's geodesic interpolation, already fully grown, and the string then has as many nodes as that path has images (gi_inputs.nimages, 10 by default). nnodes is used when GSM grows the string itself: with seeding off, or when a seeded run fails and is retried from scratch. The benchmark profile set both to 20, so the string has 20 nodes either way.

max_opt_iters#

Profile path_min_inputs.max_opt_iters; mepd's default for GSM is 500. Benchmarks: 80.

The most optimization iterations molecularGSM runs on the string. A string that reaches this number without converging (perpendicular gradient RMS below conv_tol, 0.0005 by default) is kept as it is, with a warning. It is then checked like any other path: split if it shows an intermediate, else its highest point is the TS guess.

Trade-off: the TS is optimized and checked by IRC afterwards in any case, so the string only has to give a guess or show an intermediate. mepd's default is higher because at 80 strings often stopped unconverged (gradient RMS 0.002–0.009 against 0.0005), and an unconverged string depends more on small differences in its input.

early_stop_on_minima#

Profile path_min_inputs.early_stop_on_minima; mepd's default is on. Benchmarks: off.

With it on, GSM is stopped early once a local minimum inside the string has held steady for early_stop_persistence_window iterations (10), never before early_stop_min_opt_iters (20), and only after the string has finished growing. A path with an intermediate has no single TS to converge to, so without this it runs to max_opt_iters for nothing. The stop is checked by the elementary-step check. If that finds the path elementary after all, GSM resumes from where it stopped and runs to convergence.

With it off, as in the benchmarks, every string runs to convergence or to max_opt_iters.

Other settings#

Field (CLI option) Default What it does
backend (--backend) rdkit Conformer generator: rdkit, or crest (metadynamics; much slower, also samples how molecules sit relative to each other).
n_embed (--n-embed) 0 RDKit embeddings per end; 0 = from the number of rotatable bonds.
rdkit_torsion_prefs (--rdkit-torsion-prefs) both Embed with ETKDG's torsion preferences, without, or both pooled.
rdkit_ewin (--rdkit-ewin) none Keep only embeddings within this many kcal/mol (MMFF94) of the lowest.
random_seed (--random-seed) 0 Seed for the embedding.
minimize_ends (--minimize-ends) on Minimize SMILES ends and every conformer at the profile's level.
complex_energy_tol (--complex-energy-tol) 0.5 kcal/mol Multi-molecule ends: merge conformers whose molecules each match and whose energies agree within this.
atom_mapping (--atom-mapping) on Renumber the product by SLAPMapper, and map every pair per mechanism.
max_pairs (--max-pairs) 0 (no cap) Keep only the first N conformer pairs.
pair_from (--pair-from) both start / end / start+end: build each conformer's partner in its own frame instead of pairing sampled conformers.
conformers_only (--conformers-only) off Stop before the path searches (to size a run).
start_pool, end_pool (--start-pool, --end-pool) none Reuse a saved conformer pool.

mepd channels --help lists the rest (CREST options, --pair-from options, experimental mapping options).

What a search writes#

Path What
conformers/ Each end's conformer pool (start.xyz, end.xyz; *_pool.xyz to reuse).
pair_mechanisms.json Each path search: reactant conformer, product conformer, mechanism, score.
pairs/pair_<i>_<j>/tree/ Each pair's recursive path search.
ts/ Every optimized TS and its IRC, by guess.
channels/, alternate-channels/, offtarget-exit-channels/ The sorted steps and routes.
network.json Every species and step found, as one network.
stats.json Counts and times for each stage.

mepd.read_channels(folder) reads a finished folder back.

Limits and possible improvements#

References#