Gradient Extremals

   This section of the manual, as well as the source code
to trace gradient extremals was written by Frank Jensen of 
the University of Aarhus.

   A Gradient Extremal (GE) curve consists of points where
the gradient norm on a constant energy surface is 
stationary.  This is equivalent to the condition that
the gradient is an eigenvector of the Hessian.  Such GE
curves radiate along all normal modes from a stationary
point, and the GE leaving along the lowest normal mode
from a minimum is the gentlest ascent curve.  This is not
the same as the IRC curve connecting a minimum and a TS,
but may in some cases be close.

   GEs may be divided into three groups:  those leading
to dissociation, those leading to atoms colliding, and 
those which connect stationary points.  The latter class
allows a determination of many (all?) stationary points on
a PES by tracing out all the GEs. Following GEs is thus a
semi-systematic way of mapping out stationary points.  The
disadvantages are:
   i) There are many (but finitely many!) GEs for a 
      large molecule.
  ii) Following GEs is computationally expensive.
 iii) There is no control over what type of 
      stationary point (if any) a GE will lead to. 

   Normally one is only interested in minima and TSs, but
many higher order saddle points will also be found.
Furthermore, it appears that it is necessary to follow GEs
radiating also from TSs and second (and possibly also 
higher) order saddle point to find all the TSs.

   A rather complete map of the extremals for the H2CO
potential surface is available in a paper which explains
the points just raised in greater detail:
   K.Bondensgaard, F.Jensen, 
       J.Chem.Phys. 104, 8025-8031(1996).
An earlier paper gives some of the properties of GEs:
   D.K.Hoffman, R.S.Nord, K.Ruedenberg,
       Theor. Chim. Acta 69, 265-279(1986).

   There are two GE algorithms in GAMESS, one due to Sun 
and Ruedenberg (METHOD=SR), which has been extended to
include the capability of locating bifurcation points and
turning points, and another due to Jorgensen, Jensen, and
Helgaker (METHOD=JJH):
   J. Sun, K. Ruedenberg, J.Chem.Phys. 98, 9707-9714(1993)
   P. Jorgensen, H. J. Aa. Jensen, T. Helgaker
       Theor. Chim. Acta 73, 55 (1988).

   The Sun and Ruedenberg method consist of a predictor
step taken along the tangent to the GE curve, followed by
one or more corrector steps to bring the geometry back to
the GE.  Construction of the GE tangent and the corrector
step requires elements of the third derivative of the
energy, which is obtained by a numerical differentiation 
of two Hessians.  This puts some limitations on which 
systems the GE algorithm can be used for.  First, the
numerical differentiation of the Hessian to produce third
derivatives means that the Hessian should be calculated by
analytical methods, thus only those types of wavefunctions
where this is possible can be used.  Second, each
predictor/corrector step requires at least two Hessians,
but often more.  Maybe 20-50 such steps are necessary for
tracing a GE from one stationary point to the next.  A
systematic study of all the GE radiating from a stationary
point increases the work by a factor of ~2*(3N-6).  One 
should thus be prepared to invest at least hundreds, and 
more likely thousands, of Hessian calculations.  In other
words, small systems, small basis sets, and simple wave-
functions. 

   The Jorgensen, Jensen, and Helgaker method consists of
taking a step in the direction of the chosen Hessian
eigenvector, and then a pure NR step in the perpendicular
modes.  This requires (only) one Hessian calculation for 
each step.  It is not suitable for following GEs where the
GE tangent forms a large angle with the gradient, and it 
is incapable of locating GE bifurcations.

   Although experience is limited at present, the JJH
method does not appear to be suitable for following GEs in
general (at least not in the current implementation).
Experiment with it at your own risk!

   The flow of the SR algorithm is as follows:  A 
predictor geometry is produced, either by jumping away
from a stationary point, or from a step in the tangent
direction from the previous point on the GE.  At the 
predictor geometry, we need the gradient, the Hessian, and
the third derivative in the gradient direction.  Depending
on HSDFDB, this can be done in two ways.  If .TRUE. the
gradient is calculated, and two Hessians are calculated at 
SNUMH distance to each side in the gradient direction.
The Hessian at the geometry is formed as the average of
the two displaced Hessians.  This corresponds to a double-
sided differentiation, and is the numerical most stable
method for getting the partial third derivative matrix.
If HSDFDB = .FALSE., the gradient and Hessian are 
calculated at the current geometry, and one additional 
Hessian is calculated at SNUMH distance in the gradient 
direction.  This corresponds to a single-sided differen-
tiation.  In both cases, two full Hessian calculations are 
necessary, but HSDFDB = .TRUE. require one additional
wavefunction and gradient calculation.  This is usually
a fairly small price compared to two Hessians, and the
numerically better double-sided differentiation has
therefore been made the default.

   Once the gradient, Hessian, and third derivative is 
available, the corrector step and the new GE tangent are
constructed.  If the corrector step is below a threshold,
a new predictor step is taken along the tangent vector. 
If the corrector step is larger than the threshold, the
correction step is taken, and a new micro iteration is
performed.  DELCOR thus determines how closely the GE will
be followed, and DPRED determine how closely the GE path 
will be sampled.

   The construction of the GE tangent and corrector step
involve solution of a set of linear equations, which in
matrix notation can be written as Ax=B. The A-matrix is
also the second derivative of the gradient norm on the
constant energy surface.

   After each corrector step, various things are printed
to monitor the behavior:  The projection of the gradient
along the Hessian eigenvalues (the gradient is parallel
to an eigenvector on the GE), the projection of the GE
tangent along the Hessian eigenvectors, and the overlap
of the Hessian eigenvectors with the mode being followed
from the previous (optimzed) geometry.  The sign of these
overlaps are not significant, they just refer to an
arbitrary phase of the Hessian eigenvectors.

   After the micro iterations has converged, the Hessian
eigenvector curvatures are also displayed, this is an
indication of the coupling between the normal modes.  The
number of negative eigenvalues in the A-matrix is denoted
the GE index.  If it changes, one of the eigenvalues must
have passed through zero.  Such points may either be GE 
bifurcations (where two GEs cross) or may just be "turning
points", normally when the GE switches from going uphill
in energy to downhill, or vice versa.  The distinction is
made based on the B-element corresponding to the A-matrix
eigenvalue = 0. If the B-element = 0, it is a bifurcation,
otherwise it is a turning point.

   If the GE index changes, a linear interpolation is
performed between the last two points to locate the point
where the A-matrix is singular, and the corresponding
B-element is determined.  The linear interpolation points
will in general be off the GE, and thus the evaluation of
whether the B-element is 0 is not always easy.  The
program additionally evaluates the two limiting vectors 
which are solutions to the linear sets of equations, these 
are also used for testing whether the singular point is a 
bifurcation point or turning point.

   Very close to a GE bifurcation, the corrector step 
become numerically unstable, but this is rarely a problem
in practice.  It is a priori expected that GE bifurcation
will occur only in symmetric systems, and the crossing GE
will break the symmetry.  Equivalently, a crossing GE may
be encountered when a symmetry element is formed, however
such crossings are much harder to detect since the GE
index does not change, as one of the A-matrix eigenvalues 
merely touches zero.  The program prints an message if
the absolute value of an A-matrix eigenvalue reaches a 
minimum near zero, as such points may indicate the 
passage of a bifurcation where a higher symmetry GE
crosses.  Run a movie of the geometries to see if a more 
symmetric structure is passed during the run.

   An estimate of the possible crossing GE direction is
made at all points where the A-matrix is singular, and two
perturbed geometries in the + and - direction are written 
out.  These may be used as predictor geometries for 
following a crossing GE.  If the singular geometry is a
turning point, the + and - geometries are just predictor 
geometries on the GE being followed. 

   In any case, a new predictor step can be taken to trace
a different GE from the newly discovered singular point,
using the direction determined by interpolation from the 
two end point tangents (the GE tangent cannot be uniquely
determined at a bifurcation point).  It is not possible to
determine what the sign of IFOLOW should be when starting 
off along a crossing GE at a bifurcation, one will have to 
try a step to see if it returns to the bifurcation point 
or not.

   In order to determine whether the GE index change it
is necessary to keep track of the order of the A-matrix
eigenvalues.  The overlap between successive eigenvectors
are shown as "Alpha mode overlaps".

Things to watch out for:
1) The numerical differentiation to get third derivatives
requires more accuracy than usual.  The SCF convergence 
should be at least 100 times smaller than SNUMH, and
preferably better.  With the default SNUMH of 10**(-4)
the SCF convergence should be at least 10**(-6).  Since
the last few SCF cycles are inexpensive, it is a good idea
to tighten the SCF convergence as much as possible, to
maybe 10**(-8) or better.  You may also want to increase 
the integral accuracy by reducing the cutoffs (ITOL and 
ICUT) and possibly also try more accurate integrals 
(INTTYP=HONDO).  The CUTOFF in $TRNSFM may also be reduced 
to produce more accurate Hessians.  Don't attempt to use a
value for SNUMH below 10**(-6), as you simply can't get 
enough accuracy.  Since experience is limited at present,
it is recommended that some tests runs are made to learn
the sensitivity of these factors for your system. 

2) GEs can be followed in both directions, uphill or
downhill. When stating from a stationary point, the
direction is implicitly given as away from the stationary
point.  When starting from a non-stationary point, the "+"
and "-" directions (as chosen by the sign of IFOLOW) 
refers to the gradient direction.  The "+" direction is
along the gradient (energy increases) and "-" is opposite
to the gradient (energy decreases). 

3) A switch from one GE to another may be seen when two
GE come close together.  This is especially troublesome 
near bifurcation points where two GEs actually cross.  In
such cases a switch to a GE with -higher- symmetry may
occur without any indication that this has happened, 
except possibly that a very large GE curvature suddenly 
shows up.  Avoid running the calculation with less 
symmetry than the system actually has, as this increases
the likelihood that such switches occuring.  Fix: alter
DPRED to avoid having the predictor step close to the 
crossing GE.

4) "Off track" error message:  The Hessian eigenvector
which is parallel to the gradient is not the same as
the one with the largest overlap to the previous
Hessian mode.  This usually indicate that a GE switch
has occured (note that a switch may occur without this
error message), or a wrong value for IFOLOW when starting
from a non-stationary point. Fix: check IFOLOW, if it is 
correct then reduce DPRED, and possibly also DELCOR.

5) Low overlaps of A-matrix eigenvectors.  Small overlaps
may give wrong assignment, and wrong conclusions about GE
index change. Fix: reduce DPRED.

6) The interpolation for locating a point where one of the
A-matrix eigenvalues is zero fail to converge.  Fix: 
reduce DPRED (and possibly also DELCOR) to get a shorther 
(and better) interpolation line. 

7) The GE index changes by more than 1.  A GE switch may
have occured, or more than one GE index change is located 
between the last and current point.  Fix: reduce DPRED to 
sample the GE path more closely.

8) If SNRMAX is too large the algorithm may try to locate
stationary points which are not actually on the GE being 
followed.  Since GEs often pass quite near a stationary 
point, SNRMAX should only be increased above the default 
0.10 after some consideration.