Skip to content

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:

  • nudged elastic band method 1 2
  • string method 3
  • growing string method 4

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:

Ammonia inversion animation

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.in
  • ini.in
  • fin.in
  • control.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

atom         0.00000000        0.00000000       0.00000000  N
    constrain_relaxation .true.
atom         0.81701754       -0.47169290       0.39055161  H
atom        -0.81701754       -0.47169290       0.39055161  H
atom         0.00000000        0.94340830       0.39055161  H
atom         0.00000000        0.00000000       0.00000000  N
    constrain_relaxation .true.
atom         0.81701754       -0.47169290       -0.39055161  H
atom        -0.81701754       -0.47169290       -0.39055161  H
atom         0.00000000        0.94340830       -0.39055161  H
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:

sbatch submit_aimschain.sh

After starting aimsChain, several new files and folders will be generated:

forces.log
iterations/
paths/

The file forces.log contains the residual forces of the current relaxation iteration. The content should look similar to this:

iteration0000       3.6273833086
iteration0001       1.6940221850
iteration0002       0.2063133275
iteration0003       0.0636575443
iteration0004       0.0297842837
iteration0005       0.0112770804
iteration0006       0.0022046999

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:

python plot_energy_path.py optimized -l "Optimized"

You should get a path in energy_path.png that looks similar to the following figure:

Ammonia-Inversion.png

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.

Show solutions


References