next up previous contents
Next: Optimisation of the Nudged Up: FINDING REARRANGEMENT PATHWAYS Previous: Introduction   Contents


A Double-ended Method: Nudged Elastic Band

In the present work we used the nudged elastic band [27,28] (NEB) and eigenvector-following [17,18,12,14,24,13,19,16,22,23,15] (EF) methods for locating and refining transition states. In the NEB approach the path is represented as a set of images $ \left\{\bfx _1,\bfx _2 ... \bfx _{N_i}\right\}$ that connect the endpoints $ \bfx _0$ and $ \bfx _{N_i+1}$,where $ \bfx _i$ is a vector containing the coordinates of image $ i$ (Figure 2.1) [29]. In the usual framework of double-ended methodologies [130] the endpoints are stationary points on the PES (usually minima), which are known in advance. In addition to the true potential, $ V_i$, which binds the atoms within each image, equivalent atoms in $ N_i$  adjacent images are interconnected by $ N_i+1$ springs according to a parabolic potential,

$\displaystyle \widetilde{V} = \tfrac{1}{2} k_{spr} \sum_{i=1}^{N_i+1} \vert\bfx _i - \bfx _{i-1} \vert^2.$ (4.1)

Subsequently these potentials will be referred to as the `true potential' and the `spring potential', respectively.

Figure: Graphical representation of the nudged elastic band approach. (a) The optimised nudged elastic band for a two-dimensional model surface. The band contains $ 21$ images and connects two minima $ \bfx _0$ and $ \bfx _{23}$. Image $ \bfx _9$ has the highest energy and might therefore be used to estimate transition state properties or as a starting guess for further refinement. (b) `Nudging': the NEB depicted in (a) is projected onto the $ xy$ plane and feels only the perpendicular component of the true gradient from the effective potential $ V^{\perp }$.
\begin{psfrags} \psfrag{z} [Bl][Bl]{$V\left( \bfx \right)$} \psfrag{x} [Bl][... ...terline{\includegraphics[width=.45\textheight]{dneb/surface.eps}} \end{psfrags}

The springs are intended to hold images on the path during optimisation -- otherwise they would slide down to the endpoints, or to other intermediate minima [88]. Occasionally, depending on the quality of the initial guess, we have found that some images may converge to higher index stationary points. One could imagine the whole construction as a band or rope that is stretched across the PES, which, if optimised, is capable of closely following a curve defined in terms of successive minima, transition states, and the intervening steepest-descent paths.

In practice, the above formulation encounters difficulties connected with the coupling between the `true' and `spring' components of the potential. The magnitude of the springs' interference with the true potential is system dependent and generally gives rise to corner-cutting and sliding-down problems [27]. It is convenient to discuss these difficulties in terms of the components of the true gradient, $ \gt $, and spring gradient, $ \gs $, parallel and perpendicular to the path. The parallel component of the gradient $ \gtpar $  at image $ i$ on the path is obtained by projecting out the perpendicular component $ \gtper $  using an estimate of the tangent to the path. The parallel and perpendicular components for image $ i$ are:

$\displaystyle \gtpar _i = \left( \bflp _i V_i,\bft _i \right) \bft _i, \hskip1cm \gtper _i = \bflp _i V_i - \gtpar _i,$ (4.2)

where $ V_i = V\left( \bfx _i \right)$, the unit vector $ \bft _i$ is the tangent, and $ \left( \bflp _i V_i,\bft _i \right)$ denotes the scalar product of vectors $ \bflp _i V_i$ and $ \bft _i$. Here and throughout this work we denote unit vectors by a hat.The complete gradient, $ \gt $, has $ N_i \times \eta$ components for a band of $ N_i$ images with $ \eta$  atomic degrees of freedom each.

Corner-cutting has a significant effect when a path experiences high curvature. Here $ \gsper $  is large, which prevents the images from following the path closely because the spring force necessarily has a significant component perpendicular to the tangent. The sliding-down problem occurs due to the presence of $ \gtpar $, which perturbs the distribution of images along the path, creating high-resolution regions (around the local minima) and low-resolution regions (near the transition states) [27]. Both problems significantly affect the ability of the NEB method to produce good transition state candidates. We have found that sliding-down and corner-cutting are interdependent and cannot both be remedied by adjusting the spring force constant $ k_{spr}$;increasing $ k_{spr}$ may prevent sliding-down but it will make corner-cutting worse.

The aforementioned problems can sometimes be eliminated by constructing the NEB gradient from the potential in the following way: $ \gtpar $ and $ \gsper $ are projected out, which gives the elastic band its `nudged' property [28]. Removal of $ \gtpar $ can be thought of as bringing the path into a plane or flattening the PES [Figure 2.1 (b)], while removal of $ \gsper $ is analogous to making the images heavier so that they favour the bottom of the valley at all times.

The choice of a method to estimate the tangent to the path is important for it affects the convergence of the NEB calculation. Originally, the tangent vector, $ \bft _i$, for image $ i$ was obtained by normalising the line segment between the two adjacent images, $ i+1$ and $ i-1$ [27]:

$\displaystyle \bft _i = \frac{ \bfx _{i+1}-\bfx _{i-1} }{\vert \bfx _{i+1}-\bfx _{i-1} \vert}.$ (4.3)

However, kinks can develop during optimisation of the image chain using this definition of $ \bft _i$. It has been shown [29] that kinks are likely to appear in the regions where the ratio $ g_i^{\parallel}/g_j^{\perp}$ is larger than the length of the line segment, $ \vert$$ \tau$$ \vert$, used in estimating the tangent [Figure 2.2 (a)].

Figure: Details of recent NEB implementations. (a) Conditions under which kinks appear during optimisation of the NEB using the tangent estimated from the line segment $ \bft $ connecting images $ i+1$ and $ i-1$. Displacement of image $ i$ from the path (dash-dotted line) creates forces $ {\bf F}^{\perp}_{i-1} = - \gtper _{i-1}$ and $ {\bf F}^{\perp}_{i} = - \gtper _i$. While $ {\bf F}^{\perp }_{i}$ is a restoring force that originates from $ V^{\perp }$, $ {\bf F}^{\perp }_{i-1}$ is destabilising and originates from $ V^{\parallel }$ (and is non-zero due to the fact that the tangent at image $ i-1$ has changed after displacement of image $ i$). For the case of small displacements the potential may be resolved into two contributions, $ V^{\perp }=k^{\perp }x^2/2$ and $ V^{\parallel }=-k^{\parallel }y$, and kinks will not appear if $ k^{\parallel}/k^{\perp} < \vert \bft \vert$. (b) Tangent estimate using the higher energy neighbour: image $ i+1$ is `hanging' on to image $ i$. The separation $ d$ is controlled by the lower-lying images $ (> i+1)$ but not $ V$. (c) An NEB that follows the curved region of the path: since the spring force $ {\bf F}_1$ acting on image $ i$ is compensated by projection $ {\bf F}_2$, the distribution of images becomes uneven. (d) Corner-cutting displayed on a cross-section of the curved part of the path depicted in (c): the image is displaced from the path due to the presence of $ {\bf F}_{spr}^{\perp }$.
\begin{psfrags} \psfrag{tau} [Bl][Bl]{$\bft $} \psfrag{tab} [Bl][Bl]{$\vert ... ...centerline{\includegraphics[width=.55\textheight]{dneb/tech.eps}} \end{psfrags}

Both the above ratio, the image density and $ \vert\bft \vert$ can vary depending on the system of interest, the particular pathway and other parameters of the NEB calculation. From Equation 2.3 it can be seen that $ \bft _i$, and, hence, the next step in the optimisation of image $ i$, is determined by its neighbours, which are not necessarily closer to the path than image $ i$. Therefore, a better approach in estimating the $ \bft _i$ would be to use only one neighbour, since then we only need this neighbour to be better converged than image $ i$.

There are two neighbours to select from, and it is natural to use the higher-energy one for this purpose, since steepest-descent paths are easier to follow downhill than uphill:

$\displaystyle \bft _i = \frac{ \left( j-i \right) \left( \bfx _{j}-\bfx _{i} \right) }{\vert \bfx _{j}-\bfx _{i} \vert},$ (4.4)

where $ i$ and $ j$ are two adjacent images with energies $ V_i$ and $ V_j$, and $ V_i < V_j$. In this way, an image $ i$ that has one higher-energy neighbour $ j$ behaves as if it is `hanging' on to it [Figure 2.2(b)].

The above tangent formulation requires special handling of extrema along the path, and a mechanism for switching $ \bft $ at such points was proposed [29]. It also fails to produce an even distribution of images in regions with high curvature [Figure 2.2 (c)]. We presume that Henkelman and Jónsson substitute $ \left( \gs , \bft \right) \bft $ by $ \vert \gs \vert \bft $ in Equation 2.2 to obtain a spring gradient formulation that will keep the images equispaced when the tangent from Equation 2.4 is used in the projections [28]:

$\displaystyle \gspar _i = k_{spr} \Bigl( \vert \bfx _i - \bfx _{i-1} \vert - \vert \bfx _{i+1} - \bfx _i \vert \Bigr) \bft _i.$ (4.5)

next up previous contents
Next: Optimisation of the Nudged Up: FINDING REARRANGEMENT PATHWAYS Previous: Introduction   Contents
Semen A Trygubenko 2006-04-10