A.3 Eigenvector-Following (EF) Algorithm

The work of Cerjan and Miller,[Cerjan and Miller(1981)] and later Simons and co-workers,[Simons et al.(1983)Simons, Jørgensen, Taylor, and Ozment, Banerjee et al.(1985)Banerjee, Adams, Simons, and Shepard] showed that there was a better step than simply directly following one of the Hessian eigenvectors. A simple modification to the Newton-Raphson step is capable of guiding the search away from the current region towards a stationary point with the required characteristics. This is

$\begin{equation} \label{eq:a5t} \ensuremath{\mathbf{h}} = \sum _ i \frac{-F_ i }{\left( b_ i -\lambda \right)} \ensuremath{\mathbf{u}}_ i \end{equation}$

(A.6)

in which $\lambda$ can be regarded as a shift parameter on the Hessian eigenvalue $b_ i$ . Scaling the Newton-Raphson step in this manner effectively directs the step to lie primarily, but not exclusively (unlike Poppinger’s algorithm[Poppinger(1975)]), along one of the local eigenmodes, depending on the value chosen for $\lambda$ . References Cerjan:1981,Simons:1983,Banerjee:1985 all use the same basic approach of Eq. (A.6) but differ in the means of determining the value of $\lambda$ .

The EF algorithm[Baker(1986)] uses the rational function approach presented in Refs. Banerjee:1985, yielding an eigenvalue equation of the form

$\begin{equation} \label{eq:a6} \left( {{\begin{array}{*{20}c} {\rm {\bf H}} \hfill & {\rm {\bf g}} \hfill \\ {{\rm {\bf g}}^ t} \hfill & 0 \hfill \\ \end{array} }} \right)\left( {{\begin{array}{*{20}c} {\rm {\bf h}} \hfill \\ 1 \hfill \\ \end{array} }} \right)=\lambda \left( {{\begin{array}{*{20}c} {\rm {\bf h}} \hfill \\ 1 \hfill \\ \end{array} }} \right) \end{equation}$

(A.7)

from which a suitable $\lambda$ can be obtained. Expanding Eq. (A.7) yields

$\begin{equation} \label{eq:a7a} (\ensuremath{\mathbf{H}}-\lambda )\ensuremath{\mathbf{h}} + \ensuremath{\mathbf{g}} = 0 \end{equation}$

(A.8)

and

$\begin{equation} \label{eq:a7b} \ensuremath{\mathbf{g}}^ t\ensuremath{\mathbf{h}} = \lambda \end{equation}$

(A.9)

In terms of a diagonal Hessian representation, Eq. (A.8) rearranges to Eq. (A.6), and substitution of Eq. (A.6) into the diagonal form of Eq. (A.9) gives

$\begin{equation} \sum _ i {\frac{-F_ i^2 }{\left( {b_ i -\lambda } \right)}} =\lambda \end{equation}$

(A.10)

which can be used to evaluate $\lambda$ iteratively.

The eigenvalues, $\lambda$ , of the RFO equation Eq. (A.7) have the following important properties:[Banerjee et al.(1985)Banerjee, Adams, Simons, and Shepard]

The $(n+1)$ values of $\lambda$ bracket the $n$ eigenvalues of the Hessian matrix $\lambda _ i < b_ i < \lambda _{i+1}$ .
At a stationary point, one of the eigenvalues, $\lambda$ , of Eq. (A.7) is zero and the other $n$ eigenvalues are those of the Hessian at the stationary point.
For a saddle point of order $m$ , the zero eigenvalue separates the $m$ negative and the $(n-m)$ positive Hessian eigenvalues.

This last property, the separability of the positive and negative Hessian eigenvalues, enables two shift parameters to be used, one for modes along which the energy is to be maximized and the other for which it is minimized. For a transition state (a first-order saddle point), in terms of the Hessian eigenmodes, we have the two matrix equations

$\begin{equation} \label{eq:a9a} \left( {{\begin{array}{*{20}c} {b_1 } \hfill & {F_1 } \hfill \\ {F_1 } \hfill & 0 \hfill \\ \end{array} }} \right)\left( {{\begin{array}{*{20}c} {h_1 } \hfill \\ 1 \hfill \\ \end{array} }} \right)=\lambda _ p \left( {{\begin{array}{*{20}c} {h_1 } \hfill \\ 1 \hfill \\ \end{array} }} \right) \end{equation}$

(A.11)

$\begin{equation} \label{eq:a9b} \left( {{\begin{array}{*{20}c} {b_2 } \hfill & \hfill & \hfill & {F_2 } \hfill \\ \hfill & \ddots \hfill & {\rm {\bf 0}} \hfill & \vdots \hfill \\ \hfill & {\rm {\bf 0}} \hfill & {b_ n } \hfill & {F_ n } \hfill \\ {F_2 } \hfill & \cdots \hfill & {F_ n } \hfill & 0 \hfill \\ \end{array} }} \right)\left( {{\begin{array}{*{20}c} {h_2 } \hfill \\ \vdots \hfill \\ {h_ n } \hfill \\ 1 \hfill \\ \end{array} }} \right)=\lambda _ n \left( {{\begin{array}{*{20}c} {h_2 } \hfill \\ \vdots \hfill \\ {h_ n } \hfill \\ 1 \hfill \\ \end{array} }} \right) \end{equation}$

(A.12)

where it is assumed that we are maximizing along the lowest Hessian mode $\ensuremath{\mathbf{u}}_1$ . Note that $\lambda _ p$ is the highest eigenvalue of Eq. (A.11), which is always positive and approaches zero at convergence, and $\lambda _ n$ is the lowest eigenvalue of Eq. (A.12), which it is always negative and again approaches zero at convergence.

Choosing these values of $\lambda$ gives a step that attempts to maximize along the lowest Hessian mode, while at the same time minimizing along all the other modes. It does this regardless of the Hessian eigenvalue structure (unlike the Newton-Raphson step). The two shift parameters are then used in Eq. (A.6) to give the final step

$\begin{equation} \ensuremath{\mathbf{h}} = \frac{-F_1 }{(b_1 -\lambda _ p )} \ensuremath{\mathbf{u}}_1 -\sum \limits _{i=2}^ n \frac{-F_ i }{(b_ i -\lambda _ n )}\ensuremath{\mathbf{u}}_ i \end{equation}$

(A.13)

If this step is greater than the maximum allowed, it is scaled down. For minimization only one shift parameter, $\lambda _ n$ , would be used which would act on all modes.

In Eq. (A.11) and Eq. (A.12) it was assumed that the step would maximize along the lowest Hessian mode, $b_1$ , and minimize along all the higher modes. However, it is possible to maximize along modes other than the lowest, and in this way perhaps locate transition states for alternative rearrangements/dissociations from the same initial starting point. For maximization along the $k$ th mode (instead of the lowest mode), Eq. (A.11) is replaced by

$\begin{equation} \left( {{\begin{array}{*{20}c} {b_ k } \hfill & {F_ k } \hfill \\ {F_ k } \hfill & 0 \hfill \\ \end{array} }} \right)\left( {{\begin{array}{*{20}c} {h_ k } \hfill \\ 1 \hfill \\ \end{array} }} \right)=\lambda _ p \left( {{\begin{array}{*{20}c} {h_ k } \hfill \\ 1 \hfill \\ \end{array} }} \right) \end{equation}$

(A.14)

and Eq. (A.12) would now exclude the $k$ th mode but include the lowest. Since what was originally the $k$ th mode is the mode along which the negative eigenvalue is required, then this mode will eventually become the lowest mode at some stage of the optimization. To ensure that the original mode is being followed smoothly from one cycle to the next, the mode that is actually followed is the one with the greatest overlap with the mode followed on the previous cycle. This procedure is known as mode following. For more details and some examples, see Ref. Baker:1986.