Chapter 9. Eigenvector Following

Chapter 9. Eigenvector Following
	Part I. AMPAC™ 10 User Guide

Table of Contents

Stage 1
Stage 2
Stage 3

(Note: This Chapter of the manual was adapted from the EF Manual by Frank Jensen.) The eigenvector following (EF) geometry optimization algorithm, contributed by Frank Jensen, has been added to AMPAC™. The current version of the EF optimization routine is a combination of the original EF algorithm of Simons et al.^[30] as implemented by Baker^[31] along with the QA algorithm of Culot et al.^[32] Additional features have also been added to improve stability.

Introduction

The EF geometry optimization algorithm is based on a second order Taylor expansion of the energy about the current point. At this point the energy, the gradient, and an estimate of the Hessian matrix are explicitly available. A distinction must be made between searches for minima (entered using the master keyword “ EF ”) and transition states (entered using the master keyword “ TS ”) since for a minimum on the PES, the correct Hessian has only positive eigenvalues by definition and for a TS the correct Hessian has exactly one negative eigenvalue, and the corresponding eigenvector should be in the direction of the reaction coordinate. There are three fundamental steps in determining the next geometry based on this information, and each is described below.

Stage 1

Locate the best step within or on the hypersphere using the current trust radius.

The actual geometry step is determined using a shift factor which ensures that the steplength is within or on the hypersphere. If the Hessian has the correct structure at this point, a pure Newton-Raphson step is attempted. (This corresponds to setting the shift factor to zero.) If this step is found to be longer than the trust radius, a P-RFO step is projected and if this is also too long, then the best step on the hypersphere is determined by the QA formula. This three step procedure is the default and is implemented automatically. Note that certain EF keywords can alter the default procedure (see NONR and RSCAL=n.n below).

Stage 2

Determine if this step should be retained based on various criteria.

Using the geometry increment determined from Stage 1, the new energy and gradient are evaluated. Two sets of criteria are used depending on whether the search is for a minimum or a transition state:

Minimum:: The new energy should be lower than the previous step. If this is not the case, the step is rejected and the trust radius is halved for the new step.
TS:: In this case, two criteria are used in determining whether the step is appropriate. The first is the ratio between the actual and predicted energy change which should ideally be 1.00. If it deviates substantially from this value, the second order Taylor expansion about the present point is no longer valid. Limits on how far from 1.00 the ratio can deviate before the step is rejected are governed by the keywords RMIN=n.n and RMAX=n.n (see below). If the ratio is outside the defined limits, the step is rejected, the trust radius halved, and a new step predicted. The second criteria for TS searches is that the eigenvector along which the energy is being maximized should not change substantially between iterations. The minimum overlap of the TS eigenvector with that of the previous iteration should be larger than a set value (see OMIN=n.n below), otherwise the step is rejected. Large values for RMIN and RMAX as well as setting OMIN to zero effectively inhibits step rejection.

Stage 3

Update the trust radius.

As before, Stage 3 is dependent on whether the EF search is for a minimum or for a transition state.

Minimum:: The trust radius is increased by a factor of 2 if the ratio between actual and predicted energy change is above 0.5. If the ratio is within ±0.1 of 1.0, it is multiplied by 2.
TS:: During a TS search, more restrictive criteria are used. If the ratio is within ±0.1 of 1.0, the trust radius is multiplied by 2. If the ratio is between 0.75 and 0.9, or 1.1 and 1.33, the radius is increased by a factor 2 . If the ratio is less than 0.1 or larger than 3, the trust radius is halved.

In either case above, the updating of the trust radius can be inhibited by specifying the keyword NOUPD (see below). Other trust radius parameters can be adjusted via special keywords (see DMAX=n.n, DDMAX=n.n, and DDMIN=n.n below).

Note that in all cases, as the optimization approaches the stationary point the energy changes become very small, and the actual/predicted ratio may vary wildly. For this reason, step rejection based on this ratio is only performed if both the actual and the predicted energy changes are larger than 0.02. Similarly the trust radius is only updated if both changes are larger than this value, or if the gradient norm is larger than 5. The gradient norm condition is imposed to avoid stalling the optimization when it is distant from the stationary point. Below these limits, a fixed trust radius of 0.10 is used for each step.

EF Strategies

An examination of the Hessian matrix may be obtained by specifying OPTMAX=0. This will cause the routine to process the HESS option, diagonalize the Hessian, and automatically set PRINT=4 (unless PRINT=n is given explicitly). This procedure is useful for selecting the MODE to follow in TS searches. After inspection of the output, the calculation can be restarted by specifying RESTART (to pickup the Hessian from the jobname.res file) or HESS=2. The appropriate value for MODE=n can now be specified with confidence.
Repeated use of HESS=2, varying the value of MODE, is useful for following several different TS modes from the same initial geometry. (Note that the jobname.res and jobname.den files must be copied over the name of the jobname.dat file to be properly read.)
If very highly converged geometries and wavefunctions are required, AMPAC options combined with EF options can be used to obtain them. The default value for GNORM is usually sufficient to ensure that the energy is optimized to within 0.01 kcal/mol. In general, it will be quite difficult to converge structures to a gnorm of less than 0.01. SCFCRT should be tightened and HESS=3 used instead.
The OMIN=n.n, RMIN=n.n and RMAX=n.n features (see below) have been introduced in the current version of EF to improve the stability of TS optimizations. Setting RMIN and RMAX close to 1.00 will give a very stable, albeit very slow, optimization. Wide limits on RMIN and RMAX may in certain cases give a faster convergence, but there is always the risk that very poor steps are accepted, causing the optimization to diverge to unprofitable regions of the potential surface. The default values of 0.00 and 4.00 rarely reject steps which would lead to faster convergence, but may occasionally accept poor steps. If a TS search is not converging properly, the first solution is to lower the limits to RMIN=0.50 and RMIN=2.00. Tighter limits such as these will usually slow the optimization substantially, but may be necessary in some cases.
The EF algorithm has the capability of following Hessian eigenvectors other than the one with the lowest eigenvalue toward a TS. Such higher mode following procedures are always much more difficult to converge than descent along the lowest eigenvector. Ideally, as the optimization progresses, the mode being followed by EF should become the lowest eigenvector. Care must be taken during the optimization, however, that the nature of the mode does not suddenly change, leading to optimization of a different TS than the one desired. OMIN=n.n is designed to ensure that the nature of the TS mode changes gradually by requiring some level of overlap between TS modes. While this concept at first appears quite promising, it is not without problems when the Hessian is updated. As the updated Hessian in each step is only approximately correct (not completely re-computed), there is a upper limit on a practical value for TS mode overlap. Tests have shown that 0.9 is a rough practical limit to the value of OMIN. However, a combination of RECALC=1 and OMIN=0.9 allows EF to locate very difficult transition state structures.
Following modes other than those with the lowest eigenvalue toward a TS indicates that the starting geometry is not “close” to the desired TS. It may possess multiple negative eigenvalues and be at a complex point on the potential surface. There are cases where it is very difficult to locate a starting geometry which has a Hessian with a well-defined eigenvector, and mode following may be of some use.
Using the RECALC=1 option will force computation of a full Hessian matrix at each step in the calculation. (Using the HESS=5 keyword will use the rapid LTRD algorithm for this purpose.) This is very costly in terms of time and computational effort but can be useful in locating difficult points on the PES.

Dedicated EF Keywords

DDMAX	Define the maximum size of the trust radius.
DDMIN	Define the minimum size of the trust radius.
DMAX	Define the initial trust radius.
EF	Use the eigenvector following method to locate a minimum.
GNMIN	Activate gradient test for accepting geometry steps.
HESS	Specify the source of the Hessian matrix.
IUPD	Choice of update method for the Hessian matrix.
MODE	Specify eigenvector to follow during optimization.
NONR	Specify P-RFO method for geometry projection.
NOUPD	Suppress updating of the trust radius at Stage 3.
OMIN	Specify minimum overlap between successive TS search vectors.
RECALC	Specify interval (in number of steps) for Hessian recalculation.
RMAX	Adjust maximum criterion for accepting geometry steps.
RMIN	Adjust minimum criterion for accepting geometry steps.
RSCAL	Scale the P-RFO step.
TS	Use the eigenvector following method to locate a transition state.

The following standard AMPAC keywords are also active in the EF routine, with similar functions to those they serve in the main program:

^[30]J. Simons, A. Banerjee, N. Adams, and R. Shepard. Journal of Physical Chemistry. 1985. 89. 52.

^[31]J. Baker. Journal of Computational Chemistry. 1986. 7. 385.

^[32]P. Culot, G. Dive, V.H. Nguyen, and J.M. Ghuysen. Theoretica Chimica Acta. 1992. 82. 189.


Chapter 8. CHN Methods		Chapter 10. Electrostatic Potential