Nudged elastic band (NEB) is designed to find a discrete representation of the minimum energy path (MEP) connecting a initial state minimum to a final state minimum, not just the highest point in the energy surface. In practice, some unexpected intermediate minim might be found during a NEB calculation.

• Nudging means force projection
• Elastic band refers to the harmonic spring connecting neighboring images

Trying to understand

• What should each image along the minimum energy path (MEP) looks like? (energy, force, ...)
• Why is force projection needed?
• How to define the spring force?
• How to estimate the local tangent of one image?
• How to distribute images along the path?

Initial Path

A guess of initial path affects - The number of iterations for reaching a converged MEP - Detailed transition mechanism (since there might be different pathways connecting the same set of IS and FS)

Linear interpolation

The most straightforward way to generate the initial pathway is to linearly interpolate between the Cartesian coordinates of initial and final state. $$ \textbf{r}^k_i = \textbf{r}^\text{IS} _ i + \dfrac{k}{\text{N}_\text{images}} (\textbf{r}^\text{FS}_i-\textbf{r}^\text{IS}_i) $$

• $k$ is the image index (1, 2, ..., $\text{N}_\text{images}-1$)
• $i$ is the atom index
• $\text{N}_\text{images}$ include intermediate and final images

The linear interpolation is easy to implement for both isolated and periodic systems. However, it does not guarantee a physically meaningful pathway, meaning that some atoms may stay too close/far to each other in one image. The computational cost could be high if a bad guess of pathway is used in a NEB calculation. Sometimes it is helpful to use chemical knowledge to build a intermediate image $\mu$ and do linear interpolation from $\textbf{r}^\text{IS}$ to $\textbf{r}^\mu$, then from $\textbf{r}^\mu$ to $\textbf{r}^\text{FS}$.

Image dependent pair potential (IDPP)

Inspired by Halgren and Lipscomb 1977, Smidstrup et al. proposed a method of generating initial transition pathways based on image dependent pair potential (IDPP) in J. Chem. Phys. 140(21), 214106 (2014). The method uses bond distance instead of direct Cartesian coordinates for interpolation.

$$ d_{ij}^k = r_{ij}^\text{IS} + \dfrac{k}{\text{N} _ \text{images}} (r_{ij}^\text{FS} - r_{ij}^\text{IS}) $$

$d_{ij}^k$ is the target distance between atom $i$ and atom $j$ in image $k$. There are many pairs in a image. If one pair has changed its distance, other distances will change as well. Therefore not all distances can meet $d_{ij}^k$ and a compromise needs to be made. In order to make as much as possible distances meet the ideal value $d_{ij}^k$, a object function can be defined by summing the squared deviation of pair distances $r_{ij}^k$ from the target value $d_{ij}^k$

$$ S_k^\text{IDPP}(r) = \sum_i^{\text{N} _ \text{atoms}} \sum_{j>i}^{\text{N} _ \text{atoms}} w (d_{ij}^k) (d_{ij}^k - r_{ij}^k)^2 $$

• $w$ is a weighting function that weights more on shorter distance, e.g. $w(d) = \dfrac{1}{d^4}$ since the energy rises strongly when atoms become too close to each other.
• $F_i^k(\mathbf{r}) = - \dfrac{\partial S_k^\text{IDPP}(\mathbf{r})}{\partial \mathbf{r}}$ can be used to find the optimal path on the IDPP surface.

In practice, IDPP is used in this way

1. Generating initial path from a linear interpolation
2. NEB optimization of previous path on the IDPP surface and get the IDPP path
3. Use the IDPP path as an initial path for the NEB optimization with DFT or other electronic methods

IDPP provides a pathway closer to the minium energy path than a linear interpolation of the Cartesian coordinates. - Saving computational time - Avoiding SCF issue in a nonphysical structure

Standard NEB

All images except starting and end points are moving at the same time. When forces acting on all images are converged, it's still possible that all images have lower energy than the real transition state (TS), meaning none of these images is TS. The reaction barrier is often estimated by the energy difference of the image with the highest energy and the initial state. Such a barrier ofter underestimates the real reaction barrier, because there is no guarantee that one of the image lies on the saddle point. If enough images (e.g. > 15 images, system-dependent) are inserted between the initial and final state, the image of the highest energy may relax into the saddle point. More images means smaller difference between the estimated and real barrier.

Climbing image (CI)

CI-NEB tries to relax the image with the highest energy to the saddle point. All inserted images are still relaxing simultaneously. In most cases, TS can be located easily, but the reaction pathway might be destroyed if few images are used.

force analysis

import numpy as np
import matplotlib.pyplot as plt

n_atoms = 5
n_dimension = 3
def get_data(first_seed=22, second_seed=877, force_amplitude=1):

    np.random.seed(first_seed)
    force = force_amplitude*(np.random.random((n_atoms, n_dimension)) - 0.5)
    np.random.seed(second_seed)
    tangent = np.random.random((n_atoms, n_dimension)) * np.random.random_integers(1, 5, (n_atoms, n_dimension)) - np.random.random_integers(2, 12, (n_atoms, n_dimension))
    tangent /= np.linalg.norm(tangent)
    force_CI = force - 2*np.vdot(force, tangent)*tangent

    return force, force_CI

def axes_plot(axes, force, force_CI, title):
    old_fmax_all = np.sqrt((force ** 2).sum(axis=1))
    new_fmax_all = np.sqrt((force_CI ** 2).sum(axis=1))
    y1 = force.flatten()
    y2 = force_CI.flatten()

    x = np.arange(len(y1))
    axes[0].plot(x, y1, '-o', label=f'norm(old_f) = {np.linalg.norm(y1):.3f}')
    axes[0].plot(x, y2, '-o', label=f'norm(new_f) = {np.linalg.norm(y2):.3f}')
    axes[0].legend()
    axes[0].set_title(title)

    x_force = n_dimension*np.arange(n_atoms) + 1
    axes[1].plot(x_force, old_fmax_all, '-o', label=f'old_fmax = {old_fmax_all.max():.3f}')
    axes[1].plot(x_force, new_fmax_all, '-o', label=f'new_fmax = {new_fmax_all.max():.3f}')
    axes[1].legend()
    axes[1].set_xlabel('index')

    if title == 'original':
        axes[0].set_ylabel('force component')
        axes[1].set_ylabel('atomic force magnitude')

fig, axes = plt.subplots(2, 3, sharex=True, figsize=(12, 4))

axes_plot([axes[0, 0], axes[1, 0]], *get_data(), 'original')
axes_plot([axes[0, 1], axes[1, 1]], *get_data(second_seed=1), 'different tangent direction')
axes_plot([axes[0, 2], axes[1, 2]], *get_data(force_amplitude=8), 'different DFT force magnitude')

title = f'Forces on the climbing image (CI): old(DFT), new(DFT + parallel term)\n'
title += r'$F_{CI} = F_{DFT} - 2F_{DFT}\vert_\parallel$'
plt.suptitle(title, y=1.1)
plt.savefig('fig-cGOFEE-NEB-tangent.png', dpi=200, bbox_inches='tight')
#plt.show()

AutoNEB

In standard NEB, many computational times might be wasted in relaxing two images near TS in some complex cases, because both of them are hard to converge, while the rest images are easy to converge. A possible explanation is that the "gap" between these two images is huge. More images between these two images are needed to have a better description of this region of the potential energy surface (PES), instead of wasting time trying to converge these two images. An accurate description of TS (strict force criteria, only one imaginary mode) is beneficial for description reaction mechanisms, while rough descriptions of other images (loose force criteria) are enough to have an idea of what's going on.

AutoNEB tries to identify two images near TS and insert more images between these two on the fly. The rest images are fixed after inserting new images, thus reducing computational time.