Stationary points are places on the potential energy surface (PES) with a zero gradient, i.e. zero first derivatives of the energy with respect to atomic coordinates. Two types of stationary points are of special importance to chemists. These are minima (reactants, products, intermediates) and first-order saddle points (transition states).
The two types of stationary points can be characterized by the curvature of the PES at these points. At a minimum the Hessian matrix (second derivatives of energy with respect to atomic coordinates) is positive definite, that is the curvature is positive in all directions. If there is one, and only one, negative curvature, the stationary point is a transition state (TS). Because vibrational frequencies are basically the square roots of the curvatures, a minimum has all real frequencies, and a saddle point has one imaginary vibrational “frequency”.
Structure optimizations are most effectively done by so-called quasi-Newton–Raphson methods. They require the exact gradient vector and an approximation to the Hessian matrix. The rate of convergence of the structure optimization depends on anharmonicity of the PES and of the quality of the approximation to the Hessian matrix.
The optimization procedure implemented in statpt belongs to the family of quasi-Newton–Raphsod methods [34]. It is based on the restricted second-order method, which employes Hessian shift parameter in order to control the step length and direction. This shift parameter is determined by the requirement that the step size should be equal to the actual value of the trust radius, tradius, and ensures that the shifted Hessian has the correct eigenvalue structure, all positive for a minimum search, and one negative eigenvalue for a TS search. For TS optimization there is another way of describing the same algorithm, namely as a minimization on the "image" potential. The latter is known as TRIM (Trust Radius Image Minimization) [35].
For TS optimizations the TRIM method implemented in statpt tries to maximize the energy along one of the Hessian eigenvectors, while minimizing it in all other directions. Thus, one “follows” one particular eigenvector, hereafter called the “transition” vector. After computing the Hessian for your guess structure you have to identify which vector to follow. For a good TS guess this is the eigenvector with negative eigenvalue, or imaginary frequency. A good comparison of different TS optimization methods is given in [36].
Structure optimizations using statpt are controlled by the keyword $statpt to be present in the control file. It can be set either manually or by using the stp menu of define. The type of stationary point optimization depends on the value of itrvec specified as an option within $statpt. By default itrvec is set to 0, which implies a structure minimization. A value itrvec > 0 implies a transition state optimization using the eigenvalue-following TRIM algorithm, where the index of the transition vector is specified by itrvec. Note, that statpt orders eigenvalues (and eigenvectors) of the Hessian in ascending order, shifting six (or five in the case of linear molecules) zero translation and rotation eigenvalues to the end.
Note: this order differs from that used for vibrational frequencies in the control file, where rotational and translational eigenvalues are not shifted.
By default a structure optimization is converged when all of the following criteria are met:
The default values for the convergence criteria can be changed using the stp menu of define. The necessary keywords are described in Section 20.2.23 below.
For structure optimization of minima with statpt as relaxation program just use:
jobex &
TS optimizations are performed by the jobex invokation:
jobex -trans &
The choice of the initial Hessian matrix has a great effect on the convergence of the structure optimization. At present, there are three choices for the Hessian matrix in statpt. For minimization, a diagonal matrix or approximate Hessian matrix from a forcefield calculation using uff(see Section 5.4) can be used. For transition state optimizations you have to provide either the “exact” Hessian or results from the lowest eigenvalue search (LES, see Section 14). Note also that you can calculate the Hessian with a smaller basis set and/or at a lower wavefunction level, and use it for higher level structure optimization. Usually, a Hessian matrix calculated in a minimal basis using RI-DFT is good enough for all methods implemented in TURBOMOLE.
statpt automatically takes the best choice of the Hessian from the control file. For minimizations it first looks for the exact Hessian and then for the UFF Hessian. If none of them is found it takes the scaled unit matrix. For transition state optimization the exact Hessian has a higher priority than the results of LES.
The results of LES can be used to obtain an initial Hessian matrix for transition state optimizations involving large molecules, where calculation of the full Hessian is too expensive. Note, that LES calculations for statpt, in addition to the $les keyword require the following keywords to be added manually in the control file:
The default Hessian update for minimization is bfgs, which is likely to remain positive definite. The powell update is the default for transition state optimizations, since the Hessian can develop a negative curvature as the search progresses.
Simply specify the $statpt keyword in the control file and run jobex as explained above. You can very often speedup the optimization by calculating the initial Hessian matrix using uff.
Locating minima on a PES is straightforward. In contrast, transition state optimization requires much more input. The diagonal guess Hessian will almost never work, so you must provide a computed one. The Hessian should be computed at your best guess as to what the TS should be.
The real trick here is to find a good guess for the transition state structure. The closer you are, the better. It is often difficult to guess these structures. One way to obtain a good guess is to built an approximate TS and to perform a constrained minimization by freezing internal coordinates that change most during the reaction. Alternatively, you can generate several structures intermediate to reactants and products, and compute the energy at each point. The maximum energy structure is usually a good guess for the true TS.
After obtaining a reasonable initial guess for the TS structure you have to perform a vibrational analysis (or LES calculation for a large molecule) and to identify the index of the transition vector to follow during the optimization. Ideally, this is a vector with a negative eigenvalue, or "imaginary" frequency. The best way to find the right vector is to use some graphical interface to visualize vibrations. For a reasonable guess structure there should be one vibration that resembles the reaction under study. Remember that statpt uses a different ordering of eigenvalues as compared to the aoforce output—six (five) zero eigenvalues are shifted to the end.
There is an important thing to remember at this point. Even such sophisticated optimization methods like TRIM will not replace your own chemical intuition about where transition states may be located. If you need to restart your run, do so with the coordinates which have the smallest RMS gradient. Note that the energy does not have necessarily to decrease in a transition state search (as opposed to minimizations). It is sometimes necessary to do restart several times (including a recomputation of the Hessian) before the saddle point can be located.
Assuming you do find the TS, it is always a good idea to recompute the Hessian at this structure. It is fairly common, especially when using symmetry, that at your “TS” there is a second imaginary frequency. This means that you have not found the correct TS. The proper procedure is to distort the structure along the “extra” imaginary normal mode using the tool screwer (see Section 1.5). Very often such a distortion requires also lowering the point group symmetry. The distortion must be large enough, otherwise the next run will come back to the invalid structure.