Neb Irc Activation Energy

NEB-IRC activation energy pipeline

Overview

Computes a reaction activation energy end to end — optimize reactant and product, find the minimum energy path with climbing-image NEB, refine the transition state with a Hessian, and confirm it with IRC — using GFN2-xTB through pysisyphus. Outputs a verified TS geometry, the barrier (ΔE‡; ΔG‡ after thermal corrections), an IRC energy profile (irc_energy_profile.png), and an animated TS imaginary-mode viewer (ts_imaginary_mode.html). Verification is required, not optional: a converged TS is meaningless until its single imaginary mode and IRC endpoints are checked.

When to Use

• Finding the transition state for an elementary reaction step and its activation energy
• Computing a reaction barrier (ΔE‡ or ΔG‡) to rank a series of related reactions
• Running a climbing-image NEB / minimum energy path between a reactant and product
• Verifying a candidate TS with IRC — does it connect the intended reactant and product?
• Screening barriers at a cheap semi-empirical level before committing DFT time

Reach for DFT on a multi-core node instead when you need quantitative agreement with experiment; GFN2-xTB barriers are semi-quantitative (see references/energetics.md). For a 2D scheme figure of the reaction, use the rdkit-chemdraw-cdxml skill instead.

Prerequisites

• Tools: xtb (GFN2-xTB engine), pysisyphus (pysis CLI for the path pipeline)
• Python: matplotlib for the IRC plot (the TS animation HTML needs no packages)
• Input: reactant.xyz and product.xyz with identical atom ordering (see Step 1)
• Environment: single core suffices; set OMP_NUM_THREADS to match physical cores

Work in a local scratch dir (e.g. /tmp/rxn/), not a mounted/networked workspace: pysisyphus creates and deletes symlinks and throws PermissionError mid-run on s3fs/FUSE. Copy results out at the end.

Materialize the bundled scripts into the scratch dir first. They can't be run in place from the skill directory, so use your file tools to read each one and save it into your working dir before running it. The scripts live in this skill's scripts/ folder (next to this SKILL.md):

• scripts/setup_env.sh
• scripts/pipeline.yaml
• scripts/check_result.py
• scripts/plot_irc.py

The TS imaginary-mode animation is not produced here — read the molecular-visualization-3dmol skill and use its mol_viewer.py (Step 6).

Check for the tools; install only if missing (inside pixi/conda, invoke via pixi run xtb):

cd /tmp/rxn
command -v xtb && command -v pysis || bash setup_env.sh   # xtb binary + pysisyphus, ~2-3 min
source "${ROOT:-${HOME:-/tmp}/xtbenv}/env.sh"              # re-source in every new shell

Workflow

Step 1: Prepare reactant and product geometries

Most failures originate here, not in the NEB. Build the product by editing a copy of the reactant so atom ordering is identical — a permuted order gives a path that is geometrically valid and chemically meaningless. Align non-reacting groups so a spectator conformational change does not fold into the barrier. For bimolecular reactions use a pre-reaction complex as the reactant, not separated fragments (NEB converges poorly from infinite separation, and the reference state changes the reported barrier — record it).

# Build product from a copy of the reactant, moving only the reacting atoms.
from pathlib import Path

lines = Path("reactant.xyz").read_text().splitlines()
natoms = int(lines[0])
atoms = [ln.split() for ln in lines[2:2 + natoms]]   # [symbol, x, y, z] per atom
# ... edit ONLY the coordinates of atoms that move; keep order + symbols ...
Path("product.xyz").write_text("\n".join([str(natoms), "product"] + [" ".join(a) for a in atoms]) + "\n")

Step 2: Check the endpoints sit in different basins

If the reacting groups start too close, preoptimization carries the reactant over the barrier and both endpoints relax to the same structure; the NEB then returns a flat profile and the TS search aborts. This looks like success in the log until it fails minutes later, so compare the two pre-optimized endpoints on the key reacting bond.

import numpy as np

def load_xyz(fn):
    lines = open(fn).read().splitlines()
    return np.array([[float(v) for v in ln.split()[1:4]] for ln in lines[2:2 + int(lines[0])]])

r, p = load_xyz("first_pre_opt.xyz"), load_xyz("last_pre_opt.xyz")   # written by preopt
i, j = 0, 5                                     # indices of the atoms whose bond changes
dr, dp = np.linalg.norm(r[i] - r[j]), np.linalg.norm(p[i] - p[j])
assert abs(dr - dp) > 0.3, "endpoints nearly identical: move the reacting fragment further out"

Step 3: Run the pipeline

The template chains preopt → IDPP interpolation → CI-NEB → RS-I-RFO TS optimization with Hessian → IRC both directions → endpoint reoptimization. Set charge and mult in pipeline.yaml before running — the default charge: 0 is wrong for any ion and converges silently to a meaningless TS. Add alpb: <solvent> for solution reactions. Whatever you set here must match every standalone xtb call in Step 5. Raise max_cycles to 150–200 above ~50 atoms.

pysis pipeline.yaml > pipeline.log 2>&1
tail -40 pipeline.log            # confirm it reached the endopt/IRC stage

ΔE‡ and the reaction energy come from the pipeline's | BARRIERS | block, referenced to the reactant endpoint (check_result.py in Step 4 reports the same numbers):

  Left:     0.00 kJ mol-1      # reactant endpoint
    TS:   107.84 kJ mol-1      # dE‡ = TS - Left = 107.84 kJ/mol
 Right:    10.19 kJ mol-1      # dE_rxn = Right - Left = 10.19 kJ/mol

Step 4: Verify the transition state

Never report a barrier before this passes. The checker applies all three gates (see Key Concepts) and exits 0 only if every one passes.

python3 check_result.py pipeline.log
# [PASS] imaginary frequencies  exactly 1 at -621.8 cm-1
# [PASS] IRC endpoints          forward and backward matched distinct inputs
# [PASS] NEB profile            elementary, single barrier, span 107.5 kJ/mol
#   electronic barrier dE‡ (GFN2-xTB): 107.8 kJ/mol = 25.77 kcal/mol

Step 5: Compute thermochemistry and report the barrier

Report ΔE‡ and the level of theory; for comparison against experimental rates, report ΔG‡, which needs Hessians on the TS and reactant. Every standalone xtb call must use the same --chrg, --uhf, and solvent as pipeline.yaml, on the pipeline's optimized geometries — not the raw input. A gas-phase Hessian against a solvated barrier, or a missing --chrg, yields a nonsensical (often negative) ΔG‡.

# flags must match pipeline.yaml (here: charge -1, aqueous)
xtb forward_end_final_geometry.xyz --hess --gfn 2 --alpb water --chrg -1 --uhf 0 > r_hess.log 2>&1
xtb ts_final_geometry.xyz          --hess --gfn 2 --alpb water --chrg -1 --uhf 0 > ts_hess.log 2>&1
# dG‡ = G(TS) - G(reactant), from the "TOTAL FREE ENERGY" lines
grep -i "TOTAL FREE ENERGY" r_hess.log ts_hess.log

Step 6: Deliver the two visuals

IRC energy profile — plot_irc.py reads gfn/charge/mult/solvent from pipeline.yaml and recomputes each IRC frame's energy at that exact level (the *_irc.trj comment lines carry none). Override with --charge/--mult/--alpb/--gfn only if you edited the pipeline after running.

python3 plot_irc.py               # -> irc_energy_profile.png

TS imaginary-mode animation — produced by the molecular-visualization-3dmol skill, not here. Read that skill and run its mol_viewer.py on the ts_imaginary_mode_000.trj that tsopt: do_hess: True wrote (grab the frequency from check_result.py for the label):

# after materializing mol_viewer.py per the molecular-visualization-3dmol skill
python3 mol_viewer.py ts_imaginary_mode_000.trj --mode trajectory \
    --title "Transition-state mode" --subtitle "imaginary mode -621.8 cm-1" --out ts_imaginary_mode.html

Key Parameters

Set in pipeline.yaml unless noted.