# A modified nudged elastic band algorithm with adaptive spring lengths

D. Mandelli  M. Parrinello<EMAIL_ADDRESS>Atomistic Simulations,
Italian Institute of Technology, Via Morego, 30 16163 Genova, Italy.

###### Abstract

We present a modified version of the nudged elastic band (NEB) algorithm to
find minimum energy paths connecting two known configurations. We show that
replacing the harmonic band-energy term with a discretized version of the
Onsager-Machlup action leads to a NEB algorithm with adaptive spring lengths
that automatically increase the resolution of the minimum energy path around
the saddle point of the potential energy surface. The method has the same
computational cost per optimization step of the standard NEB algorithm and
does not introduce additional parameters. We present applications to the
isomerization of alanine dipeptide, the elimination of hydrogen from ethane
and the healing of a 5-77-5 defect in graphene.

††preprint: AIP/123-QED

##  I  Introduction

An important problem in chemistry and materials science is the computation of
kinetic rates for transitions between metastable states. Typical examples are
chemical reactions, molecular conformational changes and diffusion events in
solids. All these transformations are thermally activated processes that occur
via a concerted rearrangement of atoms that bring the system from one state to
another. During the transformation, the system has to pass through high-energy
configurations that represent the kinetic bottleneck of the reaction. Although
these phenomena can be often described using classical mechanics, a direct
simulation of the dynamics of reaction events is usually impossible because in
the presence of energy barriers larger than $k_{B}T$ unfeasibly long
simulation times are needed to observe the transitions of interest. Several
advanced algorithms based on molecular dynamics have been developed to compute
kinetic rates Dellago  _et al._ (1998); van Erp, Moroni, and Bolhuis (2003);
Faradjian and Elber (2004); Tiwary and Parrinello (2013); Debnath and
Parrinello (2020); Mandelli, Hirshberg, and Parrinello (2020), however, they
can be computationally expensive, especially when the electronic degrees of
freedom play an important role and quantum mechanical approaches are needed to
compute the atomic forces. In all the cases where these direct methods become
impractical, transition state theory Eyring (1935); Wigner (1938); Keck (2007)
(TST) provides a cheaper route to compute accurate estimates of transition
rates. TST is a purely statistical approach where rates are expressed in terms
of quantities that are available directly from the potential energy surface
(PES).

When applying TST, the most challenging step is the determination of the
relevant saddle points in the multidimensional potential energy landscape.
Among the algorithms that have been designed to solve this problem one can
distinguish between local methods, like the dimer method Henkelman and Jónsson
(1999) and the activation-relaxation-technique Jay  _et al._ (2020), and the
so-called chain-of-states methods Jonsson, Mills, and Jacobsen (1998); Peters
_et al._ (2004); E and Vanden-Eijnden (2010). The first class of algorithms
uses only local information of the PES to follow the energy landscape uphill
to the saddle point and then downhill to a new minimum. These methods are
particularly useful when the final state is not known. On the other hand, when
both the initial and final states are known, it is generally more effective to
adopt chain-of-states approaches. Here, we focus on the popular nudged elastic
band (NEB) algorithm Henkelman and Jónsson (2000); Henkelman, Uberuaga, and
Jónsson (2000), which is an efficient method to find the minimum energy path
(MEP) connecting two endpoint configurations. By definition, the force acting
on the atoms is parallel to the MEP, while the energy is stationary along any
direction perpendicular to it. In the continuum limit, this is equivalent to
the requirement

$\nabla E({\bf R}_{\lambda})|_{\perp}=0$ (1)

where E is the potential energy, $R_{\lambda}$ is a path in configuration
space parameterized by $\lambda\in[0,1]$, with fixed endpoints ${\bf R}_{0}$
and ${\bf R}_{1}$, and only the component of the gradient perpendicular to the
path is considered. It follows from equation (1) that the PES saddle points
correspond to the maxima of the potential energy along the MEP.

In the NEB algorithm, an initial guess of the path connecting the two
endpoints is discretized as a sequence of $N$ images (or replicas) of the
system. The images are optimized in a concerted way so as to obtain a path
that satisfies equation (1). The configuration corresponding to the maximum of
the potential energy along the MEP identifies the saddle point of interest.
Generally, one would like to have a higher density of images in the region of
the saddle point since this can improve the accuracy with which the saddle
point is obtained by improving the resolution of the discretized MEP. However,
the standard NEB algorithm generates replicas that are equally spaced along
the path Jonsson, Mills, and Jacobsen (1998); Henkelman and Jónsson (2000). A
straightforward way to improve the sampling near the saddle is to increase the
overall number of images. Still, this approach would lead to most of the
computational time being wasted in optimizing images in irrelevant regions far
from the saddle. A more efficient solution is provided by algorithms that keep
the number of images as small as possible by focusing on increasing the
resolution only around the saddle Henkelman, Uberuaga, and Jónsson (2000);
Maragakis  _et al._ (2002); Kolsbjerg, Groves, and Hammer (2016). This is
particularly important in expensive quantum mechanical calculations, where the
number of replicas is typically limited to the range $5<N<20$.

In this work, we show how a relatively simple modification of the inter-
replica harmonic forces used in the standard NEB algorithm automatically
increases the MEP resolution around the saddle point, without the need to
increase the number of images. The method takes inspiration from a discretized
version of the Onsager-Machlup action Onsager and Machlup (1953), which leads
to a modified NEB algorithm with adaptive natural spring lengths that become
automatically shorter around the stationary points of the PES. We present
applications to the isomerization of alanine dipeptide, the elimination of
hydrogen from ethane and the healing of a 5-77-5 defect in graphene. In all
these cases, the modified algorithm reproduces the MEPs and saddle points
obtained via standard NEB calculations while improving the resolution around
the saddle point.

The article is organized as follows. In Section II.1, we review the standard
NEB algorithm. In Section II.2 we introduce our modification. In Section II.3,
we describe the model systems and the simulation protocols used to test the
method. In Section III, we present the results of simulations. In Section IV,
we present our final remarks.

##  II Materials and Methods

###  II.1 The nudged elastic band algorithm

When discretized into $N$ steps, a path is described by the ordered set
$\\{{\bf R}_{1},{\bf R}_{2},\dots,{\bf R}_{N}\\}$ of replicas of the system,
where we denote with ${\bf R}_{i}$ $i=1,\dots,N$, the system atomic
coordinates of the different images. The initial and final positions ${\bf
R}_{1}$ and ${\bf R}_{N}$ are set at the local minima of the potential energy.
The NEB algorithm optimizes the $N-2$ intermediate configurations by putting
to zero the forces

${\bf F}_{i}=-\nabla E({\bf R}_{i})|_{\perp}+{\bf F}_{i}^{s}|_{||}+{\bf
F}_{i}^{s^{\prime}}|_{\perp}$ (2)

acting on each image. The first term

$\nabla E({\bf R}_{i})|_{\perp}=\nabla E({\bf R}_{i})-\nabla E({\bf
R}_{i})\cdot{\bf\tau}_{i}{\bf\tau}_{i},$ (3)

is the component of the physical force acting perpendicular to the path at
image $i$. This term is a discretized version of equation (1) and it ensures
convergence to the MEP. In equation (3), ${\bf\tau}_{i}$ is the unit vector
tangent to the path at image $i$ and ${\bf\tau}_{i}{\bf\tau}_{i}$ indicates
the outer product. In all our calculations, ${\bf\tau}_{i}$ was computed using
the definition of Ref.Henkelman and Jónsson, 2000. The second term in equation
(2) is a fictitious spring force

${\bf F}_{i}^{s}|_{||}=k_{||}(|{\bf R}_{i}-{\bf R}_{i-1}|-|{\bf R}_{i+1}-{\bf
R}_{i}|){\bf\tau}_{i},$ (4)

acting parallel to the path. This elastic constraint is necessary when dealing
with discretized paths in order to prevent the images from sliding downhill
and to keep them equally spaced. The spring constant $k_{||}$ is a user-
defined parameter. Because only the component parallel to the path is used,
typically the results are weakly dependent on the choice of $k_{||}$. The last
term is an additional harmonic-like force given by

${\bf F}_{i}^{s^{\prime}}|_{\perp}=f(\phi_{i})({\bf F}_{i}^{s^{\prime}}-{\bf
F}_{i}^{s^{\prime}}\cdot{\bf\tau}_{i}{\bf\tau}_{i}),$ (5)

where

${\bf F}_{i}^{s^{\prime}}=k_{\perp}({\bf R}_{i+1}+{\bf R}_{i-1}-2{\bf
R}_{i}),$ (6)

and

$f(\phi_{i})=\begin{cases}\frac{1}{2}(1+\cos(\pi\cos(\phi_{i})))&\text{ if
$0<\phi_{i}<\pi/2$}\\\ 1&\text{if $\phi_{i}>\pi/2$}\end{cases}$ (7)

is a switching function that depends on the path angle at image $i$ Jonsson,
Mills, and Jacobsen (1998); Maras  _et al._ (2016):

$\cos\left(\phi_{i}\right)=\frac{({\bf R}_{i+1}-{\bf R}_{i})\cdot({\bf
R}_{i}-{\bf R}_{i-1})}{|{\bf R}_{i+1}-{\bf R}_{i}||{\bf R}_{i}-{\bf
R}_{i-1}|}$ (8)

We note that the tangent introduced in Ref.Henkelman and Jónsson, 2000 was
specifically designed to prevent oscillations of the path without the need of
including the perpendicular spring force of equation (5). However, as recently
observed in Ref.Maras  _et al._ , 2016, keeping part of the perpendicular
spring force can help improving the convergence. In this work, we have found
it useful in some cases to include this term for few initial steps of the
minimization. In all cases, to avoid systematic corner cutting problems
Jonsson, Mills, and Jacobsen (1998), the final path was eventually obtained
switching off this term by setting $k_{\perp}=0$.

Once the condition ${\bf F}_{i}=0$ is met, the sequence of images lie along
the MEP. However, usually no image relaxes exactly on top of the saddle point.
Hence, after a first relaxation, a typical NEB calculation continues with a
second “climbing image” optimization Henkelman, Uberuaga, and Jónsson (2000)
during which the force acting on the replica corresponding to the maximum
value of the potential energy is changed to

${\bf F}_{i}^{\rm climber}=-\nabla E({\bf R}_{i})+2\nabla E({\bf
R}_{i})\cdot{\bf\tau}_{i}{\bf\tau}_{i}.$ (9)

This is the physical force with the component parallel to the path changed in
sign. The effect of this force is to push the climber uphill and towards the
saddle point. The climbing image is subjected only to the force defined in
equation (9) and it does not feel any spring force, either parallel or
perpendicular. Hence, provided the estimate of the tangent to the MEP
(${\bf\tau}_{i}$) is accurate enough, the climbing image will converge to the
exact saddle.

###  II.2 Modified nudged elastic band algorithm

Equation (4) can be interpreted as the force arising from a fictitious
harmonic interaction between the replicas

$S_{\rm NEB}=\sum_{i=1}^{N-1}\frac{k_{||}}{2}|{\bf R}_{i+1}-{\bf R}_{i}|^{2}$
(10)

that keeps the images equally spaced along the path. In order to increase the
resolution of the MEP around the saddle point some variants of the NEB
algorithm have been proposed that consist in increasing automatically the
spring constant $k_{||}$ near the saddle point Henkelman, Uberuaga, and
Jónsson (2000); Ásgeirsson  _et al._ (2021) or in adding more images only in
its proximity Maragakis  _et al._ (2002); Kolsbjerg, Groves, and Hammer
(2016).

Here, we propose an alternative approach that exploits the fact that saddles
are stationary points of the PES. In designing this method, we were inspired
by the Onsager-Machlup (OM) action Onsager and Machlup (1953), $S_{\rm
OM}[{\bf R}(t)]$, that defines the probability, $p\propto e^{-S_{\rm
OM}/k_{B}T}$, of observing a stochastic trajectory ${\bf R}(t)$. In the case
of a discretized Brownian trajectory of $N$ steps, ${\bf R}_{1}\rightarrow{\bf
R}_{2}\rightarrow\dots\rightarrow{\bf R}_{N}$, connecting two endpoints ${\bf
R}_{1}$ and ${\bf R}_{N}$, the OM action is given by Mandelli, Hirshberg, and
Parrinello (2020)

$S_{\rm OM}=\sum_{i=1}^{N-1}\frac{m\nu}{4\Delta t}({\bf R}_{i+1}-{\bf
R}_{i}-{\bf L}_{i})^{2}.$ (11)

where $m$ is the diagonal mass tensor, $\nu$ is a damping coefficient and
$\Delta t$ is the time step of the discretized trajectory. Equation (11) has
the form of a harmonic interreplica interaction with variable natural spring
lengths

${\bf L}_{i}=-\frac{\Delta t}{m\nu}\nabla E({\bf R}_{i}).$ (12)

By minimizing the OM action, one generates the most probable discretized
trajectory under the assumption that the underlying dynamics is Brownian. In
this framework, the length of the path has a clear physical interpretation as
it defines its total duration, $N\Delta t$. The problem of sampling
statistically relevant dynamical trajectories based on the OM action has been
addressed elsewhere (see for example Refs. Mandelli, Hirshberg, and
Parrinello, 2020; Fujisaki, Shiga, and Kidera, 2010; Fujisaki  _et al._ ,
2013; Lee  _et al._ , 2017).

In this work, we are not interested in constructing real dynamical
trajectories, which typically requires the use of hundreds of images, but
rather in finding the MEP. Our modification to the NEB algorithm starts from
the observation that the natural spring length ${\bf L}_{i}$ vanishes at the
saddle points of the PES. This suggests that using equation (11) to define the
spring forces in a NEB-like calculation would naturally lead to an
accumulation of images around the saddle point. We therefore propose to
substitute the spring forces, ${\bf F}_{i}^{s}$ and ${\bf
F}_{i}^{s^{\prime}}$, defined in equations (4) and (6), with

${\bf F}_{i}^{\rm OM}=k_{\rm OM}({\bf R}_{i+1}+{\bf R}_{i-1}-2{\bf R}_{i}+{\bf
L}_{i-1}-{\bf L}_{i}).$ (13)

Following the NEB procedure, the MEP is found by putting to zero the forces

${\bf F}_{i}=-\nabla E({\bf R}_{i})|_{\perp}+{\bf F}_{i}^{\rm OM}|_{||}+{\bf
F}_{i}^{\rm OM}|_{\perp}.$ (14)

Where the first term is defined in equation (3), the second term is given by

${\bf F}_{i}^{\rm OM}|_{||}={\bf F}_{i}^{\rm
OM}\cdot{\bf\tau}_{i}{\bf\tau}_{i},$ (15)

and the last term is given by

${\bf F}_{i}^{\rm OM}|_{\perp}=f(\phi_{i})({\bf F}_{i}^{\rm OM}-{\bf
F}_{i}^{\rm OM}\cdot{\bf\tau}_{i}{\bf\tau}_{i}).$ (16)

At the end of the minimization, a final “climbing image” optimization is
performed as discussed in the last paragraph of Section II.1.

We note that, strictly speaking, equation (13) is not the exact gradient of
the OM action since we have neglected the dependence of the natural spring
length on the position. The exact gradient can be implemented using a finite
difference formula at the cost of two more force computations per optimization
step Mandelli, Hirshberg, and Parrinello (2020). Since our aim is to increase
the resolution of the MEP around the saddle, equation (13) is already a
reasonable choice, as demonstrated by our results. This avoids having to
compute second derivatives of the potential energy. With this implementation,
the number of force calculations per optimization step is the same as in the
standard NEB algorithm. We have checked that the results do not change when
using the exact gradient of equation (11).

Because we are not interested in the true dynamics of the system, the spring
constant $k_{\rm OM}$ loses its original dynamical interpretation and becomes
a user-defined parameter. Nevertheless, we observed that better results are
obtained defining the spring constant as $k_{\rm OM}=m\nu/2\Delta t$,
following equation (11). In all the simulations, we have used a value of
$\Delta t=1$ fs, which is a typical time step adopted in atomistic molecular
dynamics simulations. As a rule of thumb, we have found that a good choice of
$\nu$ corresponds to values satisfying $\nu\Delta t=1$.

###  II.3  Models and simulation protocols

Simulations of alanine dipeptide in vacuum were performed using the amber99-SB
Hornak  _et al._ (2006) force field. Because this force field is not
implemented in LAMMPS Plimpton (1995), we have used the convert.py python
script that is part of the InterMol Shirts  _et al._ (2017) software to
generate input files implementing an equivalent force field. In all
simulations, long-range interactions between periodic images have not been
included. Path optimizations were performed fixing the Cα carbon of all the
images in the origin, with the carbon atom of the methyl residue aligned along
the z axis and constrained to move only along it, the hydrogen atom of the Cα
carbon free to move only in the $(y,z)$ plane and all other atoms free to move
in all directions. The free energy landscape of alanine dipeptide in vacuum
was obtained using the On-the-fly Probability Enhanced Sampling (OPES) method
Invernizzi and Parrinello (2020), the details of the OPES simulation can be
found in Ref. Invernizzi and Parrinello, 2020.

Simulations of the dehydrogenation of ethane were performed using the C/H/O
ReaxFF Chenoweth, van Duin, and Goddard (2008) force field. Path optimizations
were performed fixing one carbon atom of all the images at the origin, with
the second carbon atom aligned along the $x$ axis and constrained to move only
along it and all other atoms free to move in all direction. The final
configuration adopted in the simulations corresponds to the hydrogen molecule
with its bond parallel to the carbon-carbon bond of ethylene and kept at a
distance of 2.7 Å away from it.

Simulations of graphene were performed using the LCBOPII Los  _et al._ (2005)
force field. The equilibrium carbon-carbon distance in graphene is C-C${}_{\rm
eq}=1.42$ Å. In all simulations, we considered a rectangular supercell of
sides $L_{x}=39.35$ Å and $L_{y}=34.08$ Å, containing 512 atoms, and we
applied periodic boundary conditions in $x$ and $y$. Path optimizations were
performed keeping fixed a selected atom far from the defect. We checked that
relaxations with all the atoms free to move in all directions lead to the same
results.

All simulations were performed using LAMMPS Plimpton (1995). Structural
optimizations used to generate the endpoint configurations were performed
adopting the Fast Inertial Relaxation Engine (FIRE) Bitzek  _et al._ (2006)
minimization algorithm. All standard NEB calculations were performed using the
LAMMPS implementation of the climbing-image NEB algorithm of Ref. Henkelman,
Uberuaga, and Jónsson, 2000 adopting the FIRE Bitzek  _et al._ (2006)
minimization algorithm for structural optimizations. The modified NEB
algorithm has been implemented in LAMMPS, building on our previous
implementation of the OM action Mandelli, Hirshberg, and Parrinello (2020).
Structural optimizations using our modified NEB algorithm were performed
adopting a projected velocity Verlet algorithm Jonsson, Mills, and Jacobsen
(1998). In all cases, optimization was stopped when the maximum force acting
on the atoms of each image was at least smaller than 1 meV/Å. The starting
discretized paths were constructed using a linear interpolation between the
two endpoints.

Figure 1:  (a) Ball-and-stick model of alanine dipeptide showing the two
backbone dihedral angles $\phi$ and $\psi$. Panels (b) and (c) report the four
MEPs investigated obtained using (b) the standard NEB algorithm and (c) our
modified version. As a reference, the free energy surface at $T=300$ K is
shown as a colored map. The results reported here were obtained using
$k_{||}=1$ kcal/mol/Å2, $k_{\perp}=0$ and $\nu=1$ fs-1. In panel (b), path
number 4 appears to be much coarser than the corresponding one in panel (c).
The reason is that the images near the endpoints of this path differ from each
other only in the angular orientation of the methyl groups of the blocked
terminal ends. Since these rotations do not involve the dihedral angles
$(\phi,\psi)$, these images appear to be bunched together at the two local
minima.

###  II.4 Definitions

Within the harmonic approximation to TST, the kinetic rate is given by
Vineyard (1957); Voter and Doll (1984)

$\tau^{-1}=\Omega e^{-\frac{\Delta E^{{\dagger}}}{k_{B}T}}.$ (17)

Where $\Omega=\sum_{n=1}^{D}\nu_{n}/\sum_{n=1}^{D-1}\nu^{\rm saddle}_{n}$, $D$
is the number of degrees of freedom and $\nu_{n}$, $\nu_{n}^{\rm saddle}$ are
the positive normal mode frequencies of the starting configuration and of the
configuration corresponding to the lowest first-order saddle point connecting
the potential energy basin of the initial configuration to the one of the
final state. The activation energy $\Delta E^{\dagger}=E^{\rm saddle}-E^{i}$
is the difference between the potential energy at the saddle point and the
potential energy at the minimum of the starting basin.

Throughout the text, we provide the values of the formation energy $\Delta
E=E^{f}-E^{i}$, where $E^{f}$ is the energy of the final configuration, of the
forward and backward activation energies $\Delta
E^{{\dagger},{\ddagger}}=E^{\rm saddle}-E^{i,f}$, and of the corresponding
frequency prefactors $\Omega^{{\dagger},{\ddagger}}$.

The profiles of the potential energy and of other quantities along the MEP are
plotted as a function of the reaction coordinate. The latter is defined as the
distance from the initial configuration, measured along the path, normalized
with respect to the total path length.

Table 1: The values of $\phi$ and $\psi$ after geometry optimization of alanine dipeptide starting from configurations near the three local minima of the free energy surface shown in Figure 1(b). The last column report the potential energy of the different isomers relative to the most stable $C_{\rm 7eq}$ configuration. Dihedral angles are in degrees and energies in kcal/mol. | $\phi$ | $\psi$ | $\Delta E$
---|---|---|---
$C_{\rm 7eq}$ | -77.5 | 54.1 | 0
$C_{\rm 7eq}^{{}^{\prime}}$ | -147.0 | 159.1 | 0.595
$C_{\rm 7ax}$ | 60.2 | -40.9 | 1.421
Table 2:  The difference in potential energy between the initial and final states ($\Delta E$), the forward ($\Delta E^{{\dagger}}$) and reverse ($\Delta E^{{\ddagger}}$) activation energies, the dihedral angles $(\phi_{\rm TS},\psi_{\rm TS})$ at the transition state and the frequency prefactor ($\Omega$) of the four MEPs investigated. Paths are numbered according to Figure 1(c). Energies are in kcal/mol, dihedrals in degrees, frequencies in THz. Path | $\Delta E$ | $\Delta E^{\dagger}$ | $\Delta E^{\ddagger}$ | $\phi_{\rm TS}$ | $\psi_{\rm TS}$ | $\Omega^{\dagger}$ | $\Omega^{\ddagger}$
---|---|---|---|---|---|---|---
1 | 0.595 | 1.962 | 1.365 | -81.6 | 120.2 | 2.4 | 5.9
2 | 1.421 | 8.694 | 7.272 | -2.1 | -26.4 | 2.6 | 3.6
3 | 0.826 | 14.653 | 13.827 | 121.2 | 176.7 | 4.1 | 1.2
4 | 0.595 | 8.055 | 7.460 | -149.6 | -93.6 | 3.2 | 1.3

##  III  Results

###  III.1  Isomerization of alanine dipeptide

As a first case study, we have considered the isomerization of alanine
dipeptide in vacuum. Figure 1(a) shows a ball-and-stick model of the molecule
along with the definition of the two backbone dihedral angles $(\phi,\psi)$
used to describe its different conformations. Stable configurations have been
identified looking at the free energy landscape shown in Figure 1(b), which is
characterized by the presence of three main basins, labelled $C_{\rm 7eq}$,
$C_{\rm 7eq}^{{}^{\prime}}$ and $C_{\rm 7ax}$. To generate the endpoint
configurations for the construction of the MEPs, we have used finite
temperature molecular dynamics simulations to sample configurations in each of
the three metastable basins. Starting from these configurations, we have
performed geometry optimizations and obtained fully relaxed structures of the
three isomers. Table 1 reports the values of the backbone dihedral angles of
the optimized configurations and the corresponding potential energy measured
relative to the most stable $C_{\rm 7eq}$ isomer. We have subsequently
considered four paths connecting the three basins and crossing different
energy barriers. For each of them, we have obtained the MEP using the standard
NEB algorithm and our modified version.

Figures 1(b) and (c) report the results obtained using $N=20$ images per path,
showing that our modified version systematically improves the resolution
around the saddle points. This is particularly evident for paths number 2, 3
and 4. These paths cross higher and sharper energy barriers (see also Table
2). This is the typical situation where the standard NEB algorithms leads to a
poor sampling of the transition region.

For all paths considered, both algorithms relax towards the same MEP and
locate the same saddle point. This is made clearer in Figure 2(a), where we
report the potential energy profile along the MEP corresponding to path number
2. In the region crossing the sharp energy barrier, roughly corresponding to
values of the reaction coordinate between 0.6 and 0.8, the density of images
is nearly doubled using our algorithm. In Figures S1-S3 of the Supplementary
Material we report the results for paths number 1, 3 and 4.

Figure 2:  The potential energy along path number 2, connecting the $C_{\rm
7eq}$ and the $C_{\rm 7ax}$ isomers (see Figure 1(b)), obtained using (a)
$N=20$ and (b) $N=4$ images. The insets are zoom-in of the paths in the
$(\phi,\psi)$ plane, plotted over the free energy surface. The black and the
red curves correspond to the results obtained using the standard and our
modified NEB algorithm, respectively. The green curves labeled
OM${}^{{}^{\prime}}$ are results obtained using the exact gradient of equation
(11) instead of the simplified expression of equation (13) (see discussion in
Section II.2). In panel (a), the difference between the curves at the
beginning and at the end of the path is due to rotations of the methyl groups
of the blocked terminal ends that are absent in the MEP obtained using the
modified NEB. We note that these rotations are present at the beginning of the
path when using the exact gradient of equation (11) (see green curve),
however, this does not affect the resolution around the saddle point. The
parameters used are the same as reported in the caption of Figure 1. Figure 3:
The potential energy along path number 2, obtained using the modified NEB
algorithm using different values of the parameter $\nu$, as indicated in the
legend (units are fs-1).

In order to check the dependence of the results on the number of images, we
performed additional simulations varying $N$ in the range $3\leq N\leq 20$. In
all cases, both algorithms identified the same saddle point, independent of
the number of images adopted, while our modified scheme always lead to an
improved resolution around the saddle (see Figures S4-S7 in the Supplementary
Material). As an example, Figure 2(b) shows the results obtained using $N=4$
images.

Finally, in Figure 3 we study the dependence of the results on the value of
the parameter $\nu$, defining the mass-dependent spring constant $k_{\rm
OM}=m\nu/2\Delta t$ in our modified scheme. Results show that the MEP is
independent of the value of $\nu$ across several orders of magnitude. Similar
results were obtained for all the paths considered (see Figures S8-S10 in the
Supplementary Materials).

We further checked that the results are independent of the initial condition
by performing standard NEB optimizations starting from the path optimized
using our modified algorithm. In these tests, we never observed an increase in
resolution near the saddle point and the path always relaxed to the same MEP
obtained starting from the linear interpolation.

Figure 4:  (a) The potential energy along the MEP for the dehydrogenation of
ethane obtained using $N=15$ images. The black and the red curves correspond
to the results obtained using the standard and the modified NEB algorithm,
respectively. The insets show the ethane molecule in the initial configuration
and the structure at the transition state. The hydrogen atoms that detach from
ethane during the reaction are colored in red. (b) Potential energy along the
MEP obtained using $N=5$ images. The results in both panels were obtained
using $k_{||}=1$ kcal/mol/Å2, $k_{\perp}=0$ and $\nu=1$ fs-1. Figure 5:  (a)
Top view of perfect graphene. The bond that undergoes a ninety degrees
rotation is marked in red and the two bonds that break during the
transformation are also indicated. (b) Top view of an optimized, buckled SW
defect. The rotated bond at the core of the defect is marked in red and the
two new bonds that form during the transformation are indicated. (c) Lateral
view of the sinelike buckled SW defect. (d) Lateral view of the cosinelike
buckled SW defect. In panels (c) and (d), carbon atoms are colored according
to their vertical position relative to the basal plane of graphene. The
supercells are shown with the bond that undergoes the rotation in the
foreground and in central position.

###  III.2  Dehydrogenation of ethane

As a second example, we have studied the dehydrogenation of ethane,
$C_{2}H_{6}\rightleftharpoons C_{2}H_{4}+H_{2}$. This reaction is an
endothermic elimination reaction. With the adopted force field, the difference
in energy between reactants and products is $\Delta E=43$ kcal/mol. During the
reaction, two C-H bonds break and the carbon hybridization changes from
$sp_{3}$ to $sp_{2}$. The dehydrogenation of light alkanes is of interest in
the chemical industry for the production of chemical intermediates based on
alkenes. In industrial applications, the $C-H$ bond is activated using
catalyzers. Here we consider the reaction in absence of catalyzers. In this
picture, the first step of the reaction consists in a rotation of the methyl
groups that brings two hydrogen atoms in alignment. Subsequently, the hydrogen
atoms approach, the $C-H$ bonds break, the $H_{2}$ molecule leaves and
ethylene reaches its equilibrium planar structure. In our simulations, after
the initial rotation of the methyl groups, the system maintains a symmetric
configuration with the carbon-carbon and hydrogen-hydrogen bonds parallel to
each other and centered. The carbon-carbon distance changes from 1.57 to 1.53
to 1.33 Å, respectively in ethane, at the transition state and in ethylene.
The forward and reverse activation barriers are $\Delta E^{\dagger}=76.8$
kcal/mol and $\Delta E^{\ddagger}=33.8$ kcal/mol, respectively. The frequency
prefactor in equation (17) for the dehydrogenation is $\Omega^{\dagger}=639$
THz.

In Figure 4(a) we report the MEPs obtained using the standard and modified NEB
algorithm using $N=15$ images. Our method roughly doubles the density of
images in the transition region of the path corresponding to the reaction
coordinate in the interval $[0.5,0.7]$ around the saddle point. The optimized
paths as well as the saddle points obtained from the two protocols coincide.
We checked that the results of our modified algorithm are the same in the
parameter range of $0.1<\nu<100$ fs-1 (see Figure S11 in the Supplementary
Materials).

In order to check the dependence of the results on the number of images, we
performed additional simulations varying $N$ in the range $5\leq N\leq 15$. In
all cases, both algorithms identified the same saddle point, independent of
the number of images adopted. Also in this case, our modified scheme
systematically leads to an improved resolution around the saddle. As an
example, in Figure 2(b) we report the results obtained using $N=5$ images. In
Figure S12 of the Supplementary Material we report the results for all the
values of $N$ investigated.

###  III.3  Healing of a 5-77-5 defect in graphene

An important defect in graphite is the so-called 5-77-5 topological defect,
which is formed via a rotation by $\pi/2$ of a carbon-carbon bond in the
graphene sheet (see Figure 5(a) and (b)). This rotation is also called Stone-
Wales (SW) transformation Stone and Wales (1986) and it involves the breaking
of two pristine carbon-carbon bonds (1 and 2 in panel (a)) and the subsequent
formation of two new bonds (1’ and 2’ in panel (b)). In the process, four
hexagons of the original graphene honeycomb are converted into two pentagons
and two heptagons. The SW defect plays an important role in the formation of
fullerenes and nanotubes. In fact, the planar configuration of the SW defect
is not stable and the flexible graphene sheet can easily reduce the
compressive strain at the core of the defect by acquiring a buckled
configuration. Two metastable geometries have been identified in the
literature Ma  _et al._ (2009), corresponding to the sinelike and cosinelike
structures. These are shown in Figure 5(c) and (d), as obtained after geometry
optimization using our model graphene supercell. The sinelike structure is
characterized by a larger buckling amplitude of $h=2.26$ Å as compared to the
value of $h=1.68$ Å in the cosinelike geometry. On the other hand, the length
of the rotating carbon-carbon bond is the same in both structures, being
$\sim$4 % compressed with respect to the equilibrium bond length of the
pristine sheet (see Table 3). In agreement with previous results Ma  _et al._
(2009), we found that the sinelike configuration is energetically favored. The
formation energies of the sinelike and cosinelike defects are $\Delta E=4.33$
and 4.53 eV, respectively, in fair agreement with previous density functional
theory calculations of the buckled structures Ma  _et al._ (2009) and with the
values reported in Ref. Los  _et al._ , 2005 for the unstable planar geometry
using the same force field adopted here. In the following, we focus on the
sinelike defect.

Table 3:  The values of the buckling amplitude ($h$) and of the carbon-carbon bond length of the bond that undergoes the rotations in the optimized buckled sinelike and cosinelike SW defects shown in Figure 5. Amplitudes are in Å and bond lengths are in units of the equilibrium value of the adopted force field, $C-C_{eq}=1.42$ Å. | $h$ | $C-C$
---|---|---
Sinelike | 2.26 | 0.964
Cosinelike | 1.68 | 0.965
Table 4:  The formation energy ($\Delta E$), the forward ($\Delta E^{\dagger}$) and reverse ($\Delta E^{\ddagger}$) activation energies, the angle ($\theta_{\rm TS}$ ), carbon-carbon bond length ($C-C_{\rm TS}$), and buckling amplitude ($h_{\rm TS}$) at the transition state and the frequency prefactor ($\Omega$) of the MEPs obtained for the sinelike and cosinelike SW defect. Energies are in eV, amplitudes in Å, angles in degrees, bond lengths are in units of the equilibrium value ($C-C_{\rm eq}=1.42$ Å) and frequencies are in THz. | $\Delta E$ | $\Delta E^{\dagger}$ | $\Delta E^{\ddagger}$ | $h_{\rm TS}$ | $\theta_{\rm TS}$ | $C-C_{\rm TS}$ | $\Omega^{\dagger}$ | $\Omega^{\ddagger}$
---|---|---|---|---|---|---|---|---
Sinelike | 4.43 | 4.64 | 8.97 | 1.68 | 46.41 | 0.871 | 6.5 | 0.05
Cosinelike | 4.56 | 4.40 | 8.96 | 1.17 | 39.04 | 0.872 | 5.2 | 0.07

Figure 6:  The potential energy along the MEP for the healing of the sinelike
SW defect obtained using (a) $N=20$ and (b) $N=5$ images. The black and the
red curves correspond to the results obtained using the standard and the
modified NEB algorithm, respectively. The inset in panel (a) is a zoom-in of
the region around the core of the defect at the transition state. The bond
that undergoes the ninety degrees clockwise rotation is marked in red. The
results reported here were obtained using $k_{||}=1$ eV/Å2, $k_{\perp}=0$ and
$\nu=1$ fs-1. Figure 7:  (a) The angle and (b) the carbon-carbon bond length
of the bond that undergoes the ninety degrees rotation and (c) the buckling
amplitude along the MEP for the healing of the sinelike SW defect. The black
and the red curves correspond to the results obtained using the standard and
the modified NEB algorithm, respectively, using $N=20$ images.

In Figure 6(a) we report the MEPs for the healing of the SW defect, obtained
using $N=20$ images. Also in this case, the modified algorithm leads to an
improved resolution of the MEP around the saddle point, roughly doubling the
density of images in the region straddling the sharp energy peak. This has
allowed us to investigate in details the structure of the core of the defect
along the path. In Figure 7 we report the profile of three geometrical
parameters, namely, the angle $\theta$ and the length of the rotating bond and
the overall buckling amplitude $h$. In the transition region, the angle
increases linearly (see panel (a)), while the bond length shrinks, reaching a
minimum value at the saddle point, and then increases towards the equilibrium
value (see panel (b)). The saddle point configuration corresponds to an almost
linear geometry of the four carbons at the core of the defect, suggesting an
sp-like hybridization of the two rotating carbon atoms and the formation of a
triple bond (see inset in Figure 6(a) and Table 4). Accordingly, the value of
1.24 Å measured at the saddle point coincides with the length of the triple
bond described by the adopted force field (see Table 2 in Ref. Los  _et al._ ,
2005). Along the MEP, the buckling amplitude decreases monotonously and is
characterized by an almost flat plateau around the saddle (see panel (c)). For
completeness, in Table 4 we report the energetics of the reaction and the
geometrical parameters characterizing the saddle points that one needs to
cross to heal the sinelike and cosinelike SW defects. We checked that results
obtained with our modified algorithm are the same in the range $0.1<\nu<10$
fs-1 of the parameter $\nu$ defining the spring constant $k_{\rm
OM}=m\nu/2\Delta t$ (see Figure S13 in the Supplementary Materials).

Simulations carried out using different number $N$ of images always resulted
in the identification of the same saddle point, independent of $N$ and of the
algorithm adopted. Also in this case, our modified scheme always lead to an
improved resolution around the saddle. This is exemplified here in Figure
6(b), showing results for $N=5$. In Figures S14 and S15 of the Supplementary
Material we report the analysis of the results obtained for all the values of
$N$ investigated.

##  IV  Discussion and Conclusions

In this work, we have addressed the problem of increasing the resolution
around the saddle point of discretized minimum energy paths obtained within
the framework of a nudged-elastic-band-like algorithm. Other approaches have
already been presented in the literature. Notably, in a first modification to
the original NEB algorithm Henkelman, Uberuaga, and Jónsson (2000), it was
suggested to automatically increase the spring constant in proximity of the
saddle point. Subsequently, other variants of the NEB algorithm have been
proposed, where a sequence of standard NEB calculations is performed that
iteratively increase the resolution around the image with higher energy by
locally adding new replicas Maragakis  _et al._ (2002); Kolsbjerg, Groves, and
Hammer (2016). These methods have been designed to improve computational
efficiency by relaxing only few replicas around the putative saddle point at
each iteration. Our method differs from these approaches as it modifies
directly the expression of the forces that are minimized, while the algorithm
remains the same as in the standard NEB approach Henkelman and Jónsson (2000);
Henkelman, Uberuaga, and Jónsson (2000) with fixed number of images. Unlike
Ref. Henkelman, Uberuaga, and Jónsson, 2000, rather than adjusting the spring
constant, our method exploits the stationary property of the saddle point to
define a locally adaptive natural spring length that vanishes at the saddle.

In conclusion, we have presented a modified NEB method that leads to an
automatic increase in resolution around saddle points of the potential energy
surface. We demonstrated the ability of the method in cases of practical
interest, selected out of the realm of physical chemistry and materials
science and compared it to the NEB algorithm Henkelman and Jónsson (2000);
Henkelman, Uberuaga, and Jónsson (2000). In all the cases studied, our
approach improved the resolution of the MEP around the saddle point. In light
of the results presented here, we believe that the method represents a useful
alternative tool for locating saddle points in processes with known starting
and final configurations.

Finally, we note that our approach can be combined straightforwardly with the
algorithms of Refs. Maragakis  _et al._ , 2002; Kolsbjerg, Groves, and Hammer,
2016; Zhang  _et al._ , 2016 for improved efficiency and to have more control
on the number of images as they are manually added.

## Supplementary Material

The online supplementary material includes Figures S1-S15, reporting the
results of additional simulations.

## Acknowledgments

The research was supported by the European Union Grant No.
ERC-2014-AdG-670227/VARMET. We also thank the NCCR MARVEL, funded by the Swiss
National Science Foundation.

## Data availability statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Dellago  _et al._ (1998) C. Dellago, P. G. Bolhuis, F. S. Csajka, and D. Chandler, “Transition path sampling and the calculation of rate constants,” The Journal of Chemical Physics 108, 1964–1977 (1998).
  * van Erp, Moroni, and Bolhuis (2003) T. S. van Erp, D. Moroni, and P. G. Bolhuis, “A novel path sampling method for the calculation of rate constants,” The Journal of Chemical Physics 118, 7762–7774 (2003).
  * Faradjian and Elber (2004) A. K. Faradjian and R. Elber, “Computing time scales from reaction coordinates by milestoning,” The Journal of Chemical Physics 120, 10880–10889 (2004).
  * Tiwary and Parrinello (2013) P. Tiwary and M. Parrinello, “From Metadynamics to Dynamics,” Physical Review Letters 111, 230602 (2013).
  * Debnath and Parrinello (2020) J. Debnath and M. Parrinello, “Gaussian Mixture-Based Enhanced Sampling for Statics and Dynamics,” The Journal of Physical Chemistry Letters 11, 5076–5080 (2020).
  * Mandelli, Hirshberg, and Parrinello (2020) D. Mandelli, B. Hirshberg, and M. Parrinello, “Metadynamics of Paths,” Physical Review Letters 125, 026001 (2020).
  * Eyring (1935) H. Eyring, “The Activated Complex in Chemical Reactions,” The Journal of Chemical Physics 3, 107–115 (1935).
  * Wigner (1938) E. Wigner, “The transition state method,” Transactions of the Faraday Society 34, 29 (1938).
  * Keck (2007) J. C. Keck, “Variational Theory of Reaction Rates,” (WILEY-VCH Verlag, 2007) pp. 85–121.
  * Henkelman and Jónsson (1999) G. Henkelman and H. Jónsson, “A dimer method for finding saddle points on high dimensional potential surfaces using only first derivatives,” The Journal of Chemical Physics 111, 7010–7022 (1999).
  * Jay  _et al._ (2020) A. Jay, C. Huet, N. Salles, M. Gunde, L. Martin-Samos, N. Richard, G. Landa, V. Goiffon, S. De Gironcoli, A. Hémeryck, and N. Mousseau, “Finding Reaction Pathways and Transition States: r-ARTn and d-ARTn as an Efficient and Versatile Alternative to String Approaches,” Journal of Chemical Theory and Computation 16, 6726–6734 (2020).
  * Jonsson, Mills, and Jacobsen (1998) H. Jonsson, G. Mills, and K. W. Jacobsen, “Nudged elastic band method for finding minimum energy paths of transitions,” in  _Classical and Quantum Dynamics in Condensed Phase Simulations_ (World Scientific, 1998) pp. 385–404.
  * Peters  _et al._ (2004) B. Peters, A. Heyden, A. T. Bell, and A. Chakraborty, “A growing string method for determining transition states: Comparison to the nudged elastic band and string methods,” Journal of Chemical Physics 120, 7877–7886 (2004).
  * E and Vanden-Eijnden (2010) W. E and E. Vanden-Eijnden, “Transition-Path Theory and Path-Finding Algorithms for the Study of Rare Events,” Annual Review of Physical Chemistry 61, 391–420 (2010).
  * Henkelman and Jónsson (2000) G. Henkelman and H. Jónsson, “Improved tangent estimate in the nudged elastic band method for finding minimum energy paths and saddle points,” The Journal of Chemical Physics 113, 9978–9985 (2000).
  * Henkelman, Uberuaga, and Jónsson (2000) G. Henkelman, B. P. Uberuaga, and H. Jónsson, “A climbing image nudged elastic band method for finding saddle points and minimum energy paths,” The Journal of Chemical Physics 113, 9901–9904 (2000).
  * Maragakis  _et al._ (2002) P. Maragakis, S. A. Andreev, Y. Brumer, D. R. Reichman, and E. Kaxiras, “Adaptive nudged elastic band approach for transition state calculation,” Journal of Chemical Physics 117, 4651–4658 (2002).
  * Kolsbjerg, Groves, and Hammer (2016) E. L. Kolsbjerg, M. N. Groves, and B. Hammer, “An automated nudged elastic band method,” Journal of Chemical Physics 145, 094107 (2016).
  * Onsager and Machlup (1953) L. Onsager and S. Machlup, “Fluctuations and Irreversible Processes,” Physical Review 91, 1505–1512 (1953).
  * Maras  _et al._ (2016) E. Maras, O. Trushin, A. Stukowski, T. Ala-Nissila, and H. Jónsson, “Global transition path search for dislocation formation in Ge on Si(001),” Computer Physiscs Communications 205, 13–21 (2016).
  * Ásgeirsson  _et al._ (2021) V. Ásgeirsson, B. O. Birgisson, R. Bjornsson, U. Becker, F. Neese, C. Riplinger, and H. Jónsson, “Nudged Elastic Band Method for Molecular Reactions Using Energy-Weighted Springs Combined with Eigenvector Following,” Journal of Chemical Theory and Computation 17, 4929–4945 (2021).
  * Fujisaki, Shiga, and Kidera (2010) H. Fujisaki, M. Shiga, and A. Kidera, “Onsager–Machlup action-based path sampling and its combination with replica exchange for diffusive and multiple pathways,” The Journal of Chemical Physics 132, 134101 (2010).
  * Fujisaki  _et al._ (2013) H. Fujisaki, M. Shiga, K. Moritsugu, and A. Kidera, “Multiscale enhanced path sampling based on the Onsager-Machlup action: Application to a model polymer,” The Journal of Chemical Physics 139, 054117 (2013).
  * Lee  _et al._ (2017) J. Lee, I.-H. Lee, I. Joung, J. Lee, and B. R. Brooks, “Finding multiple reaction pathways via global optimization of action,” Nature Communications 8, 15443 (2017).
  * Hornak  _et al._ (2006) V. Hornak, R. Abel, A. Okur, B. Strockbine, A. Roitberg, and C. Simmerling, “Comparison of multiple Amber force fields and development of improved protein backbone parameters,” Proteins: Structure, Function, and Bioinformatics 65, 712–725 (2006).
  * Plimpton (1995) S. Plimpton, “Fast Parallel Algorithms for Short-Range Molecular Dynamics,” Journal of Computational Physics 117, 1–19 (1995).
  * Shirts  _et al._ (2017) M. R. Shirts, C. Klein, J. M. Swails, J. Yin, M. K. Gilson, D. L. Mobley, D. A. Case, and E. D. Zhong, “Lessons learned from comparing molecular dynamics engines on the SAMPL5 dataset,” Journal of Computer-Aided Molecular Design 31, 147–161 (2017).
  * Invernizzi and Parrinello (2020) M. Invernizzi and M. Parrinello, “Rethinking Metadynamics: From Bias Potentials to Probability Distributions,” The Journal of Physical Chemistry Letters 11, 2731–2736 (2020).
  * Chenoweth, van Duin, and Goddard (2008) K. Chenoweth, A. C. T. van Duin, and W. A. Goddard, “ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation,” The Journal of Physical Chemistry A 112, 1040–1053 (2008).
  * Los  _et al._ (2005) J. H. Los, L. M. Ghiringhelli, E. J. Meijer, and A. Fasolino, “Improved long-range reactive bond-order potential for carbon. I. Construction,” Physical Review B 72, 214102 (2005).
  * Bitzek  _et al._ (2006) E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, and P. Gumbsch, “Structural Relaxation Made Simple,” Physical Review Letters 97, 170201 (2006).
  * Vineyard (1957) G. H. Vineyard, “Frequency factors and isotope effects in solid state rate processes,” Journal of Physics and Chemistry of Solids 3, 121–127 (1957).
  * Voter and Doll (1984) A. F. Voter and J. D. Doll, “Transition state theory description of surface self-diffusion: Comparison with classical trajectory results,” The Journal of Chemical Physics 80, 5832–5838 (1984).
  * Stone and Wales (1986) A. J. Stone and D. J. Wales, “Theoretical studies of icosahedral C60 and some related species,” Chemical Physics Letters 128, 501–503 (1986).
  * Ma  _et al._ (2009) J. Ma, D. Alfè, A. Michaelides, and E. Wang, “Stone-Wales defects in graphene and other planar sp2-bonded materials,” Physical Review B 80, 033407 (2009).
  * Zhang  _et al._ (2016) J. Zhang, H. Zhang, H. Ye, and Y. Zheng, “Free-end adaptive nudged elastic band method for locating transition states in minimum energy path calculation,” The Journal of Chemical Physics 145, 094104 (2016).