Nudged Elastic Band (NEB)¶

vibe-qc is growing a native Nudged Elastic Band driver for finding minimum-energy paths between reactant and product geometries. It is the right transition-state finder for systems where the reactive event is delocalised — most notably for surface catalysis (the flagship target is N₂ dissociation on Fe(100)), where a cluster model is a poor approximation of the actual catalytic surface and single-point TS searches struggle.

Status¶

NEB is landing across five increments on main. This page tracks what is shippable today; the rest is signposted as roadmap.

Increment	Status	Public API surface
1. Interpolators + dataclasses	Available	`interpolate_linear`, `interpolate_idpp`, `NEBImage`, `NEBPath`
2. NEB driver (improved tangent + spring)	Available	`run_neb(...)`, `NEBResult`
3. Climbing image (CI-NEB)	Available	`run_neb(..., climbing_image=True)`
4. Periodic dispatch (BIPOLE + FD gradient)	Available	`run_neb(periodic_r, periodic_p, ..., kpoints=...)`
5. QVF `reaction.path` emitter	Available	`result.write_qvf(stem)`
+ Density warm-start (HANDOVER_PERIODIC_NEB M4)	Available	`run_neb(..., warm_start=True)`

Increment 1 — path-construction primitives (available today)¶

Two interpolators construct an initial chain of images between a reactant and a product geometry. They do not run any SCF — they operate purely on atomic positions.

Linear Cartesian¶

from vibeqc import Atom, Molecule, interpolate_linear

reactant = Molecule([Atom(1, [0.0, 0.0, 0.0]),
                     Atom(1, [0.0, 0.0, 1.4])])
product  = Molecule([Atom(1, [0.0, 0.0, 0.0]),
                     Atom(1, [0.0, 0.0, 4.0])])

path = interpolate_linear(reactant, product, n_images=5)
#   len(path) == 7   (5 intermediate + the two endpoints)

Endpoints are returned as the original objects (no copy), so path[0] is reactant and path[-1] is product.

IDPP (Image-Dependent Pair Potential)¶

Linear interpolation between two bonded structures often places atoms inside each other in the intermediate images. IDPP (Smidstrup, Pedersen, Stokbro, Jónsson 2014) constructs a smooth chain by minimising

S^(i)(R) = sum_{j<k}  (d^(i)_{jk} - r_{jk}(R))^2 / r_{jk}(R)^4

for each intermediate image i, where d^(i)_{jk} is the linear interpolation between reactant and product pair-distance matrices. The 1/r⁴ weighting drives images away from atom-atom clashes while still tracking the interpolated distance manifold.

from vibeqc import interpolate_idpp

path = interpolate_idpp(reactant, product, n_images=7)

IDPP is the recommended default starting path for NEB whenever reactant and product differ in bonding connectivity.

Dataclasses¶

from vibeqc import NEBImage, NEBPath

images = [NEBImage(system=s) for s in path]
neb = NEBPath(images=images, spring_constant=0.1)
neb.n_images        # 7 (intermediate + endpoints)
neb.n_intermediate  # 5
neb.energies()      # np.ndarray of NaN (energies set by the driver)

NEBImage carries (system, energy, gradient, tangent) — the energy/gradient/tangent slots stay None after interpolation; they are populated by the NEB driver (Increment 2+).

Periodic systems¶

Both interpolators accept PeriodicSystem. The lattice must match between reactant and product (variable-cell NEB is out of scope — fix the cell to the reactant’s). Today the periodic IDPP uses Cartesian (non-minimum-image) pair distances, which is correct for the slab + adsorbate target where no atom-pair wraps the cell, but will produce nonsense paths if reactant and product differ by an atom hopping across the PBC. Minimum-image distances will land together with the periodic NEB driver.

Increment 2 — NEB driver (available today, molecular only)¶

run_neb runs an improved-tangent NEB end to end: it computes per-image SCFs (in parallel via joblib), assembles the NEB force (tangent + spring + perpendicular-true-force projection), and optimises with a damped quick-min outer loop (the same scheme ASE’s MDMin uses).

from vibeqc import Atom, Molecule, run_neb

# H + H₂ → H₂ + H, collinear.
reactant = Molecule(
    [Atom(1, [0.0, 0.0, 0.0]),
     Atom(1, [0.0, 0.0, 1.4]),
     Atom(1, [0.0, 0.0, 4.4])],
    0, 2,  # neutral doublet
)
product = Molecule(
    [Atom(1, [0.0, 0.0, 0.0]),
     Atom(1, [0.0, 0.0, 3.0]),
     Atom(1, [0.0, 0.0, 4.4])],
    0, 2,
)
result = run_neb(
    reactant, product,
    basis="sto-3g",
    n_images=5,
    method="UHF",
    spring_constant=0.1,
    interpolation="idpp",
    max_iter=80,
    conv_tol_force=2e-3,
)
print(result.converged, result.n_iter, result.max_force)
print("TS energy:", result.energies[result.transition_state_index])
ts = result.path.images[result.transition_state_index].system

NEBResult carries the final path, the per-image energies, a boolean converged, the index of the highest-energy intermediate image (transition_state_index), the number of outer iterations, and the final maximum NEB force magnitude.

Parallelism¶

Per-image SCFs run inside a joblib.Parallel block per outer iteration. n_jobs=-1 (default) uses every available core; pass n_jobs=1 for serial. The endpoints are evaluated once and cached — they don’t move while free_endpoints=False (the default).

Frozen atoms¶

Pass freeze_indices=[i, j, …] to zero the NEB force on the listed atoms; their geometry stays at the endpoint position throughout the optimisation. This is implemented at the NEB layer (the per-image SCF still sees the full atom set, we just zero the force components before the quick-min step). Useful for keeping slab-substrate atoms fixed during surface NEB once Increment 4 lands.

Method support¶

method="RHF" | "UHF" | "RKS" | "UKS" — same set as optimize_molecule. Pass the matching *_options (e.g. uhf_options=UHFOptions()) to tune the SCF; default options work for the textbook H + H₂ benchmark but more demanding systems will likely want a higher max_iter on the SCF + maybe a tighter conv_tol_grad.

Algorithmic notes¶

Improved tangent. The tangent at each intermediate image is determined by the energy ordering of its neighbours (eqs. 8–11 of Henkelman & Jónsson 2000). When the central image is the local energy extremum, the tangent is a weighted combination of the forward and backward differences with the magnitude of the larger energy difference upfront.
Outer loop. Quick-min (damped MD) with adaptive step. Each intermediate image carries a velocity; per step the velocity is projected onto the NEB force and zeroed if the projection is negative. Step size grows by 1.1 when aligned and resets on direction flips. Equivalent in spirit to ASE’s MDMin. L-BFGS-B on the concatenated coordinate would also work but assumes the NEB force is conservative — it isn’t (the spring contribution depends only on the parallel projection of the displacement), so quick-min is the standard choice.
Slow convergence tail. Near convergence the geometry can be essentially locked in (TS energy stable to mHa) while the spring-component max-force decays slowly. This is intrinsic to the quick-min / MDMin scheme; for tighter convergence increase max_iter, or pre-relax the band with a coarser conv_tol_force and finish the saddle with a quasi-Newton TS search.

Increment 3 — climbing image (available today)¶

Pass climbing_image=True to run_neb to promote the highest- energy intermediate image to a “climbing” image partway through the optimisation:

result = vibeqc.run_neb(
    reactant, product, basis="sto-3g",
    n_images=5, method="UHF",
    climbing_image=True,
    climbing_image_start_fraction=0.3,  # default; warm-up fraction
    conv_tol_force=1e-3, max_iter=80,
)
ts_image = result.path.images[result.path.climbing_image_index]

The CI image’s spring contribution is removed and the tangent- parallel component of its true force is inverted, so it climbs uphill along the path to the saddle while still relaxing perpendicular. Other images keep standard NEB dynamics, which keeps the path itself anchored at the right shape. Henkelman, Uberuaga, Jónsson 2000 is the original reference.

Warm-up. The highest-energy intermediate image is selected after climbing_image_start_fraction * max_iter iterations of plain NEB and then fixed for the remainder of the run. The warm-up keeps the climbing selection from flipping between iterations before the band has found its rough shape; the selection lock keeps the climber from losing momentum to a later re-selection. Default fraction is 0.3.
Accuracy. On the textbook H + H₂ → H₂ + H benchmark the climbing image lands at the symmetric saddle to ≤ 1e-4 bohr — six orders of magnitude tighter than the plain-NEB highest- energy image (which Inc 2’s test bounds at 0.05 bohr).
NEBResult.path.climbing_image_index records which path image was promoted; transition_state_index then equals that index. Both are None on a plain (non-climbing) run.

The climbing image is the recommended saddle estimator from run_neb: it’s accurate enough to seed a follow-up quasi-Newton TS search (or in many cases to just report directly as the saddle).

Increment 4 — periodic dispatch (available today, BIPOLE + FD gradient)¶

run_neb accepts PeriodicSystem endpoints. The path is laid out exactly as in the molecular case (improved-tangent forces, quick-min outer loop, optional climbing image); per-image SCFs go through run_pbc_bipole_rhf/uhf/rks/uks and per-image gradients through a central-difference fallback (the J^LR reciprocal-Ewald contribution is still missing from the analytic BIPOLE gradient — see HANDOVER_PERIODIC_NEB.md).

import numpy as np
from vibeqc import (
    Atom, PeriodicSystem, run_neb,
    PeriodicRHFOptions, LatticeSumOptions,
)

L = np.diag([10.0, 10.0, 10.0])
reactant = PeriodicSystem(3, L, [
    Atom(1, [0.0, 0.0, 0.0]),
    Atom(1, [0.0, 0.0, 1.4]),
])
product = PeriodicSystem(3, L, [
    Atom(1, [0.0, 0.0, 0.0]),
    Atom(1, [0.0, 0.0, 1.8]),
])

# Periodic SCF options. Use the same options object the
# non-NEB BIPOLE drivers take.
rhf_opts = PeriodicRHFOptions()
rhf_opts.max_iter = 50
result = run_neb(
    reactant, product,
    basis="sto-3g",
    n_images=5,
    method="RHF",
    rhf_options=rhf_opts,
    kpoints=(2, 2, 2),       # Monkhorst-Pack mesh; passed to monkhorst_pack
    fd_step_bohr=1e-3,       # central-difference half-step
    interpolation="idpp",
    max_iter=20,
    conv_tol_force=2e-3,
)
ts = result.path.images[result.transition_state_index].system

The FD-gradient surface. Per-image gradient cost is 6N + 1 BIPOLE SCFs (one per Cartesian degree of freedom, ± displaced, plus the reference). For N = 10 atoms × 5 intermediate images × 50 outer iterations that’s ~15 000 BIPOLE SCFs total — slow, but correct in the limit fd_step_bohr → 0. When the J^LR derivative lands in the analytic BIPOLE gradient, the per-image cost drops to 2 SCFs (one for the energy + one for the analytic gradient).

kpoints accepts either a 3-tuple of ints (converted internally via vibeqc.monkhorst_pack) or a pre-built BlochKMesh. None ⇒ Γ-only mesh (only useful as a sanity-check; pick a real mesh for production).

Same lattice for both endpoints. Variable-cell NEB is out of scope — fix the cell to the reactant’s lattice. Endpoint lattices that don’t match raise.

Frozen substrate atoms. The same freeze_indices mechanism works for surface NEB: pass the indices of the substrate atoms you want to keep fixed. The frozen-atom forces are zeroed at the NEB layer; the BIPOLE SCF still sees the full system.

Mixing molecular + periodic endpoints raises at the type-dispatch guard inside run_neb. The endpoints must be the same type.

Increment 5 — QVF reaction.path emitter (available today)¶

NEBResult carries a write_qvf(stem) method that emits a vibe-view-renderable archive of the converged (or last-evaluated) path:

result = vibeqc.run_neb(reactant, product, basis="sto-3g", ...)
qvf_path = result.write_qvf("my_neb_run")
#   → Path('my_neb_run.qvf')

The archive contains three sections:

structure — the reactant geometry (single-frame fallback for viewers that don’t yet animate reaction.path).
reaction.path — per-image coords + energies + reaction coordinate (cumulative arc length normalised to 0–1) + waypoint annotations (reactant / product / transition_state). The TS waypoint points at the climbing image when CI-NEB ran, else at the highest-energy intermediate.
citations — BibTeX assembled with uses_neb=True (plus uses_ci_neb=True when the run was climbing-image). Always includes Henkelman+Jónsson 2000 (improved tangent) + Smidstrup 2014 (IDPP); CI-NEB adds Henkelman+Uberuaga+Jónsson 2000.

For periodic NEB the archive ships as QVF v2 automatically: the writer detects periodic frames and emits the per-frame lattice

dim on the reaction.path section so vibe-view can draw the unit cell. See the “Periodic reaction paths (QVF v2)” subsection of the vibe-view guide for the schema details. Inc B (vibe-view renderer track) will add the cell-rendering + atom-wrapping logic to the renderer itself.

Density warm-start (available today)¶

run_neb reuses each image’s converged SCF density across outer iterations: the density at outer iter N is fed in as the SCF initial guess for the same image at outer iter N+1. The geometry change between iterations is typically small, so the SCF converges in many fewer iterations than from a cold SAD/Hcore guess.

result = vq.run_neb(
    reactant, product, basis="def2-svp",
    method="UKS", functional="pbe",
    n_images=7,
    warm_start=True,   # default; pass False to force cold SCFs (benchmarking)
)

The behaviour is bit-exact vs. cold-start — warm-start changes only the initial guess, not the converged density. The NEB trajectory (energies + forces + per-image positions) is identical at any fixed max_iter budget regardless of whether warm_start is on.

Coverage today¶

Method	Molecular	Periodic (BIPOLE)
RHF	yes	yes
UHF	yes	yes
RKS	yes	API-ready *
UKS	yes	API-ready *

* Periodic RKS / UKS warm-start kwargs (initial_density=, init_alpha=/init_beta=) are wired through run_pbc_bipole_rks / run_pbc_bipole_uks, but end-to-end testing is blocked by an unrelated upstream build_xc_periodic arg-order bug in those drivers’ SCF loops. When that bug is fixed (separate chat) the NEB dispatch automatically picks them up — no further changes needed in run_neb. Until then, periodic NEB on KS-DFT methods is limited to the existing cold-start behaviour (the warm-start kwarg is a harmless no-op on those methods today).

Observed speedups¶

Molecular UKS on the H + H₂ → H₂ + H benchmark (3 intermediate images, 10 outer iters): ~2× wall-time reduction (~84 s → ~42 s). Open-shell DFT SCFs are the biggest win — expensive per iteration AND slow to converge from SAD/Hcore on a doublet manifold.
Periodic UHF on H₃-in-cubic-box (single SCF call): 76 SCF iters → 25 SCF iters, ~3× speedup.
Molecular RHF / RKS / periodic RHF: speedups in the flat-to-modest range on small systems — the cold SCF already converges in 3-5 iters at the test scale. Real surface NEB workloads should see 2-4× the maintainer originally scoped in HANDOVER_PERIODIC_NEB.md § M4.

Bonus: FD-gradient inner warm-start (periodic only)¶

Periodic NEB computes per-image gradients by central differences (6N + 1 SCFs per image per outer iteration; see Increment 4). Each of those 6N displaced SCFs additionally warm-starts from the reference SCF’s converged density at the same image — geometrically small perturbations (default fd_step_bohr=1e-3) are essentially the same electronic structure, so the displaced SCFs converge in a couple of iterations instead of from cold. This stacking multiplies the warm-start savings for periodic NEB specifically.

When to turn warm-start off¶

Benchmarking the cold-start cost (warm_start=False).
Paranoid bisection when an unrelated convergence problem is suspected and you want to remove warm-start as a variable.

There’s no correctness reason to disable warm-start at production time — it’s bit-exact in both directions.

What is not in vibe-qc today¶

The dedicated periodic-NEB chat (HANDOVER_PERIODIC_NEB.md) is the place where the J^LR analytic-gradient + periodic-NEB- specific optimisations get the focused attention they need. The run_neb(PeriodicSystem, ...) path landed in Inc 4 + density warm-start lands the most impactful pieces; the dedicated chat will close the remaining gap (FD → analytic gradient).
The vibe-view renderer’s lattice-box + atom-wrap logic for periodic reaction paths (QVF “Inc B”). The archive already carries everything the renderer needs; the rendering pass is scheduled separately.
Frame-to-frame anchored atom wrap in vibe-view’s renderer. The current modulo-1 wrap can produce a visual “jump” when an atom crosses a cell boundary mid-animation. Rare in practice for chemistry NEBs; cosmetic.

For TS finders that ship today see the geometry-optimization + Hessian tooling. The relaxed-scan workflow (relaxed_scan / ScanResult) can also provide an interpolated starting path that NEB consumes.

Citations¶

Both NEB papers fire automatically whenever run_neb runs (via the routes.drivers.neb route registered against uses_neb=True):

Henkelman, G.; Jónsson, H. Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points. J. Chem. Phys. 113, 9978 (2000). doi:10.1063/1.1323224

Smidstrup, S.; Pedersen, A.; Stokbro, K.; Jónsson, H. Improved initial guess for minimum energy path calculations. J. Chem. Phys. 140, 214106 (2014). doi:10.1063/1.4878664

For climbing_image=True, additionally:

Henkelman, G.; Uberuaga, B. P.; Jónsson, H. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J. Chem. Phys. 113, 9901 (2000). doi:10.1063/1.1329672

The matching routes.drivers.ci_neb route fires in addition to the base routes.drivers.neb route when CI-NEB is enabled.