NEBCert

Quality checks for nudged-elastic-band (NEB) reaction paths, transition states and tunnelling-corrected rate constants. NEBCert reads an ORCA NEB / NEB-CI / NEB-TS output (or a table of image energies), checks the band and the transition state, and computes Eyring rate constants with Wigner and Eckart tunnelling factors.

A PASS means that the checks that were run found nothing wrong. It does not show that the path is the minimum-energy path, that the transition state connects the intended minima, or that the barrier is accurate at the chosen level of theory.

What it checks

With an electronic barrier (no zero-point, thermal or entropic terms, as from a band) the Eyring number is reported as NOT_APPLICABLE: it is not a TST rate constant, and it does not count towards the overall status, which is the worst of the individual statuses. NEBCert does not run or read IRC calculations; --irc records the user's statement that an IRC connected the intended minima.

For ORCA outputs NEBCert reads the last PATH SUMMARY table, the climbing image, the NEB and TS convergence banners, the PATH SUMMARY FOR NEB-TS transition state and the final frequencies. When an optimised NEB-TS transition state is present its energy is the barrier; otherwise the energy of a converged climbing image, and otherwise the maximum of a natural cubic spline through the images. A table printed before the images were computed (E = 0, as left by a job whose image calculations failed) is rejected.

The Eckart factor is κ(T) = exp(V₁/kT)/kT ∫ P(E) exp(−E/kT) dE with the transmission probability of the asymmetric Eckart barrier (Eckart 1930; Johnston and Heicklen 1962), the barrier width fixed by the imaginary frequency, and the integral evaluated by adaptive quadrature (scipy quad).

Validation (validation/)

All scripts and results are in validation/; every number below is regenerated by them.

• Tunnelling and rate constants (compare_overreact.py, eckart_highprecision.py): on 720 combinations of imaginary frequency (200–3000 cm⁻¹), forward and reverse barriers (2–40 and 2–60 kcal/mol) and temperature (200–1000 K), κ_Eckart agrees with a 40-digit mpmath reference to 1.7 × 10⁻⁹ (relative). overreact deviates from the same reference by up to 1.3 × 10⁻², above 10⁻³ at 73 points, all with κ < 1.9. κ_Wigner and the Eyring rate agree with overreact to 3 × 10⁻¹¹ and 1.3 × 10⁻⁹.
• Version 1.0.0 (the only release published before 1.2.0; legacy100_grid.py): its κ_Eckart differs from the reference by more than a factor of 2 at 529 of the 720 points (median ratio 6.8, range 3 × 10⁻⁸ to 2 × 10⁴³).
• Real ORCA 6.1.1 calculations (make_neb_inputs.py, benchmark_neb.py; B3LYP-D3(BJ)/def2-SVP; results in validation/results/neb_benchmark.csv, figure validation/figures/fig3_profiles.png). Eight reactions with NEB-TS and frequencies (H + HCl, OH + H₂, CH₃ + H₂, Cl⁻ + CH₃Cl, HCOH → H₂CO, H₂CO → H₂ + CO, HCN → HNC, NH₃ inversion), three degraded runs and one failed job.
  • Values read from the .out file agree with ORCA's full-precision files: band energies within 0.006 kcal/mol (.final.interp; the output table has five decimals), TS energies within 0.003 kcal/mol (_NEB-TS_converged.xyz), frequencies within 0.005 cm⁻¹ (.hess). κ_Eckart for the eight transition states agrees with the mpmath reference to 3 × 10⁻⁹.
  • Verdicts: PASS for six reactions; WARNING for HCOH → H₂CO and H₂CO → H₂ + CO, whose imaginary frequencies (2060 and 1880 cm⁻¹) put 298.15 K below the crossover temperature (κ_Eckart = 1.8 × 10⁵ and 7.6 × 10⁶). FAIL for a band stopped after three iterations (its highest image, 45.8 kcal/mol, against a converged barrier of 2.3 kcal/mol) and for a band whose end points collapsed into the same structure; the failed job is rejected. A 3-image NEB-CI passes, with the climbing-image barrier equal to the 8-image NEB-TS barrier (84.36 kcal/mol; the spline maximum is 85.62).
  • Version 1.1.0 warned on every reaction (the --irc bug), passed the failed job and crashed on the collapsed band. Version 1.0.0 passed the failed job, crashed on the collapsed band, and its κ_Eckart was 1.6 to 6 × 10⁴⁶ times the reference for the eight transition states.

Not validated

• The VASP neb.dat and Gaussian IRC parsers have not been tested on real output; the Gaussian IRC parser is not used by the command-line interface.
• Paths with several steps are flagged, not split.

Installation

pip install "git+https://github.com/sircalch/nebcert@v1.2.0"

From source:

git clone https://github.com/sircalch/nebcert.git
cd nebcert
pip install -e .[dev]

Command line

# ORCA NEB / NEB-CI / NEB-TS output (frequencies are read from a NEB-TS ... Freq run)
nebcert assess -i reaction.out --reaction "H + HCl -> H2 + Cl" --functional "B3LYP-D3(BJ)/def2-SVP" -o report/

# Gibbs energy of activation for the Eyring rate constant
nebcert assess -i reaction.out --gibbs-barrier-kcal 12.4 -o report/

# Table of image energies plus frequencies
nebcert assess -i neb.dat --frequencies "-1250,150,300,800" -o report/

# Synthetic demonstration (not a real calculation)
nebcert demo -o demo/

The report directory contains report.html, figures, a summary table (CSV and LaTeX), a draft methods paragraph that states the results of the checks, and citation.bib.

Python

from nebcert import (calculate_neb_profile_analysis, verify_ts_frequency_and_irc,
                     calculate_eyring_tst_rates, calculate_quantum_tunneling_corrections,
                     assess_reaction_pathway_quality)
from nebcert.parsers.orca_neb import parse_orca_neb_output

d = parse_orca_neb_output("reaction.out")
ts = d["ts"]
neb = calculate_neb_profile_analysis(d["energies_ev"], d["coordinates_s_ang"], d["tangent_forces"],
                                     band_converged=d["neb_converged"],
                                     ts_energy_rel_ev=ts["barrier_fwd_ev"] if ts else None,
                                     ts_optimisation_converged=ts["converged"] if ts else None)
freq = verify_ts_frequency_and_irc(d["frequencies"])
tst = calculate_eyring_tst_rates(neb.e_forward_barrier_ev, barrier_type="electronic")
tun = calculate_quantum_tunneling_corrections(abs(freq.imaginary_frequency_cm1),
                                              neb.e_forward_barrier_ev, neb.e_reverse_barrier_ev)
report = assess_reaction_pathway_quality({"reaction": "H + HCl -> H2 + Cl"}, neb, freq, tst, tun)
print(report.overall_status, report.validation_score)
for message in report.recommendations:
    print("-", message)

Citation

nebcert cite

prints the APA and BibTeX entries for the installed version (concept DOI 10.5281/zenodo.22217586).

License

MIT, see LICENSE.