ORCA IRC Calculation Skill

This skill handles the full IRC workflow: generating an ORCA input file from a confirmed transition state geometry, and post-processing the output to extract the reaction-path energy profile and endpoint geometries.

Scope: Input generation, local post-processing, and (optionally) HPC submission via the CatGo workflow engine. The "Submitting to HPC" section below covers the proven Expanse flow. If the user is running on their own non-CatGo infrastructure, just generate the input file from the template in Stage 1 and skip the submission section.

Target version: ORCA 6.x. The parser and defaults below are written against the ORCA 6 output layout used by CatGo's own parser (server/catgo/utils/orca_output.py::OrcaIrcOutput).

Prerequisites

IRC requires a transition state geometry that has been confirmed by a frequency calculation showing exactly one imaginary frequency corresponding to the reaction coordinate.

If the user has not run a frequency calculation on the candidate TS, suggest the orca-freq skill before generating an IRC input. Running IRC on a geometry that isn't a true first-order saddle point produces meaningless paths. This is a suggestion, not a hard stop — if the user confirms they know what they're doing and still wants the input, generate it.

Stage 1: Input Generation

Before generating the input, gather from the user:

• TS geometry — file path or pasted XYZ coordinates
• Method and basis set — must match whatever was used for the TS optimization and frequency calculation. Method/basis consistency across Opt → Freq → IRC is critical; mixing levels means the TS may not be a stationary point on the new surface.
• Charge and multiplicity
• Number of cores (%pal nprocs)
• Memory per core (%maxcore, in MB)

Optional (apply only if user specifies):

• Solvent model (CPCM). Gas-phase IRC is standard practice for confirming connectivity; solvation mainly affects barrier heights, not topology.
• Dispersion correction — again, match whatever was used for the TS search.

Hessian Handling (important — read all of this)

IRC moves downhill from the TS along the imaginary-frequency mode. To identify that mode, ORCA needs a Hessian at the TS geometry. Getting this step right is the single biggest factor in whether the IRC gives a sensible path.

What the Hessian is used for, concretely:

1. Diagonalize to find the mode with a negative eigenvalue (the imaginary frequency). That eigenvector is the initial descent direction.
2. Build a mass-weighted quadratic model around the TS. The very first displacement along the mode is scaled by Scale_Displ_SD, chosen (via Init_Displ_DE) to correspond to a target drop in energy.
3. Subsequent steps use an adaptive integrator in mass-weighted coordinates; the Hessian is not re-diagonalized every step (that's what SD vs. LQA flavors are about), but the first-step quality dominates the final path.

Consequences:

• A Hessian with more than one negative eigenvalue means your geometry is not a clean TS. IRC will pick one of the negative modes and the result is arbitrary. Fix the TS first.
• A numerically noisy Hessian (bad SCF convergence, too-loose grid) can flip the sign of a small imaginary mode or rotate the eigenvector into a nearby low-frequency mode. The IRC then follows the wrong coordinate — often tumbling, rotation, or a methyl rock instead of the bond-making/breaking mode. calc_numfreq or InitHess read from a converged AnFreq job is the fix.
• Reusing a .hess file from a calculation at a different method/basis/grid/solvation is a silent footgun. The eigenmode numbering can differ, so "the imaginary mode" may not be the same mode on the IRC surface. Only use InitHess read if the .hess was produced at the same level as this IRC job.

The three options:

Option A — InitHess read (preferred when applicable): Read a pre-computed .hess file from the frequency calculation at the same level of theory. Fastest and most accurate, because the Hessian is already well-converged and you skip recomputing it. Ask the user for the filename.

%irc
  InitHess read
  Hess_Filename "previous_freq.hess"
end

Option B — InitHess calc_anfreq: Compute analytical frequencies at the start of the IRC. Use when no .hess file is available and the method supports analytical Hessians (common DFT functionals, HF). Cleanest and still fairly fast for small-to-medium systems.

%irc
  InitHess calc_anfreq
end

Option C — InitHess calc_numfreq: Compute numerical frequencies. Use when analytical Hessians aren't available or aren't reliable for the chosen method (some double hybrids, multireference, custom-parameter methods). Much slower (scales with 6N SCFs) but works everywhere.

%irc
  InitHess calc_numfreq
end

Decision flow for the skill: ask "Do you have a .hess file from your frequency calculation at this level of theory?" If yes → Option A with filename. If no and the method supports analytical Hessians → Option B. If no and it doesn't → Option C.

IRC Parameters and ORCA Defaults

ORCA 6's documented defaults:

• MaxIter = 20 iterations per arm
• Step = adaptive, governed by Scale_Displ_SD (initial value derived from Init_Displ_DE) and dynamically rescaled — bounded to 1/16× – 4× of the initial value. There is no single fixed step-size default to quote.

The default 20 iterations is often too few to reach a clear minimum, so the recommended override is:

%irc
  MaxIter 70
  PrintLevel 1
  Direction both
end

Do not set Step unless the user specifically asks for a fixed step — the adaptive algorithm is usually better than any constant. If the path behaves badly, the right lever is usually MaxIter, InitHess, or the TS quality, not a fixed step size.

Rationale for the overrides:

• MaxIter 70 — ORCA's default of 20 is almost always too few for a path that reaches a clear minimum. Bump to 100–150 for flat PES regions.
• PrintLevel 1 — prints energy and gradient at each step so the run can be monitored and the path summary is written in full.
• Direction both — follows the path in both directions from the TS. Almost always what you want when confirming connectivity.

Composite Methods

For composite methods like PBEh-3c or r²SCAN-3c, do NOT specify a separate basis set — the basis is built into the method keyword.

! PBEh-3c IRC

This automatically satisfies method/basis consistency since the composite method has a fixed internal basis.

Template

# ORCA IRC Calculation
# Method: {method} {basis}
# TS confirmed by frequency calculation: {num_imaginary} imaginary mode(s)

%pal nprocs {nprocs} end
%maxcore {maxcore}

! {method} {basis} {dispersion} IRC

%irc
  InitHess {inithess_option}
  {hess_filename_line}
  MaxIter 70
  PrintLevel 1
  Direction both
end

%output jsongbwfile True jsonpropfile True end

* xyzfile {charge} {multiplicity} {ts_geometry_file}

When generating the input:

• If the user provides inline coordinates, use * xyz {charge} {mult} followed by the coordinates and closing *.
• Omit {dispersion} if already included in a composite method.
• Omit {basis} for composite methods.
• Omit the Hess_Filename line unless using InitHess read.
• Keep the %output jsongbwfile True jsonpropfile True end line — it makes ORCA emit the JSON files OPI's Output.parse() consumes during post-processing (Stage 2).

About OPI input builders

OPI (pip install orca-pi) ships typed BlockIrc builders that validate keys at construction. However, the catgo backend's irc node already emits its own %irc block from node params (irc_max_iter, irc_direction, init_hess, etc.). Pasting an OPI-built %irc block via extra_blocks would produce two %irc blocks in the same .inp, which is undefined behavior.

For this skill, stick with node params for %irc content and use extra_blocks only for %output. The OPI parsing wins (Stage 2) still apply. If you need a knob BlockIrc exposes that the node params don't (scale_init_displ, tolmaxg, monitor_internals, etc.), open that as a node-def gap rather than dual-emitting blocks.

Submitting to HPC (Expanse) — proven flow

Use this when the user wants the CatGo workflow engine to run the IRC on Expanse. Skip if they only want the input file.

Use catgo_workflow (graph-based), NOT catgo_workflow_engine (task-based). The graph-based tool auto-captures the viewer structure on create and lets you wire opt → freq → IRC explicitly. Param keys differ from the task-based API: graph-based uses method/basis, task-based uses orca_method/orca_basis.

1. Confirm TS structure is loaded and find session_id

catgo_view(action: "get_state")

curl -s http://localhost:8000/api/hpc/connections

Copy the session_id for host: login.expanse.sdsc.edu.

2. Create the workflow (auto-captures the TS structure)

catgo_workflow(action: "create", name: "IRC from TS")

3. Add the IRC node

catgo_workflow(action: "batch", workflow_id: "<wf_id>", operations: [
  {"op": "add_node", "node_type": "irc", "label": "irc1",
   "params": {
     "software": "orca",
     "method": "B3LYP",
     "basis": "def2-SVP",
     "charge": 0,
     "multiplicity": 1,
     "dispersion": "D3BJ",
     "max_iterations": 70
   }},
  {"op": "connect", "from_id": "<structure_input_id>", "to_id": "irc1",
   "from_handle": "structure", "to_handle": "structure"}
])

Canonical IRC node params (what the engine actually reads)

Parameter	Default	Description
method	r2SCAN-3c	DFT functional
basis	6-31G	Basis set (omit for composite methods)
charge / multiplicity	0 / 1
dispersion	(none)	D4 | D3BJ | D3 — use this field, NOT extra_keywords
grid	DefGrid2	DefGrid1/2/3
max_iterations	30	per-arm iteration cap (legacy alias: max_irc_iterations)
num_cores / max_core_mb	4 / 4000	%pal nprocs / %maxcore

⚠️ Hessian-source / step-control params (init_hess, irc_direction, Init_Displ_DE, etc.) are not currently read by the IRC engine code path (server/workflow/engines/orca.py:218). The engine emits a minimal %irc MaxIter ... end block from max_iterations only. To control Hessian source, direction, or initial displacement, either (a) extend the node-def, or (b) skip the workflow engine and run ORCA directly on a hand-edited input file using the Stage 1 template above.

4. Run with the full HPC run_config

IRC chains many SCFs along the path — bump walltime and prefer shared/compute over debug. Read server/templates/orca_generic.sh and pass its contents as default_template.

catgo_workflow(action: "run", workflow_id: "<wf_id>", run_config: {
  "execution_mode": "hpc",
  "default_session_id": "<expanse_session_id>",
  "base_work_dir": "/expanse/lustre/projects/sdp126/jyang25/ORCA/catgo",
  "default_job_params": {
    "nodes": 1, "ntasks": 8, "cpus_per_task": 1,
    "walltime": "04:00:00", "partition": "shared"
  },
  "cluster_configs": {
    "<expanse_session_id>": {
      "account": "sdp126",
      "partition": "shared",
      "module_loads": "module load cpu/0.17.3b\nmodule load gcc/10.2.0/npcyll4\nexport PATH=$HOME/openmpi-4.1.8/bin:$PATH\nexport LD_LIBRARY_PATH=$HOME/openmpi-4.1.8/lib:$LD_LIBRARY_PATH",
      "orca_dir": "/home/jyang25/orca_6_1_1_RRP8",
      "default_template": "<contents of server/templates/orca_generic.sh>",
      "default_job_params": {
        "nodes": 1, "ntasks": 8, "cpus_per_task": 1,
        "walltime": "04:00:00", "partition": "shared"
      }
    }
  }
})

The local-scratch template stages I/O to $TMPDIR/orca_$SLURM_JOB_ID and copies results back. Required on Expanse — Lustre kills ORCA's many-small-file I/O during per-step SCFs.

5. Monitor

catgo_workflow(action: "status", workflow_id: "<wf_id>")

6. Pull files for post-processing

When status is COMPLETED, pull the output and trajectory files referenced in Stage 2 below:

mkdir -p ./local_run
for f in ORCA.out ORCA_IRC.dat ORCA_IRC.xyz ORCA_IRC_F.xyz ORCA_IRC_B.xyz \
         ORCA.property.json ORCA.json; do
  curl -s -X POST http://localhost:8000/api/hpc/files/read-content \
    -H 'Content-Type: application/json' \
    -d "{\"session_id\":\"<expanse_session_id>\",\"file_path\":\"<work_dir>/$f\"}" \
    > ./local_run/$f
done

Then run the parser/plotter from Stage 2 against the local copy.

Submission gotchas

• catgo_workflow_engine.add_task doesn't auto-attach the viewer structure → "No input structure provided".
• partition=workq (Shaheen default) is invalid on Expanse → use debug/shared/compute.
• partition=debug capped at 30 min — IRC almost always needs more.
• Missing account=sdp126 → "Invalid account or account/partition combination".
• Missing module_loads + orca_dir → orca not on PATH; per-step SCFs silently produce nothing.
• After re-connecting to Expanse, the session_id changes — re-discover via /api/hpc/connections and update both default_session_id and the cluster_configs key.
• A .hess from a different method/basis/grid is a silent footgun for init_hess: "read" — only reuse a .hess produced at exactly this IRC's level of theory.

Stage 2: Post-Processing

After the IRC completes, ORCA produces:

File	Contents
{basename}.out	Main output with per-step energies/gradients and the IRC PATH SUMMARY table
{basename}_IRC.dat	Clean tabular file: step, energy (Eh), path length — easiest programmatic source when present
{basename}_IRC.xyz	Full trajectory: concatenated geometries from backward end → TS → forward end
{basename}_IRC_F_trj.xyz	Forward-arm trajectory
{basename}_IRC_B_trj.xyz	Backward-arm trajectory
{basename}_IRC_F.xyz	Final endpoint geometry (forward direction)
{basename}_IRC_B.xyz	Final endpoint geometry (backward direction)

The forward/backward labeling is arbitrary — it depends on the sign of the initial displacement. The user identifies which endpoint is reactant vs. product by looking at the geometries.

Parsing the Energy Profile

Parse the IRC PATH SUMMARY table in the main output file. OPI does not model this block as a typed object, but its grepper recipe replaces the bespoke "find marker, slice 10 KB, walk lines" boundary handling with a single call. Requires pip install orca-pi.

Each row has the shape:

<step>   <E_Eh>   <dE_kcal/mol>   <max_grad>   <rms_grad>   [<= TS]

Reference parser using OPI's grepper:

import re
import sys
sys.path.insert(0, ".claude/skills")  # for the _shared helper
from _shared.orca_opi import grep_block

ROW_PATTERN = re.compile(
    r"(\d+)\s+([-\d.]+)\s+([-\d.]+)\s+([\d.]+)\s+([\d.]+)\s*(.*)"
)


def parse_irc_path(output_file):
    # offset=2 skips the header + ruler; count is a safe upper bound.
    lines = grep_block(output_file, "IRC PATH SUMMARY", offset=2, count=400)

    steps = []
    ts_raw_step = None
    for line in lines:
        line = line.strip()
        if not line:
            continue
        m = ROW_PATTERN.match(line)
        if not m:
            continue
        raw_step = int(m.group(1))
        is_ts = "<= TS" in m.group(6)
        if is_ts:
            ts_raw_step = raw_step
        steps.append({
            "step":          raw_step,
            "energy_eh":     float(m.group(2)),
            "dE_kcal_mol":   float(m.group(3)),
            "max_gradient":  float(m.group(4)),
            "rms_gradient":  float(m.group(5)),
            "is_ts":         is_ts,
        })

    # Re-index so backward steps are negative, TS = 0, forward positive
    if ts_raw_step is not None:
        for s in steps:
            s["step"] -= ts_raw_step
    return steps


# Sanity check before parsing — replaces hand-grep of "ORCA TERMINATED NORMALLY"
from opi.output.grepper.recipes import has_terminated_normally
assert has_terminated_normally(output_file)

Plotting the Energy Profile

Plotting matches the style used by CatGo's IrcPathPlot.svelte: ΔE in kcal/mol on the y-axis, IRC step number (re-indexed relative to the TS) on the x-axis, backward arm in purple, forward arm in green, TS marked with a red point and an amber dashed vertical line.

import matplotlib.pyplot as plt


BACKWARD_COLOR = "#8b5cf6"
FORWARD_COLOR = "#10b981"
TS_COLOR = "#ef4444"
TS_GUIDE = "#f59e0b"


def plot_irc_profile(steps, output_png="irc_profile.png"):
    xs = [s["step"] for s in steps]
    ys = [s["dE_kcal_mol"] for s in steps]

    ts_idx = next((i for i, s in enumerate(steps) if s["is_ts"]), None)

    fig, ax = plt.subplots(figsize=(8, 5))

    if ts_idx is not None:
        ax.plot(xs[: ts_idx + 1], ys[: ts_idx + 1], color=BACKWARD_COLOR,
                linewidth=1.5, marker="o", markersize=3, label="Backward")
        ax.plot(xs[ts_idx:], ys[ts_idx:], color=FORWARD_COLOR,
                linewidth=1.5, marker="o", markersize=3, label="Forward")
        ax.scatter([xs[ts_idx]], [ys[ts_idx]], color=TS_COLOR, s=60,
                   zorder=5, label="TS")
        ax.axvline(xs[ts_idx], color=TS_GUIDE, linestyle="--",
                   linewidth=1, alpha=0.7)
    else:
        ax.plot(xs, ys, color=FORWARD_COLOR, linewidth=1.5, marker="o",
                markersize=3)

    ax.set_xlabel("IRC Step (relative to TS)")
    ax.set_ylabel("ΔE (kcal/mol)")
    ax.legend(loc="best", frameon=False)
    fig.tight_layout()
    fig.savefig(output_png, dpi=150)
    print(f"Profile saved to {output_png}")


if __name__ == "__main__":
    import sys
    steps = parse_irc_path(sys.argv[1])
    out = sys.argv[2] if len(sys.argv) > 2 else "irc_profile.png"
    plot_irc_profile(steps, out)

Viewing the IRC profile in the IDE

After plot_irc_profile(...) writes the PNG, surface it inline with the shared helper:

sys.path.insert(0, ".claude/skills")
from _shared.orca_opi import show_png
show_png("irc_profile.png", "IRC energy profile")
# prints `![IRC energy profile](irc_profile.png)`

Then reply to the user with that markdown link so Claude Code renders the figure inline in chat.

Notes:

• dE_kcal_mol in the PATH SUMMARY is already relative to the TS, so TS sits at ΔE = 0 and both arms descend to negative values. No extra conversion needed.
• If the user wants path length on the x-axis instead of step number, use {basename}_IRC.dat as the source — it's cleaner than regex-parsing the main output.

Extracting Endpoint Geometries

Read {basename}_IRC_F.xyz and {basename}_IRC_B.xyz to show the user the endpoint structures. Remind the user that these endpoint geometries are not fully optimized — they should run a follow-up geometry optimization on each endpoint to get proper reactant and product structures.

Common Issues and Troubleshooting

IRC terminates early / doesn't reach a minimum: Increase MaxIter — try 100 or 150 for flat PES regions. Recompute with InitHess calc_numfreq if the initial direction looked wrong.

IRC follows the wrong mode: The initial Hessian is bad or the TS has extra low-lying imaginary frequencies. Use InitHess read with a .hess file from a converged frequency job at the same level, or switch from calc_anfreq to calc_numfreq. If the TS has more than one imaginary frequency, the geometry is not a true TS — fix it first.

IRC connects to unexpected structures: The TS probably corresponds to a different reaction than intended. Common with NEB-TS when the algorithm finds a stepwise mechanism instead of a concerted one. Re-run the TS search with a different guess geometry or a tighter NEB.

Energy goes up instead of down: Sign of a bad Hessian or the geometry not being a true TS. Recompute the Hessian with NumFreq and verify exactly one imaginary mode before rerunning IRC.

Solvation: Gas-phase IRC is standard for connectivity confirmation. If the user needs a solvated energy profile, single-point CPCM calculations on the gas-phase IRC geometries are more practical than running the full IRC with CPCM.