3.11. Transition States

---

download:pdf

---

The overall objective of the methods encompassed by the Transition State Search module is the calculation of the rate of transitions such as chemical reactions, surface processes or diffusion in solids. These transitions are many orders of magnitude slower than vibrations of atoms and are therefore rather rare events. Consequently, a direct simulation of the dynamics of these processes is not feasible within an accessible timeframe.

Whereas the time scale problem rules out direct dynamical simulation it opens up the use of statistical methods for estimating rates of transitions, i.e. transition state theory [1]. Two basic assumptions have to be made in transition state theory:

The rate is slow enough that Boltzmann statistics are applicable for the initial state, and

an energy divide (a n-1 dimensional energy surface, where n is the number of degrees of freedom) can be identified such that the reaction trajectory crosses only once.

The transition state (TS) energy divide represents the bottleneck of the transition. The rate constant \(\Gamma\) is then given by

(1)\[\Gamma = P'\overline{v}\]

(2)\[P' = \dfrac{P_{TS}}{P_{0}} = \dfrac{Q_{TS}}{Q_{0}} = \dfrac{e^{\Delta E}}{k_{B}T}\]

\(\overline{v}\) mean velocity of crossing the transition state,

\(P\) probability ratio of being at the transition state \(P_{TS}\) or initial state (\(P_{0}\)),

\(Q\) partition functions of the system in transition state (\(Q_{TS}\)) and initial state (\(Q_{0}\)),

\(\Delta E\) energy difference between the transition state and the initial state

Transition state theory gives a rate of escape from the initial state irrespective of the final state and is known for always overestimating the reactivity of the system. Final states can be found by starting from the energy divide and sliding down the barrier using dynamical simulations.

In crystals, the harmonic transition state theory is a good approximation simplifying the rate estimation for e.g. diffusion processes and reaction in solids [2]. The approach requires the knowledge of lowest lying saddle points on a pathway between initial and final state, the so-called minimum energy path (MEP). The rate constant for transitions through each saddle point can be derived from the energy and the frequency of the normal modes at the corresponding initial state and the transition state. The transition state is characterized by one unstable normal mode with an imaginary frequency pointing in the direction of the transition. Whereas the Transition State Search module enables one to find the minimum energy path with all its saddle points and intermediate local minima, the frequencies of normal modes and thermodynamic functions at minima and saddle points can be obtained from MedeA Phonon. Combining these two computational tools provides access to kinetic information, such as rate constants of reactions and temperature-dependent diffusion coefficients [3].

In the following, methods for searching minimum energy paths and transition states are discussed in more detail.

3.11.1. Drag Method

Although not implemented in MedeA TSS, the most simple and intuitive approach for mapping minimum energy paths is briefly discussed. One degree of freedom (the drag coordinate) is selected and kept fixed whereas all other n-1 degrees of freedom are relaxed in each step of the simulation, which corresponds to an optimization of the system in a n-1 dimensional hyperplane. In small steps, the drag coordinate is increased, incrementally shifting the system from the initial towards the final state. The method relies on an ad hoc but good guess for a drag coordinate, which is one of the drawbacks of this approach. Furthermore, if the minimum energy path is curved, as is shown as dotted line in the Figure to the left, the system may snap (along the long dashed line) from the valley of slowest ascent (solid line) directly into the product valley (solid line) without ever approaching structurally the saddle point [4]. The problem is caused by the fact that the drag coordinate (short dashed line) is at a quite large angle to the direction of the unstable normal mode at the saddle point.

If the search succeeds the drag method tends to be quite efficient in terms of number of force evaluations, however, there are many cases for which the method fails to find the saddle point.

3.11.2. Elastic Band Method

In the elastic band method, a string of images \(R_{i}\) are connected by springs to form a path from reactants \(R_{0}\) to products \(R_{N}\), where \(R_{i}\) is a vector of all positions of atoms in image \(i\). The springs connect each atom of one image \(i\) to the corresponding atom in neighboring images \(i-1\) and \(i+1\). The resulting chains of springs between atoms of the initial and final image keep the intermediate images on track. The strength of the spring connections is given by a spring constant \(k\).

The initial path of images is created as close as possible to the minimum energy path. A linear interpolation between initial and final configuration is usually applied as the initial path. However, there are cases where linear interpolation fails such as strongly non-linear pathways or situations where linear interpolation moves atoms too close to each other. For these processes, the string of initial images may require a guideline about the expected pathway.

The string of images \(\left[ R_{0},R_{1},R_{2},...,R_{N} \right]\) has fixed endpoints \(R_{0}\) and \(R_{N}\), whereas the \(N-1\) intermediate images are adjusted by an optimization procedure. In the elastic band method an object function can be constructed as

(3)\[S \left( R_{1},R_{2},...,R_{N} \right) = \sum_{i=1}^{N-1} E (R_{i}) + \sum_{i=1}^{N} \dfrac{k}{2} \left( R_{i} - R_{i-1} \right)^2\]

and the set of intermediate images \(\left[ R_{0},R_{1},R_{2},...,R_{N} \right]\) can be relaxed to the minimum energy path by minimization of the object function. In practice, however, the TSS module does not make use of the object function for minimization purposes. The string of images is rather optimized by the minimization of all forces on all atoms including spring force components. Nevertheless, due to the mere existence of an object function the convergence behavior is usually much better than for the nudged elastic band approach.

3.11.2.1. Problems of the Elastic Band Method

Despite these advantages, two main issues are observed with the elastic band method, which can be overcome by the so-called nudging: corner cutting and down sliding.

3.11.2.1.1. Corner Cutting

As illustrated by the Figure above, in regions where the minimum energy path is curved there is a tendency to pull off the images from the minimum energy path into higher regions of the potential well [5]. As a consequence, saddle points may be overestimated by the plain elastic band method in such cases.

This effect is enhanced for increasing spring constants and is caused by spring force components perpendicular to the minimum energy path.

3.11.2.1.2. Down Sliding

Another unfavorable effect of elastic band calculations is the tendency of images to slide down from the saddle point towards the endpoints. As a result, images tend to avoid the barrier region providing the lowest resolution to the most important part, the saddle point. This may cause an underestimation of transition state energies.

This effect is enhanced for smaller spring constants and is caused by the true force components parallel to the minimum energy path.

Although these problems of the plain elastic band method can be separately addressed by the choice for the spring constant, they cannot be resolved simultaneously because of the opposite influence of this parameter.

3.11.3. The Nudged Elastic Band Method (NEB)

Both problems can be simultaneously resolved by a force projection scheme, i.e. the so-called nudging [4],[5].

Those force components causing corner-cutting and down sliding are selectively projected out of the total force of the elastic band method. Specifically, the component of the spring force perpendicular to the minimum energy path causing corner-cutting and the component of the true force parallel to the minimum energy path being responsible for down sliding are subtracted from the original elastic band total force. The projections are calculated via an estimated tangent to the minimum energy path at each image.

The total force of the nudged elastic band method acting on image i is then given by the sum of the spring force parallel to the local tangent and the true force perpendicular to the local tangent

(4)\[F_{i}^{NEB} = F_{i}^{s\parallel} - F_{i}^{\perp}\]

with the true force component

(5)\[F_{i}^{\perp} = \nabla E (R_{i})^{\perp} = \nabla E (R_{i}) - \nabla E (R_{i}) \tau_{i} \tau_{i}\]

and the spring force component

(6)\[F_{i}^{s\parallel} = k \left( \vert R_{i+1} - R_{i} \vert - \vert R_{i} - R_{i-1} \vert \right) \tau_{i}\]

Herein, \(\tau_{i}\) is the normalized local tangent to the minimum energy path at image \(i\). The above definition of the spring force ensures equal spacing of the images even in regions of high curvature.

3.11.3.1. Tangent Estimation in NEB

The nudged elastic band method relies on a good estimate for the local tangent to the minimum energy path. The originally proposed tangent estimation by bisection of the unit vectors

(7)\[\tau_{i} = \dfrac{R_{i} - R_{i-1}}{\vert R_{i} - R_{i-1} \vert} + \dfrac{R_{i+1} - R_{i}}{\vert R_{i+1} - R_{i} \vert}\]

can develop kinks when the energy along the path changes rapidly, but the restoring force perpendicular to the path is weak (see left panel of Figure below). This may cause problems with the convergence of the optimization process [4].

Improvement of this behavior has been achieved by using only the adjacent image with the higher energy instead of both neighboring images [6].

\begin{equation} \tau_{i} = \Bigg \lbrace \begin{matrix} R_{i+1} - R_{i} \quad if \quad E_{i+1} > E_{i} > E_{i-1} \\ R_{i} - R_{i-1} \quad if \quad E_{i+1} < E_{i} < E_{i-1} \\ \end{matrix} \end{equation}

If the energies of both adjacent images \(E_{i+1}\) and \(E_{i-1}\) are either both higher or both lower than the energy \(E_{i}\) of image \(i\), the tangent is constructed as a weighted average of the vectors to the neighboring two images

\begin{equation} \tau_{i} = \Bigg \lbrace \begin{matrix} (R_{i+1} - R_{i}) \Delta E_{i}^{max} + (R_{i} - R_{i-1}) \Delta E_{i}^{max} \quad if \quad E_{i+1} > E_{i-1} \\ (R_{i+1} - R_{i}) \Delta E_{i}^{min} + (R_{i} - R_{i-1}) \Delta E_{i}^{max} \quad if \quad E_{i+1} < E_{i-1} \\ \end{matrix} \end{equation}

with

(10)\[\begin{split}\Delta E_{i}^{max} = max \big( (E_{i+1} - E_{i}),(E_{i-1} - E_{i}) \big) \\ \Delta E_{i}^{min} = min \big( (E_{i+1} - E_{i}),(E_{i-1} - E_{i}) \big)\end{split}\]

The weighted average is applied only at the extrema of the minimum energy path to smoothly switch between the two possible tangents. In the right panel of the Figure, the improvement due to this modified tangent construction is illustrated [4].

3.11.3.2. The Climbing Image Nudged Elastic Band Method

The nudged elastic band method is able to map out the minimum energy path by a chain of mostly equidistant images. The saddle point, however, may be located somewhat between the images and neither the structure nor the energy of the transition state can be derived. One solution is to set up a second nudged elastic band calculation between the converged images both sides of the expected saddle point thereby refining resolution at the energy barrier.

(11)\[F_{i_{max}} = - \nabla E(R_{i_{max}}) + 2 \nabla E(R_{i_{max}}) \tau_{i_{max}} \tau_{i_{max}}\]

A further rather efficient option is a small modification of the projection scheme for one selected image closest to the saddle point. For the image of highest energy, the spring force is omitted and the inverted parallel component of the true force is included [4].

Due to this modification, the selected image climbs up along the elastic band and converges rigorously to the transition state in an automatic fashion. The equidistance requirement of images close to the saddle point is lifted by this procedure.

3.11.4. Double Nudging

A modification of the nudged elastic band method was proposed to stabilize convergence within a Broyden-Fletcher-Goldfarb-Shanno optimization also implemented in the Transition State Search module, the so-called doubly nudged elastic band method [7]. In plain nudged elastic band method the spring forces act only along the band and have no tendency to shorten the band. Double nudging includes an additional component of the spring force perpendicular to the path. This may improve the stability of the convergence to the MEP, however, it actually prevents from accurately converging to the MEP.

As shown in the Figure above from Sheppard [8], the component of the spring force perpendicular to the tangent

(12)\[F_{i}^{S} = k \left( \vert R_{i+1} - R_{i} \vert - \vert R_{i} - R_{i-1} \vert \right)\]

(13)\[F_{i}^{S \perp} = F_{i}^{S} - F_{i}^{S} \tau_{i} \tau_{i}\]

is constructed. The figure indicates perpendicular and parallel components of the spring force, the right panel shows the plane normal to the tangent. The double nudging force is the component of the perpendicular spring force which is orthogonal to the perpendicular true force

(14)\[F_{i}^{DNEB} = F_{i}^{S} - F_{i}^{S \perp} F_{i}^{\perp} F_{i}^{\perp}\]

[1]	Henry Eyring, “The Activated Complex in Chemical Reactions,” Journal of Physical Chemistry 3, no. 2 (1934): 107.

[2]	C Wert and C Zener, “Interstitial Atomic Diffusion Coefficients,” Physical Review 76, no. 8 (October 15, 1949): 1169.

[3]	Erich Wimmer, Walter Wolf, J rgen Sticht, Saxe, Clint Geller, et al., “Temperature-Dependent Diffusion Coefficients From Ab Initio Computations: Hydrogen, Deuterium, and Tritium in Nickel,” Physical Review B 77, no. 13 (April 2008): 134305.

[4]	(1, 2, 3, 4, 5) Graeme Henkelman, BP Uberuaga, and Hannes Jónsson, “A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths,” Journal of Chemical Physics 113, no. 22 (2000): 9901.

[5]	(1, 2) Hannes Jónsson, Greg Mills, and K W Jacobsen, “Nudged Elastic Band Method for Finding Minimum Energy Paths of Transitions,” Classical and Quantum Dynamics in Condensed Phase Simulations (1998): 385.

[6]	Graeme Henkelman and Hannes Jónsson, “A Dimer Method for Finding Saddle Points on High Dimensional Potential Surfaces Using Only First Derivatives,” Journal of Chemical Physics 111, no. 15 (1999): 7010.

[7]	Semen A Trygubenko and David J Wales, “A Doubly Nudged Elastic Band Method for Finding Transition States,” Journal of Chemical Physics 120, no. 5 (2004): 2082.

[8]	Daniel Sheppard, Rye Terrell, and Graeme Henkelman, “Optimization Methods for Finding Minimum Energy Paths,” Journal of Chemical Physics 128, no. 13 (2008): 134106.

download:pdf