"""Nudged Elastic Band (NEB) — minimum energy path finder.
Public surface
==============
Path construction (Increment 1):
* :class:`NEBImage`, :class:`NEBPath`
* :func:`interpolate_linear`, :func:`interpolate_idpp`
Driver (Increment 2 — molecular only):
* :func:`run_neb` — improved-tangent NEB with parallel per-image
SCFs and a quick-min outer loop on the concatenated NEB force.
Climbing image is Increment 3; periodic dispatch is Increment 4.
* :class:`NEBResult` — converged path + energies + transition-state
index + iteration count.
Both interpolators accept either :class:`vibeqc.Molecule` or
:class:`vibeqc.PeriodicSystem`. Periodic interpolation today uses
Cartesian (non-minimum-image) pair distances; that is correct for
the slab + adsorbate flagship 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. MIC support will land
with the periodic NEB driver.
The driver currently rejects :class:`PeriodicSystem` endpoints —
periodic NEB is Increment 4, gated on the J^LR (reciprocal Ewald)
contribution being added to the BIPOLE analytic gradient (today's
periodic analytic gradient is incomplete; converging a NEB against
it would mask saddle-point errors). See
``HANDOVER_PERIODIC_NEB.md`` for the dedicated periodic-NEB chat.
References
==========
* Henkelman & Jónsson, "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, Pedersen, Stokbro, Jónsson,
"Improved initial guess for minimum energy path calculations",
J. Chem. Phys. 140, 214106 (2014). doi:10.1063/1.4878664.
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import Any, List, Optional, Sequence, Union
import numpy as np
from scipy.optimize import minimize
from ._vibeqc_core import Atom, Molecule, PeriodicSystem
System = Union[Molecule, PeriodicSystem]
# ---------------------------------------------------------------------------
# Dataclasses
# ---------------------------------------------------------------------------
@dataclass
class NEBImage:
"""One image on a NEB path.
``energy``, ``gradient``, and ``tangent`` are populated by the
NEB driver during optimisation. In Increment 1 they are all
``None`` for freshly interpolated images.
"""
system: System
energy: Optional[float] = None
gradient: Optional[np.ndarray] = None # (n_atoms, 3) Ha/bohr
tangent: Optional[np.ndarray] = None # (n_atoms, 3) unit-norm
@dataclass
class NEBPath:
"""A NEB path as an ordered list of images.
First and last entries are the reactant and product respectively
(fixed by default during optimisation). ``spring_constant`` is in
Ha/bohr². ``climbing_image_index`` is None for standard NEB; the
CI-NEB driver sets it to the highest-energy intermediate after
the warm-up phase.
"""
images: List[NEBImage]
spring_constant: float = 0.1
climbing_image_index: Optional[int] = None
@property
def n_images(self) -> int:
return len(self.images)
@property
def n_intermediate(self) -> int:
return max(0, len(self.images) - 2)
def energies(self) -> np.ndarray:
"""Per-image energies as an array; NaN where unset."""
return np.array(
[
img.energy if img.energy is not None else np.nan
for img in self.images
],
dtype=float,
)
# ---------------------------------------------------------------------------
# Geometry helpers
# ---------------------------------------------------------------------------
def _atom_iter(system: System):
return system.atoms if isinstance(system, Molecule) else system.unit_cell
def _positions_of(system: System) -> np.ndarray:
return np.array([list(a.xyz) for a in _atom_iter(system)], dtype=float)
def _atomic_numbers_of(system: System) -> np.ndarray:
return np.array([int(a.Z) for a in _atom_iter(system)], dtype=int)
def _rebuild_with_positions(template: System, positions: np.ndarray) -> System:
"""Return a copy of ``template`` with Cartesian ``positions`` (bohr)."""
if isinstance(template, Molecule):
new_atoms = [
Atom(int(a.Z), list(p))
for a, p in zip(template.atoms, positions)
]
return Molecule(new_atoms, template.charge, template.multiplicity)
new_atoms = [
Atom(int(a.Z), list(p))
for a, p in zip(template.unit_cell, positions)
]
return PeriodicSystem(
template.dim,
np.asarray(template.lattice, dtype=float),
new_atoms,
charge=template.charge,
multiplicity=template.multiplicity,
)
def _check_compatible(reactant: System, product: System) -> None:
if type(reactant) is not type(product):
raise ValueError(
"reactant and product must be the same system type; got "
f"{type(reactant).__name__} and {type(product).__name__}"
)
zr = _atomic_numbers_of(reactant)
zp = _atomic_numbers_of(product)
if zr.shape != zp.shape or not np.array_equal(zr, zp):
raise ValueError(
"reactant and product must have matching atomic-number "
"sequences in the same order (NEB does not reorder atoms; "
"pre-align if needed)."
)
if isinstance(reactant, PeriodicSystem):
lr = np.asarray(reactant.lattice, dtype=float)
lp = np.asarray(product.lattice, dtype=float)
if not np.allclose(lr, lp):
raise ValueError(
"reactant and product must share the same lattice — "
"variable-cell NEB is out of scope. Fix the cell to the "
"reactant's lattice before constructing endpoints."
)
# ---------------------------------------------------------------------------
# Linear interpolation
# ---------------------------------------------------------------------------
def interpolate_linear(
reactant: System,
product: System,
n_images: int,
) -> List[System]:
"""Linear Cartesian interpolation between two endpoints.
Parameters
----------
reactant, product
Same-type Molecule or PeriodicSystem with matching atom
ordering. For PeriodicSystem the lattices must agree.
n_images
Number of *intermediate* images. The returned list has length
``n_images + 2`` (endpoints included).
Returns
-------
list[System]
``[reactant, img_1, ..., img_n, product]``. Endpoints are
returned as the original objects (not copies).
"""
if n_images < 0:
raise ValueError(f"n_images must be >= 0; got {n_images}")
_check_compatible(reactant, product)
r0 = _positions_of(reactant)
r1 = _positions_of(product)
out: List[System] = [reactant]
for k in range(1, n_images + 1):
t = k / (n_images + 1)
pos = (1.0 - t) * r0 + t * r1
out.append(_rebuild_with_positions(reactant, pos))
out.append(product)
return out
# ---------------------------------------------------------------------------
# IDPP interpolation — Smidstrup et al. 2014
# ---------------------------------------------------------------------------
#
# For each intermediate image i = 1..n_images we define a target
# pair-distance matrix
#
# d^(i)_{jk} = (1 - t_i) * d^(R)_{jk} + t_i * d^(P)_{jk}
#
# (linear in pair-distance space between reactant and product). We
# then minimise the image-dependent pair-potential objective
#
# S^(i)(R) = sum_{j<k} (d^(i)_{jk} - r_{jk}(R))^2 / r_{jk}(R)^4
#
# over each image's Cartesian coordinates independently. The 1/r^4
# weighting drives images strongly to relieve close atom-atom
# contacts while still tracking the interpolated distance manifold.
# Analytic gradient is provided to L-BFGS-B for speed.
def _pair_distance_matrix(positions: np.ndarray) -> np.ndarray:
diff = positions[:, None, :] - positions[None, :, :]
return np.sqrt(np.einsum("ijc,ijc->ij", diff, diff))
def _idpp_value_and_grad(
flat: np.ndarray,
target: np.ndarray,
n_atoms: int,
) -> tuple[float, np.ndarray]:
positions = flat.reshape(n_atoms, 3)
diff = positions[:, None, :] - positions[None, :, :]
r2 = np.einsum("ijc,ijc->ij", diff, diff)
# Off-diagonal mask; we never read the diagonal because it's masked
# out of every aggregation below.
mask = ~np.eye(n_atoms, dtype=bool)
# Replace diagonal r2=0 with 1.0 so divisions are finite — masked
# back to zero before any sum. Also floor off-diagonal r2 at a
# small positive value so an L-BFGS-B step that briefly drives two
# atoms onto each other produces a finite (very large, repulsive)
# gradient instead of NaN — the optimiser can then escape.
_R_FLOOR2 = 1e-6 # bohr²
r2_off = np.maximum(r2, _R_FLOOR2)
safe_r2 = np.where(mask, r2_off, 1.0)
r = np.sqrt(safe_r2)
inv_r2 = 1.0 / safe_r2
inv_r4 = inv_r2 * inv_r2
inv_r5 = inv_r4 / r
delta = target - r # (d - r)
# Energy: 0.5 * sum over all (i,j) of (d-r)^2 / r^4 (symmetric →
# the 1/2 converts to the j<k sum).
pair_energy = delta * delta * inv_r4
energy = 0.5 * float(np.sum(np.where(mask, pair_energy, 0.0)))
# d/dr [(d-r)^2 / r^4] = -2 (d-r) / r^4 - 4 (d-r)^2 / r^5
dS_dr = -2.0 * delta * inv_r4 - 4.0 * delta * delta * inv_r5
dS_dr = np.where(mask, dS_dr, 0.0)
# Chain rule: dr_{jk}/dR_j = diff[j,k] / r[j,k], dr/dR_k = - of that.
# grad[k] = sum_j dS_dr[k,j] / r[k,j] * diff[k,j]
inv_r = np.where(mask, 1.0 / r, 0.0)
factor = dS_dr * inv_r
grad = np.einsum("ij,ijc->ic", factor, diff)
return energy, grad.ravel()
def interpolate_idpp(
reactant: System,
product: System,
n_images: int,
*,
max_iter: int = 1000,
tol: float = 1e-5,
) -> List[System]:
"""Image-Dependent Pair Potential interpolation (Smidstrup 2014).
Builds a linear-Cartesian starting path, then for each
intermediate image independently minimises the IDPP objective
``S^(i)(R) = sum_{j<k} (d^(i)_{jk} - r_{jk}(R))^2 / r_{jk}(R)^4``
with target distances interpolated linearly between reactant and
product pair-distance matrices. Endpoints are returned unchanged.
Cites: Smidstrup, Pedersen, Stokbro, Jónsson,
J. Chem. Phys. 140, 214106 (2014). doi:10.1063/1.4878664.
"""
if n_images < 0:
raise ValueError(f"n_images must be >= 0; got {n_images}")
_check_compatible(reactant, product)
if n_images == 0:
return [reactant, product]
n_atoms = len(_positions_of(reactant))
d_R = _pair_distance_matrix(_positions_of(reactant))
d_P = _pair_distance_matrix(_positions_of(product))
linear_path = interpolate_linear(reactant, product, n_images)
out: List[System] = [reactant]
for k in range(1, n_images + 1):
t = k / (n_images + 1)
target = (1.0 - t) * d_R + t * d_P
x0 = _positions_of(linear_path[k]).ravel()
res = minimize(
_idpp_value_and_grad,
x0,
args=(target, n_atoms),
jac=True,
method="L-BFGS-B",
options={"maxiter": max_iter, "gtol": tol, "ftol": tol},
)
out.append(
_rebuild_with_positions(reactant, res.x.reshape(n_atoms, 3))
)
out.append(product)
return out
# ===========================================================================
# Driver (Increment 2) — improved-tangent NEB, molecular only.
# ===========================================================================
#
# The driver below implements the textbook formulation of NEB:
#
# F_total_i = F_spring_∥_i + F_true_⊥_i
#
# with the improved-tangent estimator of Henkelman & Jónsson 2000
# replacing the original "central-difference" tangent of Mills/
# Jónsson/Schenter 1995 — energy ordering of neighbours decides
# whether τ_i points uphill or downhill, with a transitional weighted
# mix when the central image is the local extremum (eqs. 8-11 of JCP
# 113, 9978, 2000). Spring forces are projected onto τ_i; true
# nuclear gradients are projected *off* τ_i. The result is a
# discretised elastic band that relaxes to the minimum energy path.
#
# Outer loop: damped quick-min (also Henkelman+Jónsson). Each image
# carries a velocity; per step it is projected onto F_total before
# advancing — the projection keeps the band from drifting away from
# the MEP when v has accumulated tangential momentum. L-BFGS-B on
# the concatenated coordinate vector would *also* work, but only if
# the NEB force were a true gradient. It isn't: F_spring_∥ depends on
# the parallel projection of the difference vector, which is not the
# derivative of a scalar potential. Quick-min is the standard choice
# in ASE / VASP / Quantum ESPRESSO for the same reason.
#
# Parallelism: joblib.Parallel over images per outer iteration.
# Endpoints are evaluated once and cached. n_jobs defaults to
# joblib's loky-backend pool size (≈ logical cores); pass n_jobs=1
# for serial.
#
# Out of Increment 2 scope:
# * Climbing image — Increment 3 ("climbing_image=True").
# * Periodic dispatch — Increment 4 (PeriodicSystem endpoints are
# rejected here with a ValueError pointing at HANDOVER_PERIODIC_NEB).
# * SCF density warm-start across NEB steps — Increment 4+ tuning.
#
# References:
# Henkelman, Jónsson, J. Chem. Phys. 113, 9978 (2000).
# doi:10.1063/1.1323224.
[docs]
@dataclass
class NEBResult:
"""Outcome of a :func:`run_neb` run.
Attributes
----------
path
The converged (or last-evaluated) :class:`NEBPath` — endpoints
plus intermediate images with their final geometry, energy,
gradient, and tangent.
energies
Per-image energy array (Hartree), length ``n_images + 2``.
Endpoints are at indices 0 and -1.
converged
True iff the max-norm of the NEB force fell below
``conv_tol_force`` within ``max_iter`` outer iterations.
transition_state_index
Index into ``path.images`` of the highest-energy image — the
best non-climbing estimate of the saddle. ``None`` if energies
are unset.
n_iter
Number of completed outer iterations.
max_force
Final maximum-norm NEB force (Ha/bohr) over intermediate
images. Compare against ``conv_tol_force`` to gauge how close
to converged a non-converged run was.
"""
path: NEBPath
energies: np.ndarray
converged: bool
transition_state_index: Optional[int]
n_iter: int
max_force: float
# Captured at run_neb call time so write_qvf can build the
# citation surface + manifest provenance without the user
# having to repeat the SCF flavour. All four can be None if
# NEBResult is constructed by hand (e.g. in tests).
method: Optional[str] = None
basis: Optional[str] = None
functional: Optional[str] = None
is_periodic: bool = False
[docs]
def write_qvf(
self,
stem: Any,
*,
compression: Optional[int] = None,
) -> Any:
"""Emit a vibe-view ``reaction.path`` QVF archive.
Builds an ``OutputPlan`` + context dict from this result and
delegates to :func:`vibeqc.output.formats.qvf.write_qvf`. The
archive contains a ``structure`` section (reactant geometry),
a ``reaction.path`` section (every image's coords + energies
+ waypoints), and a ``citations`` section (BibTeX assembled
with ``uses_neb=True`` — plus ``uses_ci_neb=True`` when this
result came from a climbing-image run).
For periodic NEB results the archive ships as QVF v2 — the
``reaction.path`` section additionally carries the per-frame
lattice + dim (see ``docs/user_guide/vibe_view.md`` §
"Periodic reaction paths"). The writer detects periodic
frames automatically.
Parameters
----------
stem
Path stem; the writer appends ``.qvf``.
compression
Optional ``zipfile`` compression constant; if ``None``
the writer uses its default (``ZIP_DEFLATED`` or
``ZIP_ZSTANDARD`` when ``zipfile-zstd`` is installed).
Returns
-------
pathlib.Path
The written archive path.
"""
from .output.formats.qvf import write_reaction_path_qvf
n_total = len(self.path.images)
# Waypoints: reactant, product, TS (climbing image when
# available, else the highest-energy intermediate).
waypoints: list[dict[str, Any]] = [
{
"frame_index": 0,
"label": "reactant",
"kind": "reactant",
"energy_eh": float(self.energies[0]),
},
{
"frame_index": n_total - 1,
"label": "product",
"kind": "product",
"energy_eh": float(self.energies[-1]),
},
]
ts_idx = (
self.path.climbing_image_index
if self.path.climbing_image_index is not None
else self.transition_state_index
)
if ts_idx is not None and 0 < ts_idx < n_total - 1:
waypoints.append(
{
"frame_index": int(ts_idx),
"label": "TS",
"kind": "transition_state",
"energy_eh": float(self.energies[ts_idx]),
}
)
# Reaction coordinate: cumulative arc length over frame
# geometries, normalised to [0, 1].
positions = [
_positions_of(img.system) for img in self.path.images
]
arc = [0.0]
for i in range(1, n_total):
arc.append(
arc[-1]
+ float(np.linalg.norm(positions[i] - positions[i - 1]))
)
rc = (
[a / arc[-1] for a in arc]
if arc[-1] > 0.0
else [0.0] * n_total
)
return write_reaction_path_qvf(
stem,
frames=[img.system for img in self.path.images],
energies=[float(e) for e in self.energies],
waypoints=waypoints,
reaction_coordinate=rc,
method=(self.method or "RHF").upper(),
basis=self.basis or "sto-3g",
functional=self.functional,
extra_assemble_kwargs={
"uses_neb": True,
"uses_ci_neb": self.path.climbing_image_index is not None,
},
compression=compression,
)
# --- improved tangent (Henkelman + Jónsson 2000) ---------------------------
def _improved_tangent(
R_prev: np.ndarray,
R_curr: np.ndarray,
R_next: np.ndarray,
E_prev: float,
E_curr: float,
E_next: float,
) -> np.ndarray:
"""Improved tangent τ_i per Henkelman+Jónsson 2000 eq. 8-11.
Returns a unit-norm Cartesian vector with the same shape as the
input positions. Geometry-only fallback (no energy data) is the
original central-difference tangent ``R_next - R_prev``.
"""
tau_plus = R_next - R_curr
tau_minus = R_curr - R_prev
if E_next > E_curr > E_prev:
tau = tau_plus
elif E_next < E_curr < E_prev:
tau = tau_minus
else:
dE_max = max(abs(E_next - E_curr), abs(E_prev - E_curr))
dE_min = min(abs(E_next - E_curr), abs(E_prev - E_curr))
if E_next > E_prev:
tau = tau_plus * dE_max + tau_minus * dE_min
else:
tau = tau_plus * dE_min + tau_minus * dE_max
norm = float(np.linalg.norm(tau))
if norm < 1e-12:
# Degenerate band (all three images coincident); fall back to
# the central-difference tangent. If that's also zero, return
# zero — the outer loop is at a fixed point.
tau = R_next - R_prev
norm = float(np.linalg.norm(tau))
if norm < 1e-12:
return np.zeros_like(tau)
return tau / norm
# --- per-image SCF + gradient (the worker called inside joblib) ------------
def _nuclear_repulsion_molecular(mol: Molecule) -> float:
"""Sum Z_i Z_j / r_ij for a molecule (bohr in, Ha out).
The high-level SCF entry points compute this internally; the
low-level ``run_*_scf_with_jk`` path wants ``E_nuc`` as a scalar
so the NEB driver's warm-start helper computes it here.
"""
atoms = list(mol.atoms)
n = len(atoms)
e = 0.0
for i in range(n):
zi = int(atoms[i].Z)
xi = np.asarray(atoms[i].xyz, dtype=float)
for j in range(i + 1, n):
zj = int(atoms[j].Z)
xj = np.asarray(atoms[j].xyz, dtype=float)
r = float(np.linalg.norm(xi - xj))
if r > 0:
e += zi * zj / r
return e
def _build_scf_common_pieces(
mol: Molecule,
basis: Any,
) -> tuple[np.ndarray, np.ndarray, float, Any]:
"""Shared S / Hcore / E_nuc / JKBuilder construction for the
warm-start path. All four methods (RHF / UHF / RKS / UKS) need
the same building blocks before calling their low-level
``run_*_scf_with_jk`` entry point.
"""
from ._vibeqc_core import (
compute_kinetic,
compute_nuclear,
compute_overlap,
make_direct_jk_builder,
)
S = compute_overlap(basis)
Hcore = compute_kinetic(basis) + compute_nuclear(basis, mol)
e_nuc = _nuclear_repulsion_molecular(mol)
jk = make_direct_jk_builder(basis)
return S, Hcore, e_nuc, jk
def _empty_density() -> np.ndarray:
return np.zeros((0, 0), dtype=np.float64)
def _run_rhf_warm_start(
mol: Molecule,
basis: Any,
options: Any,
initial_density: Optional[np.ndarray],
) -> Any:
"""RHF SCF via the low-level ``run_rhf_scf_with_jk`` entry point.
Used by :func:`_evaluate_image` to seed the SCF from a previous
outer-iteration's converged density (within-image warm-start —
HANDOVER_PERIODIC_NEB.md M4). The high-level ``run_rhf`` doesn't
expose ``initial_density``, so the NEB driver routes RHF through
the lower-level entry that does.
When ``initial_density is None`` the SCF behaves identically to
``run_rhf`` (the C++ binding falls back to a diagonalisation of
Hcore for the initial guess). The cost difference vs. a cold
``run_rhf`` is the Python-side construction of S/Hcore/E_nuc/JK
— a single-pass set of compute_* calls, much cheaper than the
SCF iterations themselves.
"""
from ._vibeqc_core import RHFOptions, run_rhf_scf_with_jk
opts = options if options is not None else RHFOptions()
S, Hcore, e_nuc, jk = _build_scf_common_pieces(mol, basis)
return run_rhf_scf_with_jk(
basis,
mol.n_electrons(),
S,
Hcore,
e_nuc,
jk,
opts,
initial_density=initial_density
if initial_density is not None
else _empty_density(),
)
def _run_uhf_warm_start(
mol: Molecule,
basis: Any,
options: Any,
initial_density: Optional[tuple[np.ndarray, np.ndarray]],
) -> Any:
"""UHF SCF via ``run_uhf_scf_with_jk``. Density cache is the
``(alpha, beta)`` tuple of per-spin density matrices."""
from ._vibeqc_core import UHFOptions, run_uhf_scf_with_jk
opts = options if options is not None else UHFOptions()
S, Hcore, e_nuc, jk = _build_scf_common_pieces(mol, basis)
n_total = mol.n_electrons()
mult = mol.multiplicity
# Standard alpha/beta partition: n_alpha = (N + 2S) / 2,
# n_beta = N - n_alpha. Matches the molecular runner's
# convention.
n_alpha = (n_total + (mult - 1)) // 2
n_beta = n_total - n_alpha
init_a, init_b = (
(_empty_density(), _empty_density())
if initial_density is None
else initial_density
)
return run_uhf_scf_with_jk(
basis,
n_alpha,
n_beta,
S,
Hcore,
e_nuc,
jk,
opts,
init_alpha=init_a,
init_beta=init_b,
)
def _resolve_xc_grid(
mol: Molecule,
grid_options: Any,
) -> Any:
"""Build (or pass through) the XC integration grid for KS-DFT."""
from ._vibeqc_core import GridOptions, build_grid
gopt = grid_options if grid_options is not None else GridOptions()
return build_grid(mol, gopt)
def _run_rks_warm_start(
mol: Molecule,
basis: Any,
options: Any,
functional: Optional[str],
grid_options: Any,
initial_density: Optional[np.ndarray],
) -> Any:
"""RKS SCF via ``run_rks_scf_with_jk``. Same closed-shell density
convention as RHF; additionally needs an XC integration grid."""
from ._vibeqc_core import RKSOptions, run_rks_scf_with_jk
opts = options if options is not None else RKSOptions()
if functional is not None and not opts.functional:
opts.functional = functional
S, Hcore, e_nuc, jk = _build_scf_common_pieces(mol, basis)
grid = _resolve_xc_grid(mol, grid_options)
return run_rks_scf_with_jk(
basis,
mol.n_electrons(),
S,
Hcore,
e_nuc,
jk,
grid,
opts,
initial_density=initial_density
if initial_density is not None
else _empty_density(),
)
def _run_uks_warm_start(
mol: Molecule,
basis: Any,
options: Any,
functional: Optional[str],
grid_options: Any,
initial_density: Optional[tuple[np.ndarray, np.ndarray]],
) -> Any:
"""UKS SCF via ``run_uks_scf_with_jk``. Open-shell α/β densities +
XC grid."""
from ._vibeqc_core import UKSOptions, run_uks_scf_with_jk
opts = options if options is not None else UKSOptions()
if functional is not None and not opts.functional:
opts.functional = functional
S, Hcore, e_nuc, jk = _build_scf_common_pieces(mol, basis)
grid = _resolve_xc_grid(mol, grid_options)
n_total = mol.n_electrons()
mult = mol.multiplicity
n_alpha = (n_total + (mult - 1)) // 2
n_beta = n_total - n_alpha
init_a, init_b = (
(_empty_density(), _empty_density())
if initial_density is None
else initial_density
)
return run_uks_scf_with_jk(
basis,
n_alpha,
n_beta,
S,
Hcore,
e_nuc,
jk,
grid,
opts,
init_alpha=init_a,
init_beta=init_b,
)
def _evaluate_image(
positions: np.ndarray,
template: Molecule,
basis_name: str,
method: str,
*,
functional: Optional[str],
rhf_options: Any,
uhf_options: Any,
rks_options: Any,
uks_options: Any,
gradient_options: Any,
grid_options: Any,
dispersion_params: Any,
initial_density: Optional[np.ndarray] = None,
) -> tuple[float, np.ndarray, Optional[np.ndarray]]:
"""Run SCF + gradient at one geometry.
Returns ``(energy, gradient, converged_density)``. Energies in Ha;
gradients in Ha/bohr, shape (n_atoms, 3). ``converged_density`` is
the converged closed-shell density matrix as a numpy array for
method=``"rhf"``; ``None`` for any other method (warm-start not
yet wired for UHF / RKS / UKS — same low-level
``run_*_scf_with_jk`` entry points exist; landing under the same
pattern is a follow-up).
The molecule is reconstructed from ``positions`` (bohr) using
``template`` for atomic numbers + charge + multiplicity. The basis
set is rebuilt per geometry — vibe-qc's BasisSet is bound to the
nuclei it was constructed with.
When ``initial_density`` is provided (RHF only), the SCF starts
from that density matrix instead of SAD/Hcore — within-image
density warm-start across NEB outer iterations.
"""
from ._vibeqc_core import BasisSet
from .molecular_optimize import _compute_molecular_gradient, _run_molecular_scf
mol = _rebuild_with_positions(template, positions)
basis = BasisSet(mol, basis_name)
method_lower = method.lower()
# The density cache slot is one of:
# - None (no warm-start density available)
# - np.ndarray for closed-shell (RHF / RKS)
# - (alpha, beta) tuple of np.ndarrays for open-shell (UHF / UKS)
# _evaluate_image returns the cache in whichever shape matches
# the method, and the outer loop carries it back per image.
converged_density: Optional[Any] = None
if method_lower == "rhf":
scf_result = _run_rhf_warm_start(
mol, basis, rhf_options, initial_density
)
energy = float(scf_result.energy)
converged_density = np.asarray(
scf_result.density, dtype=np.float64
).copy()
elif method_lower == "uhf":
scf_result = _run_uhf_warm_start(
mol, basis, uhf_options, initial_density
)
energy = float(scf_result.energy)
converged_density = (
np.asarray(scf_result.density_alpha, dtype=np.float64).copy(),
np.asarray(scf_result.density_beta, dtype=np.float64).copy(),
)
elif method_lower == "rks":
scf_result = _run_rks_warm_start(
mol, basis, rks_options, functional, grid_options, initial_density
)
energy = float(scf_result.energy)
converged_density = np.asarray(
scf_result.density, dtype=np.float64
).copy()
elif method_lower == "uks":
scf_result = _run_uks_warm_start(
mol, basis, uks_options, functional, grid_options, initial_density
)
energy = float(scf_result.energy)
converged_density = (
np.asarray(scf_result.density_alpha, dtype=np.float64).copy(),
np.asarray(scf_result.density_beta, dtype=np.float64).copy(),
)
else:
energy, scf_result = _run_molecular_scf(
mol,
basis,
method_lower,
functional=functional,
rhf_options=rhf_options,
uhf_options=uhf_options,
rks_options=rks_options,
uks_options=uks_options,
)
grad = _compute_molecular_gradient(
mol,
basis,
scf_result,
method_lower,
gradient_options=gradient_options,
grid_options=grid_options,
dispersion_params=dispersion_params,
)
if dispersion_params is not None:
from .dispersion import compute_d3bj
disp = compute_d3bj(mol, dispersion_params)
energy = energy + float(disp.energy)
return float(energy), np.asarray(grad, dtype=float), converged_density
def _evaluate_image_periodic(
positions: np.ndarray,
template: PeriodicSystem,
basis_name: str,
method: str,
*,
kmesh: Any,
functional: Optional[str],
rhf_options: Any,
uhf_options: Any,
rks_options: Any,
uks_options: Any,
fd_step_bohr: float,
initial_density: Optional[np.ndarray] = None,
) -> tuple[float, np.ndarray, Optional[np.ndarray]]:
"""Run a periodic BIPOLE SCF + finite-difference gradient.
Returns (energy, gradient). Gradient shape (n_atoms, 3), Ha/bohr,
in Cartesian coordinates — *not* fractional, even though the
underlying SCF is BIPOLE. The NEB outer loop works in Cartesian
space throughout (positions, tangents, spring distances).
The gradient is computed by central-differencing the SCF energy
along each Cartesian degree of freedom. This is the path the
periodic-NEB chat's handover recommends until the J^LR
reciprocal-Ewald contribution lands in the analytic BIPOLE
gradient (see ``python/vibeqc/bipole_gradient.py`` and
``HANDOVER_PERIODIC_NEB.md``). Cost: 6N + 1 BIPOLE SCFs per
image per outer iteration. Correctness is exact in the limit
``fd_step_bohr → 0``.
"""
from ._vibeqc_core import BasisSet
from .bipole_optimize import _run_scf
sys_current = _rebuild_with_positions(template, positions)
basis_current = BasisSet(sys_current.unit_cell_molecule(), basis_name)
method_upper = method.upper()
opts = {
"rhf": rhf_options,
"uhf": uhf_options,
"rks": rks_options,
"uks": uks_options,
}[method.lower()]
# Periodic SCF density warm-start. The outer-loop cache feeds in
# the previous outer iter's converged density at this geometry's
# image, and the 6N FD-displaced SCFs below additionally
# warm-start from the reference SCF's converged density
# (within-image FD speedup). Density shape varies by method:
#
# * RHF / RKS — single closed-shell density (list of cell
# blocks); the SCF driver takes ``initial_density=blocks``.
# * UHF / UKS — open-shell (alpha_blocks, beta_blocks) tuple;
# the SCF driver takes ``init_alpha=`` + ``init_beta=``.
is_closed_shell = method_upper in ("RHF", "RKS")
is_open_shell = method_upper in ("UHF", "UKS")
def _warm_kwargs(density: Any) -> dict[str, Any]:
if density is None:
return {}
if is_closed_shell:
return {"initial_density": density}
if is_open_shell:
alpha_blocks, beta_blocks = density
return {"init_alpha": alpha_blocks, "init_beta": beta_blocks}
return {}
energy, scf_result = _run_scf(
sys_current,
basis_current,
kmesh,
opts,
method_upper,
functional,
**_warm_kwargs(initial_density),
)
# Reference density: the converged density from the just-run
# reference SCF, in the shape this method requires. Used to
# warm-start the FD-displaced SCFs *and* returned to the outer
# loop's per-image cache for the next iteration.
reference_density: Optional[Any] = None
if is_closed_shell:
reference_density = [
np.asarray(b, dtype=float).copy()
for b in scf_result.density.blocks
]
elif is_open_shell:
reference_density = (
[np.asarray(b, dtype=float).copy() for b in scf_result.density_alpha.blocks],
[np.asarray(b, dtype=float).copy() for b in scf_result.density_beta.blocks],
)
n_atoms = positions.shape[0]
grad = np.zeros((n_atoms, 3), dtype=float)
for a in range(n_atoms):
for c in range(3):
disp_plus = positions.copy()
disp_plus[a, c] += fd_step_bohr
sys_p = _rebuild_with_positions(template, disp_plus)
basis_p = BasisSet(sys_p.unit_cell_molecule(), basis_name)
disp_kwargs = _warm_kwargs(reference_density)
e_p, _ = _run_scf(
sys_p,
basis_p,
kmesh,
opts,
method_upper,
functional,
**disp_kwargs,
)
disp_minus = positions.copy()
disp_minus[a, c] -= fd_step_bohr
sys_m = _rebuild_with_positions(template, disp_minus)
basis_m = BasisSet(sys_m.unit_cell_molecule(), basis_name)
e_m, _ = _run_scf(
sys_m,
basis_m,
kmesh,
opts,
method_upper,
functional,
**disp_kwargs,
)
grad[a, c] = (e_p - e_m) / (2.0 * fd_step_bohr)
return float(energy), grad, reference_density
# --- NEB force kernel ------------------------------------------------------
def _neb_forces(
positions: list[np.ndarray],
energies: list[float],
gradients: list[np.ndarray],
spring_constant: float,
frozen_mask: Optional[np.ndarray],
climbing_index: Optional[int] = None,
) -> tuple[list[np.ndarray], list[np.ndarray]]:
"""Compute NEB total forces + tangents for every intermediate image.
Returns (forces, tangents). Both are length n_intermediate
(excluding endpoints). Frozen atoms (per ``frozen_mask`` of shape
(n_atoms,) bool, True = frozen) have their force components set
to zero before return.
``climbing_index`` is the index (into ``positions``, i.e. global
image index — 0 is reactant, ``len(positions) - 1`` is product)
of the climbing image. For that image the standard
``F_spring_∥ + F_true_⊥`` is replaced by the climbing-image force
``F_climb = -∇E + 2 (∇E · τ) τ`` (Henkelman, Uberuaga, Jónsson
2000): the true force with its tangent-parallel component
inverted, no spring contribution. The climbing image then
relaxes uphill *along* τ while still relaxing perpendicular,
landing on the saddle. ``None`` ⇒ standard NEB on every
intermediate image.
"""
n_total = len(positions)
forces: list[np.ndarray] = []
tangents: list[np.ndarray] = []
for i in range(1, n_total - 1):
tau = _improved_tangent(
positions[i - 1],
positions[i],
positions[i + 1],
energies[i - 1],
energies[i],
energies[i + 1],
)
true_force = -gradients[i]
parallel = float(np.sum(true_force * tau))
if climbing_index is not None and i == climbing_index:
# Climbing image: invert the parallel component of the
# true force and drop the spring contribution.
# F_climb = -∇E + 2 (∇E · τ) τ
# = (true_force - parallel τ) + (-parallel τ)
# = F_true_⊥ − parallel·τ
f_total = true_force - 2.0 * parallel * tau
else:
# Spring force projected onto τ (signed).
d_next = float(np.linalg.norm(positions[i + 1] - positions[i]))
d_prev = float(np.linalg.norm(positions[i] - positions[i - 1]))
f_spring = spring_constant * (d_next - d_prev) * tau
# True force perpendicular to τ.
f_perp = true_force - parallel * tau
f_total = f_spring + f_perp
if frozen_mask is not None:
f_total = f_total.copy()
f_total[frozen_mask] = 0.0
forces.append(f_total)
tangents.append(tau)
return forces, tangents
# --- public entry point ----------------------------------------------------
def run_neb(
reactant: System,
product: System,
basis: str,
n_images: int = 7,
*,
method: str = "RKS",
functional: Optional[str] = "pbe",
spring_constant: float = 0.1,
interpolation: str = "idpp",
max_iter: int = 100,
conv_tol_force: float = 1e-3,
freeze_indices: Optional[Sequence[int]] = None,
dispersion_params: Any = None,
rhf_options: Any = None,
uhf_options: Any = None,
rks_options: Any = None,
uks_options: Any = None,
gradient_options: Any = None,
grid_options: Any = None,
n_jobs: int = -1,
initial_step: float = 0.05,
max_step: float = 0.2,
progress: bool = False,
climbing_image: bool = False,
climbing_image_start_fraction: float = 0.3,
kpoints: Optional[Any] = None,
fd_step_bohr: float = 1e-3,
warm_start: bool = True,
) -> NEBResult:
"""Find a minimum-energy path with improved-tangent NEB.
Parameters
----------
reactant, product
:class:`Molecule` endpoints. Same atom-number sequence and
length; the order is preserved through the path. Periodic
systems are rejected (Increment 4).
basis
Basis-set name passed to :class:`BasisSet`. Rebuilt per
image per outer iteration (the basis is bound to its nuclei).
n_images
Number of *intermediate* images. The returned path has
``n_images + 2`` images (endpoints included). Endpoints are
fixed by default.
method
``"RHF"`` / ``"UHF"`` / ``"RKS"`` / ``"UKS"`` (case-insensitive).
Same dispatch as :func:`vibeqc.run_job` and
:func:`vibeqc.molecular_optimize.optimize_molecule`.
functional
XC functional for KS-DFT (e.g. ``"pbe"``, ``"b3lyp"``).
Ignored for HF methods.
spring_constant
``k`` in Ha/bohr² for the spring term
``F_spring_i = k (|R_{i+1} - R_i| − |R_i − R_{i-1}|) τ_i``.
0.1 is the canonical default; turn down if the path is
chemically smooth, up if it kinks.
interpolation
``"idpp"`` (default, recommended for bonded paths) or
``"linear"``.
max_iter
Hard cap on outer iterations. Quick-min iterations are cheap
in this loop's terms (one parallel SCF batch each).
conv_tol_force
Convergence threshold on the max-norm of the NEB force over
all intermediate atoms (Ha/bohr). 1e-3 ≈ 0.05 eV/Å — the
Kolsbjerg 2016 recommendation.
freeze_indices
Atom indices to freeze (NEB-local; the SCF + gradient still
sees them, but their force components are zeroed before the
step). Useful for slab substrate atoms in surface NEB.
dispersion_params
Optional D3-BJ parameters; the dispersion energy and gradient
are folded into each image's energy + gradient.
n_jobs
joblib.Parallel ``n_jobs``; -1 ⇒ all cores. 1 ⇒ serial.
initial_step, max_step
Quick-min step sizes (bohr). The integrator starts at
``initial_step`` and scales up (capped at ``max_step``) when
the velocity is aligned with the force.
progress
If True, prints one line per outer iteration summarising
max-force, total path length, and TS-image energy.
climbing_image
Enable climbing-image NEB (CI-NEB, Henkelman+Uberuaga+Jónsson
2000). After a warm-up phase the highest-energy intermediate
image is promoted to "climbing": its spring contribution is
dropped and the tangent-parallel component of its true force
is inverted, so it climbs uphill along τ to the saddle while
still relaxing perpendicular. Other images keep standard NEB
dynamics. Default ``False`` (plain improved-tangent NEB).
climbing_image_start_fraction
Fraction of ``max_iter`` to spend in the plain-NEB warm-up
phase before promoting an image to climbing. Default 0.3 —
the band has typically found its rough shape by this point,
so the identity of the highest-energy image is reliable.
Ignored when ``climbing_image=False``.
kpoints
Periodic-only. Either a ``BlochKMesh`` (passed straight
through to ``run_pbc_bipole_*``) or a 3-tuple of ints (the
Monkhorst-Pack mesh sizes; converted with
:func:`vibeqc.monkhorst_pack`). Ignored when
``reactant``/``product`` are ``Molecule``. ``None`` ⇒
Γ-only mesh (``(1, 1, 1)``) for a sanity-check periodic
run; pick a real k-mesh for production.
fd_step_bohr
Periodic-only. Half-step for the central-difference
per-image BIPOLE gradient (the J^LR reciprocal contribution
is still missing from the analytic BIPOLE gradient, so the
periodic NEB driver uses the FD fallback to keep things
bit-exact). Cost per image is 6N + 1 BIPOLE SCFs;
``HANDOVER_PERIODIC_NEB.md`` tracks the work that will let
us switch to the analytic gradient. Default 1e-3 bohr —
same value the BIPOLE-FD gradient unit tests use.
warm_start
Within-image density warm-start across outer iterations
(HANDOVER_PERIODIC_NEB.md M4 — the 2-4× SCF-cost reduction
the original handover scoped). When True (default) the
converged density from outer iter N is fed in as the SCF
initial guess at iter N+1 for the same image — the geometry
change between outer iters is small relative to a SAD/Hcore
guess, so the SCF converges in fewer iterations. Currently
active for ``method="RHF"`` only; UHF / RKS / UKS and the
periodic path will land on the same hook as follow-ups (the
evaluator accepts ``initial_density=`` uniformly; only the
molecular RHF leg threads it down to the SCF entry today).
Pass ``False`` to force every SCF to cold-start (useful for
benchmarking the speedup).
Returns
-------
:class:`NEBResult`
Notes
-----
The outer loop is the quick-min (damped MD) integrator of
Henkelman+Jónsson 2000. Each intermediate image carries a
velocity; at each step the velocity is projected onto the NEB
force, zeroed if the projection is negative (preventing climbing
against the force). Step length grows by a factor of 1.1 when
aligned and is capped at ``max_step``; it is reset to
``initial_step`` on direction flips. This is the textbook
"FIRE-lite" cousin used by ASE's MDMin.
"""
is_periodic = isinstance(reactant, PeriodicSystem)
if is_periodic and not isinstance(product, PeriodicSystem):
raise ValueError(
"run_neb: reactant and product must be the same system "
"type — got PeriodicSystem and Molecule."
)
if not is_periodic and not isinstance(reactant, Molecule):
raise NotImplementedError(
"run_neb: only Molecule and PeriodicSystem endpoints are "
"supported."
)
if n_images < 1:
raise ValueError(f"n_images must be >= 1; got {n_images}")
_check_compatible(reactant, product)
method_lower = method.lower()
if method_lower not in ("rhf", "uhf", "rks", "uks"):
raise ValueError(
f"run_neb: unsupported method {method!r}. "
"Use one of 'RHF', 'UHF', 'RKS', 'UKS'."
)
# --- 1. Initial path ---------------------------------------------------
if interpolation == "idpp":
initial_systems = interpolate_idpp(reactant, product, n_images)
elif interpolation == "linear":
initial_systems = interpolate_linear(reactant, product, n_images)
else:
raise ValueError(
f"interpolation must be 'idpp' or 'linear'; got {interpolation!r}"
)
n_atoms = len(_positions_of(reactant))
frozen_mask: Optional[np.ndarray] = None
if freeze_indices is not None:
fi = {int(i) for i in freeze_indices}
bad = [i for i in fi if i < 0 or i >= n_atoms]
if bad:
raise ValueError(
f"run_neb: freeze_indices {bad} out of range "
f"[0, {n_atoms})"
)
frozen_mask = np.zeros(n_atoms, dtype=bool)
for i in fi:
frozen_mask[i] = True
positions: list[np.ndarray] = [
_positions_of(s) for s in initial_systems
]
# --- 2. Periodic-only setup: resolve the k-point mesh ------------------
kmesh = None
if is_periodic:
from ._vibeqc_core import monkhorst_pack as _mp
if kpoints is None:
kmesh = _mp(reactant, (1, 1, 1))
elif hasattr(kpoints, "kpoints"):
# Already a BlochKMesh-shaped object.
kmesh = kpoints
else:
kmesh = _mp(reactant, tuple(int(k) for k in kpoints))
# --- 3. Evaluate endpoints once + cache --------------------------------
if is_periodic:
eval_kwargs: dict[str, Any] = {
"template": reactant,
"basis_name": basis,
"method": method_lower,
"kmesh": kmesh,
"functional": functional,
"rhf_options": rhf_options,
"uhf_options": uhf_options,
"rks_options": rks_options,
"uks_options": uks_options,
"fd_step_bohr": fd_step_bohr,
}
evaluator = _evaluate_image_periodic
else:
eval_kwargs = {
"template": reactant,
"basis_name": basis,
"method": method_lower,
"functional": functional,
"rhf_options": rhf_options,
"uhf_options": uhf_options,
"rks_options": rks_options,
"uks_options": uks_options,
"gradient_options": gradient_options,
"grid_options": grid_options,
"dispersion_params": dispersion_params,
}
evaluator = _evaluate_image
e_R, g_R, _d_R = evaluator(positions[0], **eval_kwargs)
e_P, g_P, _d_P = evaluator(positions[-1], **eval_kwargs)
energies: list[float] = [e_R] + [0.0] * n_images + [e_P]
gradients: list[np.ndarray] = (
[g_R] + [np.zeros((n_atoms, 3))] * n_images + [g_P]
)
# Per-image density cache for within-image SCF warm-start
# (HANDOVER_PERIODIC_NEB.md M4). Only the RHF path on the
# molecular evaluator returns a density today; UHF / RKS / UKS
# and the periodic evaluator return None, so warm-start is a
# no-op for them. ``warm_start`` (kwarg on run_neb) defaults to
# True; pass False to force every SCF to cold-start (useful for
# benchmarking + for paranoid bisection).
image_densities: list[Optional[np.ndarray]] = [None] * (n_images + 2)
# --- 3. Outer loop (quick-min) -----------------------------------------
from joblib import Parallel, delayed
velocities: list[np.ndarray] = [
np.zeros((n_atoms, 3)) for _ in range(n_images)
]
step = initial_step
max_force = float("inf")
converged = False
n_iter = 0
last_forces: list[np.ndarray] = [
np.zeros((n_atoms, 3)) for _ in range(n_images)
]
# CI-NEB warm-up: the band needs to be roughly settled before
# we promote its highest-energy image to climbing; otherwise the
# climbing-image selection can flip between outer iterations and
# destabilise the loop.
if not 0.0 <= climbing_image_start_fraction <= 1.0:
raise ValueError(
"climbing_image_start_fraction must be in [0, 1]; got "
f"{climbing_image_start_fraction!r}"
)
climbing_warmup_iters = (
int(round(max_iter * climbing_image_start_fraction))
if climbing_image
else 0
)
climbing_index: Optional[int] = None
for outer in range(max_iter):
n_iter = outer + 1
# Per-image warm-start density: pass last iter's converged
# density when ``warm_start`` is on; otherwise pass None
# which falls back to the cold initial guess (SAD/Hcore).
results = Parallel(n_jobs=n_jobs, prefer="processes")(
delayed(evaluator)(
positions[i],
**eval_kwargs,
initial_density=image_densities[i] if warm_start else None,
)
for i in range(1, n_images + 1)
)
for k, (e_k, g_k, d_k) in enumerate(results):
energies[k + 1] = e_k
gradients[k + 1] = g_k
image_densities[k + 1] = d_k # None for non-RHF; harmless.
# Promote a climbing image once we are out of the warm-up
# phase. Henkelman+Uberuaga+Jónsson 2000: pick the highest-
# energy intermediate image at the moment of promotion;
# keep that selection fixed for the remainder of the run so
# the climber doesn't lose its momentum to a re-selection.
if (
climbing_image
and climbing_index is None
and outer >= climbing_warmup_iters
):
climbing_index = int(np.argmax(energies[1:-1])) + 1
forces, tangents = _neb_forces(
positions,
energies,
gradients,
spring_constant,
frozen_mask,
climbing_index=climbing_index,
)
max_force = max(float(np.max(np.abs(f))) for f in forces)
if progress:
ts_idx = int(np.argmax(energies[1:-1])) + 1
path_len = sum(
float(np.linalg.norm(positions[i + 1] - positions[i]))
for i in range(len(positions) - 1)
)
climb_tag = (
f" [CI={climbing_index}]" if climbing_index is not None else ""
)
print(
f"neb iter {n_iter:3d}{climb_tag}: max|F| = {max_force:.4e} "
f"Ha/bohr, E_TS = {energies[ts_idx]:.6f} Ha, "
f"path length = {path_len:.3f} bohr"
)
if max_force < conv_tol_force:
converged = True
break
# Quick-min velocity update + step. Each intermediate image:
# v ← v + Δt F
# if v·F < 0: v ← 0
# else: v ← (v·F̂) F̂ (project onto F direction)
# R ← R + Δt v
# Adaptive Δt: grow by 1.1 when aligned, reset on flips. This
# is the standard MDMin / quick-min recipe (Henkelman+Jónsson
# 2000 § III.C and ASE's MDMin optimiser).
aligned = True
for k in range(n_images):
f_k = forces[k]
f_flat = f_k.ravel()
f_norm = float(np.linalg.norm(f_flat))
if f_norm < 1e-15:
velocities[k] = np.zeros_like(f_k)
continue
v_k = velocities[k] + step * f_k
dot = float(np.sum(v_k * f_k))
if dot < 0.0:
v_k = np.zeros_like(v_k)
aligned = False
else:
# Project velocity onto force direction.
f_hat = f_k / f_norm
v_k = float(np.sum(v_k * f_hat)) * f_hat
# Cap displacement: |Δx| ≤ max_step per atom.
disp = step * v_k
disp_norm = float(np.max(np.abs(disp)))
if disp_norm > max_step:
disp = disp * (max_step / disp_norm)
v_k = np.zeros_like(v_k)
velocities[k] = v_k
positions[k + 1] = positions[k + 1] + disp
last_forces = forces
# Grow step when every image stayed aligned; reset on any
# flip. Cap at 10× initial_step to keep the integrator stable
# near convergence.
if aligned:
step = min(step * 1.1, 10.0 * initial_step)
else:
step = initial_step
# --- 4. Build NEBResult ------------------------------------------------
images: list[NEBImage] = []
final_systems = [
_rebuild_with_positions(reactant, p) for p in positions
]
final_systems[0] = reactant
final_systems[-1] = product
# Tangents for endpoints aren't defined; populate intermediates only.
forces_for_tangents, tangents = _neb_forces(
positions,
energies,
gradients,
spring_constant,
frozen_mask,
climbing_index=climbing_index,
)
for i, s in enumerate(final_systems):
img = NEBImage(system=s)
img.energy = energies[i]
img.gradient = gradients[i].copy()
if 1 <= i <= n_images:
img.tangent = tangents[i - 1]
images.append(img)
path = NEBPath(
images=images,
spring_constant=spring_constant,
climbing_image_index=climbing_index,
)
energies_arr = np.array(energies, dtype=float)
ts_index: Optional[int] = None
if n_images >= 1:
# Highest-energy intermediate image (excluding endpoints).
inner = energies_arr[1:-1]
ts_index = int(np.argmax(inner)) + 1
return NEBResult(
path=path,
energies=energies_arr,
converged=converged,
transition_state_index=ts_index,
n_iter=n_iter,
max_force=max_force,
method=method,
basis=basis,
functional=functional,
is_periodic=is_periodic,
)
__all__ = [
"NEBImage",
"NEBPath",
"NEBResult",
"interpolate_linear",
"interpolate_idpp",
"run_neb",
]