Hessian matrix in FORCE calculations

The Hessian matrix is the matrix of second derivatives of the energy with respect to geometry. The most important Hessian is that used in the FORCE calculation. Normal modes are expressed as Cartesian displacements, consequently the Hessian is based on Cartesian rather than internal coordinates.

Although first derivatives are relatively easy to calculate, second derivatives are not. The simplest, although not an elegant, way to calculate [73] second derivatives is to calculate first derivatives for a given geometry, then perturb the geometry, do an SCF calculation on the new geometry, and re-calculate the derivatives. The second derivatives can then be calculated from the difference of the two first derivatives divided by the step size. This method, which is used in the EigenFollowing routine, is called `single-sided' derivatives.

The Hessian is quite sensitive to geometry, and should only be evaluated at stationary points. Because of this sensitivity, "double-sided" derivatives   are used:

begin{displaymath}H_{i,j} = frac{g_i^{+delta_j}-g_i^{-delta_j}}{2delta}.end{displaymath}

Note the asymmetry in the treatment of the Cartesian coordinates i and j. It can be shown that

begin{displaymath}frac{g_j^{+delta_i}-g_j^{-delta_i}}{2delta} = frac{g_i^{+delta_j}-g_i^{-delta_j}}{2delta}.end{displaymath}

To help improve precision, the Hessian is calculated from

begin{displaymath}H_{i,j} = frac{1}{2}left ( frac{g_j^{+delta_i}-g_j^{-del......a} + frac{g_i^{+delta_j}-g_i^{-delta_j}}{2delta} right ).end{displaymath}


Calculation of the Hessian matrix

The elements of the Hessian are defined as:

\begin{displaymath}H_{i,j} = \frac{\delta^2E}{\delta x_i\delta x_j} \end{displaymath}

and are generated by use of finite displacements, that is, for each atomic coordinate xi, the coordinate is first incremented by a small amount, the gradients calculated, then the coordinate is decremented and the gradients re-calculated. The second derivative is then obtained from the difference of the two derivatives and the step size:

\begin{displaymath}H_{i,j} = \frac{(\frac{\delta E}{\delta x_i})_{_{+0.5\Delta x... ...rac{\delta E}{\delta x_i})_{_{-0.5\Delta x_j}}} {\Delta x_j}. \end{displaymath}

This is done for all 3N Cartesian coordinates. Because the Hessian is symmetric, that is Hi,j=Hj,i, the random errors that occur in the gradient calculation can be reduced (by a factor of (1/ 2)1/2) by re-defining the Hessian as:

\begin{displaymath}H_{i,j} = \frac{1}{2}\left(\frac{(\frac{\delta E}{\delta x_i}... ...ta E}{\delta x_j})_{_{-0.5\Delta x_i}}} {\Delta x_i}\right). \end{displaymath}

A call to the energy - gradient function COMPFG will generate the gradients in kcal/mol/Ångstrom at a given geometry. These can then be converted into millidynes/Ångstrom (or 105 dynes/cm) as follows:

\begin{displaymath}H_{i,j} (\mathrm{millidynes/\AA}) = 10^5\frac{({\mathrm{Kcal\... ...athrm{Mole\ to\ molecule})}H_{i,j}(\mathrm{kcal/mol/\AA ^2 }) \end{displaymath}

or

\begin{displaymath}H_{i,j} (\mathrm{millidynes/\AA}) = 10^5\frac{4.184*10^3*10^7}{(10^{-8*2}) (6.022*10^{23})}H_{i,j}(\mathrm{kcal/mol/\AA ^2}). \end{displaymath}

Diagonalization of this matrix yields the force constants of the system.

In order to calculate the vibrational frequencies, the Hessian matrix is first mass-weighted:

\begin{displaymath}H^m_{i,j} = \frac{H_{i,j}}{\sqrt{M_i*M_j}}. \end{displaymath}

Then the Hessian is converted from millidynes per Ångstrom to dynes per centimeter by multiplying by 105.

Diagonalization of this matrix yields eigenvalues, e, which represent the quantities (k/m)1/2 , from which the vibrational frequencies can be calculated:

\begin{displaymath}\bar{\nu}_i = \frac{1}{2\pi c}\sqrt{\epsilon_i}. \end{displaymath}