First Project

This tutorial walks you through a complete workflow: from raw NEB calculation output to a publication-quality energy landscape figure. We’ll use a surface diffusion problem as our example.

Prerequisites

Installation

# Core package with all pip-available extras
pip install "rgpycrumbs[all]"

# Or for development
git clone https://github.com/HaoZeke/rgpycrumbs
cd rgpycrumbs
uv sync --all-extras

Sample Data

For this tutorial, we assume you have:

  1. NEB calculation output files (neb_*.dat) from eOn

  2. A trajectory file (path.con) with image geometries

Download sample data:

wget https://example.com/sample-neb-data.tar.gz
tar xzf sample-neb-data.tar.gz
cd sample-neb-data

Step 1: Visualize the NEB Path

First, let’s see how the NEB optimization converged:

python -m rgpycrumbs.cli eon plt-neb --start 0 --end 50 -o neb_convergence.pdf

This produces a plot showing energy vs. iteration for each image. Look for:

  • Smooth convergence (lines should flatten)

  • Consistent spacing between images

  • No images crossing over each other

Troubleshooting: If the plot is empty, check that your neb_*.dat files exist and contain energy data in column 2.

Step 2: Extract the Final Path

Now extract the energy profile from the final iteration:

# SKIP-DOCTEST -- requires sample NEB data files
from rgpycrumbs.basetypes import nebpath
from rgpycrumbs.eon.plt_neb import extract_neb_data
import matplotlib.pyplot as plt

# Load the final iteration data
data = extract_neb_data("neb_49.dat")  # Last iteration

# Access the path
path = data.nebpath
print(f"Number of images: {len(path.energy)}")
print(f"Barrier: {max(path.energy) - min(path.energy):.4f} eV")

Step 3: Fit a Smooth Surface

For publication figures, we want a smooth curve through the data points:

# SKIP-DOCTEST -- requires sample NEB data files
from rgpycrumbs.surfaces import get_surface_model
import numpy as np

# Prepare 1D data for surface fitting
x_obs = path.arc_dist.reshape(-1, 1)  # Arc distance as 2D array
y_obs = path.energy.to("hartree").magnitude  # Convert to hartree

# Fit a gradient-enhanced Matérn surface
# (Use standard Matern if you don't have force data)
model = get_surface_model("matern")(
    x_obs, y_obs,
    ls=0.5,  # Initial length scale guess
    smoothing=1e-4  # Noise regularization
)

# Evaluate on a fine grid for smooth plotting
x_fine = np.linspace(x_obs.min(), x_obs.max(), 300).reshape(-1, 1)
y_fine = model(x_fine)

# Get uncertainty estimates
y_var = model.predict_var(x_fine)

Troubleshooting: If optimization fails, try:

  1. Different kernel: get_surface_model("tps") is more stable

  2. Better length scale: increase ls if data is smooth, decrease if rough

  3. More smoothing: increase smoothing to 1e-3

Step 4: Create the Publication Figure

Now combine everything into a publication-ready plot:

# SKIP-DOCTEST -- requires sample NEB data files
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator

fig, ax = plt.subplots(figsize=(4, 3), dpi=300)

# Plot raw data points
ax.plot(
    path.arc_dist.to("angstrom").magnitude,
    path.energy.to("eV").magnitude,
    'o',
    color='#004D40',  # Teal
    markersize=6,
    label='NEB images',
    zorder=10
)

# Plot fitted curve
ax.plot(
    x_fine.flatten(),
    y_fine * 27.2114,  # Hartree to eV
    '-',
    color='#FF655D',  # Coral
    linewidth=2,
    label='Matérn fit',
    zorder=5
)

# Add uncertainty band
ax.fill_between(
    x_fine.flatten(),
    (y_fine - 2*np.sqrt(y_var)) * 27.2114,
    (y_fine + 2*np.sqrt(y_var)) * 27.2114,
    alpha=0.2,
    color='#FF655D',
    label='95% CI'
)

# Formatting
ax.set_xlabel('Reaction Coordinate (Å)')
ax.set_ylabel('Energy (eV)')
ax.legend(frameon=False, loc='upper right')
ax.xaxis.set_major_locator(MaxNLocator(5))
ax.yaxis.set_major_locator(MaxNLocator(5))

# Tight layout for publication
fig.tight_layout()
fig.savefig('energy_profile.pdf', bbox_inches='tight')
print("Saved energy_profile.pdf")

Step 5: Analyze the Transition State

Find and characterize the saddle point:

# SKIP-DOCTEST -- requires sample NEB data files
# Find maximum (transition state)
ts_idx = np.argmax(y_fine)
ts_x = x_fine[ts_idx][0]
ts_energy = y_fine[ts_idx] * 27.2114  # eV

# Find reactant and product minima
reactant_energy = y_fine[0] * 27.2114
product_energy = y_fine[-1] * 27.2114

# Calculate barriers
forward_barrier = ts_energy - reactant_energy
reverse_barrier = ts_energy - product_energy

print(f"Transition state at: {ts_x:.2f} Å")
print(f"Forward barrier: {forward_barrier:.3f} eV")
print(f"Reverse barrier: {reverse_barrier:.3f} eV")
print(f"Reaction energy: {product_energy - reactant_energy:.3f} eV")

Next Steps

Now that you have the basics, explore:

  1. **2D Landscapes**: Use plot_landscape_surface for RMSD-based projections

  2. **Gradient-Enhanced Fitting**: Include force data for better accuracy

  3. **Batch Processing**: Process multiple NEB calculations automatically

See also:

  • tools/eon/plt_neb for NEB plotting options

  • tools/surfaces for surface fitting details

  • tools/eon/con_splitter for trajectory processing

Complete Script

Here’s the complete workflow as a single script:

# SKIP-DOCTEST -- requires sample NEB data files
#!/usr/bin/env python
"""
First Project: NEB Energy Landscape Analysis

This script processes eOn NEB output and creates a publication-quality figure.
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator

from rgpycrumbs.basetypes import nebpath
from rgpycrumbs.eon.plt_neb import extract_neb_data
from rgpycrumbs.surfaces import get_surface_model

# 1. Load data
data = extract_neb_data("neb_49.dat")
path = data.nebpath

# 2. Fit surface
x_obs = path.arc_dist.reshape(-1, 1)
y_obs = path.energy.to("hartree").magnitude

model = get_surface_model("matern")(x_obs, y_obs, ls=0.5, smoothing=1e-4)

# 3. Evaluate
x_fine = np.linspace(x_obs.min(), x_obs.max(), 300).reshape(-1, 1)
y_fine = model(x_fine)
y_var = model.predict_var(x_fine)

# 4. Plot
fig, ax = plt.subplots(figsize=(4, 3), dpi=300)

ax.plot(
    path.arc_dist.to("angstrom").magnitude,
    path.energy.to("eV").magnitude,
    'o', color='#004D40', markersize=6, label='NEB images', zorder=10
)

ax.plot(
    x_fine.flatten(), y_fine * 27.2114,
    '-', color='#FF655D', linewidth=2, label='Matérn fit', zorder=5
)

ax.fill_between(
    x_fine.flatten(),
    (y_fine - 2*np.sqrt(y_var)) * 27.2114,
    (y_fine + 2*np.sqrt(y_var)) * 27.2114,
    alpha=0.2, color='#FF655D', label='95% CI'
)

ax.set_xlabel('Reaction Coordinate (Å)')
ax.set_ylabel('Energy (eV)')
ax.legend(frameon=False, loc='upper right')
fig.tight_layout()
fig.savefig('energy_profile.pdf', bbox_inches='tight')

# 5. Analyze
ts_idx = np.argmax(y_fine)
print(f"Barrier: {(y_fine[ts_idx] - y_fine[0]) * 27.2114:.3f} eV")

Run with:

python first_project.py