"""Screen every Reaction-QM B3LYP reaction for validity at the g-xTB level. python -m benchmark.screen_gxtb --sample 200 # a random subset first python -m benchmark.screen_gxtb # all 199,890 B3LYP-RXN reactions python -m benchmark.screen_gxtb --report # summarize the results so far Run from the folder above src/, with PYTHONPATH=src. Needs mepd, g-xTB as `gxtb-xtb`, and B3LYPD3_TZVP.h5 and B3LYPD3_TZVP_reaction_info.csv from doi:10.5281/zenodo.18551029 in data/reaction_qm/. The question this answers: when `mepd channels` fails to recover a Reaction-QM reaction, is that the search's fault or is the reaction simply not there on the g-xTB surface? Those need separating before any channel statistic means anything. The test takes the *reference* B3LYP saddle as a gift -- no searching at all -- and asks only whether g-xTB agrees it is a saddle connecting the same two species: 1. Re-converge the B3LYP TS geometry to a g-xTB saddle. 2. Run an IRC from it. 3. Require the IRC's two endpoints to match the reference reactant and product by bond connectivity, in either direction. A reaction that fails this cannot be found by any search at this level of theory, because the thing being searched for does not exist there. A reaction that passes but that `channels` misses is a genuine search failure. The pass set is also the honest denominator for every later benchmark number, and the converged saddles are themselves a g-xTB-level reference TS set. Endpoint graphs come from perceiving the reference *geometries*, not from Reaction-QM's SMILES strings -- the same convention mepd uses, and the only one that can be compared against a perceived IRC endpoint. Results stream to JSONL so the run is resumable; re-invoking skips reaction ids already recorded. """ from __future__ import annotations import argparse import json import os import random import signal import time from pathlib import Path from typing import Optional import numpy as np from benchmark import reaction_qm as rqm _REPO_ROOT = Path(__file__).resolve().parents[2] DEFAULT_INPUTS = _REPO_ROOT / "mepd_profile.toml" DEFAULT_OUT = _REPO_ROOT / "gxtb_screen" / "results.jsonl" # Most reactions finish in single-digit seconds; this only exists so one # pathological case cannot hold a worker forever. DEFAULT_TIMEOUT_S = 300 _WORKER: dict = {} class _Timeout(Exception): pass def _alarm(signum, frame): raise _Timeout() def _init_worker(inputs_path: str, timeout: int) -> None: """One engine and one open h5 handle per worker, built once.""" from mepd.inputs import RunInputs os.environ.setdefault("OMP_NUM_THREADS", "1") import h5py _WORKER["run_inputs"] = RunInputs.open(Path(inputs_path)) _WORKER["h5"] = h5py.File(rqm.H5_PATH, "r") _WORKER["index"] = rqm.build_index() _WORKER["timeout"] = timeout _WORKER["frag_cache"] = {} signal.signal(signal.SIGALRM, _alarm) def _node(symbols, coords, charge: int = 0, multiplicity: int = 1): from qcconst.constants import ANGSTROM_TO_BOHR from qcdata.models.structure import Structure from mepd.nodes.node import StructureNode geometry = np.asarray(coords, dtype=float) * ANGSTROM_TO_BOHR return StructureNode( structure=Structure( symbols=list(symbols), geometry=geometry, charge=charge, multiplicity=multiplicity, ) ) def _separated_energy(engine, species, cache: dict) -> Optional[float]: """Sum of the g-xTB energies of each fragment optimized *on its own*. This is the reference point Reaction-QM's own `dE_dagger` uses: for a bimolecular reaction the barrier is measured from infinitely separated reactants, not from a pre-reaction complex. Taking the IRC's reactant-side endpoint instead would measure from the complex, which on the calibration sample inflated multi-fragment errors to +7.9 kcal/mol against +4.1 for unimolecular ones -- an artifact of the reference point, not a real disagreement, and one that would have biased the agreement filter against every bimolecular reaction. Fragments repeat across the dataset (2.2x on average), so results are memoized per worker by fragment SMILES.""" total = 0.0 for sp in species: key = sp.smiles if key in cache: if cache[key] is None: return None total += cache[key] continue try: traj = engine.compute_geometry_optimization( _node(sp.symbols, sp.coordinates, sp.charge, sp.multiplicity) ) node = traj[-1] if hasattr(traj, "__len__") else traj energy = float(node.energy) except Exception: cache[key] = None return None cache[key] = energy total += energy return total def _smiles(node) -> str: try: return getattr(getattr(node, "graph", None), "smiles", "") or "" except Exception: return "" def screen_one(reaction_id: str) -> dict: """Re-converge one reference saddle at g-xTB and check where its IRC goes.""" import mepd.cli as cli rec = {"reaction_id": reaction_id, "status": "error"} started = time.time() signal.setitimer(signal.ITIMER_REAL, _WORKER["timeout"]) try: rxn = rqm.load_reaction(_WORKER["h5"], reaction_id, _WORKER["index"]) if rxn.ts is None: rec["status"] = "no_ts" return rec charge, mult = rxn.charge, rxn.reactants[0].multiplicity ts_guess = _node(rxn.ts.symbols, rxn.ts.coordinates, charge, mult) ref_r = _node(*zip(*rqm._separate_fragments(rxn.reactants)), charge, mult) ref_p = _node(*zip(*rqm._separate_fragments(rxn.products)), charge, mult) rec["ref_reactant"] = _smiles(ref_r) rec["ref_product"] = _smiles(ref_p) rec["n_atoms"] = rxn.ts.n_atoms rec["n_heavy"] = rxn.ts.n_heavy rec["ref_barrier_kcal"] = round(rxn.barrier_kcal(), 4) if rxn.barrier_kcal() else None engine = _WORKER["run_inputs"].engine t0 = time.time() try: ts = engine.compute_transition_state(node=ts_guess) except Exception as exc: rec.update(status="ts_failed", detail=f"{type(exc).__name__}: {exc}"[:200]) return rec rec["t_ts"] = round(time.time() - t0, 2) from mepd.nodes.node import StructureNode if not isinstance(ts, StructureNode): rec["status"] = "ts_failed" return rec t0 = time.time() try: irc = engine.compute_irc_chain(ts) except Exception as exc: rec.update(status="irc_failed", detail=f"{type(exc).__name__}: {exc}"[:200]) return rec rec["t_irc"] = round(time.time() - t0, 2) if irc is None or len(irc) < 2: rec["status"] = "irc_failed" return rec a, b = irc[0], irc[-1] rec["irc_a"] = _smiles(a) rec["irc_b"] = _smiles(b) forward = cli._connectivity_matches(a, ref_r) and cli._connectivity_matches(b, ref_p) reverse = cli._connectivity_matches(b, ref_r) and cli._connectivity_matches(a, ref_p) if forward or reverse: rec["status"] = "match" rec["direction"] = "forward" if forward else "reverse" # Two barriers, because they answer different questions. # `complex` is measured from the IRC's own reactant-side well; # `sep` is measured from separately optimized fragments and is the # only one comparable to Reaction-QM's B3LYP dE_dagger. try: reactant_side = a if forward else b rec["gxtb_barrier_complex_kcal"] = round( (ts.energy - reactant_side.energy) * rqm.HARTREE_TO_KCAL, 4 ) except Exception: pass cache = _WORKER["frag_cache"] e_sep_r = _separated_energy(engine, rxn.reactants, cache) if e_sep_r is not None: rec["gxtb_barrier_sep_kcal"] = round( (ts.energy - e_sep_r) * rqm.HARTREE_TO_KCAL, 4 ) e_sep_p = _separated_energy(engine, rxn.products, cache) if e_sep_p is not None: rec["gxtb_dE_kcal"] = round( (e_sep_p - e_sep_r) * rqm.HARTREE_TO_KCAL, 4 ) # Keep the converged saddle: this is the g-xTB-level reference TS. rec["ts_xyz"] = "\n".join( f"{s} {c[0]:.6f} {c[1]:.6f} {c[2]:.6f}" for s, c in zip( ts.structure.symbols, np.asarray(ts.structure.geometry) / 1.8897261254578281, ) ) elif cli._connectivity_matches(a, b): # IRC relaxed into the same well both ways -- the optimizer found # a saddle, but not one that separates two species. rec["status"] = "degenerate_irc" else: rec["status"] = "mismatch" return rec except _Timeout: rec["status"] = "timeout" return rec except Exception as exc: rec["detail"] = f"{type(exc).__name__}: {exc}"[:200] return rec finally: signal.setitimer(signal.ITIMER_REAL, 0) rec["t_total"] = round(time.time() - started, 2) def _already_done(out_path: Path) -> set: done = set() if out_path.is_file(): with out_path.open() as fh: for line in fh: try: done.add(json.loads(line)["reaction_id"]) except Exception: continue return done def run_screen( out_path: Path = DEFAULT_OUT, *, inputs: Path = DEFAULT_INPUTS, jobs: int = 22, sample: Optional[int] = None, seed: int = 0, timeout: int = DEFAULT_TIMEOUT_S, chunksize: int = 4, ) -> None: import multiprocessing as mp out_path = Path(out_path) out_path.parent.mkdir(parents=True, exist_ok=True) ids = list(rqm.build_index().keys()) if sample: random.Random(seed).shuffle(ids) ids = ids[:sample] done = _already_done(out_path) todo = [i for i in ids if i not in done] print(f"{len(ids)} target reaction(s); {len(done)} already recorded; {len(todo)} to run") if not todo: return started = time.time() n_done = 0 counts: dict = {} ctx = mp.get_context("fork") with ctx.Pool( jobs, initializer=_init_worker, initargs=(str(inputs), timeout) ) as pool, out_path.open("a") as fh: for rec in pool.imap_unordered(screen_one, todo, chunksize=chunksize): fh.write(json.dumps(rec) + "\n") n_done += 1 counts[rec["status"]] = counts.get(rec["status"], 0) + 1 if n_done % 200 == 0 or n_done == len(todo): fh.flush() rate = n_done / (time.time() - started) eta = (len(todo) - n_done) / rate / 3600 if rate else float("nan") summary = " ".join(f"{k}={v}" for k, v in sorted(counts.items())) print(f" {n_done}/{len(todo)} {rate:.1f}/s ETA {eta:.1f}h {summary}", flush=True) print(f"finished in {(time.time()-started)/3600:.2f} h") def report(out_path: Path = DEFAULT_OUT) -> None: import collections import statistics out_path = Path(out_path) if not out_path.is_file(): print(f"no results at {out_path}") return recs = [] with out_path.open() as fh: for line in fh: try: recs.append(json.loads(line)) except Exception: continue n = len(recs) counts = collections.Counter(r["status"] for r in recs) print(f"{n} reaction(s) screened at g-xTB") for k, v in counts.most_common(): print(f" {v:>7} ({100*v/n:5.1f}%) {k}") matched = [r for r in recs if r["status"] == "match"] if matched: pairs = [ (r["ref_barrier_kcal"], r["gxtb_barrier_sep_kcal"]) for r in matched if r.get("ref_barrier_kcal") is not None and r.get("gxtb_barrier_sep_kcal") is not None ] if pairs: errs = [g - b for b, g in pairs] mae = statistics.mean(abs(e) for e in errs) print(f"\nbarrier agreement on {len(pairs)} matched reaction(s):") print(f" g-xTB - B3LYP: MAE {mae:.1f} mean {statistics.mean(errs):+.1f} " f"median {statistics.median(errs):+.1f} kcal/mol") for thresh in (5, 10, 20): within = sum(1 for e in errs if abs(e) <= thresh) print(f" within {thresh:>2} kcal/mol: {within} ({100*within/len(errs):.0f}%)") t = [r["t_total"] for r in recs if "t_total" in r] if t: print(f"\nper-reaction wall: median {statistics.median(t):.1f}s mean {statistics.mean(t):.1f}s " f"max {max(t):.0f}s total {sum(t)/3600:.1f} CPU-h") def main(argv=None) -> int: p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) p.add_argument("--out", type=Path, default=DEFAULT_OUT) p.add_argument("--inputs", type=Path, default=DEFAULT_INPUTS) p.add_argument("--jobs", type=int, default=os.cpu_count() or 1, help="reactions screened in parallel") p.add_argument("--sample", type=int, help="screen only a random subset (calibration)") p.add_argument("--seed", type=int, default=0) p.add_argument("--timeout", type=int, default=DEFAULT_TIMEOUT_S) p.add_argument("--chunksize", type=int, default=4) p.add_argument("--report", action="store_true") args = p.parse_args(argv) if args.report: report(args.out) return 0 run_screen(args.out, inputs=args.inputs, jobs=args.jobs, sample=args.sample, seed=args.seed, timeout=args.timeout, chunksize=args.chunksize) report(args.out) return 0 if __name__ == "__main__": raise SystemExit(main())