next up previous contents
Next: Results Up: FINDING REARRANGEMENT PATHWAYS Previous: L-BFGS Minimiser   Contents

A Single-ended Method: Eigenvector-following

Single-ended methods use only the function and its derivatives to search for a transition state from the starting point (the initial guess). Since a transition state is a local maximum in one direction but a minimum in all the others it is not generally possible to use standard minimisation methods for this purpose.

Newton-type optimisation methods are based on approximating the objective function locally by a quadratic model and then minimising that function approximately using, for example, the Newton-Raphson (NR) approach [31]. However, the NR algorithm can converge to a stationary point of any index [130]. To converge to a transition state the algorithm may need to start from a geometry at which the Hessian has exactly one negative eigenvalue and the corresponding eigenvector is at least roughly parallel to the reaction coordinate. Locating such a point can be a difficult task and several methods have been developed that help to increase the basins of attraction of transition states, which gives more freedom in choosing the starting geometry [9,22,24,19].

The most widely used single-ended transition state search method is eigenvector-following (EF). In it simplest form it requires diagonalisation of the Hessian matrix [12,14,13,15]. For larger systems it is advantageous to use hybrid EF methods that avoid either calculating the Hessian or diagonalising it [22,133,23,8].

We start by reviewing the theory behind the Newton-Raphson methods. Solving the eigenvalue problem

$\displaystyle {\bf H}{\bf C}={\bf C}{\bf\Lambda}$ (4.21)

yields the matrix $ {\bf C}$, the columns of which are the eigenvectors $ {\bf c}_i$, and the diagonal matrix $ {\bf\Lambda}$, with eigenvalues $ \lambda_i$ on the diagonal. The matrix $ {\bf C}$ defines a unitary transformation to a new orthogonal basis $ \{{\bf c}_i\}$, in which the original Hessian matrix $ {\bf H}$ is diagonal

$\displaystyle {\bf C}^{-1} {\bf H} {\bf C} = {\bf\Lambda},$ (4.22)

so that Equation 2.13 is simplified to

$\displaystyle p'_i = - \dfrac{g'_i}{\lambda_i},$ (4.23)

where $ {\bf g}'$ is the gradient vector in the new basis, which is related to the original gradient vector $ {\bf g}$ as

$\displaystyle {\bf g}'={\bf C}^T {\bf g},$ (4.24)

because of the unitarity of $ {\bf C}$

$\displaystyle {\bf C}^{-1}={\bf C}^{T}.$ (4.25)

The unnormalised search direction $ {\bf p}'$ is defined similarly.

After taking a step of length $ \varpi=\vert{\bf p}\vert$ in the direction $ {\bf p}$, as prescribed by Equation 2.23, the energy changes by

$\displaystyle \Delta V' = - \dfrac{1}{2} {\bf g}' {\bf\Lambda}^{-1} {\bf g}'$ (4.26)

provided our quadratic approximation to $ V({\bf X})$ is perfect. In general $ \lambda_i\in\mathbb{R}$ and terms in Equation 2.26 with $ \lambda_i < 0$ and $ \lambda_i > 0$ increase and decrease the energy, respectively. For isolated molecules and bulk models with periodic boundary conditions the Hessian matrix is singular: there are three zero eigenvalues that correspond to overall translations. Non-linear isolated molecules have three additional zero eigenvalues due to rotations. Luckily, in all cases for which $ \lambda_i=0$ the analytic form of the corresponding eigenvectors $ {\bf c}_i$ is known so the eigenvalue shifting procedure could be applied:

$\displaystyle {\bf H}'={\bf H}+\sum_i \tilde{\lambda}_i {\bf c}_i\otimes{\bf c}_i^T,$ (4.27)

where $ \tilde{\lambda}_i$ are the new eigenvalues. In this work all the zero eigenvalues were shifted by the same amount ($ 10^6$). As a result, $ \det {\bf H}'\neq0$ and problems with applying Equation 2.23 do not arise.

The analytic forms of the $ {\bf c}_i$'s corresponding to translations along the $ x$, $ y$ and $ z$ directions are

$\displaystyle \hat{\bf c}_1=\sqrt{\dfrac{1}{N}}\left( \begin{array}{c} 1 \\ 0 \... ...eft( \begin{array}{c} 0 \\ 0 \\ 1 \\ 0 \\ 0 \\ 1 \\ \vdots \end{array} \right),$ (4.28)

respectively. The eigenvector that corresponds to rotation about the $ x$ axis by a small angle $ \phi$ can be obtained as

$\displaystyle \hat{\bf c}_4=\tilde{\bf R}_x \hat{\bf x} - \hat{\bf x},$ (4.29)

where $ \hat{\bf x}=(x_1,y_1,z_1,x_2,\dots)^T$ is a normalised $ 3N$-dimensional vector describing the position before the rotation, and $ \tilde{\bf R}_x$ is Maclaurin series expansion of the $ 3N$-dimensional rotation matrix [134] $ {\bf R}_x$ with respect to a small angle $ \phi$ truncated to the second order. The displacement vectors due to infinitesimal rotation about the $ x$, $ y$ and $ z$ axes are therefore

$\displaystyle \hat{\bf c}_4=\left( \begin{array}{c} 0 \\ -\dfrac{\phi^2}{2} y_1... ...hi y_2 \\ -\phi x_2 - \dfrac{\phi^2}{2} y_2 \\ 0 \\ \vdots \end{array} \right),$ (4.30)

respectively.

If started from a point at which the Hessian matrix has $ n$ negative eigenvalues, the NR method is expected to converge to a stationary point of index $ n$. However, the step taking procedure may change the index as optimisation proceeds.

Introducing Lagrange multipliers allows a stationary point of specified index to be located. Consider a Lagrangian function

$\displaystyle L=-V'({\bf x}'_0)-\sum_{\alpha=1}^{3N}\left[ g'_\alpha({\bf x}'_0... ...alpha\, {p'}_\alpha^2 - {1\over2} \mu_\alpha({p'}_\alpha^2-c_\alpha^2) \right],$ (4.31)

where a separate Lagrange multiplier $ \mu_\alpha$ is used for every eigendirection $ \alpha $, $ c_\alpha$ is a constraint on the step $ p'_\alpha$, $ \lambda_\alpha$ is the $ \alpha $th eigenvalue of $ {\bf H}'$, which is defined in Equation 2.27, and $ {\bf x}'_0$ is the point about which the potential energy function $ V'({\bf x})$ was expanded. Primes, as before, denote that components are expressed in the orthogonal basis. Differentiating Equation 2.31 with respect to $ {\bf p}'$ and using the condition for a stationary point the optimal step can be obtained:

$\displaystyle p'_\alpha={g'_\alpha({\bf x}'_0)\over \mu_\alpha-\lambda_\alpha}.$ (4.32)

The predicted energy change corresponding to this step is

$\displaystyle \Delta V'= \sum_{\alpha=1}^{3N} {(\mu_\alpha-\lambda_\alpha/2) \over (\mu_\alpha-\lambda_\alpha)^2} g'_\alpha({\bf x}'_0)^2.$ (4.33)

We have some freedom in choosing the Lagrange multipliers $ \mu_\alpha$ as long as for the eigendirections for which an uphill (downhill) step is to be taken $ \mu_\alpha > \lambda_\alpha/2$ ( $ \mu_\alpha < \lambda_\alpha/2$). The following choice of Lagrange multipliers allows the recovery of the Newton-Raphson step in the vicinity of a stationary point [9]:

$\displaystyle \mu_\alpha=\lambda_\alpha\pm{1\over2}\vert\lambda_\alpha\vert\left(1+\sqrt{1+4g'_\alpha({\bf x}'_0)^2/\lambda^2_\alpha}\right),$ (4.34)

with the plus sign for an uphill step, and the minus sign for a downhill step.

In cases when Hessian diagonalisation is impossible or undesirable the smallest eigenvalue and a corresponding eigenvector can be found using an iterative method [135].

next up previous contents
Next: Results Up: FINDING REARRANGEMENT PATHWAYS Previous: L-BFGS Minimiser   Contents
Semen A Trygubenko 2006-04-10