Introduction to aimsChain
The general problem that aimsChain is trying to solve is described by a chain of states that consists of a set of distinct configurations of structures (referred to as images) usually spanned between two minima on a potential energy surface. The goal is to relax the chain of states in a way that it coincides with the minimum energy path.
In aimsChain, several methods are implemented to perform the relaxation of the chain of states to find a good approximation for the minimum energy path:
To relax each image onto the minimal energy path, the computation of forces for each atom are needed. aimsChain uses FHI-aims to perform these calculations. As for FHI-aims, aimsChain can handle periodic and non-periodic structures.
To learn about aimsChain, we will look into the inversion of the ammonia molecule as an example:
Initial setup
Info
This is a concise tutorial created for the FHI-aims webinar series given in January 2025. It introduces participants to finding minimum energy pathways with aimsChain and covers the most important aspects, but it is far from exhaustive. For learning more about these kinds of calculations, we strongly recommend to also check the detailed description of aimsChain in the FHI-aims manual.
aimsChain requires several input files:
chain.inini.infin.incontrol.in
The file chain.in is the main input file for aimsChain and specifies the runtime settings for aimsChain. The files ini.in and fin.in contain the initial and final geometry of your structure, respectively, in the FHI-aims input format of geometry.in. However, you can give them individual names (see below). As usual, the file control.in sets the runtime parameters for FHI-aims. The same control.in file will be used for all images of your aimsChain calculation.
Let's take a closer look into the syntax of chain.in. Here is the one we will use for our example:
####################################################
## Mandatory settings
####################################################
#specify the command used for single point calculation
run_aims mpirun -n 8 aims.x
#initial geometry
initial_file ini.in
#final geometry
final_file fin.in
####################################################
## General settings for task
####################################################
#Method to use, neb/string
method string
#number of images to use
n_images 9
#convergence criterion for forces, in eV/A
force_thres 0.01
The option run_aims specifies the run command for FHI-aims to be used. initial_file and final_file specify the name for the initial and the final geometry.in input file. The option method allows selection of the chain of states method. Here, the string method is selected. n_images specifies the number of images to be used between the initial and final configuration. The option force_thres sets the threshold for the maximum allowed residual force component for the entire chain of states.
Below we also provide the input file ini.in and fin.in, as well as the control.in
xc pbe
vdw_correction_hirshfeld
relativistic atomic_zora scalar
charge_mix_param 0.4
compute_forces .true.
################################################################################
#
# FHI-aims code project
# VB, Fritz-Haber Institut, 2009
#
# Suggested "light" defaults for N atom (to be pasted into control.in file)
# Be sure to double-check any results obtained with these settings for post-processing,
# e.g., with the "tight" defaults and larger basis sets.
#
################################################################################
species N
# global species definitions
nucleus 7
mass 14.0067
#
l_hartree 4
#
cut_pot 3.5 1.5 1.0
basis_dep_cutoff 1e-4
#
radial_base 35 5.0
radial_multiplier 1
angular_grids specified
division 0.2599 50
division 0.4601 110
division 0.5885 194
division 0.6503 302
# division 0.6939 434
# division 0.7396 590
# division 0.7632 770
# division 0.8122 974
# division 1.1604 1202
# outer_grid 974
outer_grid 302
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 2 s 2.
valence 2 p 3.
# ion occupancy
ion_occ 2 s 1.
ion_occ 2 p 2.
################################################################################
#
# Suggested additional basis functions. For production calculations,
# uncomment them one after another (the most important basis functions are
# listed first).
#
# Constructed for dimers: 1.0 A, 1.1 A, 1.5 A, 2.0 A, 3.0 A
#
################################################################################
# "First tier" - improvements: -1193.42 meV to -220.60 meV
hydro 2 p 1.8
hydro 3 d 6.8
hydro 3 s 5.8
# "Second tier" - improvements: -80.21 meV to -6.86 meV
# hydro 4 f 10.8
# hydro 3 p 5.8
# hydro 1 s 0.8
# hydro 5 g 16
# hydro 3 d 4.9
# "Third tier" - improvements: -4.29 meV to -0.53 meV
# hydro 3 s 16
# ionic 2 p auto
# hydro 3 d 6.6
# hydro 4 f 11.6
# "Fourth tier" - improvements: -0.75 meV to -0.25 meV
# hydro 2 p 4.5
# hydro 2 s 2.4
# hydro 5 g 14.4
# hydro 4 d 14.4
# hydro 4 f 16.8
# Further basis functions - -0.21 meV and below
# hydro 3 p 14.8
# hydro 3 s 4.4
# hydro 3 d 19.6
# hydro 5 g 12.8
################################################################################
#
# FHI-aims code project
# VB, Fritz-Haber Institut, 2009
#
# Suggested "light" defaults for H atom (to be pasted into control.in file)
# Be sure to double-check any results obtained with these settings for post-processing,
# e.g., with the "tight" defaults and larger basis sets.
#
################################################################################
species H
# global species definitions
nucleus 1
mass 1.00794
#
l_hartree 4
#
cut_pot 3.5 1.5 1.0
basis_dep_cutoff 1e-4
#
radial_base 24 5.0
radial_multiplier 1
angular_grids specified
division 0.2421 50
division 0.3822 110
division 0.4799 194
division 0.5341 302
# division 0.5626 434
# division 0.5922 590
# division 0.6542 770
# division 0.6868 1202
# outer_grid 770
outer_grid 302
################################################################################
#
# Definition of "minimal" basis
#
################################################################################
# valence basis states
valence 1 s 1.
# ion occupancy
ion_occ 1 s 0.5
################################################################################
#
# Suggested additional basis functions. For production calculations,
# uncomment them one after another (the most important basis functions are
# listed first).
#
# Basis constructed for dimers: 0.5 A, 0.7 A, 1.0 A, 1.5 A, 2.5 A
#
################################################################################
# "First tier" - improvements: -1014.90 meV to -62.69 meV
hydro 2 s 2.1
hydro 2 p 3.5
# "Second tier" - improvements: -12.89 meV to -1.83 meV
# hydro 1 s 0.85
# hydro 2 p 3.7
# hydro 2 s 1.2
# hydro 3 d 7
# "Third tier" - improvements: -0.25 meV to -0.12 meV
# hydro 4 f 11.2
# hydro 3 p 4.8
# hydro 4 d 9
# hydro 3 s 3.2
Due to the symmetry of the inversion process of the ammonia molecule, we can constrain the Nitrogen atom. Strictly, this is not necessary to get an estimation for the barrier height, but it will help to converge in less iterations towards the minimum energy path. By default, aimsChain linearly interpolates between start initial and final configuration (however it is also possible to provide the whole chain of states).
Start aimsChain
To start aimsChain on the provided resources during the hands-on, write a new submission script submit_aimschain.sh with the following:
#!/bin/bash -l
# Standard output and error:
#SBATCH -o ./Output
#SBATCH -e ./Error
# Initial working directory:
#SBATCH -D ./
# Job Name:
#SBATCH -J Jobscript
# Number of nodes and MPI tasks per node:
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=8
#
# Wall clock limit:
#SBATCH --time=00:10:00
module purge
module load fhi-aims
ulimit -s unlimited
export OMP_NUM_THREADS=1
export MKL_NUM_THREADS=1
runchain
Running calculations on other resources
If you are running the calculations on a different machine, make sure to adjust the submission script specific to your system and submission system.
Submit the calculation with:
After starting aimsChain, several new files and folders will be generated:
The file forces.log contains the residual forces of the current relaxation iteration. The content should look similar to this:
The unit of the residual forces in the second column are given in eV/Angstrom.
The folder iterations/ contains the FHI-aims raw data for each iteration and for each node that has been computed. Note that the initial and final structures, ini.in and fin.in, are only computed once in the first iteration, since they are kept fixed and not optimized.
The folder paths/ contains all the structure files per path and iteration as well as the parsed total energy from the FHI-aims output file stored in ener.lst.
By default, you can (abort and) restart the aimsChain run command at any time by just retyping runchain and aimsChain will continue from the latest unfinished calculation.
Analyzing the results
Once aimsChain has finished, the final results are stored in the folder optimized, which contains all FHI-aims output files per image, as well as a summarized list of energies and structures per image. To analyze the minimum energy path create the following script plot_energy_path.py:
"""
Plot energy profiles from NEB (Nudged Elastic Band) or string method calculations from AimsChain.
This script reads energy profiles from one or more calculation directories and creates
a comparative plot. For each directory, it reads:
- ener.lst: Contains energy values and point types (FIXED/NORMAL)
- path.lst: Contains reaction coordinates for each point
- path.xyz: Contains atomic coordinates for calculating reaction coordinates
Usage:
python plot_energy_path.py path/to/calculation1 [path/to/calculation2 ...]
python plot_energy_path.py path/to/calculation1 -o custom_name.png
python plot_energy_path.py calc1 calc2 calc3 -l "Path 1" "Path 2" "Path 3"
Options:
-o, --output Specify custom output filename (default: energy_path.png)
-l, --labels Optional labels for the paths (must match number of folders)
"""
import numpy as np
import matplotlib.pyplot as plt
import argparse
import os
from scipy.interpolate import CubicSpline
def parse_arguments():
parser = argparse.ArgumentParser(description="Plot NEB/String energy profile")
parser.add_argument(
"folders",
type=str,
nargs="+",
help="Paths to folders containing ener.lst and path.lst",
)
parser.add_argument(
"--output",
"-o",
type=str,
default="energy_path.png",
help="Output filename (default: energy_path.png)",
)
parser.add_argument(
"--labels",
"-l",
type=str,
nargs="+",
help="Optional labels for the paths in the legend",
)
args = parser.parse_args()
if args.labels and len(args.labels) != len(args.folders):
parser.error(
f"Number of labels ({len(args.labels)}) must match number of folders ({len(args.folders)})"
)
return args
def read_xyz_trajectory(filename):
"""Read atomic coordinates from an xyz trajectory file."""
trajectory = []
n_atoms = None
current_frame = []
with open(filename, "r") as f:
for line in f:
if n_atoms is None:
n_atoms = int(line.strip())
continue
if "Energy:" in line:
if current_frame:
trajectory.append(np.array(current_frame))
current_frame = []
continue
if len(line.split()) == 4: # Atom line: symbol x y z
x, y, z = map(float, line.split()[1:])
current_frame.append([x, y, z])
if current_frame:
trajectory.append(np.array(current_frame))
return trajectory
def calculate_reaction_coordinate(trajectory):
"""Calculate normalized reaction coordinates based on RMSD between images."""
n_frames = len(trajectory)
distances = np.zeros(n_frames)
# Calculate RMSD between consecutive images
for i in range(1, n_frames):
diff = trajectory[i] - trajectory[i - 1]
rmsd = np.sqrt(np.mean(np.sum(diff**2, axis=1)))
distances[i] = rmsd
# Calculate cumulative distance
cumulative = np.cumsum(distances)
# Normalize to [0,1] range
reaction_coordinates = cumulative / cumulative[-1]
return reaction_coordinates
def read_folder_data(folder_path):
"""Read energy data and fixed points from a folder."""
ener_file = os.path.join(folder_path, "ener.lst")
if not os.path.exists(ener_file):
raise FileNotFoundError(f"Could not find ener.lst in {folder_path}")
# Read the energy data
energies = []
fixed_points = []
with open(ener_file, "r") as f:
next(f) # Skip header
for line in f:
columns = line.split()
energies.append(float(columns[1]))
fixed_points.append(columns[2].strip() == "FIXED")
energies = np.array(energies)
fixed_points = np.array(fixed_points)
# Align energies to minimum
min_energy = np.min(energies)
relative_energies = energies - min_energy
return relative_energies, fixed_points, min_energy
def plot_energy_profiles(folders, output_file, labels=None):
"""Create energy profile plot with smooth interpolation."""
plt.figure(figsize=(10, 6))
colors = plt.cm.tab10(np.linspace(0, 1, len(folders)))
all_coordinates = []
all_energies = []
for i, (folder_path, color) in enumerate(zip(folders, colors)):
try:
# Read energy data
energies, fixed, min_energy = read_folder_data(folder_path)
# Read trajectory and calculate reaction coordinate
xyz_file = os.path.join(folder_path, "path.xyz")
trajectory = read_xyz_trajectory(xyz_file)
reaction_coordinate = calculate_reaction_coordinate(trajectory)
# Create spline interpolation
cs = CubicSpline(reaction_coordinate, energies)
smooth_coords = np.linspace(0, 1, 200)
smooth_energies = cs(smooth_coords)
plt.plot(
smooth_coords,
smooth_energies,
"-",
color=color,
linewidth=2,
label=labels[i] if labels and i < len(labels) else folder_path,
)
plt.plot(
reaction_coordinate,
energies,
"o",
color=color,
markersize=6,
markerfacecolor=color,
)
# Add fixed points with different marker
fixed_coords = reaction_coordinate[fixed]
fixed_energies = energies[fixed]
plt.plot(
fixed_coords,
fixed_energies,
"s",
color=color,
markersize=8,
markeredgecolor="black",
markerfacecolor=color,
)
all_coordinates.extend(reaction_coordinate)
all_energies.extend(energies)
print(f"\nEnergy profile summary for {folder_path}:")
print(f"Minimum energy: {min_energy:.6f} eV")
print(f"Maximum barrier: {np.max(energies):.6f} eV")
except Exception as e:
print(f"Error processing {folder_path}: {e}")
continue
plt.xlabel("Reaction Coordinate", fontsize=15)
plt.ylabel("Relative Energy [eV]", fontsize=15)
plt.grid(True, linestyle="--", alpha=0.7)
plt.legend(loc="upper right", fontsize=15)
# Set axis limits with padding
x_padding = 0.05
y_padding = (max(all_energies) - min(all_energies)) * 0.1
plt.xlim(-x_padding, 1 + x_padding)
plt.ylim(min(all_energies) - y_padding, max(all_energies) + y_padding)
plt.savefig(output_file, dpi=300, bbox_inches="tight")
print(f"\nPlot saved as: {output_file}")
plt.close()
def main():
args = parse_arguments()
try:
plot_energy_profiles(args.folders, args.output, args.labels)
except Exception as e:
print(f"Error: {e}")
return 1
return 0
if __name__ == "__main__":
main()
and execute:
You should get a path in energy_path.png that looks similar to the following figure:
Here, we find an energy barrier for the ammonia inversion of about 0.236 eV, which compares well with the reference value of 0.25 eV 5. You can also use GIMS (and other software of your choice) to visualize the individual structures of each image.
Solutions
You find the solutions to the above exercise by clicking on the button below.
References
-
G. Mills K. W. Jacobsen H. Jonsson. Nudged Elastic Band Method for Finding Minimum Energy Paths of Transitions. World Scientific, 1998. ↩
-
G. Henkelman, B. P. Uberuagga, and H. Jonsson. A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J. Chem. Phys., 113:9901, 2000. ↩
-
E Weinan, Weiqing Ren, and Eric Vanden-Eijnden. Simplified and improved string method for computing the minimum energy paths in barrier-crossing events. J. Chem. Phys., 126:164103, 2007. ↩
-
Baron Peters, Andreas Heyden, Alexis T. Bell, and Arup Chakraborty. A growing string method for determining transition states: Comparison to the nudged elastic band and string methods. The Journal of Chemical Physics, 120(17):7877–7886, 2004. ↩
-
A. Rauk and L. C. Allen. Electronic Structure and Inversion Barrier of Ammonia. J. Chem. Phys. 52(8):4133-4144, 1970. ↩