How-to: Create Publication NEB Figures

Problem

You have parsed NEB data (eOn or ORCA) and want publication-quality figures for energy profiles, 2D landscapes, and structure galleries.

1D Energy Profiles

Hermite Spline (Default)

Uses cubic Hermite splines with Savitzky-Golay smoothed force derivatives:

import matplotlib.pyplot as plt
import numpy as np
from chemparseplot.plot.neb import plot_energy_path, SmoothingParams

fig, ax = plt.subplots(figsize=(5.37, 5.37), dpi=200)

plot_energy_path(
    ax,
    rc=rmsd_r,             # reaction coordinate (RMSD from reactant)
    energy=energies,       # energy array
    f_para=f_parallel,     # parallel force component
    color="#1f77b4",
    alpha=1.0,
    zorder=10,
    method="hermite",
    smoothing=SmoothingParams(window_length=5, polyorder=2),
)

ax.set_xlabel(r"RMSD from Reactant ($\AA$)")
ax.set_ylabel("Energy (eV)")
fig.savefig("profile.pdf", dpi=300, bbox_inches="tight")

Standard Cubic Spline

Use method”spline”= for B-spline interpolation without force information:

plot_energy_path(
    ax, rc, energy, f_para,
    color="#1f77b4", alpha=1.0, zorder=10,
    method="spline",
)

Eigenvalue Profile

Plot the lowest eigenvalue along the path (useful for saddle point characterization):

from chemparseplot.plot.neb import plot_eigenvalue_path

plot_eigenvalue_path(
    ax,
    rc=rmsd_r,
    eigenvalue=eigenvalues,
    color="#d62728",
    alpha=0.8,
    zorder=10,
    grid_color="white",    # horizontal zero-line color
)

2D Landscape Surfaces

        flowchart LR
  RMSD[RMSD-A / RMSD-B] --> PROJ{project_path?}
  PROJ -->|true| SD["(s, d) valley"]
  PROJ -->|false| RAW["raw RMSD axes"]
  SD --> FIT[Surface fit]
  RAW --> FIT
  CFG[SurfaceFitConfig] --> FIT
  FIT --> CONT[contourf + variance]
  CONT --> STRIP[optional structure strip]
    

Tip

Prefer SurfaceFitConfig / TOML keys over ad-hoc kwargs when wiring suite plots. auto_thin stays off unless dense movies force it.

Gradient-Enhanced Matern (Default)

from chemparseplot.plot.neb import plot_landscape_surface, plot_landscape_path_overlay

fig, ax = plt.subplots(figsize=(5.37, 5.37), dpi=200)

plot_landscape_surface(
    ax,
    rmsd_r=df["r"].to_numpy(),
    rmsd_p=df["p"].to_numpy(),
    grad_r=df["grad_r"].to_numpy(),
    grad_p=df["grad_p"].to_numpy(),
    z_data=df["z"].to_numpy(),
    step_data=df["step"].to_numpy(),
    method="grad_matern",
    cmap="viridis",
    show_pts=True,
    variance_threshold=0.05,
    project_path=True,       # reaction valley projection
)

# Overlay the colored NEB path
plot_landscape_path_overlay(
    ax,
    r=latest_r, p=latest_p, z=latest_z,
    cmap="viridis",
    z_label="Energy (eV)",
    project_path=True,
)

Gradient-Enhanced IMQ

For very large NEB histories, grad_imq can switch to Nystrom approximation above the suite threshold. Dense single-ended force-eval movies are a different problem: prefer opt-in thinning (below) rather than always thinning.

plot_landscape_surface(
    ax, rmsd_r, rmsd_p, grad_r, grad_p, z_data,
    method="grad_imq",
    n_inducing=200,          # number of inducing points for Nystrom
)

Dense clouds: SurfaceFitConfig / auto_thin (default off)

eOn write_movies can dump every force evaluation. On 100+ nearly collinear points, GradientIMQ may raise ValueError: Surface prediction produced no finite values for contourf. Thinning is opt-in so historical calls stay unchanged.

from chemparseplot.plot.neb import SurfaceFitConfig, plot_landscape_surface

# TOML-friendly mapping (same keys as rgpycrumbs plot.toml)
cfg = SurfaceFitConfig.from_mapping({
    "auto_thin": True,
    "max_surface_points": 64,
})

plot_landscape_surface(
    ax, rmsd_r, rmsd_p, grad_r, grad_p, z_data,
    method="grad_imq",
    surface_fit=cfg,
)
# Equivalent kwargs: auto_thin=True, max_surface_points=64

When enabled, only the GP training set is subsampled (first/last + evenly spaced intermediates). Path scatter and viewport still use the full cloud.

From the CLI suite, set the same keys in plot TOML (not CLI flags):

[shared]
auto_thin = true
max_surface_points = 64

Standard RMSD Axes (No Projection)

Set project_path=False to plot in raw (RMSD-R, RMSD-P) coordinates:

plot_landscape_surface(
    ax, rmsd_r, rmsd_p, grad_r, grad_p, z_data,
    project_path=False,
)
ax.set_xlabel(r"RMSD from Reactant ($\AA$)")
ax.set_ylabel(r"RMSD from Product ($\AA$)")

Extra Points on the Surface

Plot additional structures (saddle points, local minima) on the landscape:

import numpy as np

extra = np.array([
    [sp_rmsd_r, sp_rmsd_p],
    [min_rmsd_r, min_rmsd_p],
])

plot_landscape_surface(
    ax, rmsd_r, rmsd_p, grad_r, grad_p, z_data,
    extra_points=extra,
)

Structure Rendering

Structure Strip (Gallery)

Renders a horizontal gallery of atomic structures below the energy profile:

from chemparseplot.plot.neb import plot_structure_strip

fig, (ax_energy, ax_struct) = plt.subplots(
    2, 1, figsize=(7, 6), height_ratios=[3, 1]
)

plot_structure_strip(
    ax_struct,
    atoms_list=[atoms_list[0], atoms_list[3], atoms_list[-1]],
    labels=["Reactant", "TS", "Product"],
    zoom=0.3,
    rotation="0x,90y,0z",
    theme_color="black",
    max_cols=6,
    renderer="ase",          # or "xyzrender" for ray-traced images
)

Inset Structure Annotation

Place a structure image as an annotation on an energy profile:

from chemparseplot.plot.neb import plot_structure_inset

plot_structure_inset(
    ax,
    atoms=atoms_list[3],
    x=rmsd_r[3],             # data x coordinate
    y=energies[3],           # data y coordinate
    xybox=(40, 40),          # offset in points
    rad=0.3,                 # arrow curvature
    zoom=0.4,
    rotation="0x,90y,0z",
    renderer="ase",          # or "xyzrender"
)

xyzrender Backend

For ray-traced structure images, use the xyzrender renderer:

# xyzrender stages via rgpycrumbs ensure_import when AUTO_DEPS is on
# (or: pip install 'xyzrender>=0.1.3'). Uses the Python API, not a PATH binary.
plot_structure_strip(
    ax, atoms_list, labels,
    renderer="xyzrender",
)

Theme System

Built-in Themes

from chemparseplot.plot.theme import get_theme, setup_global_theme

# Available themes: "ruhi", "cmc.batlow"
theme = get_theme("ruhi")
setup_global_theme(theme)

Custom Theme Overrides

theme = get_theme("ruhi", font_size=14, facecolor="#f5f5f5")
setup_global_theme(theme)

Publication Preset

Adds line width, spine removal, and DPI defaults on top of the base theme:

from chemparseplot.plot.theme import setup_publication_theme

theme = get_theme("ruhi")
setup_publication_theme(theme)

Per-Axis Theming

from chemparseplot.plot.theme import apply_axis_theme

apply_axis_theme(ax, theme)

Journal Export

DPI and Size for Common Journals

# ACS journals: single column = 3.25 in, double = 7.0 in
fig, ax = plt.subplots(figsize=(3.25, 3.25), dpi=300)

# Nature: single column = 89 mm = 3.5 in
fig, ax = plt.subplots(figsize=(3.5, 3.5), dpi=300)

Save as Vector (PDF/EPS)

fig.savefig("figure.pdf", dpi=300, bbox_inches="tight", pad_inches=0.05)
fig.savefig("figure.eps", dpi=300, bbox_inches="tight")

Save as Raster (PNG/TIFF)

fig.savefig("figure.png", dpi=600, bbox_inches="tight")
fig.savefig("figure.tiff", dpi=600, bbox_inches="tight")

Single-Ended Method Visualization

The same 2D reaction valley projection works for dimer/saddle search and minimization trajectories. The interpretation changes: s becomes optimization progress and d becomes lateral deviation.

Dimer/Saddle Search Landscape

Requires eOn output with write_movies=true and a climb movie file. Metadata-rich full CON movies are treated as the canonical structured source; climb.dat remains available as a fallback/export table when the embedded metadata is incomplete.

from chemparseplot.parse.eon.dimer_trajectory import load_dimer_trajectory
from chemparseplot.parse.neb_utils import calculate_landscape_coords, compute_synthetic_gradients
from chemparseplot.plot.optimization import plot_optimization_landscape
import numpy as np
import matplotlib.pyplot as plt

traj = load_dimer_trajectory(Path("saddle_job/"))

# Use IRA for RMSD (if available)
try:
    from rgpycrumbs._aux import _import_from_parent_env
    ira_mod = _import_from_parent_env("ira_mod")
    ira_instance = ira_mod.IRA()
except ImportError:
    ira_instance = None

rmsd_a, rmsd_b = calculate_landscape_coords(
    traj.atoms_list, ira_instance, ira_kmax=1.8,
    ref_a=traj.initial_atoms,     # initial structure
    ref_b=traj.saddle_atoms,      # saddle point
)

energies = traj.dat_df["delta_e"].to_numpy()
f_para = -np.gradient(energies)
grad_a, grad_b = compute_synthetic_gradients(rmsd_a, rmsd_b, f_para)

fig, ax = plt.subplots(figsize=(5.37, 5.37), dpi=200)
plot_optimization_landscape(
    ax, rmsd_a, rmsd_b, grad_a, grad_b, energies,
    project_path=True,
    label_mode="optimization",
)
fig.savefig("saddle_landscape.pdf", dpi=300, bbox_inches="tight")

Minimization Landscape

from chemparseplot.parse.eon.min_trajectory import load_min_trajectory

traj = load_min_trajectory(Path("min_job/"))
rmsd_a, rmsd_b = calculate_landscape_coords(
    traj.atoms_list, ira_instance, ira_kmax=1.8,
    ref_a=traj.initial_atoms,
    ref_b=traj.final_atoms,
)
# ... same pattern as above with plot_optimization_landscape

Energy/Convergence Profiles

from chemparseplot.plot.optimization import plot_optimization_profile, plot_convergence_panel

fig, (ax_e, ax_ev) = plt.subplots(1, 2, figsize=(10, 4))
dat = traj.dat_df
plot_optimization_profile(
    ax_e,
    dat["iteration"].to_numpy(),
    dat["delta_e"].to_numpy(),
    eigenvalues=dat["eigenvalue"].to_numpy(),
    ax_eigen=ax_ev,
)

OCI-NEB/RONEB: MMF Peak Overlay

When OCI-NEB produces peak{NN}_pos.con files, overlay them on the landscape. If a metadata-rich climb / climb.con refinement movie is also present in the same directory, higher-level tools can distinguish the dimer/MMF samples from the base NEB band on the same 2D projection:

from chemparseplot.plot.neb import plot_mmf_peaks_overlay

plot_mmf_peaks_overlay(
    ax, peak_rmsd_r, peak_rmsd_p, peak_energies,
    project_path=True,
)

CLI Shortcuts

# Saddle search landscape
rgpycrumbs eon plt-saddle --job-dir saddle_job/ --plot-type landscape

# Minimization profile
rgpycrumbs eon plt-min --job-dir min_job/ --plot-type profile

# NEB with MMF peaks
rgpycrumbs eon plt-neb --mmf-peaks --peak-dir . --plot-type landscape

See Also

• How-to: Parse eOn NEB Calculations

• How-to: Parse ORCA NEB Calculations

• Plotting API Reference

• NEB Parsing Design