Table of Contents

Nudged elastic band (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):

    force = force_amplitude*(np.random.random((n_atoms, n_dimension)) - 0.5)
    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}')

    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}')

    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')


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.