9.1.2 Theoretical Background

Geometry optimization refers to the determination of stationary points, principally minima and transition states, on molecular potential energy surfaces. It is an iterative process, requiring the calculation of the energy, gradient and (possibly) Hessian at each optimization cycle. The energy, gradient and Hessian information is used to compute an optimization step, $𝐡$, which displaces the geometry to be closer to the target stationary point. This process is repeated until convergence is attained. The art of geometry optimization lies in calculating the step, $𝐡$, so as to converge in as few cycles as possible.

Consider the energy, $E({𝐱}_0)$ at some point ${𝐱}_0$ on a potential energy surface. We can express the energy at a nearby point $𝐱={𝐱}_0+𝐡$ by means of the Taylor series

If we knew the exact form of the energy functional $E(𝐱)$ and all its derivatives, we could move from the current point ${𝐱}_0$ directly to a stationary point, (i.e., we would know exactly what the step $𝐡$ ought to be). Since we typically know only the lower derivatives of $E(𝐱)$ at best, then we can estimate the step $𝐡$ by differentiating the Taylor series with respect to $𝐡$, keeping only the first few terms on the right hand side, and setting the left hand side, $dE({𝐱}_0+𝐡)/d𝐡$, to zero, which is the value it would have at a genuine stationary point. Thus

from which

where

Equation (9.3) is known as the Newton-Raphson step. It is the major component of almost all geometry optimization algorithms in quantum chemistry.

The above derivation assumed exact first (gradient) and second (Hessian) derivative information. Analytical gradients are available for most methodologies supported in Q-Chem, however, analytical second derivatives are less likely to be available, see Table 9.2. Furthermore, even if they were available, it would not necessarily be advantageous to use them as their evaluation is computationally demanding, and efficient optimizations can be performed without an exact Hessian. An excellent compromise in practice is to begin with an approximate Hessian matrix, and update this using gradient and displacement information generated as the optimization progresses. In this way the starting Hessian can be “improved” at essentially no cost. Using Eq. (9.3) with an approximate Hessian is called the quasi Newton-Raphson step.

The nature of the Hessian matrix (in particular its eigenvalue structure) plays a crucial role in a successful optimization. All stationary points on a potential energy surface have a zero gradient vector, however, the character of the stationary point (i.e., what type of structure it corresponds to) is determined by the Hessian matrix. Diagonalization of the Hessian yields a set of mutually orthogonal directions on the energy surface (the eigenvectors) together with the curvature along these directions (the eigenvalues). At a local minimum (corresponding to a well in the potential energy surface) the curvature along all of these directions must be positive, reflecting the fact that a small displacement along any of these directions causes the energy to rise. At a transition state, the curvature is negative (i.e., the energy is a maximum) along one direction, but positive along all other directions. Thus, for a stationary point to be a transition state, the Hessian matrix at that point must have one, and only one, negative eigenvalue, while for a minimum the Hessian must have all positive eigenvalues. In the latter case the Hessian is called positive definite. If searching for a minimum it is important that the Hessian matrix be positive definite; in fact, unless the Hessian is positive definite there is no guarantee that the step predicted by Eq. (9.3) is will actually lower the energy. Similarly, for a transition state search, the Hessian must have one negative eigenvalue. Maintaining the correct Hessian eigenvalue structure is not difficult for minimization, but can be problematic when searching for a transition state.

In a diagonal Hessian representation the Newton-Raphson step can be written

where ${𝐮}_i$ and $b_i$ are the eigenvectors and eigenvalues of the Hessian matrix $𝐇$ and $F_i={𝐮}_i^t𝐠$ is the component of $𝐠$ along the local direction (eigenmode) ${𝐮}_i$. As discussed by Simons et al., ¹⁰⁹⁸ Simons J. et al.
J. Phys. Chem.
(1983), 87, pp. 2745. Link the Newton-Raphson step can be considered as minimizing along directions ${𝐮}_i$, which have positive eigenvalues, and maximizing along directions with negative eigenvalues. Thus, if the user is searching for a minimum and the Hessian matrix is positive definite, then the Newton-Raphson step is appropriate since it is attempting to minimize along all directions simultaneously. However, if the Hessian has one or more negative eigenvalues, then the basic Newton-Raphson step is not appropriate for a minimum search, since it will be maximizing and not minimizing along one or more directions. Exactly the same arguments apply during a transition state search except that the Hessian must have one negative eigenvalue, because the user has to maximize along one direction. However, there must be only one negative eigenvalue. A positive definite Hessian is a disaster for a transition state search because the Newton-Raphson step will head towards a minimum.

If firmly in a region of the potential energy surface with the right Hessian character, then a careful search (based on the Newton-Raphson step) will almost always lead to a stationary point of the correct type. However, this is only true if the Hessian is exact. If the Hessian is being updated approximately, then there is no guarantee that the Hessian eigenvalue structure will be preserved from one cycle to the next unless one is very careful during the update. Updating procedures that “guarantee” conservation of a positive definite Hessian do exist (or at least warn the user if the update is likely to introduce negative eigenvalues). This can be very useful during a minimum search; but there are no such guarantees for preserving the Hessian character required for a transition state.

In addition to the difficulties in retaining the correct Hessian character, there is the matter of obtaining a “correct” Hessian in the first instance. This is particularly acute for a transition state search. For a minimum search it is possible to “guess” a reasonable, positive-definite starting Hessian (for example, by carrying out a molecular mechanics minimization and calculating the Hessian, also using molecular mechanics, at the minimum) but this option is usually not available for transition states. Even if the user calculates the Hessian exactly at the starting geometry, the guess for the structure may not be sufficiently accurate, and the expensive, exact Hessian may not have the desired eigenvalue structure.

Consequently, particularly for a transition state search, an alternative to the basic Newton-Raphson step is clearly needed, especially when the Hessian matrix is inappropriate for the stationary point being sought.

One of the first algorithms that was capable of taking corrective action during a transition state search if the Hessian had the wrong eigenvalue structure, was developed by Poppinger, ⁹⁷⁶ Poppinger D.
Chem. Phys. Lett.
(1975), 35, pp. 550. Link who suggested that, instead of taking the Newton-Raphson step, if the Hessian had all positive eigenvalues, the lowest Hessian mode be followed uphill; whereas, if there were two or more negative eigenvalues, the mode corresponding to the least negative eigenvalue be followed downhill. While this step should lead the user back into the right region of the energy surface, it has the disadvantage that the user is maximizing or minimizing along one mode only, unlike the Newton-Raphson step which maximizes/minimizes along all modes simultaneously. Another drawback is that successive such steps tend to become linearly dependent, which degrades most of the commonly used Hessian updates.