#!/usr/bin/env python3
"""Plots Nudged Elastic Band (NEB) reaction paths and landscapes.
.. versionadded:: 0.0.2
This script provides a command-line interface (CLI) to visualize data
generated from NEB calculations. It can plot:
1. **Energy/Eigenvalue Profiles:** Shows the evolution of the energy or
lowest eigenvalue along the reaction coordinate. It can overlay multiple
paths (e.g., from different optimization steps) and use a
physically-motivated Hermite spline interpolation using force data.
2. **2D Reaction Landscapes:** Plots the path on a 2D coordinate system
defined by the Root Mean Square Deviation (RMSD) from the reactant
and product structures. This requires the 'ira_mod' library.
It can also interpolate and display the 2D energy/eigenvalue surface.
The script can also render atomic structures from a .con file as insets
on the plots for key points (reactant, saddle, product).
This script follows the guidelines laid out here:
https://realpython.com/python-script-structure/
"""
# /// script
# requires-python = ">=3.11"
# dependencies = [
# "click",
# "matplotlib",
# "numpy",
# "scipy",
# "jax",
# "adjustText>=1.0",
# "cmcrameri",
# "rich",
# "ase",
# "polars",
# "h5py",
# "chemparseplot[neb,plot]>=1.9.17,<2",
# "xyzrender>=0.1.3",
# "readcon>=0.13.1",
# "rgpycrumbs>=1.10.4",
# ]
# ///
#
# Dispatch model (prefer ``rgpycrumbs eon plt-neb`` / ``python -m rgpycrumbs.cli``;
# not raw ``uv run plt_neb.py`` as the primary path):
# - ``uv run`` (default when the active env lacks the readcon/plot stack, or
# RGPYCRUMBS_FORCE_UV=1): PEP 723 deps above + optional SBOM constraints.
# - In-env / ``--dev``: ensure_import + AUTO_DEPS (CLI defaults AUTO_DEPS=1).
# - Optional CycloneDX SBOM: ``--sbom`` / RGPYCRUMBS_SBOM (PyPI pins only).
import logging
import sys
from pathlib import Path
from typing import Any
import click
import matplotlib as mpl
import matplotlib.patheffects as path_effects
import matplotlib.pyplot as plt
import numpy as np
import polars as pl
try:
from rgpycrumbs._aux import ensure_import, warn_on_direct_script_import
except ImportError: # pragma: no cover - direct script execution without package root
warn_on_direct_script_import = None
if warn_on_direct_script_import is not None:
warn_on_direct_script_import(__name__, "rgpycrumbs eon plt-neb")
# Optional heavy: resolve at use site so library import does not require
# adjustText / AUTO_DEPS until a plot actually needs label adjustment.
[docs]
def adjust_text(*args, **kwargs):
"""Lazy adjustText (ensure_import when AUTO_DEPS is on)."""
if ensure_import is not None:
from rgpycrumbs._aux import enable_library_auto_deps
enable_library_auto_deps()
return ensure_import("adjustText").adjust_text(*args, **kwargs)
from adjustText import adjust_text as _adjust_text # pragma: no cover
return _adjust_text(*args, **kwargs)
try:
from ._render_cli import add_config_option, add_render_options
from .plot_config import library_plot, run_from_click
except ImportError: # pragma: no cover - direct script execution
from rgpycrumbs.eon._render_cli import add_config_option, add_render_options
from rgpycrumbs.eon.plot_config import library_plot, run_from_click
# Lazy plot stack: AUTO_DEPS + ensure_import (same as jax / adjustText)
try:
from rgpycrumbs._aux import enable_library_auto_deps
from rgpycrumbs._aux import ensure_import as _ei
enable_library_auto_deps()
_ei("chemparseplot")
except ImportError:
pass
from chemparseplot.parse.eon.dimer_trajectory import load_dimer_trajectory
from chemparseplot.parse.eon.neb import (
aggregate_neb_landscape_data,
compute_profile_rmsd,
estimate_rbf_smoothing,
load_structures_and_calculate_additional_rmsd,
)
# --- Library Imports ---
from chemparseplot.parse.file_ import find_file_paths
from chemparseplot.parse.trajectory.hdf5 import (
history_to_landscape_df as hdf5_history_to_landscape_df,
)
from chemparseplot.parse.trajectory.hdf5 import (
history_to_profile_dats,
result_to_atoms_list,
)
from chemparseplot.parse.trajectory.hdf5 import (
result_to_profile_dat as hdf5_result_to_profile_dat,
)
from chemparseplot.parse.trajectory.neb import (
load_trajectory,
trajectory_to_landscape_df,
trajectory_to_profile_dat,
)
from matplotlib.gridspec import GridSpec
from matplotlib.patches import ArrowStyle
from rich.logging import RichHandler
try:
from chemparseplot.parse.projection import (
compute_projection_basis,
project_to_sd,
)
except ImportError:
[docs]
compute_projection_basis = None
project_to_sd = None
from chemparseplot.plot.neb import (
convert_neb_values,
default_neb_ylabel,
landscape_half_span,
landscape_projection_basis,
mark_saddle_point,
plot_energy_path,
plot_landscape_path_overlay,
plot_landscape_surface,
plot_neb_evolution,
plot_phase_points_overlay,
plot_structure_inset,
plot_structure_strip,
profile_strip_payload,
profile_structure_indices,
save_plot,
)
from chemparseplot.plot.structs import (
StructurePlacement,
convert_energy,
)
from chemparseplot.plot.theme import (
apply_axis_theme,
get_theme,
setup_global_theme,
)
# --- Logging Setup ---
logging.basicConfig(
level=logging.INFO,
format="%(levelname)s - %(message)s",
handlers=[RichHandler(rich_tracebacks=True, show_path=False, markup=True)],
)
[docs]
log = logging.getLogger("rich")
# --- Constants ---
[docs]
DEFAULT_PATH_PATTERN = "neb_path_*.con"
[docs]
IRA_KMAX_DEFAULT = 14.0
[docs]
NEB_LANDSCAPE_STRIP_ZOOM_MULT = 3.15
[docs]
NEB_PROFILE_STRIP_ZOOM_MULT = 3.4
# Landscape "all" strip: two rows of six so each molecule is larger.
[docs]
NEB_LANDSCAPE_STRIP_MAX_COLS = 6
# Profile strip needs real inches so R/SP/P stay legible under the curve.
[docs]
NEB_PROFILE_STRIP_HEIGHT_IN = 1.95
[docs]
NEB_PROFILE_MAIN_HEIGHT_IN = 3.55
# --- CLI ---
[docs]
def plot_neb_from_settings(settings: dict[str, Any]) -> Path | None:
"""Run the eOn NEB plot pipeline from a resolved settings mapping.
Prefer :func:`plot_neb` for library callers. This entry is the shared
implementation used by the Click CLI after ``resolve_from_click``.
Parameters
----------
settings:
Mapping from :func:`rgpycrumbs.eon.plot_config.merge_plot_settings`
(or ``resolve_from_click``).
Returns
-------
pathlib.Path | None
Output path when a figure is written.
.. versionadded:: 1.10.2
"""
from rgpycrumbs._aux import enable_library_auto_deps
enable_library_auto_deps()
input_dat_pattern = settings["input_dat_pattern"]
input_path_pattern = settings["input_path_pattern"]
con_file = settings.get("con_file")
additional_con = settings.get("additional_con")
source = settings["source"]
input_traj = settings.get("input_traj")
input_h5 = settings.get("input_h5")
plot_type = settings["plot_type"]
landscape_mode = settings["landscape_mode"]
landscape_path = settings["landscape_path"]
project_path = settings["project_path"]
rc_mode = settings["rc_mode"]
plot_structures = settings["plot_structures"]
rbf_smoothing = settings.get("rbf_smoothing")
show_pts = settings["show_pts"]
plot_mode = settings["plot_mode"]
surface_type = settings["surface_type"]
auto_thin = bool(settings.get("auto_thin", False))
max_surface_points = int(settings.get("max_surface_points", 64))
n_inducing = settings.get("n_inducing")
output_file = settings.get("output_file")
start = settings.get("start")
end = settings.get("end")
normalize_rc = settings["normalize_rc"]
title = settings["title"]
xlabel = settings.get("xlabel")
ylabel = settings.get("ylabel")
energy_unit = settings["energy_unit"]
highlight_last = settings["highlight_last"]
theme = settings["theme"]
cmap_profile = settings.get("cmap_profile")
cmap_landscape = settings.get("cmap_landscape")
facecolor = settings.get("facecolor")
fontsize_base = settings.get("fontsize_base")
figsize = settings["figsize"]
fig_height = settings.get("fig_height")
aspect_ratio = settings.get("aspect_ratio")
dpi = settings["dpi"]
zoom_ratio = settings["zoom_ratio"]
rotation = settings["rotation"]
perspective_tilt = settings["perspective_tilt"]
strip_renderer = settings["strip_renderer"]
xyzrender_config = settings["xyzrender_config"]
strip_spacing = settings["strip_spacing"]
strip_dividers = settings["strip_dividers"]
arrow_head_length = settings.get("arrow_head_length", 0.2)
arrow_head_width = settings.get("arrow_head_width", 0.3)
arrow_tail_width = settings.get("arrow_tail_width", 0.1)
spline_method = settings["spline_method"]
# Hermite uses -f_para as dE/d(rc). With image index that derivative is
# not dE/dindex, so Hermite overshoots and looks jerky. Use a plain cubic
# spline through the energies instead (smooth curve, no force slope).
# Explicit none|spline from the user is left alone.
if rc_mode == "index" and spline_method == "hermite":
log.info(
"rc-mode=index: switching spline-method hermite → spline "
"(force-based Hermite is meaningless on image index; cubic "
"through energies keeps a smooth curve)"
)
spline_method = "spline"
draw_reactant = settings.get("draw_reactant", (15, 60, 0.1))
draw_saddle = settings.get("draw_saddle", (15, 60, 0.1))
draw_product = settings.get("draw_product", (15, 60, 0.1))
mmf_peaks = settings.get("mmf_peaks")
peak_dir = settings.get("peak_dir")
show_evolution = settings["show_evolution"]
show_legend = settings["show_legend"]
cache_file = settings["cache_file"]
force_recompute = settings["force_recompute"]
ira_kmax = settings["ira_kmax"]
sp_file = settings["sp_file"]
augment_dat = settings.get("augment_dat")
augment_con = settings.get("augment_con")
# 1. Setup Theme
active_theme = get_theme(
theme,
cmap_profile=cmap_profile,
cmap_landscape=cmap_landscape,
font_size=fontsize_base,
facecolor=facecolor,
)
setup_global_theme(active_theme)
if fig_height and aspect_ratio:
figsize = (fig_height * aspect_ratio, fig_height)
elif fig_height or aspect_ratio:
log.error(
"Both --fig-height and --aspect-ratio must be provided together. Using default figsize."
)
fig = plt.figure(figsize=figsize, dpi=dpi)
# Layout Logic
has_strip = plot_structures in ["all", "crit_points"] and plot_type in {
"landscape",
"profile",
}
if has_strip:
# Single-row strips need little vertical gap; multi-row needs more.
n_expected = (3 if plot_structures == "crit_points" else 12) + len(
additional_con or []
)
# Prefer one strip row (matches plot_structure_strip single-row preference).
n_rows = 1 if n_expected <= 16 else 2
if plot_type == "profile":
# Temporary GridSpec; profile uses content-sized layout later so the
# strip is re-rendered at final pixel height (not crushed by a cap).
calc_hspace = 0.22 if n_rows > 1 else 0.12
height_ratios = [1, 0.62 if n_rows == 1 else 0.75]
else:
# Landscape: keep strip close; avoid a tall empty canvas under the map.
calc_hspace = 0.28 if n_rows > 1 else 0.12
height_ratios = [1, 0.32 if n_rows == 1 else 0.45]
gs = GridSpec(2, 1, height_ratios=height_ratios, hspace=calc_hspace, figure=fig)
ax = fig.add_subplot(gs[0])
ax_strip = fig.add_subplot(gs[1])
apply_axis_theme(ax_strip, active_theme)
else:
ax = fig.add_subplot(111)
ax_strip = None
apply_axis_theme(ax, active_theme)
atoms_list = None
additional_atoms_data = []
sp_data = None
# Filled so strips re-render after content-sized layout (profile + landscape).
landscape_strip_render = None
profile_strip_render = None
# Only attempt to load structures if specifically requested or needed for the plot type
if con_file:
try:
overlay_bundle = load_structures_and_calculate_additional_rmsd(
con_file, additional_con, ira_kmax, sp_file
)
atoms_list = overlay_bundle.atoms_list
additional_atoms_data = overlay_bundle.additional_structures
sp_data = overlay_bundle.saddle_point
except Exception as e:
log.error(f"Error loading structures: {e}")
# Critical failure for landscape/RMSD modes
if plot_type == "landscape" or rc_mode == "rmsd":
log.critical("Cannot proceed without structures.")
raise RuntimeError(
"Cannot proceed without structures "
f"(plot_type={plot_type!r}, rc_mode={rc_mode!r}): {e}"
) from e
# --- Trajectory source: load once if applicable ---
# In-memory ConFrame path (path_frames / plot(frames, kind="neb")): treat as traj.
frames = settings.get("frames")
traj_atoms_list = settings.get("atoms_list")
if traj_atoms_list is None and frames is not None:
from chemparseplot.parse.eon.frame_series import atoms_list_from_frames
traj_atoms_list = atoms_list_from_frames(frames)
source = "traj"
if source == "traj" and traj_atoms_list is None:
if not input_traj:
msg = "--input-traj is required when --source traj is used"
log.critical(msg)
raise RuntimeError(msg)
traj_atoms_list = load_trajectory(str(input_traj))
if plot_type == "landscape":
# --- Landscape Plot ---
z_label = default_neb_ylabel(plot_mode, energy_unit)
if source == "traj":
df = trajectory_to_landscape_df(traj_atoms_list, ira_kmax=ira_kmax)
# Use traj structures for con_file features when not provided
if atoms_list is None:
atoms_list = traj_atoms_list
elif source == "hdf5":
if not input_h5:
msg = "--input-h5 is required when --source hdf5 is used"
log.critical(msg)
raise RuntimeError(msg)
h5_str = str(input_h5)
# Prefer history file for multi-step landscape
try:
df = hdf5_history_to_landscape_df(h5_str, ira_kmax=ira_kmax)
except Exception:
log.warning("History read failed, falling back to single-step result.")
from chemparseplot.parse.trajectory.hdf5 import (
result_to_atoms_list as _r2a,
)
from chemparseplot.parse.trajectory.neb import (
trajectory_to_landscape_df as _traj_ldf,
)
hdf5_atoms = _r2a(h5_str)
df = _traj_ldf(hdf5_atoms, ira_kmax=ira_kmax)
if atoms_list is None:
atoms_list = result_to_atoms_list(h5_str)
else:
dat_paths = find_file_paths(input_dat_pattern)
con_paths = find_file_paths(str(input_path_pattern))
if not dat_paths:
msg = f"No data files found for pattern: {input_dat_pattern}"
log.critical(msg)
raise RuntimeError(msg)
# Fallback if no path files found but main file exists
if not con_paths and con_file:
con_paths = [con_file]
y_col = 2 if plot_mode == "energy" else 4
df = aggregate_neb_landscape_data(
dat_paths,
con_paths,
y_col,
None,
cache_file=cache_file,
force_recompute=force_recompute,
ira_kmax=ira_kmax,
augment_dat=augment_dat,
augment_con=augment_con,
ref_atoms=atoms_list[0] if atoms_list else None, # main reactant
prod_atoms=atoms_list[-1] if atoms_list else None, # main product
)
# Compute a SINGLE projection basis from the full dataset's endpoints.
# This basis is used consistently for: surface grid, path overlay,
# additional-con overlay, and viewport calculation.
r_full = df["r"].to_numpy()
p_full = df["p"].to_numpy()
global_basis = compute_projection_basis(r_full, p_full) if project_path else None
# Surface Generation
if landscape_mode == "surface":
if landscape_path == "last":
max_step = df["step"].max()
df_surface = df.filter(pl.col("step") == max_step)
else:
df_surface = df
# Prepare arrays
r_all = df_surface["r"].to_numpy()
p_all = df_surface["p"].to_numpy()
z_all = convert_neb_values(df_surface["z"].to_numpy(), plot_mode, energy_unit)
gr_all = convert_neb_values(
df_surface["grad_r"].to_numpy(), plot_mode, energy_unit
)
gp_all = convert_neb_values(
df_surface["grad_p"].to_numpy(), plot_mode, energy_unit
)
step_all = df_surface["step"].to_numpy()
# Heuristic for RBF smoothing if missing
if rbf_smoothing is None:
rbf_smoothing = estimate_rbf_smoothing(df)
log.info(f"Calculated heuristic RBF smoothing: {rbf_smoothing:.4f}")
extra_pts = []
if sp_data:
extra_pts.append([sp_data.r, sp_data.p])
for overlay in additional_atoms_data:
extra_pts.append([overlay.r, overlay.p])
extra_pts_arr = np.array(extra_pts) if extra_pts else None
# Pre-compute viewport from FULL data (not filtered surface data)
vp_xlim = vp_ylim = None
if project_path and global_basis is not None:
_s, _d = project_to_sd(r_full, p_full, global_basis)
# Modest s pad; d half-span matches s so Δs=Δd (true 1:1 Å).
_s_pad = max((_s.max() - _s.min()) * 0.04, 0.02)
vp_xlim = (float(_s.min() - _s_pad), float(_s.max() + _s_pad))
_s_half = 0.5 * (vp_xlim[1] - vp_xlim[0])
_half = max(
_s_half,
abs(float(_d.max())) * 1.12,
abs(float(_d.min())) * 1.12,
0.02,
)
if sp_data is not None:
_, _sd = project_to_sd(
np.array([sp_data.r]), np.array([sp_data.p]), global_basis
)
_half = max(_half, abs(float(_sd[0])) * 1.12)
for overlay in additional_atoms_data:
_, _ad = project_to_sd(
np.array([overlay.r]), np.array([overlay.p]), global_basis
)
_half = max(_half, abs(float(_ad[0])) * 1.12)
vp_ylim = (-_half, _half)
plot_landscape_surface(
ax,
r_all,
p_all,
gr_all,
gp_all,
z_all,
step_data=step_all,
method=surface_type,
rbf_smooth=rbf_smoothing,
cmap=active_theme.cmap_landscape,
show_pts=show_pts,
# so we always show 5% and 95%, this is the user defined additional one
# TODO(rg): just be a user parameter..
variance_threshold=0.5, # 50% uncertainty
project_path=project_path,
extra_points=extra_pts_arr,
n_inducing=n_inducing,
xlim=vp_xlim,
ylim=vp_ylim,
basis=global_basis,
auto_thin=auto_thin,
max_surface_points=max_surface_points,
)
# Path Overlay (Final Step)
max_step = df["step"].max()
df_final = df.filter(pl.col("step") == max_step)
final_r = df_final["r"].to_numpy()
final_p = df_final["p"].to_numpy()
final_z = convert_neb_values(df_final["z"].to_numpy(), plot_mode, energy_unit)
# Pass all-iteration data for triangulated background when no GP surface
bg_kwargs = {}
if landscape_mode != "surface":
bg_kwargs = {
"all_r": df["r"].to_numpy(),
"all_p": df["p"].to_numpy(),
"all_z": convert_neb_values(df["z"].to_numpy(), plot_mode, energy_unit),
}
plot_landscape_path_overlay(
ax,
final_r,
final_p,
final_z,
active_theme.cmap_landscape,
z_label,
project_path=project_path,
basis=global_basis,
**bg_kwargs,
)
# --- OCI-NEB/RONEB: refinement-sample overlay ---
_show_mmf = mmf_peaks
_peak_search_dir = peak_dir or Path(".")
if _show_mmf is None:
_show_mmf = any(
(_peak_search_dir / name).exists() for name in ("climb", "climb.con")
)
if _show_mmf:
from rgpycrumbs.geom.api.alignment import calculate_rmsd_from_ref
try:
from rgpycrumbs._aux import _import_from_parent_env
_ira_mod = _import_from_parent_env("ira_mod")
ira_instance = _ira_mod.IRA()
except (ImportError, AttributeError):
ira_instance = None
# Use same references as the main path
ref_r = atoms_list[0] if atoms_list else None
ref_p = atoms_list[-1] if atoms_list else None
if ref_r is not None and ref_p is not None:
try:
dimer_traj = load_dimer_trajectory(_peak_search_dir)
except (FileNotFoundError, ValueError):
dimer_traj = None
if dimer_traj is not None and dimer_traj.atoms_list:
try:
dimer_rmsd_r = calculate_rmsd_from_ref(
dimer_traj.atoms_list,
ira_instance,
ref_atom=ref_r,
ira_kmax=ira_kmax,
)
dimer_rmsd_p = calculate_rmsd_from_ref(
dimer_traj.atoms_list,
ira_instance,
ref_atom=ref_p,
ira_kmax=ira_kmax,
)
except ValueError as exc:
log.warning(
"Skipping MMF refinement overlay from %s: %s",
_peak_search_dir,
exc,
)
else:
plot_phase_points_overlay(
ax,
dimer_rmsd_r,
dimer_rmsd_p,
project_path=project_path,
path_rmsd_r=final_r,
path_rmsd_p=final_p,
)
log.info(
"Plotted %d MMF refinement frame(s)",
len(dimer_traj.atoms_list),
)
# --- OCI-NEB/RONEB: Band Evolution ---
if show_evolution:
unique_steps = sorted(df["step"].unique().to_list())
if len(unique_steps) > 1:
step_r_list = []
step_p_list = []
for step in unique_steps:
step_df = df.filter(pl.col("step") == step)
step_r_list.append(step_df["r"].to_numpy())
step_p_list.append(step_df["p"].to_numpy())
plot_neb_evolution(
ax,
step_r_list,
step_p_list,
project_path=project_path,
)
log.info("Plotted band evolution (%d steps)", len(unique_steps))
# Saddle Point Marker
if sp_data:
# Use explicit SP coordinates
sp_x_raw, sp_y_raw = sp_data.r, sp_data.p
log.info(f"Plotting explicit SP at R={sp_x_raw:.3f}, P={sp_y_raw:.3f}")
else:
# Fallback to heuristic
if plot_mode == "energy":
saddle_idx = np.argmax(final_z[1:-1]) + 1
else:
saddle_idx = np.argmin(final_z)
sp_x_raw, sp_y_raw = final_r[saddle_idx], final_p[saddle_idx]
# Apply projection to saddle point if enabled
if project_path:
_sp_basis = landscape_projection_basis(global_basis, final_r, final_p)
sp_sd = project_to_sd(np.array([sp_x_raw]), np.array([sp_y_raw]), _sp_basis)
sp_x, sp_y = float(sp_sd[0][0]), float(sp_sd[1][0])
else:
sp_x, sp_y = sp_x_raw, sp_y_raw
# Star only on the landscape — R/SP/P text already comes from the
# strip-driven main-axes labels (avoids a second floating SP box).
mark_saddle_point(
ax,
sp_x,
sp_y,
font_size=active_theme.font_size,
vline=False,
annotate=False,
)
if additional_atoms_data:
marker_cmap = mpl.colormaps.get_cmap("tab10")
for i, overlay in enumerate(additional_atoms_data):
color = marker_cmap(i % 10)
if project_path:
_add_basis = landscape_projection_basis(
global_basis, final_r, final_p
)
_s, _d = project_to_sd(
np.array([overlay.r]), np.array([overlay.p]), _add_basis
)
plot_add_r, plot_add_p = float(_s[0]), float(_d[0])
else:
plot_add_r, plot_add_p = overlay.r, overlay.p
ax.plot(
plot_add_r,
plot_add_p,
marker="*",
markersize=int(active_theme.font_size * 1.1),
color=color,
markeredgecolor="white",
markeredgewidth=1.0,
linestyle="None",
zorder=102,
label=overlay.label,
)
if has_strip and atoms_list:
strip_payload = []
# Helper to calculate projected coordinates for labels
def get_projected_coords(r_val, p_val):
if project_path:
_pc_basis = landscape_projection_basis(global_basis, final_r, final_p)
_s, _d = project_to_sd(
np.array([r_val]), np.array([p_val]), _pc_basis
)
return float(_s[0]), float(_d[0])
return r_val, p_val
# Saddle index for labels (explicit SP may still be a separate geometry).
if plot_mode == "energy":
s_idx = int(np.argmax(final_z[1:-1]) + 1)
else:
s_idx = int(np.argmin(final_z))
if plot_structures == "all":
# One entry per band image (no R/SP/P *and* numeric duplicates at
# the same x, which stacked strip captions on top of each other).
for i, atoms_i in enumerate(atoms_list):
ix, iy = get_projected_coords(final_r[i], final_p[i])
if i == 0:
label = "R"
elif i == len(atoms_list) - 1:
label = "P"
elif i == s_idx:
label = "SP"
else:
label = str(i)
strip_payload.append(
StructurePlacement(
atoms=atoms_i,
x=ix,
y=iy,
label=label,
)
)
# Explicit SP is marked on the main axes; the band already has an
# SP-labeled image so we do not double it in the strip.
else:
# Crit points only: R, SP, P.
rx, ry = get_projected_coords(final_r[0], final_p[0])
strip_payload.append(
StructurePlacement(atoms=atoms_list[0], x=rx, y=ry, label="R")
)
if sp_data:
sx, sy = get_projected_coords(sp_data.r, sp_data.p)
strip_payload.append(
StructurePlacement(
atoms=sp_data.atoms,
x=sx,
y=sy,
label=sp_data.label,
)
)
else:
sx, sy = get_projected_coords(final_r[s_idx], final_p[s_idx])
strip_payload.append(
StructurePlacement(
atoms=atoms_list[s_idx],
x=sx,
y=sy,
label="SP",
)
)
px, py = get_projected_coords(final_r[-1], final_p[-1])
strip_payload.append(
StructurePlacement(atoms=atoms_list[-1], x=px, y=py, label="P")
)
# Add additional structures
for overlay in additional_atoms_data:
ax_r, ax_p = get_projected_coords(overlay.r, overlay.p)
strip_payload.append(
StructurePlacement(
atoms=overlay.atoms,
x=ax_r,
y=ax_p,
label=overlay.label,
)
)
strip_payload.sort(key=lambda entry: entry.x)
# Defer strip drawing until after equal-aspect figure layout so
# structures are rendered at the final (large) strip pixel size.
# A first draw on the small GridSpec cell left tiny molecules.
landscape_strip_render = {
"payload": strip_payload,
"zoom": zoom_ratio * NEB_LANDSCAPE_STRIP_ZOOM_MULT,
"rotation": rotation,
"theme_color": active_theme.textcolor,
"renderer": strip_renderer,
"xyzrender_config": xyzrender_config,
"col_spacing": strip_spacing,
"show_dividers": strip_dividers,
"perspective_tilt": perspective_tilt,
"width_fill_fraction": 0.92,
# Two rows (e.g. 12 → 6+6) for larger per-cell molecules.
"max_cols": NEB_LANDSCAPE_STRIP_MAX_COLS,
"prefer_single_row": False,
}
# Annotate Main Plot -- only label R, SP, P (not additional con;
# those are identified by the legend markers instead)
main_plot_texts = []
main_labels = {"R", "SP", "P"}
for d in strip_payload:
if d.label not in main_labels:
continue
t = ax.text(
d.x,
d.y,
d.label,
fontsize=11,
fontweight="bold",
color="white",
ha="center",
va="bottom",
zorder=102,
)
t.set_path_effects(
[path_effects.withStroke(linewidth=2.5, foreground="black")]
)
main_plot_texts.append(t)
if main_plot_texts:
adjust_text(
main_plot_texts,
ax=ax,
arrowprops={"arrowstyle": "-", "color": "white", "lw": 1.0},
expand_points=(2.0, 2.0),
force_text=(1.0, 2.0),
force_points=(1.0, 2.0),
)
# Labels
if project_path:
final_xlabel = xlabel or r"Reaction progress $s$ ($\AA$)"
final_ylabel = ylabel or r"Orthogonal deviation $d$ ($\AA$)"
final_title = "Reaction Valley Projection" if title == "NEB Path" else title
else:
final_xlabel = xlabel or r"RMSD from Reactant ($\AA$)"
final_ylabel = ylabel or r"RMSD from Product ($\AA$)"
final_title = "RMSD(R,P) projection" if title == "NEB Path" else title
else:
# --- Profile Plot ---
strip_payload = []
if source == "hdf5":
if not input_h5:
msg = "--input-h5 is required when --source hdf5 is used"
log.critical(msg)
raise RuntimeError(msg)
h5_str = str(input_h5)
# Use history final step if available, else result
try:
dats = history_to_profile_dats(h5_str)
data = dats[-1]
except Exception:
data = hdf5_result_to_profile_dat(h5_str)
if atoms_list is None:
atoms_list = result_to_atoms_list(h5_str)
if rc_mode == "index":
data[1] = np.arange(data.shape[1])
elif normalize_rc:
data[1] = data[1] / data[1].max() if data[1].max() > 0 else data[1]
y_col = 2 if plot_mode == "energy" else 4
data[y_col] = convert_neb_values(data[y_col], plot_mode, energy_unit)
if plot_mode == "energy":
data[3] = convert_energy(data[3], energy_unit)
color = active_theme.highlight_color
plot_energy_path(
ax,
data[1],
data[y_col],
data[3],
color,
1.0,
20,
method=spline_method,
)
elif source == "traj":
# Trajectory source: single extxyz file -> one profile
data = trajectory_to_profile_dat(traj_atoms_list)
if atoms_list is None:
atoms_list = traj_atoms_list
if rc_mode == "index":
data[1] = np.arange(data.shape[1])
elif normalize_rc:
data[1] = data[1] / data[1].max() if data[1].max() > 0 else data[1]
y_col = 2 if plot_mode == "energy" else 4
data[y_col] = convert_neb_values(data[y_col], plot_mode, energy_unit)
if plot_mode == "energy":
data[3] = convert_energy(data[3], energy_unit)
color = active_theme.highlight_color
plot_energy_path(
ax,
data[1],
data[y_col],
data[3],
color,
1.0,
20,
method=spline_method,
)
if atoms_list and plot_structures != "none":
if has_strip:
strip_payload.extend(
profile_strip_payload(
atoms_list,
data[1],
data[y_col],
plot_structures,
plot_mode,
)
)
else:
indices = profile_structure_indices(
atoms_list, data[y_col], plot_structures, plot_mode
)
for i in indices:
if i == 0:
xybox, rad = draw_reactant[:2], draw_reactant[2]
elif i == len(atoms_list) - 1:
xybox, rad = draw_product[:2], draw_product[2]
else:
xybox, rad = draw_saddle[:2], draw_saddle[2]
if plot_structures == "all":
xybox = (15.0, 60.0 if i % 2 == 0 else -60.0)
rad = 0.1 if i % 2 == 0 else -0.1
plot_structure_inset(
ax,
atoms_list[i],
data[1][i],
data[y_col][i],
xybox,
rad,
zoom=zoom_ratio,
rotation=rotation,
renderer=strip_renderer,
xyzrender_config=xyzrender_config,
perspective_tilt=perspective_tilt,
)
else:
# eOn source: multiple .dat files
dat_paths = find_file_paths(input_dat_pattern)
file_paths_to_plot = dat_paths[start:end]
if not file_paths_to_plot:
msg = "No profile data files found in the requested start:end range"
log.error(msg)
raise RuntimeError(msg)
# Optional: Load RMSD for X-axis
rmsd_rc = None
if rc_mode == "rmsd" and atoms_list:
df_rmsd = compute_profile_rmsd(
atoms_list,
cache_file=cache_file,
force_recompute=force_recompute,
ira_kmax=ira_kmax,
)
rmsd_rc = df_rmsd["r"].to_numpy()
# Plot Loop
cm = plt.get_cmap(active_theme.cmap_profile)
color_divisor = (
len(file_paths_to_plot) - 1 if len(file_paths_to_plot) > 1 else 1.0
)
y_col = 2 if plot_mode == "energy" else 4
last_profile_rc = None
last_profile_y = None
for idx, fpath in enumerate(file_paths_to_plot):
try:
data = np.loadtxt(fpath, skiprows=1).T
except Exception as ex:
log.error(ex)
continue
# X-Axis Logic
if rc_mode == "rmsd" and rmsd_rc is not None:
if len(rmsd_rc) == data.shape[1]:
data[1] = rmsd_rc
elif rc_mode == "index":
data[1] = np.arange(data.shape[1])
elif normalize_rc:
data[1] = data[1] / data[1].max() if data[1].max() > 0 else data[1]
# Style Logic
is_last = idx == len(file_paths_to_plot) - 1
if highlight_last and is_last:
color, alpha, zorder = active_theme.highlight_color, 1.0, 20
step_label = f"Step {idx + 1} (final)"
else:
color = cm(idx / color_divisor)
alpha = 1.0 if idx == 0 else 0.5
zorder = 10 if idx == 0 else 5
step_label = f"Step {idx + 1}" if idx == 0 else None
# Plot
data[y_col] = convert_neb_values(data[y_col], plot_mode, energy_unit)
if plot_mode == "energy":
data[3] = convert_energy(data[3], energy_unit)
plot_energy_path(
ax,
data[1],
data[y_col],
data[3], # Forces
color,
alpha,
zorder,
method=spline_method,
label=step_label,
)
if is_last or last_profile_rc is None:
last_profile_rc = np.asarray(data[1], dtype=float)
last_profile_y = np.asarray(data[y_col], dtype=float)
if (
highlight_last
and is_last
and atoms_list
and plot_structures != "none"
):
if has_strip:
strip_payload.extend(
profile_strip_payload(
atoms_list,
data[1],
data[y_col],
plot_structures,
plot_mode,
)
)
else:
indices = profile_structure_indices(
atoms_list, data[y_col], plot_structures, plot_mode
)
for i in indices:
if i == 0:
xybox, rad = draw_reactant[:2], draw_reactant[2]
elif i == len(atoms_list) - 1:
xybox, rad = draw_product[:2], draw_product[2]
else:
xybox, rad = draw_saddle[:2], draw_saddle[2]
if plot_structures == "all":
xybox = (15.0, 60.0 if i % 2 == 0 else -60.0)
rad = 0.1 if i % 2 == 0 else -0.1
# Call library function
plot_structure_inset(
ax,
atoms_list[i],
data[1][i],
data[y_col][i],
xybox,
rad,
zoom=zoom_ratio,
rotation=rotation,
renderer=strip_renderer,
xyzrender_config=xyzrender_config,
)
# Saddle marker on the (final) 1D profile — gold star + vertical guide.
if (
last_profile_rc is not None
and last_profile_y is not None
and len(last_profile_y) > 2
):
if plot_mode == "energy":
s_idx = int(np.argmax(last_profile_y[1:-1]) + 1)
else:
s_idx = int(np.argmin(last_profile_y))
mark_saddle_point(
ax,
float(last_profile_rc[s_idx]),
float(last_profile_y[s_idx]),
font_size=active_theme.font_size,
vline=True,
# Strip already labels SP; skip the floating box on the curve.
annotate=not has_strip,
)
# --- Profile Additional Structures ---
if additional_atoms_data and rc_mode == "rmsd":
for i, overlay in enumerate(additional_atoms_data):
ax.axvline(
overlay.r,
color=active_theme.gridcolor,
linestyle=":",
linewidth=2,
zorder=90,
)
if has_strip:
strip_payload.append(
StructurePlacement(
atoms=overlay.atoms,
x=float(overlay.r),
label=overlay.label,
)
)
elif plot_structures != "none":
y_span = ax.get_ylim()[1] - ax.get_ylim()[0]
y_pos = ax.get_ylim()[0] + 0.9 * y_span
plot_structure_inset(
ax,
overlay.atoms,
overlay.r,
y_pos,
xybox=(
draw_saddle[0] + (i * 15),
draw_saddle[1],
), # Stagger slightly
rad=draw_saddle[2],
zoom=zoom_ratio,
rotation=rotation,
renderer=strip_renderer,
xyzrender_config=xyzrender_config,
arrow_props={
"arrowstyle": ArrowStyle.Fancy(
head_length=arrow_head_length,
head_width=arrow_head_width,
tail_width=arrow_tail_width,
),
"connectionstyle": f"arc3,rad={rad}",
"linestyle": "-",
"alpha": 0.8,
"color": "black",
"linewidth": 1.2,
},
)
if has_strip and strip_payload:
deduped_payload = []
seen = set()
for entry in sorted(strip_payload, key=lambda d: d.x):
key = (entry.label, round(entry.x, 8))
if key in seen:
continue
seen.add(key)
deduped_payload.append(entry)
# Defer strip draw until content-sized layout (same as landscape).
profile_strip_render = {
"payload": deduped_payload,
"zoom": zoom_ratio * NEB_PROFILE_STRIP_ZOOM_MULT,
"rotation": rotation,
"theme_color": active_theme.textcolor,
"renderer": strip_renderer,
"xyzrender_config": xyzrender_config,
"col_spacing": strip_spacing,
"show_dividers": strip_dividers,
"perspective_tilt": perspective_tilt,
"width_fill_fraction": 0.72,
}
# Profile Labels
if xlabel:
final_xlabel = xlabel
elif rc_mode == "rmsd":
final_xlabel = r"RMSD from reactant ($\AA$)"
elif rc_mode == "index":
final_xlabel = "Image index"
else:
final_xlabel = r"Path length ($\AA$)"
final_ylabel = ylabel or default_neb_ylabel(plot_mode, energy_unit)
if title == "NEB Path":
final_title = {
"path": "Energy against path length",
"index": "Energy against image index",
"rmsd": "Energy against reactant RMSD",
}.get(rc_mode, title)
else:
final_title = title
# Final Aesthetics
ax.set_xlabel(final_xlabel, weight="bold")
ax.set_ylabel(final_ylabel, weight="bold")
ax.set_title(final_title)
ax.minorticks_on()
projected_layout_done = False
if plot_type == "landscape" and not aspect_ratio:
if project_path:
# True 1:1 metric: Δd window = Δs window so set_aspect('equal') is a
# square panel where 1 Å of s matches 1 Å of d (same RMSD unit).
x_min, x_max = ax.get_xlim()
half_span = landscape_half_span(
(x_min, x_max),
final_r,
final_p,
additional_atoms_data,
global_basis,
equal_metric=True,
)
if sp_data is not None and global_basis is not None:
_sp_basis = landscape_projection_basis(global_basis, final_r, final_p)
_, _spd = project_to_sd(
np.array([sp_data.r]), np.array([sp_data.p]), _sp_basis
)
half_span = max(half_span, abs(float(_spd[0])) * 1.12)
ax.set_ylim(-half_span, half_span)
# Center s so the data window is exactly square (Δs == Δd == 2*half).
s_mid = 0.5 * (x_min + x_max)
ax.set_xlim(s_mid - half_span, s_mid + half_span)
x_min, x_max = ax.get_xlim()
y_min, y_max = -half_span, half_span
# Square data window (Δs == Δd == 2*half) + square axes box so
# 1 Å of s equals 1 Å of d on screen (true 1:1 RMSD metric).
# Note: d ticks show ±half (e.g. ±0.67), not 0→1.34 — the *span*
# matches s (e.g. 0→1.34); both windows are the same length in Å.
map_in = 6.2 # square: map_w_in == map_h_in
# Two-row strip needs real height so molecules stay large.
n_strip = (
len(landscape_strip_render["payload"])
if landscape_strip_render is not None
else 0
)
n_strip_rows = (
max(
1,
(n_strip + NEB_LANDSCAPE_STRIP_MAX_COLS - 1)
// NEB_LANDSCAPE_STRIP_MAX_COLS,
)
if n_strip
and landscape_strip_render is not None
and not landscape_strip_render.get("prefer_single_row", True)
else (1 if n_strip else 0)
)
strip_h_in = (3.35 if n_strip_rows >= 2 else 2.2) if has_strip else 0.0
y_label_in = 0.95
cbar_in = 1.20
top_in = 0.55
gap_map_strip_in = 0.20 if has_strip else 0.0
bottom_in = 0.18
fig_w = y_label_in + map_in + cbar_in
fig_h = top_in + map_in + gap_map_strip_in + strip_h_in + bottom_in
fig.set_size_inches(fig_w, fig_h, forward=True)
left = y_label_in / fig_w
map_w_frac = map_in / fig_w
map_h_frac = map_in / fig_h
map_bottom = (bottom_in + strip_h_in + gap_map_strip_in) / fig_h
ax.set_position([left, map_bottom, map_w_frac, map_h_frac])
ax.set_aspect("equal", adjustable="box", anchor="C")
cbar_left = left + map_w_frac + 0.012
cbar_w_frac = max(0.015, 0.20 / fig_w)
for other in list(fig.axes):
if other is ax or other is ax_strip:
continue
other.set_position([cbar_left, map_bottom, cbar_w_frac, map_h_frac])
if ax_strip is not None:
strip_bottom = bottom_in / fig_h
strip_h_frac = strip_h_in / fig_h
ax_strip.set_position([left, strip_bottom, map_w_frac, strip_h_frac])
ax_strip.set_navigate(False)
# Re-render structure strip at the final large axes size.
if landscape_strip_render is not None and ax_strip is not None:
ax_strip.clear()
apply_axis_theme(ax_strip, active_theme)
ax_strip.axis("off")
fig.canvas.draw()
plot_structure_strip(
ax_strip,
landscape_strip_render["payload"],
zoom=landscape_strip_render["zoom"],
rotation=landscape_strip_render["rotation"],
theme_color=landscape_strip_render["theme_color"],
renderer=landscape_strip_render["renderer"],
xyzrender_config=landscape_strip_render["xyzrender_config"],
col_spacing=landscape_strip_render["col_spacing"],
show_dividers=landscape_strip_render["show_dividers"],
perspective_tilt=landscape_strip_render["perspective_tilt"],
width_fill_fraction=landscape_strip_render["width_fill_fraction"],
max_cols=landscape_strip_render.get(
"max_cols", NEB_LANDSCAPE_STRIP_MAX_COLS
),
prefer_single_row=landscape_strip_render.get(
"prefer_single_row", False
),
)
ax_strip.set_position([left, strip_bottom, map_w_frac, strip_h_frac])
projected_layout_done = True
log.info(
"Set 1:1 (s,d) square panel: Δs=Δd=%.3f Å "
"(s=[%.3f, %.3f], d=[%.3f, %.3f]); strip_rows=%d; "
"figsize=(%.2f, %.2f) in",
2.0 * half_span,
x_min,
x_max,
y_min,
y_max,
n_strip_rows,
fig_w,
fig_h,
)
else:
# Raw RMSD(R) vs RMSD(P) is also Å–Å.
ax.set_aspect("equal", adjustable="box", anchor="C")
if show_legend:
ax.legend(
# Path usually sits at d ≳ 0; lower-left negative-d corner is free.
# Avoid loc='best' on wide equal-aspect maps (drops on the band).
loc="lower left",
borderaxespad=0.3,
frameon=True,
framealpha=1.0,
facecolor="white",
edgecolor="black",
fontsize=int(active_theme.font_size * 0.75),
).set_zorder(101)
if not projected_layout_done:
if plot_type == "profile" and ax_strip is not None and profile_strip_render:
# Content-sized profile: main panel + dedicated strip inches so R/SP/P
# are large enough to read. Re-render the strip at the final pixel size
# (the old min(height, 0.20) cap crushed molecules after a first draw).
y_label_in = 0.85
right_in = 0.35
top_in = 0.55
bottom_in = 0.28
# Room for the main-axes xlabel between the curve and R/SP/P strip.
gap_main_strip_in = 0.58
main_h_in = NEB_PROFILE_MAIN_HEIGHT_IN
strip_h_in = NEB_PROFILE_STRIP_HEIGHT_IN
# Prefer the requested width when provided; grow only if too narrow.
fig_w = max(float(figsize[0]), 6.4)
fig_h = top_in + main_h_in + gap_main_strip_in + strip_h_in + bottom_in
fig.set_size_inches(fig_w, fig_h, forward=True)
left = y_label_in / fig_w
main_w_frac = max(0.55, 1.0 - (y_label_in + right_in) / fig_w)
main_h_frac = main_h_in / fig_h
strip_h_frac = strip_h_in / fig_h
strip_bottom = bottom_in / fig_h
main_bottom = (bottom_in + strip_h_in + gap_main_strip_in) / fig_h
ax.set_position([left, main_bottom, main_w_frac, main_h_frac])
ax.xaxis.labelpad = 6
ax_strip.set_position([left, strip_bottom, main_w_frac, strip_h_frac])
ax_strip.set_navigate(False)
ax_strip.clear()
apply_axis_theme(ax_strip, active_theme)
ax_strip.axis("off")
fig.canvas.draw()
plot_structure_strip(
ax_strip,
profile_strip_render["payload"],
zoom=profile_strip_render["zoom"],
rotation=profile_strip_render["rotation"],
theme_color=profile_strip_render["theme_color"],
renderer=profile_strip_render["renderer"],
xyzrender_config=profile_strip_render["xyzrender_config"],
col_spacing=profile_strip_render["col_spacing"],
show_dividers=profile_strip_render["show_dividers"],
perspective_tilt=profile_strip_render["perspective_tilt"],
width_fill_fraction=profile_strip_render["width_fill_fraction"],
)
ax_strip.set_position([left, strip_bottom, main_w_frac, strip_h_frac])
projected_layout_done = True
log.info(
"Profile content layout: main=%.2f in, strip=%.2f in, figsize=(%.2f, %.2f)",
main_h_in,
strip_h_in,
fig_w,
fig_h,
)
elif not ax_strip:
plt.tight_layout(pad=0.5)
elif ax_strip:
# Fallback for non-profile strips without a deferred payload.
fig.canvas.draw()
pos_main = ax.get_position()
pos_strip = ax_strip.get_position()
gap = 0.04
strip_h = max(pos_strip.height, 0.28)
strip_y = max(0.02, pos_main.y0 - gap - strip_h)
ax_strip.set_position([pos_main.x0, strip_y, pos_main.width, strip_h])
from matplotlib.offsetbox import AnnotationBbox
for artist in ax_strip.get_children():
if isinstance(artist, AnnotationBbox):
artist.set_clip_on(True)
if output_file:
# Content-sized layouts already set figure inches; tight crop is still fine.
if projected_layout_done:
fig.canvas.draw()
fig.savefig(
output_file,
transparent=False,
dpi=dpi,
bbox_inches="tight",
pad_inches=0.12,
facecolor=fig.get_facecolor(),
)
else:
save_plot(output_file, dpi, has_strip=ax_strip is not None)
else:
plt.show()
return Path(output_file) if output_file else None
[docs]
plot_neb = library_plot("neb", plot_neb_from_settings)
@click.command()
@click.pass_context
@add_config_option
@click.option(
"--input-dat-pattern",
default=DEFAULT_INPUT_PATTERN,
help="Glob pattern for input data files.",
)
@click.option(
"--input-path-pattern",
default=DEFAULT_PATH_PATTERN,
help="Glob pattern for input path files.",
)
@click.option(
"--con-file",
type=click.Path(exists=True, dir_okay=False, path_type=Path),
default=None,
help="Path to .con trajectory file.",
)
@click.option(
"--additional-con",
type=(
click.Path(exists=True, dir_okay=False, path_type=Path),
str,
), # Takes (Path, Label)
multiple=True,
default=None,
help="Path(s) to additional .con file(s) and label.",
)
@click.option(
"--augment-dat",
type=str,
default=None,
help="Glob pattern for extra .dat files for surface fitting.",
)
@click.option(
"--augment-con",
type=str,
default=None,
help="Glob pattern for extra .con files for surface fitting.",
)
@click.option(
"--sp-file",
type=click.Path(exists=False, dir_okay=False, path_type=Path),
default=Path("sp.con"),
help="Path to explicit saddle point file (eOn sp.con).",
)
@click.option(
"--source",
type=click.Choice(["eon", "traj", "hdf5"]),
default="eon",
help="Data source: 'eon' for .dat/.con pairs, 'traj' for extxyz, 'hdf5' for ChemGP HDF5.",
)
@click.option(
"--input-h5",
type=click.Path(exists=True, dir_okay=False, path_type=Path),
default=None,
help="Path to ChemGP NEB HDF5 file (result or history).",
)
@click.option(
"--input-traj",
type=click.Path(exists=True, dir_okay=False, path_type=Path),
default=None,
help="Path to extxyz trajectory file (used with --source traj).",
)
@click.option(
"--plot-type",
type=click.Choice(["profile", "landscape"]),
default="profile",
help="Type of plot to generate.",
)
@click.option(
"--rbf-smoothing",
type=float,
default=None,
show_default=True,
help="Smoothing term for 2D RBF.",
)
@click.option(
"--landscape-mode",
type=click.Choice(["path", "surface"]),
default="surface",
help="For landscape plot: 'path' or 'surface'.",
)
@click.option(
"--landscape-path",
type=click.Choice(["last", "all"]),
default="all",
help="Last uses an interpolation only on the last path, otherwise use all points.",
)
@click.option(
"--project-path/--no-project-path",
is_flag=True,
default=True,
help="Project landscape coordinates into the reaction valley (s, d).",
)
@click.option(
"--rc-mode",
type=click.Choice(["path", "rmsd", "index"]),
default="path",
help="Reaction coordinate for profile plot.",
)
@click.option(
"--plot-structures",
type=click.Choice(["none", "all", "crit_points"]),
default="none",
help="Structures to render on the path. Requires --con-file.",
)
@click.option(
"--surface-type",
type=click.Choice(
[
"grid",
"rbf",
"grad_matern",
"grad_imq",
"grad_imq_ny",
"matern",
"imq",
"grad_rq",
"grad_se",
]
),
default="rbf",
help="Interpolation method for the 2D surface.",
)
@click.option(
"--n-inducing",
type=int,
default=None,
help="Number of inducing points for Nystrom or RFF features. Defaults to 300 (Nystrom) or 500 (RFF).",
)
@click.option(
"--show-pts/--no-show-pts",
default=True,
help="Show all paths from the optimization on the RMSD 2D plot.",
)
@click.option(
"--plot-mode",
type=click.Choice(["energy", "eigenvalue"]),
default="energy",
help="Quantity to plot.",
)
@click.option(
"-o",
"--output-file",
type=click.Path(path_type=Path),
default=None,
help="Output image filename.",
)
@click.option(
"--start", type=int, default=None, help="Start file index for profile plot."
)
@click.option("--end", type=int, default=None, help="End file index for profile plot.")
@click.option(
"--normalize-rc", is_flag=True, default=False, help="Normalize reaction coordinate."
)
@click.option("--title", default="NEB Path", help="Plot title.")
@click.option("--xlabel", default=None, help="X-axis label.")
@click.option("--ylabel", default=None, help="Y-axis label.")
@click.option(
"--energy-unit",
type=click.Choice(["eV", "kcal/mol", "kJ/mol"]),
default="eV",
show_default=True,
help="Presentation unit for energy-like axes and color scales.",
)
# --- Theme and Override Options ---
@click.option(
"--theme",
default="ruhi",
help="The plotting theme to use.",
)
@click.option("--cmap-profile", default=None, help="Colormap for profile plot.")
@click.option("--cmap-landscape", default=None, help="Colormap for landscape plot.")
@click.option("--facecolor", type=str, default=None, help="Background color.")
@click.option("--fontsize-base", type=int, default=None, help="Base font size.")
# --- Figure and Inset Options ---
@click.option(
"--figsize",
nargs=2,
type=(float, float),
default=(5.37, 5.37),
show_default=True,
help="Figure width, height in inches.",
)
@click.option(
"--fig-height",
type=float,
default=None,
help="Figure height in inches.",
)
@click.option(
"--aspect-ratio",
type=float,
default=None,
help="Figure aspect ratio.",
)
@click.option(
"--dpi",
type=int,
default=200,
show_default=True,
help="Resolution in Dots Per Inch.",
)
@click.option(
"--zoom-ratio",
type=float,
default=0.5,
show_default=True,
help="Scale the inset image.",
)
@add_render_options
@click.option(
"--arrow-head-length",
type=float,
default=0.2,
show_default=True,
help="Arrow head length.",
)
@click.option(
"--arrow-head-width",
type=float,
default=0.3,
show_default=True,
help="Arrow head width.",
)
@click.option(
"--arrow-tail-width",
type=float,
default=0.1,
show_default=True,
help="Arrow tail width.",
)
# --- Path/Spline Options ---
@click.option(
"--highlight-last/--no-highlight-last",
is_flag=True,
default=True,
help="Highlight last path.",
)
@click.option(
"--spline-method",
type=click.Choice(["hermite", "spline", "none"]),
default="hermite",
help=(
"Profile interpolant: hermite (uses -f_para as dE/drc; path RC only), "
"spline (cubic through energies), none (markers + straight segments). "
"Image-index RC defaults hermite→spline (force slopes are not dE/dindex)."
),
)
# --- Inset Position Options ---
@click.option(
"--draw-reactant",
type=(float, float, float),
nargs=3,
default=(15, 60, 0.1),
show_default=True,
help="Reactant inset pos (x, y, rad).",
)
@click.option(
"--draw-saddle",
type=(float, float, float),
nargs=3,
default=(15, 60, 0.1),
show_default=True,
help="Saddle inset pos (x, y, rad).",
)
@click.option(
"--draw-product",
type=(float, float, float),
nargs=3,
default=(15, 60, 0.1),
show_default=True,
help="Product inset pos (x, y, rad).",
)
@click.option(
"--cache-file",
type=click.Path(path_type=Path),
default=Path(".neb_landscape.parquet"),
help="Parquet cache file.",
)
@click.option(
"--force-recompute",
is_flag=True,
default=False,
help="Force re-calculation of RMSD.",
)
@click.option(
"--mmf-peaks/--no-mmf-peaks",
is_flag=True,
default=None,
help="Overlay OCI/MMF refinement samples on landscape (auto-detected from climb movies).",
)
@click.option(
"--peak-dir",
type=click.Path(exists=True, file_okay=False, path_type=Path),
default=None,
help="Directory containing OCI/MMF refinement outputs (for example climb/climb.con).",
)
@click.option(
"--show-evolution",
is_flag=True,
default=False,
help="Show band evolution across iterations (requires write_movies data).",
)
@click.option(
"--show-legend/--no-legend",
default=True,
help="Show the legends.",
)
@click.option(
[docs]
"--ira-kmax",
default=IRA_KMAX_DEFAULT,
help="kmax factor for IRA.",
)
def main(ctx, config, **params):
"""CLI entry: merge flags/config then run plot_neb_from_settings."""
return run_from_click("neb", plot_neb_from_settings, ctx, config=config, **params)
if __name__ == "__main__":
main()
