Note
Go to the end to download the full example code.
Finding Reaction Paths with eOn and a Metatomic Potential¶
- Authors:
Rohit Goswami @HaoZeke, Hanna Tuerk @HannaTuerk, Arslan Mazitov @abmazitov, Michele Ceriotti @ceriottim
This example describes how to find the reaction pathway for oxadiazole formation from N₂O and ethylene. We will use the PET-MAD metatomic model to calculate the potential energy and forces.
The primary goal is to contrast a standard Nudged Elastic Band (NEB) calculation using the atomic simulation environment (ASE) with more sophisticated methods available in the eOn package. For even a relatively simple reaction like this, a basic NEB implementation can struggle to converge or may time out. We will show how eOn’s advanced features, such as energy-weighted springs and mixing in single-ended dimer search steps, can efficiently locate and refine the transition state along the path.
Our approach will be:
Set up the PET-MAD metatomic calculator.
Use ASE to generate an initial IDPP reaction path.
Illustrate the limitations of a standard NEB calculation in ASE.
Refine the path and locate the transition state saddle point using eOn’s optimizers, including energy-weighted springs and the dimer method.
Visualize the final converged pathway.
Demonstrate endpoint relaxation with eOn
Importing Required Packages¶
First, we import all the necessary python packages for this task.
By convention, all import
statements are at the top of the file.
import os
import subprocess
import sys
from contextlib import chdir
from pathlib import Path
import ase.io as aseio
import ira_mod
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
import numpy as np
from ase.mep import NEB
from ase.optimize import LBFGS
from ase.visualize import view
from ase.visualize.plot import plot_atoms
from metatomic.torch.ase_calculator import MetatomicCalculator
from rgpycrumbs.eon.helpers import write_eon_config
from rgpycrumbs.run.jupyter import run_command_or_exit
# sphinx_gallery_thumbnail_number = 4
Obtaining the Foundation Model - PET-MAD¶
PET-MAD is an instance of a point edge transformer model trained on
the diverse MAD dataset
which learns equivariance through data driven measures
instead of having equivariance baked in [1]. In turn, this enables
the PET model to have greater design space to learn over. Integration in
Python and the C++ eOn client occurs through the metatomic software [2],
which in turn relies on the atomistic machine learning toolkit build
over metatensor. Essentially using any of the metatomic models involves
grabbing weights off of HuggingFace and loading them with
metatomic before running the
engine of choice.
repo_id = "lab-cosmo/upet"
tag = "v1.5.0"
url_path = f"models/pet-mad-xs-{tag}.ckpt"
fname = Path(f"models/pet-mad-xs-{tag}.pt")
url = f"https://huggingface.co/{repo_id}/resolve/main/{url_path}"
fname.parent.mkdir(parents=True, exist_ok=True)
subprocess.run(
[
"mtt",
"export",
url,
"-o",
str(fname), # Convert Path to string
],
check=True,
)
print(f"Successfully exported {fname}.")
Successfully exported models/pet-mad-xs-v1.5.0.pt.
Nudged Elastic Band (NEB)¶
Given two known configurations on a potential energy surface (PES), often one wishes to determine the path of highest probability between the two. Under the harmonic approximation to transition state theory, connecting the configurations (each point representing a full molecular structure) by a discrete set of images allows one to evolve the path under an optimization algorithm, and allows approximating the reaction to three states: the reactant, product, and transition state.
The location of this transition state (≈ the point with the highest energy along this path) determines the barrier height of the reaction. This saddle point can be found by transforming the second derivatives (Hessian) to step along the softest mode. However, an approximation which is free from explicitly finding this mode involves moving the highest image of a NEB path: the “climbing” image.
Mathematically, the saddle point has zero first derivatives and a single negative eigenvalue. The climbing image technique moves the highest energy image along the reversed NEB tangent force, avoiding the cost of full Hessian diagonalization used in single-ended methods [3].
Concretely, in this example, we will consider a reactant and product state, for oxadiazole formation, namely N₂O and ethylene.
reactant = aseio.read("data/min_reactant.con")
product = aseio.read("data/min_product.con")
We can visualize these structures using ASE.
fig, (ax1, ax2) = plt.subplots(1, 2)
plot_atoms(reactant, ax1, rotation=("-90x,0y,0z"))
plot_atoms(product, ax2, rotation=("-90x,0y,0z"))
ax1.text(0.3, -1, "reactant")
ax2.text(0.3, -1, "product")
ax1.set_axis_off()
ax2.set_axis_off()
Endpoint minimization¶
For finding reaction pathways, the endpoints should be minimized. We provide initial configurations which are already minimized, but in order to see how to relax endpoints with eOn, please have a look at the end of this tutorial.
Initial path generation¶
To begin an NEB method, an initial path is required, the optimal construction of which still forms an active area of research. The simplest approximation to an initial path for NEB methods linearly interpolate between the two known configurations building on intuition developed from “drag coordinate” methods. This may break bonds or otherwise also unphysically pass atoms through each other, similar to the effect of incorrect permutations. To ameliorate this effect, the NEB algorithm is often started from the linear interpolation but then the path is optimized on a surrogate potential energy surface, commonly something cheap and analytic, like the IDPP (Image dependent pair potential, [5]) which provides a surface based on bond distances, and thus preventing atom-in-atom collisions.
Here we use the IDPP from ASE to setup the initial path. You can find more information about this method in the corresponding ASE tutorial or in the original IDPP publication [5]. A brief pedagogical discussion of the transition state methods may be found on the Rowan blog, though the software is proprietary there.
N_INTERMEDIATE_IMGS = 10
# total includes the endpoints
TOTAL_IMGS = N_INTERMEDIATE_IMGS + 2
images = [reactant]
images += [reactant.copy() for _ in range(N_INTERMEDIATE_IMGS)]
images += [product]
neb = NEB(images)
neb.interpolate("idpp")
We don’t cover subtleties in setting the number of images, typically too many intermediate images may cause kinks but too few will be unable to resolve the tangent to any reasonable quality.
For eOn, we write the initial path to a file called idppPath.dat.
output_dir = Path("path")
output_dir.mkdir(exist_ok=True)
output_files = [output_dir / f"{num:02d}.con" for num in range(TOTAL_IMGS)]
for outfile, img in zip(output_files, images):
aseio.write(outfile, img)
print(f"Wrote {len(output_files)} IDPP images to '{output_dir}/'.")
summary_file = Path("idppPath.dat")
summary_file.write_text("\n".join(str(f.resolve()) for f in output_files) + "\n")
print(f"Wrote absolute paths to '{summary_file}'.")
Wrote 12 IDPP images to 'path/'.
Wrote absolute paths to 'idppPath.dat'.
Running NEBs¶
We will now consider actually running the Nudged Elastic Band with different codes.
ASE and Metatomic¶
We first consider using metatomic with the ASE calculator.
# define the calculator
def mk_mta_calc():
return MetatomicCalculator(
fname,
device="cpu",
non_conservative=False,
uncertainty_threshold=0.001,
)
# set calculators for images
ipath = [reactant] + [reactant.copy() for _ in range(10)] + [product]
for img in ipath:
img.calc = mk_mta_calc()
print(img.calc._model.capabilities().outputs)
neb = NEB(ipath, climb=True, k=5, method="improvedtangent")
neb.interpolate("idpp")
# store initial path guess for plotting
initial_energies = np.array([img.get_potential_energy() for img in ipath])
# setup the NEB calculation
optimizer = LBFGS(neb, trajectory="A2B.traj", logfile="opt.log")
conv = optimizer.run(fmax=0.01, steps=100)
print("Check if calculation has converged:", conv)
if conv:
print(neb)
final_energies = np.array([img.get_potential_energy() for img in ipath])
# Plot initial and final path
plt.figure(figsize=(8, 6))
# Initial Path (Blue)
plt.plot(
initial_energies - initial_energies[0],
"o-",
label="Initial Path (IDPP)",
color="xkcd:blue",
)
# Final Path (Orange)
plt.plot(
final_energies - initial_energies[0],
"o-",
label="Optimized Path (LBFGS)",
color="xkcd:orange",
)
# Metadata
plt.xlabel("Image number")
plt.ylabel("Potential Energy (eV)")
plt.legend()
plt.grid(True, alpha=0.3)
plt.title("NEB Path Evolution")
plt.show()
{'energy': <torch.ScriptObject object at 0x560d0b0ef860>, 'energy_ensemble': <torch.ScriptObject object at 0x560d0b72c2c0>, 'energy_uncertainty': <torch.ScriptObject object at 0x560d0b85d0e0>, 'features': <torch.ScriptObject object at 0x560d0b72c600>, 'mtt::aux::energy_last_layer_features': <torch.ScriptObject object at 0x560d0b85d2a0>, 'mtt::aux::non_conservative_forces_last_layer_features': <torch.ScriptObject object at 0x560d069eace0>, 'mtt::aux::non_conservative_forces_uncertainty': <torch.ScriptObject object at 0x560d09ac92c0>, 'mtt::aux::non_conservative_stress_last_layer_features': <torch.ScriptObject object at 0x560d09a7c0d0>, 'mtt::aux::non_conservative_stress_uncertainty': <torch.ScriptObject object at 0x560d0b85b150>, 'non_conservative_forces': <torch.ScriptObject object at 0x560d0aff5c50>, 'non_conservative_stress': <torch.ScriptObject object at 0x560d0b85c6f0>}
/home/runner/work/atomistic-cookbook/atomistic-cookbook/.nox/eon-pet-neb/lib/python3.13/site-packages/ase/calculators/calculator.py:519: UserWarning: Some of the atomic energy uncertainties are larger than the threshold of 0.001 eV. The prediction is above the threshold for atoms [0 1 2 3 4 5 6 7 8].
self.calculate(atoms, [name], system_changes)
Check if calculation has converged: False
In the 100 NEB steps we took, the structure did unfortunately not converge. The metatomic calculator for PET-MAD v1.5.0 provides LLPR based energy uncertainties. As we obtain a warning that the uncertainty of the path structure sampled is very high, we stop after 100 steps. The ASE algorithm with LBFGS optimizer does not find good intermediate structures and does not converge at all. Our test showed that the FIRE optimizer works better in this context, but still takes over 500 steps to converge, and since second order methods are faster, we consider the LBFGS routine throughout this notebook.
We thus want to look at a different code, which manages to compute a NEB for this simple system more efficiently.
eOn and Metatomic¶
eOn has two improvements to accurately locate the saddle point.
Energy weighting for improving tangent resolution near the climbing image
The Off-path climbing image (6) which involves iteratively switching to the dimer method for faster convergence by the climbing image.
To use eOn, we setup a function that writes the desired eOn input for us and
runs the eonclient binary. Since we are in a notebook environment, we will
use several abstractions over raw subprocess calls. In practice, writing
and using eOn involves a configuration file, which we define as a dictionary
to be used with a helper to generate the final output.
# Define configuration as a dictionary for clarity
neb_settings = {
"Main": {"job": "nudged_elastic_band", "random_seed": 706253457},
"Potential": {"potential": "Metatomic"},
"Metatomic": {"model_path": fname.absolute()},
"Nudged Elastic Band": {
"images": N_INTERMEDIATE_IMGS,
# initialization
"initializer": "file",
"initial_path_in": "idppPath.dat",
"minimize_endpoints": "false",
# CI-NEB settings
"climbing_image_method": "true",
"climbing_image_converged_only": "true",
"ci_after": 0.5,
"ci_after_rel": 0.8,
# energy weighing
"energy_weighted": "true",
"ew_ksp_min": 0.972,
"ew_ksp_max": 9.72,
# OCI-NEB settings
"ci_mmf": "true",
"ci_mmf_after": 0.1,
"ci_mmf_after_rel": 0.5,
"ci_mmf_penalty_strength": 1.5,
"ci_mmf_penalty_base": 0.4,
"ci_mmf_angle": 0.9,
"ci_mmf_nsteps": 1000,
},
"Optimizer": {
"max_iterations": 1000,
"opt_method": "lbfgs",
"max_move": 0.1,
"converged_force": 0.01,
},
"Debug": {"write_movies": "true"},
}
Which now let’s us write out the final triplet of reactant, product, and configuration of the eOn-NEB.
write_eon_config(Path("."), neb_settings)
aseio.write("reactant.con", reactant)
aseio.write("product.con", product)
Wrote eOn config to 'config.ini'
Run the main C++ client¶
This runs ‘eonclient’ and streams output live. If it fails, the notebook execution stops here.
run_command_or_exit(["eonclient"], capture=True, timeout=300)
eOn Client
2.12.0 (d586950)
compiled Sun Mar 08 01:19:26 PM GMT 2026
OS: linux
Arch: x86_64
Hostname: runnervm46oaq
PID: 2739
DIR: /home/runner/work/atomistic-cookbook/atomistic-cookbook/examples/eon-pet-neb
Loading parameter file config.ini
* [Main] job: nudged_elastic_band
* [Main] random_seed: 706253457
* [Potential] potential: Metatomic
* [Debug] write_movies: true
* [Optimizer] opt_method: lbfgs
* [Optimizer] converged_force: 0.01
* [Optimizer] max_iterations: 1000
* [Optimizer] max_move: 0.1
* [Metatomic] model_path: /home/runner/work/atomistic-cookbook/atomistic-cookbook/examples/eon-pet-neb/models/pet-mad-xs-v1.5.0.pt
* [Nudged Elastic Band] images: 10
* [Nudged Elastic Band] energy_weighted: true
* [Nudged Elastic Band] ew_ksp_min: 0.972
* [Nudged Elastic Band] ew_ksp_max: 9.72
* [Nudged Elastic Band] climbing_image_method: true
* [Nudged Elastic Band] climbing_image_converged_only: true
* [Nudged Elastic Band] ci_after: 0.5
* [Nudged Elastic Band] ci_after_rel: 0.8
* [Nudged Elastic Band] ci_mmf: true
* [Nudged Elastic Band] ci_mmf_after: 0.1
* [Nudged Elastic Band] ci_mmf_after_rel: 0.5
* [Nudged Elastic Band] ci_mmf_nsteps: 1000
* [Nudged Elastic Band] ci_mmf_angle: 0.9
* [Nudged Elastic Band] ci_mmf_penalty_strength: 1.5
* [Nudged Elastic Band] ci_mmf_penalty_base: 0.4
* [Nudged Elastic Band] initializer: file
* [Nudged Elastic Band] initial_path_in: idppPath.dat
* [Nudged Elastic Band] minimize_endpoints: false
[MetatomicPotential] Initializing...
[MetatomicPotential] Loading model from '/home/runner/work/atomistic-cookbook/atomistic-cookbook/examples/eon-pet-neb/models/pet-mad-xs-v1.5.0.pt'
[MetatomicPotential] Using device: cpu
[MetatomicPotential] Using dtype: float32
[MetatomicPotential] Initialization complete.
minimize_endpoints == false: not minimizing endpoints.
[MetatomicPotential] Atomic numbers changed, re-creating types tensor.
Nudged elastic band calculation started.
===============================================================
NEB Optimization Configuration
===============================================================
Baseline Force : 4.8417
Climbing Image (CI) : ENABLED
- Relative Trigger : 3.8733 (Factor: 0.80)
- Absolute Trigger : 0.5000
- Converged Only : true
Hybrid MMF (RONEB) : ENABLED
- Initial Threshold : 2.4208 (Factor: 0.50)
- Absolute Floor : 0.1000
- Penalty Scheme : 0.40 (Base: 0.40, Str: 1.50)
- Angle Tolerance : 0.9000
---------------------------------------------------------------
iteration step size ||Force|| max image max energy
---------------------------------------------------------------
1 0.0000e+00 4.8417e+00 7 1.725
2 5.0376e-02 4.5821e+00 7 1.516
3 3.7332e-02 4.2743e+00 7 1.36
4 3.2983e-02 3.8123e+00 8 1.313
5 2.5497e-02 2.5528e+00 8 1.274
6 1.7107e-02 1.9778e+00 8 1.261
7 2.3434e-02 1.9507e+00 8 1.239
8 9.2884e-02 3.6056e+00 8 1.103
9 6.7319e-02 2.1734e+00 8 0.9666
10 4.8547e-02 1.8217e+00 8 1.012
11 7.8434e-03 1.6218e+00 8 1.009
Triggering MMF. Force: 1.6218, Threshold: 2.4208 (0.50x baseline)
Saddle point search started from reactant with energy -55.538291931152344 eV.
[Dimer] Step Step Size Delta E ||Force|| Curvature Torque Angle Rots Align
[IDimerRot] ----- --------- ---------- ------------------ -13.0514 5.092 3.255 1 1.000
[Dimer] 1 0.0782749 0.0288 3.01974e+00 -13.0514 5.092 3.255 1
[IDimerRot] ----- --------- ---------- ------------------ -14.6870 3.436 3.452 1 0.998
[Dimer] 2 0.0611317 -0.0292 1.59326e+00 -14.6870 3.436 3.452 1
[IDimerRot] ----- --------- ---------- ------------------ -14.8460 8.122 4.636 1 0.995
[Dimer] 3 0.0261520 -0.0538 1.09470e+00 -14.8460 8.122 4.636 1
[IDimerRot] ----- --------- ---------- ------------------ -15.5198 6.618 4.383 1 0.995
[Dimer] 4 0.0153358 -0.0627 1.02465e+00 -15.5198 6.618 4.383 1
[IDimerRot] ----- --------- ---------- ------------------ -15.4939 2.925 1.425 1 0.990
[Dimer] 5 0.1146375 -0.0911 1.21043e+00 -15.4939 2.925 1.425 1
[IDimerRot] ----- --------- ---------- ------------------ -14.9935 4.169 2.243 1 0.990
[Dimer] 6 0.0443296 -0.0994 4.89101e-01 -14.9935 4.169 2.243 1
[IDimerRot] ----- --------- ---------- ------------------ -14.8470 3.323 2.498 1 0.989
[Dimer] 7 0.0124607 -0.1016 2.70053e-01 -14.8470 3.323 2.498 1
[IDimerRot] ----- --------- ---------- ------------------ -14.5060 5.041 ------ ---- 0.988
[Dimer] 8 0.0071106 -0.1027 3.16591e-01 -14.5060 2.520 4.828 0
[IDimerRot] ----- --------- ---------- ------------------ -14.5580 2.667 1.749 1 0.988
[Dimer] 9 0.0251421 -0.1057 3.57330e-01 -14.5580 2.667 1.749 1
[IDimerRot] ----- --------- ---------- ------------------ -14.5460 4.484 ------ ---- 0.988
[Dimer] 10 0.0486354 -0.0987 1.32748e+00 -14.5460 2.242 3.734 0
[IDimerRot] ----- --------- ---------- ------------------ -14.6727 4.946 ------ ---- 0.988
[Dimer] 11 0.0200490 -0.1075 1.91103e-01 -14.6727 2.473 4.783 0
[IDimerRot] ----- --------- ---------- ------------------ -14.7318 4.515 ------ ---- 0.988
[Dimer] 12 0.0052334 -0.1077 9.15985e-02 -14.7318 2.258 4.356 0
[IDimerRot] ----- --------- ---------- ------------------ -14.8018 4.564 ------ ---- 0.988
[Dimer] 13 0.0062779 -0.1078 4.24296e-02 -14.8018 2.282 4.382 0
[IDimerRot] ----- --------- ---------- ------------------ -14.8764 4.637 ------ ---- 0.988
[Dimer] 14 0.0019485 -0.1079 3.78475e-02 -14.8764 2.318 4.429 0
[IDimerRot] ----- --------- ---------- ------------------ -14.8935 4.661 ------ ---- 0.988
[Dimer] 15 0.0042606 -0.1079 7.10907e-02 -14.8935 2.330 4.446 0
[IDimerRot] ----- --------- ---------- ------------------ -14.8688 4.770 ------ ---- 0.988
[Dimer] 16 0.0016818 -0.1079 2.81388e-02 -14.8688 2.385 4.556 0
[IDimerRot] ----- --------- ---------- ------------------ -14.8596 4.732 ------ ---- 0.988
[Dimer] 17 0.0007758 -0.1079 1.24376e-02 -14.8596 2.366 4.524 0
[IDimerRot] ----- --------- ---------- ------------------ -14.8529 4.766 ------ ---- 0.988
[Dimer] 18 0.0004586 -0.1079 9.68664e-03 -14.8529 2.383 4.558 0
NEB converged after MMF. Force: 0.0097
0 0.000000 0.000000 0.207697
1 0.282181 -0.006584 -0.164238
2 0.535350 0.034321 -0.111273
3 0.793048 0.054756 -0.054319
4 1.066606 0.061256 0.014155
5 1.368841 0.075085 -0.079275
6 1.689444 0.201614 -1.015680
7 1.953833 0.636887 -1.119796
8 2.308917 0.901016 -0.000616
9 2.658222 0.836441 0.421163
10 2.878520 0.091713 3.579631
11 3.229941 -0.731091 -0.335196
Found 5 extrema
Energy reference: -56.54723358154297
extrema #1 at image position 0.5647820345116134 with energy -0.016631910798466265 and curvature 0.10530467667611913
extrema #2 at image position 3.8372236726847824 with energy 0.061567213020090605 and curvature -0.02280495003613408
extrema #3 at image position 4.0853545814596135 with energy 0.06107631389293999 and curvature 0.048078473574813936
extrema #4 at image position 8.002283294316609 with energy 0.901016481020001 and curvature -0.09432545742591755
extrema #5 at image position 10.958360874662544 with energy -0.7335247394473186 and curvature 2.7658094341382315
Final state:
Nudged elastic band, successful.
Generated MMF peak 00 at position 3.837 (Energy: 0.062 eV)
Generated MMF peak 01 at position 8.002 (Energy: 0.901 eV)
Timing Information:
Real time: 3.952 seconds
User time: 6.818 seconds
System time: 0.428 seconds
CompletedProcess(args=['eonclient'], returncode=0, stdout="eOn Client\n2.12.0 (d586950)\ncompiled Sun Mar 08 01:19:26 PM GMT 2026\nOS: linux\nArch: x86_64\nHostname: runnervm46oaq\nPID: 2739\nDIR: /home/runner/work/atomistic-cookbook/atomistic-cookbook/examples/eon-pet-neb\nLoading parameter file config.ini\n* [Main] job: nudged_elastic_band\n* [Main] random_seed: 706253457\n* [Potential] potential: Metatomic\n* [Debug] write_movies: true\n* [Optimizer] opt_method: lbfgs\n* [Optimizer] converged_force: 0.01\n* [Optimizer] max_iterations: 1000\n* [Optimizer] max_move: 0.1\n* [Metatomic] model_path: /home/runner/work/atomistic-cookbook/atomistic-cookbook/examples/eon-pet-neb/models/pet-mad-xs-v1.5.0.pt\n* [Nudged Elastic Band] images: 10\n* [Nudged Elastic Band] energy_weighted: true\n* [Nudged Elastic Band] ew_ksp_min: 0.972\n* [Nudged Elastic Band] ew_ksp_max: 9.72\n* [Nudged Elastic Band] climbing_image_method: true\n* [Nudged Elastic Band] climbing_image_converged_only: true\n* [Nudged Elastic Band] ci_after: 0.5\n* [Nudged Elastic Band] ci_after_rel: 0.8\n* [Nudged Elastic Band] ci_mmf: true\n* [Nudged Elastic Band] ci_mmf_after: 0.1\n* [Nudged Elastic Band] ci_mmf_after_rel: 0.5\n* [Nudged Elastic Band] ci_mmf_nsteps: 1000\n* [Nudged Elastic Band] ci_mmf_angle: 0.9\n* [Nudged Elastic Band] ci_mmf_penalty_strength: 1.5\n* [Nudged Elastic Band] ci_mmf_penalty_base: 0.4\n* [Nudged Elastic Band] initializer: file\n* [Nudged Elastic Band] initial_path_in: idppPath.dat\n* [Nudged Elastic Band] minimize_endpoints: false\n[MetatomicPotential] Initializing...\n[MetatomicPotential] Loading model from '/home/runner/work/atomistic-cookbook/atomistic-cookbook/examples/eon-pet-neb/models/pet-mad-xs-v1.5.0.pt'\n[MetatomicPotential] Using device: cpu\n[MetatomicPotential] Using dtype: float32\n[MetatomicPotential] Initialization complete.\nminimize_endpoints == false: not minimizing endpoints.\n[MetatomicPotential] Atomic numbers changed, re-creating types tensor.\nNudged elastic band calculation started.\n===============================================================\n NEB Optimization Configuration\n===============================================================\n Baseline Force : 4.8417\n Climbing Image (CI) : ENABLED\n - Relative Trigger : 3.8733 (Factor: 0.80)\n - Absolute Trigger : 0.5000\n - Converged Only : true\n Hybrid MMF (RONEB) : ENABLED\n - Initial Threshold : 2.4208 (Factor: 0.50)\n - Absolute Floor : 0.1000\n - Penalty Scheme : 0.40 (Base: 0.40, Str: 1.50)\n - Angle Tolerance : 0.9000\n---------------------------------------------------------------\n iteration step size ||Force|| max image max energy\n---------------------------------------------------------------\n 1 0.0000e+00 4.8417e+00 7 1.725\n 2 5.0376e-02 4.5821e+00 7 1.516\n 3 3.7332e-02 4.2743e+00 7 1.36\n 4 3.2983e-02 3.8123e+00 8 1.313\n 5 2.5497e-02 2.5528e+00 8 1.274\n 6 1.7107e-02 1.9778e+00 8 1.261\n 7 2.3434e-02 1.9507e+00 8 1.239\n 8 9.2884e-02 3.6056e+00 8 1.103\n 9 6.7319e-02 2.1734e+00 8 0.9666\n 10 4.8547e-02 1.8217e+00 8 1.012\n 11 7.8434e-03 1.6218e+00 8 1.009\nTriggering MMF. Force: 1.6218, Threshold: 2.4208 (0.50x baseline)\nSaddle point search started from reactant with energy -55.538291931152344 eV.\n[Dimer] Step Step Size Delta E ||Force|| Curvature Torque Angle Rots Align\n[IDimerRot] ----- --------- ---------- ------------------ -13.0514 5.092 3.255 1 1.000\n[Dimer] 1 0.0782749 0.0288 3.01974e+00 -13.0514 5.092 3.255 1\n[IDimerRot] ----- --------- ---------- ------------------ -14.6870 3.436 3.452 1 0.998\n[Dimer] 2 0.0611317 -0.0292 1.59326e+00 -14.6870 3.436 3.452 1\n[IDimerRot] ----- --------- ---------- ------------------ -14.8460 8.122 4.636 1 0.995\n[Dimer] 3 0.0261520 -0.0538 1.09470e+00 -14.8460 8.122 4.636 1\n[IDimerRot] ----- --------- ---------- ------------------ -15.5198 6.618 4.383 1 0.995\n[Dimer] 4 0.0153358 -0.0627 1.02465e+00 -15.5198 6.618 4.383 1\n[IDimerRot] ----- --------- ---------- ------------------ -15.4939 2.925 1.425 1 0.990\n[Dimer] 5 0.1146375 -0.0911 1.21043e+00 -15.4939 2.925 1.425 1\n[IDimerRot] ----- --------- ---------- ------------------ -14.9935 4.169 2.243 1 0.990\n[Dimer] 6 0.0443296 -0.0994 4.89101e-01 -14.9935 4.169 2.243 1\n[IDimerRot] ----- --------- ---------- ------------------ -14.8470 3.323 2.498 1 0.989\n[Dimer] 7 0.0124607 -0.1016 2.70053e-01 -14.8470 3.323 2.498 1\n[IDimerRot] ----- --------- ---------- ------------------ -14.5060 5.041 ------ ---- 0.988\n[Dimer] 8 0.0071106 -0.1027 3.16591e-01 -14.5060 2.520 4.828 0\n[IDimerRot] ----- --------- ---------- ------------------ -14.5580 2.667 1.749 1 0.988\n[Dimer] 9 0.0251421 -0.1057 3.57330e-01 -14.5580 2.667 1.749 1\n[IDimerRot] ----- --------- ---------- ------------------ -14.5460 4.484 ------ ---- 0.988\n[Dimer] 10 0.0486354 -0.0987 1.32748e+00 -14.5460 2.242 3.734 0\n[IDimerRot] ----- --------- ---------- ------------------ -14.6727 4.946 ------ ---- 0.988\n[Dimer] 11 0.0200490 -0.1075 1.91103e-01 -14.6727 2.473 4.783 0\n[IDimerRot] ----- --------- ---------- ------------------ -14.7318 4.515 ------ ---- 0.988\n[Dimer] 12 0.0052334 -0.1077 9.15985e-02 -14.7318 2.258 4.356 0\n[IDimerRot] ----- --------- ---------- ------------------ -14.8018 4.564 ------ ---- 0.988\n[Dimer] 13 0.0062779 -0.1078 4.24296e-02 -14.8018 2.282 4.382 0\n[IDimerRot] ----- --------- ---------- ------------------ -14.8764 4.637 ------ ---- 0.988\n[Dimer] 14 0.0019485 -0.1079 3.78475e-02 -14.8764 2.318 4.429 0\n[IDimerRot] ----- --------- ---------- ------------------ -14.8935 4.661 ------ ---- 0.988\n[Dimer] 15 0.0042606 -0.1079 7.10907e-02 -14.8935 2.330 4.446 0\n[IDimerRot] ----- --------- ---------- ------------------ -14.8688 4.770 ------ ---- 0.988\n[Dimer] 16 0.0016818 -0.1079 2.81388e-02 -14.8688 2.385 4.556 0\n[IDimerRot] ----- --------- ---------- ------------------ -14.8596 4.732 ------ ---- 0.988\n[Dimer] 17 0.0007758 -0.1079 1.24376e-02 -14.8596 2.366 4.524 0\n[IDimerRot] ----- --------- ---------- ------------------ -14.8529 4.766 ------ ---- 0.988\n[Dimer] 18 0.0004586 -0.1079 9.68664e-03 -14.8529 2.383 4.558 0\nNEB converged after MMF. Force: 0.0097\n 0 0.000000 0.000000 0.207697\n 1 0.282181 -0.006584 -0.164238\n 2 0.535350 0.034321 -0.111273\n 3 0.793048 0.054756 -0.054319\n 4 1.066606 0.061256 0.014155\n 5 1.368841 0.075085 -0.079275\n 6 1.689444 0.201614 -1.015680\n 7 1.953833 0.636887 -1.119796\n 8 2.308917 0.901016 -0.000616\n 9 2.658222 0.836441 0.421163\n 10 2.878520 0.091713 3.579631\n 11 3.229941 -0.731091 -0.335196\nFound 5 extrema\nEnergy reference: -56.54723358154297\nextrema #1 at image position 0.5647820345116134 with energy -0.016631910798466265 and curvature 0.10530467667611913\nextrema #2 at image position 3.8372236726847824 with energy 0.061567213020090605 and curvature -0.02280495003613408\nextrema #3 at image position 4.0853545814596135 with energy 0.06107631389293999 and curvature 0.048078473574813936\nextrema #4 at image position 8.002283294316609 with energy 0.901016481020001 and curvature -0.09432545742591755\nextrema #5 at image position 10.958360874662544 with energy -0.7335247394473186 and curvature 2.7658094341382315\nFinal state: \nNudged elastic band, successful.\nGenerated MMF peak 00 at position 3.837 (Energy: 0.062 eV)\nGenerated MMF peak 01 at position 8.002 (Energy: 0.901 eV)\nTiming Information:\n Real time: 3.952 seconds\n User time: 6.818 seconds\n System time: 0.428 seconds\n")
Visual interpretation¶
rgpycrumbs is a visualization toolkit designed to bridge the gap between raw computational output and physical intuition, mapping high-dimensional NEB trajectories onto interpretable 1D energy profiles and 2D RMSD landscapes. As it is normally used from the command-line, here we define a helper.
def run_neb_plot(
mode: str,
con_file: str = "neb.con",
output_file: str = "plot.png",
title: str = "",
rotation: str = "90x,0y,0z",
) -> None:
"""
Constructs the CLI command for rgpycrumbs plotting to avoid clutter in notebooks.
mode: 'profile' (1D) or 'landscape' (2D)
"""
base_cmd = [
sys.executable,
"-m",
"rgpycrumbs.cli",
"eon",
"plt-neb",
"--con-file",
con_file,
"--output-file",
output_file,
"--ase-rotation",
rotation,
"--facecolor",
"white",
"--figsize",
"5.37",
"5.37",
"--fontsize-base",
"16",
"--show-pts",
"--zoom-ratio",
"0.25",
"--plot-structures",
"crit_points",
]
if title:
base_cmd.extend(["--title", title])
if mode == "profile":
base_cmd.extend(
[
"--plot-type",
"profile",
]
)
elif mode == "landscape":
base_cmd.extend(
[
"--plot-type",
"landscape",
"--rc-mode",
"path",
"--landscape-mode",
"surface",
"--landscape-path",
"all",
"--surface-type",
"grad_imq",
]
)
else:
raise ValueError(f"Unknown plot mode: {mode}")
# Run the generated command
run_command_or_exit(base_cmd, capture=False, timeout=60)
We check both the standard 1D profile against the path reaction coordinate, or the distance between intermediate images:
# Clean env to prevent backend conflicts in notebooks
os.environ.pop("MPLBACKEND", None)
# Run the 1D plotting command using the helper
run_neb_plot("profile", title="NEB Path Optimization", output_file="1D_oxad.png")
# Display the result
img = mpimg.imread("1D_oxad.png")
plt.figure(figsize=(5.37, 5.37))
plt.imshow(img)
plt.axis("off")
plt.show()
--> Dispatching to: uv run /home/runner/work/atomistic-cookbook/atomistic-cookbook/.nox/eon-pet-neb/lib/python3.13/site-packages/rgpycrumbs/eon/plt_neb.py --con-file neb.con --output-file 1D_oxad.png --ase-rotation 90x,0y,0z --facecolor white --figsize 5.37 5.37 --fontsize-base 16 --show-pts --zoom-ratio 0.25 --plot-structures crit_points --title NEB Path Optimization --plot-type profile
Downloading matplotlib (8.3MiB)
Downloading polars-runtime-32 (44.7MiB)
Downloading ase (2.8MiB)
Downloading kiwisolver (1.4MiB)
Downloading numpy (15.8MiB)
Downloading pillow (6.7MiB)
Downloading pygments (1.2MiB)
Downloading ml-dtypes (4.8MiB)
Downloading scipy (33.6MiB)
Downloading fonttools (4.7MiB)
Downloading jaxlib (78.7MiB)
Downloading h5py (5.2MiB)
Downloading jax (2.9MiB)
Downloaded kiwisolver
Downloaded pygments
Downloaded ml-dtypes
Downloaded h5py
Downloaded fonttools
Downloaded pillow
Downloaded jax
Downloaded matplotlib
Downloaded polars-runtime-32
Downloaded ase
Downloaded numpy
Downloaded jaxlib
Downloaded scipy
Installed 34 packages in 74ms
[03/16/26 14:55:56] INFO INFO - Setting global rcParams for ruhi theme
WARNING WARNING - Font 'Atkinson Hyperlegible' not found.
Falling back to 'sans-serif'.
INFO INFO - Reading structures from neb.con
[03/16/26 14:55:57] INFO INFO - Loaded 12 structures.
INFO INFO - Loading explicit saddle point from sp.con
INFO INFO - Searching for files with pattern:
'neb_*.dat'
INFO INFO - Found 12 file(s).
The 2D PES landscape is projected onto reaction-valley coordinates [3, 7]: progress along the path and orthogonal deviation, computed from permutation-corrected RMSD distances to the reactant and product. The energy surface is interpolated using a gradient-enhanced inverse multiquadric (IMQ) Gaussian process that incorporates both energies and projected tangential forces from the NEB optimization history.
# Run the 2D plotting command using the helper
run_neb_plot("landscape", title="NEB-RMSD Surface", output_file="2D_oxad.png")
# Display the result
img = mpimg.imread("2D_oxad.png")
plt.figure(figsize=(5.37, 5.37))
plt.imshow(img)
plt.axis("off")
plt.show()
--> Dispatching to: uv run /home/runner/work/atomistic-cookbook/atomistic-cookbook/.nox/eon-pet-neb/lib/python3.13/site-packages/rgpycrumbs/eon/plt_neb.py --con-file neb.con --output-file 2D_oxad.png --ase-rotation 90x,0y,0z --facecolor white --figsize 5.37 5.37 --fontsize-base 16 --show-pts --zoom-ratio 0.25 --plot-structures crit_points --title NEB-RMSD Surface --plot-type landscape --rc-mode path --landscape-mode surface --landscape-path all --surface-type grad_imq
[03/16/26 14:55:59] INFO INFO - Setting global rcParams for ruhi theme
WARNING WARNING - Font 'Atkinson Hyperlegible' not found.
Falling back to 'sans-serif'.
INFO INFO - Reading structures from neb.con
INFO INFO - Loaded 12 structures.
INFO INFO - Loading explicit saddle point from sp.con
INFO INFO - Searching for files with pattern:
'neb_*.dat'
INFO INFO - Found 12 file(s).
INFO INFO - Searching for files with pattern:
'neb_path_*.con'
INFO INFO - Found 12 file(s).
INFO INFO - Computing Landscape data...
INFO INFO - Saving Landscape cache to
.neb_landscape.parquet...
INFO INFO - Calculated heuristic RBF smoothing: 0.0973
INFO INFO - Generating 2D surface using grad_imq
(Projected: True)...
INFO:2026-03-16 14:55:59,627:jax._src.xla_bridge:822: Unable to initialize backend 'tpu': INTERNAL: Failed to open libtpu.so: libtpu.so: cannot open shared object file: No such file or directory
INFO INFO - Unable to initialize backend 'tpu':
INTERNAL: Failed to open libtpu.so: libtpu.so:
cannot open shared object file: No such file or
directory
[03/16/26 14:56:06] INFO INFO - Plotting explicit SP at R=0.723, P=0.599
WARNING WARNING - Looks like you are using a tranform that
doesn't support FancyArrowPatch, using ax.annotate
instead. The arrows might strike through texts.
Increasing shrinkA in arrowprops might help.
INFO INFO - Set symmetric Y-axis limits: [-0.74, 0.74]
/home/runner/work/atomistic-cookbook/atomistic-cookbook/.nox/eon-pet-neb/lib/python3.13/site-packages/rgpycrumbs/eon/plt_neb.py:1132: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
plt.tight_layout(pad=0.5)
Each black dot is a configuration evaluated during NEB optimization [7]. The horizontal axis measures progress along the converged path; the vertical axis measures perpendicular displacement. Both coordinates derive from permutation-corrected RMSD (via IRA [4]) to the reactant and product. The energy surface is interpolated by a gradient-enhanced inverse multiquadric GP that uses both the energy and the tangential NEB force at each evaluated configuration. See [3, Chapter 4] for details.
Relaxing the endpoints with eOn¶
In this final part we come back to an essential point of performing NEB calculations, and that is the relaxation of the initial states. In the tutorials above we used directly relaxed structures, and here we are demonstrating how these can be relaxed. We first load structures which are not relaxed.
reactant = aseio.read("data/reactant.con")
product = aseio.read("data/product.con")
# For compatibility with eOn, we also need to provide
# a unit cell
def center_cell(atoms):
atoms.set_cell([20, 20, 20])
atoms.pbc = True
atoms.center()
return atoms
reactant = center_cell(reactant)
product = center_cell(product)
The resulting reactant has a larger box:
fig, (ax1, ax2) = plt.subplots(1, 2)
plot_atoms(reactant, ax1, rotation=("-90x,0y,0z"))
plot_atoms(product, ax2, rotation=("-90x,0y,0z"))
ax1.text(0.3, -1, "reactant")
ax2.text(0.3, -1, "product")
ax1.set_axis_off()
ax2.set_axis_off()
# Reactant setup
dir_reactant = Path("min_reactant")
dir_reactant.mkdir(exist_ok=True)
aseio.write(dir_reactant / "pos.con", reactant)
# Product setup
dir_product = Path("min_product")
dir_product.mkdir(exist_ok=True)
aseio.write(dir_product / "pos.con", product)
# Shared minimization settings
min_settings = {
"Main": {"job": "minimization", "random_seed": 706253457},
"Potential": {"potential": "Metatomic"},
"Metatomic": {"model_path": fname.absolute()},
"Optimizer": {
"max_iterations": 2000,
"opt_method": "lbfgs",
"max_move": 0.1,
"converged_force": 0.01,
},
}
write_eon_config(dir_reactant, min_settings)
write_eon_config(dir_product, min_settings)
Wrote eOn config to 'min_reactant/config.ini'
Wrote eOn config to 'min_product/config.ini'
Run the minimization¶
The ‘eonclient’ will use the correct configuration within the folder.
with chdir(dir_reactant):
run_command_or_exit(["eonclient"], capture=True, timeout=300)
with chdir(dir_product):
run_command_or_exit(["eonclient"], capture=True, timeout=300)
eOn Client
2.12.0 (d586950)
compiled Sun Mar 08 01:19:26 PM GMT 2026
OS: linux
Arch: x86_64
Hostname: runnervm46oaq
PID: 2968
DIR: /home/runner/work/atomistic-cookbook/atomistic-cookbook/examples/eon-pet-neb/min_reactant
Loading parameter file config.ini
* [Main] job: minimization
* [Main] random_seed: 706253457
* [Potential] potential: Metatomic
* [Optimizer] opt_method: lbfgs
* [Optimizer] converged_force: 0.01
* [Optimizer] max_iterations: 2000
* [Optimizer] max_move: 0.1
* [Metatomic] model_path: /home/runner/work/atomistic-cookbook/atomistic-cookbook/examples/eon-pet-neb/models/pet-mad-xs-v1.5.0.pt
[MetatomicPotential] Initializing...
[MetatomicPotential] Loading model from '/home/runner/work/atomistic-cookbook/atomistic-cookbook/examples/eon-pet-neb/models/pet-mad-xs-v1.5.0.pt'
[MetatomicPotential] Using device: cpu
[MetatomicPotential] Using dtype: float32
[MetatomicPotential] Initialization complete.
Beginning minimization of pos.con
[Matter] Iter Step size ||Force|| Energy
[MetatomicPotential] Atomic numbers changed, re-creating types tensor.
[Matter] 0 0.00000e+00 5.54443e+00 -56.42982
[Matter] 1 2.97071e-02 3.12868e+00 -56.52656
[Matter] 2 7.57456e-03 1.55711e+00 -56.55402
[Matter] 3 1.35664e-02 9.31052e-01 -56.57292
[Matter] 4 5.29642e-03 6.00603e-01 -56.57933
[Matter] 5 8.08060e-03 3.57719e-01 -56.58324
[Matter] 6 3.40677e-03 4.18714e-01 -56.58333
[Matter] 7 1.46115e-03 1.26751e-01 -56.58393
[Matter] 8 7.16582e-04 9.85748e-02 -56.58410
[Matter] 9 3.95169e-03 9.95970e-02 -56.58454
[Matter] 10 2.01624e-03 5.61905e-02 -56.58468
[Matter] 11 1.95003e-03 3.97002e-02 -56.58477
[Matter] 12 2.30683e-03 6.36212e-02 -56.58485
[Matter] 13 4.66550e-03 9.43427e-02 -56.58500
[Matter] 14 1.22800e-02 1.83059e-01 -56.58521
[Matter] 15 9.89259e-03 1.20386e-01 -56.58556
[Matter] 16 1.61361e-02 8.66203e-02 -56.58602
[Matter] 17 2.11452e-02 1.09236e-01 -56.58640
[Matter] 18 4.16577e-02 1.74626e-01 -56.58646
[Matter] 19 1.91988e-02 7.57295e-02 -56.58683
[Matter] 20 2.29560e-03 5.05628e-02 -56.58693
[Matter] 21 7.75299e-03 7.10412e-02 -56.58712
[Matter] 22 2.21505e-03 6.77133e-02 -56.58721
[Matter] 23 6.80430e-03 1.49062e-01 -56.58727
[Matter] 24 4.00760e-03 6.41103e-02 -56.58733
[Matter] 25 1.17224e-03 4.57010e-02 -56.58735
[Matter] 26 7.80825e-04 2.77179e-02 -56.58738
[Matter] 27 1.87048e-03 3.39392e-02 -56.58740
[Matter] 28 1.53786e-03 3.90554e-02 -56.58744
[Matter] 29 1.23242e-02 1.60178e-01 -56.58748
[Matter] 30 3.79748e-03 6.96497e-02 -56.58764
[Matter] 31 5.26884e-03 5.55254e-02 -56.58779
[Matter] 32 1.25438e-02 1.00909e-01 -56.58803
[Matter] 33 1.44681e-02 9.65765e-02 -56.58823
[Matter] 34 2.02486e-02 2.82730e-01 -56.58789
[Matter] 35 6.10932e-03 4.63355e-02 -56.58842
[Matter] 36 4.10213e-03 3.67786e-02 -56.58846
[Matter] 37 1.34446e-02 5.52373e-02 -56.58856
[Matter] 38 1.23548e-02 8.28853e-02 -56.58861
[Matter] 39 1.75677e-03 3.46862e-02 -56.58866
[Matter] 40 2.11208e-03 2.11309e-02 -56.58868
[Matter] 41 1.32127e-03 1.93667e-02 -56.58870
[Matter] 42 2.12929e-03 2.21846e-02 -56.58872
[Matter] 43 5.18183e-03 6.48984e-02 -56.58876
[Matter] 44 4.81726e-03 5.28584e-02 -56.58882
[Matter] 45 2.32792e-02 1.04619e-01 -56.58904
[Matter] 46 4.82330e-02 2.07603e-01 -56.58913
[Matter] 47 1.77333e-02 4.73333e-01 -56.58892
[Matter] 48 1.71859e-02 9.60735e-02 -56.58949
[Matter] 49 1.94490e-02 8.20720e-02 -56.58958
[Matter] 50 7.67642e-03 4.66023e-02 -56.58968
[Matter] 51 1.79312e-02 3.59808e-02 -56.58978
[Matter] 52 7.19134e-03 5.99919e-02 -56.58979
[Matter] 53 2.30165e-03 2.68956e-02 -56.58981
[Matter] 54 1.02601e-03 1.38612e-02 -56.58982
[Matter] 55 2.07252e-03 1.18351e-02 -56.58982
[Matter] 56 2.36135e-03 1.53048e-02 -56.58984
[Matter] 57 6.13896e-03 3.14988e-02 -56.58984
[Matter] 58 1.52674e-03 1.01827e-02 -56.58985
[Matter] 59 1.27621e-03 9.09600e-03 -56.58986
Minimization converged within tolerence
Saving result to min.con
Final Energy: -56.5898551940918
Timing Information:
Real time: 3.534 seconds
User time: 6.030 seconds
System time: 0.364 seconds
eOn Client
2.12.0 (d586950)
compiled Sun Mar 08 01:19:26 PM GMT 2026
OS: linux
Arch: x86_64
Hostname: runnervm46oaq
PID: 2971
DIR: /home/runner/work/atomistic-cookbook/atomistic-cookbook/examples/eon-pet-neb/min_product
Loading parameter file config.ini
* [Main] job: minimization
* [Main] random_seed: 706253457
* [Potential] potential: Metatomic
* [Optimizer] opt_method: lbfgs
* [Optimizer] converged_force: 0.01
* [Optimizer] max_iterations: 2000
* [Optimizer] max_move: 0.1
* [Metatomic] model_path: /home/runner/work/atomistic-cookbook/atomistic-cookbook/examples/eon-pet-neb/models/pet-mad-xs-v1.5.0.pt
[MetatomicPotential] Initializing...
[MetatomicPotential] Loading model from '/home/runner/work/atomistic-cookbook/atomistic-cookbook/examples/eon-pet-neb/models/pet-mad-xs-v1.5.0.pt'
[MetatomicPotential] Using device: cpu
[MetatomicPotential] Using dtype: float32
[MetatomicPotential] Initialization complete.
Beginning minimization of pos.con
[Matter] Iter Step size ||Force|| Energy
[MetatomicPotential] Atomic numbers changed, re-creating types tensor.
[Matter] 0 0.00000e+00 3.06158e+00 -57.21474
[Matter] 1 2.39016e-02 1.77853e+00 -57.26960
[Matter] 2 6.46355e-03 1.06725e+00 -57.28843
[Matter] 3 1.52115e-02 3.87136e-01 -57.30855
[Matter] 4 6.73302e-03 3.07369e-01 -57.30998
[Matter] 5 1.56556e-03 1.61397e-01 -57.31073
[Matter] 6 3.81971e-03 1.37591e-01 -57.31147
[Matter] 7 4.66163e-03 1.11214e-01 -57.31198
[Matter] 8 4.21086e-03 6.02403e-02 -57.31223
[Matter] 9 1.41723e-03 6.77039e-02 -57.31226
[Matter] 10 4.58927e-04 2.92407e-02 -57.31230
[Matter] 11 6.08592e-04 2.39671e-02 -57.31232
[Matter] 12 1.38596e-03 3.20175e-02 -57.31234
[Matter] 13 2.56312e-03 4.41857e-02 -57.31240
[Matter] 14 8.27994e-03 9.74039e-02 -57.31249
[Matter] 15 8.20212e-03 7.47819e-02 -57.31265
[Matter] 16 1.07967e-02 7.45609e-02 -57.31285
[Matter] 17 1.69010e-02 6.13897e-02 -57.31311
[Matter] 18 8.78472e-03 6.84289e-02 -57.31321
[Matter] 19 3.27224e-03 3.58288e-02 -57.31327
[Matter] 20 8.90267e-04 2.36595e-02 -57.31329
[Matter] 21 1.07312e-03 2.46566e-02 -57.31331
[Matter] 22 2.84623e-03 3.60104e-02 -57.31335
[Matter] 23 6.35455e-03 6.20421e-02 -57.31339
[Matter] 24 3.80874e-03 4.38551e-02 -57.31342
[Matter] 25 1.85012e-03 3.53511e-02 -57.31346
[Matter] 26 4.79637e-03 3.53574e-02 -57.31352
[Matter] 27 9.43590e-03 4.98727e-02 -57.31364
[Matter] 28 2.34210e-02 1.02105e-01 -57.31388
[Matter] 29 4.09174e-02 1.57146e-01 -57.31428
[Matter] 30 4.02217e-02 1.78874e-01 -57.31490
[Matter] 31 6.60623e-02 1.48779e-01 -57.31609
[Matter] 32 4.83925e-02 2.59226e-01 -57.31651
[Matter] 33 4.23671e-02 2.95475e-01 -57.31561
[Matter] 34 2.91962e-02 7.46349e-02 -57.31726
[Matter] 35 2.37041e-02 8.97803e-02 -57.31749
[Matter] 36 2.09002e-02 1.20788e-01 -57.31773
[Matter] 37 4.26262e-02 1.87694e-01 -57.31808
[Matter] 38 2.18288e-02 1.28266e-01 -57.31835
[Matter] 39 1.08579e-02 7.55075e-02 -57.31850
[Matter] 40 4.49118e-02 1.03403e-01 -57.31871
[Matter] 41 3.34855e-02 1.47937e-01 -57.31876
[Matter] 42 3.44553e-02 9.13546e-02 -57.31900
[Matter] 43 2.30040e-02 5.38286e-02 -57.31907
[Matter] 44 1.35690e-02 2.94944e-02 -57.31915
[Matter] 45 3.89008e-02 8.58173e-02 -57.31919
[Matter] 46 4.81503e-03 6.72985e-02 -57.31923
[Matter] 47 3.30005e-03 2.05615e-02 -57.31926
[Matter] 48 8.13624e-03 1.37025e-02 -57.31927
[Matter] 49 9.43128e-03 1.83992e-02 -57.31928
[Matter] 50 5.98552e-03 1.31633e-02 -57.31928
[Matter] 51 5.09006e-03 1.25965e-02 -57.31928
[Matter] 52 3.57024e-04 4.93817e-03 -57.31929
Minimization converged within tolerence
Saving result to min.con
Final Energy: -57.31928634643555
Timing Information:
Real time: 3.097 seconds
User time: 5.096 seconds
System time: 0.315 seconds
Additionally, the relative ordering must be preserved, for which we use IRA [4].
reactant = aseio.read(dir_reactant / "min.con")
product = aseio.read(dir_product / "min.con")
ira = ira_mod.IRA()
# Default value
kmax_factor = 1.8
nat1 = len(reactant)
typ1 = reactant.get_atomic_numbers()
coords1 = reactant.get_positions()
nat2 = len(product)
typ2 = product.get_atomic_numbers()
coords2 = product.get_positions()
print("Running ira.match to find rotation, translation, AND permutation...")
# r = rotation, t = translation, p = permutation, hd = Hausdorff distance
r, t, p, hd = ira.match(nat1, typ1, coords1, nat2, typ2, coords2, kmax_factor)
print(f"Matching complete. Hausdorff Distance (hd) = {hd:.6f} Angstrom")
# Apply rotation (r) and translation (t) to the original product coordinates
# This aligns the product's orientation to the reactant's
coords2_aligned = (coords2 @ r.T) + t
# Apply the permutation (p)
# This re-orders the aligned product atoms to match the reactant's atom order
# p[i] = j means reactant atom 'i' matches product atom 'j'
# So, the new coordinate array's i-th element should be coords2_aligned[j]
coords2_aligned_permuted = coords2_aligned[p]
# Save the new aligned-and-permuted structure
# CRUCIAL: Use chemical symbols from the reactant,
# because we have now permuted the product coordinates to match the reactant order.
product = reactant.copy()
product.positions = coords2_aligned_permuted
Running ira.match to find rotation, translation, AND permutation...
Matching complete. Hausdorff Distance (hd) = 1.001631 Angstrom
Finally we can visualize these with ASE.
view(reactant, viewer="x3d")
view(product, viewer="x3d")
fig, (ax1, ax2) = plt.subplots(1, 2)
plot_atoms(reactant, ax1, rotation=("-90x,0y,0z"))
plot_atoms(product, ax2, rotation=("-90x,0y,0z"))
ax1.text(0.3, -1, "reactant")
ax2.text(0.3, -1, "product")
ax1.set_axis_off()
ax2.set_axis_off()
References¶
Mazitov, A.; Bigi, F.; Kellner, M.; Pegolo, P.; Tisi, D.; Fraux, G.; Pozdnyakov, S.; Loche, P.; Ceriotti, M. PET-MAD, a Universal Interatomic Potential for Advanced Materials Modeling. arXiv March 18, 2025. https://doi.org/10.48550/arXiv.2503.14118.
Bigi, F.; Abbott, J. W.; Loche, P.; Mazitov, A.; Tisi, D.; Langer, M. F.; Goscinski, A.; Pegolo, P.; Chong, S.; Goswami, R.; Chorna, S.; Kellner, M.; Ceriotti, M.; Fraux, G. Metatensor and Metatomic: Foundational Libraries for Interoperable Atomistic Machine Learning. arXiv August 21, 2025. https://doi.org/10.48550/arXiv.2508.15704.
Goswami, R. Efficient Exploration of Chemical Kinetics. arXiv October 24, 2025. https://doi.org/10.48550/arXiv.2510.21368.
Gunde, M.; Salles, N.; Hémeryck, A.; Martin-Samos, L. IRA: A Shape Matching Approach for Recognition and Comparison of Generic Atomic Patterns. J. Chem. Inf. Model. 2021, 61 (11), 5446–5457. https://doi.org/10.1021/acs.jcim.1c00567.
Smidstrup, S.; Pedersen, A.; Stokbro, K.; Jónsson, H. Improved Initial Guess for Minimum Energy Path Calculations. J. Chem. Phys. 2014, 140 (21), 214106. https://doi.org/10.1063/1.4878664.
Goswami, R; Gunde, M; Jónsson, H. Enhanced climbing image nudged elastic band method with hessian eigenmode alignment, Jan. 22, 2026, arXiv: arXiv:2601.12630. doi: 10.48550/arXiv.2601.12630.
R. Goswami, Two-dimensional RMSD projections for reaction path visualization and validation, MethodsX, p. 103851, Mar. 2026, doi: 10.1016/j.mex.2026.103851.
Total running time of the script: (1 minutes 10.764 seconds)