I am trying to understand the construction of Hessian matrix using numerical differentiation (NUMFREQ option in ORCA for example). I understand that Hessian is a square matrix containing second derivatives of energy with respect to 3N atomic coordinates. If I try to perform a NUMFREQ calculation in ORCA, I see that the atoms are displaced in +x,+y,+z,-x,-y,-z directions and energies are evaluated for each of those structures. I figured out the calculation of first derivatives from the corresponding SCF energies. But I am unclear on how to get second mixed partial derivatives from the first derivatives. I tried asking this question in ChatGPT, and I am attaching the answer I got below (after Andrea Pellegrini's polishing of the formatting):

$$ \frac{d^2E}{(dR_idR_j)} = \frac{\frac{dE}{dR_i} \frac{dE}{dR_j}}{(2 \delta R)} $$

I am not sure if this expression is right. Any comments on how this calculation is actually done is welcome.


1 Answer 1


Finite differences

Let's take a general $2n+1$ point central finite difference rule with uniform grid spacing $$ f'(x) \approx h^{-1} \sum_{i=-n}^{n} w_i f(x+ ih), $$ such as two-point central finite differences with $n=1$ and ${\bf w}=(-1/2, 0, 1/2)$ which yields the approximation $$ f'(x) \approx \frac {f(x+h)-f(x-h)} {2h}. $$

Geometry optimization

For geometrical gradients, our function is the energy $E(\{\bf R\})$ where ${\bf R}_I = (X_I, Y_I, Z_I)$ are the coordinates of the $I$-th nucleus. (If we have $N$ nuclei, $\{\bf R\}$ is a vector with $3N$ elements, consisting of the $x$, $y$, and $z$ coordinates of each $N$ nuclei.)

The $k$-th component ($x$, $y$ or $z$) of the force acting on nucleus $I$ is $$ F_{Ik} = -\frac {\partial E} {\partial R_{Ik}}, $$ and it can be evaluated efficiently with analytical methods for a lot of quantum chemical methods. Note that here I use two indices: one for the nucleus and another for the cartesian component.

Simplification of the notation with compound indices

I will now simplify the notation $Ik \to i$. If you number the nuclei from 0 to N-1, then the x, y, and z coordinates of the Ith nucleus will be stored at indices $i=3I$, $i=3I+1$ and $i=3I+2$, respectively. I will also denote the vector of nuclear coordinates with ${\bf x}$ and then the force is $$ F_{i} = -\frac {\partial E} {\partial x_i}. $$


The Hessian is the second derivative, and like you said, it is a $3N \times 3N$ square matrix. Analytical Hessians are less often available than analytical gradients, and they can be evaluated either seminumerically, if analytical gradients are available, or fully numerically by applying two rounds of finite differences.

The Hessian is $$ H_{ij} = \frac {\partial^2 E} {\partial x_i \partial x_j} $$ where I now use the compound indexing.

Seminumerical scheme

Note that if I know the gradients $g_i = \partial E / \partial x_i = -F_i$, I can evaluate the Hessian with finite differences as

$$ H_{ij} = \frac {\partial^2 E} {\partial x_i \partial x_j} =\frac {\partial g_i} {\partial x_j} $$

This form is explicit, since it tells you that the Hessian measures how much the $i$th gradient changes when the $j$th coordinate is modified. How is this evaluated? Let's take the two-point equation above for simplicity

$$ H_{ij} =\frac {\partial g_i} {\partial x_j} \approx \frac {g_i({\bf x}+h\hat{e}_j)-g_i({\bf x}-h\hat{e}_j)} {2h}. $$

That is, for each column $j$ in the Hessian, you get the row $i$ by taking the difference of gradients of all the atoms where the gradients are evaluated with coordinates where the $j$-th coordinate has been modified. This means that for this rule, you need to do 2 gradient evaluations for every $3N$ columns, or $6N$ gradient calculations.

Note in the above that if the $i$ derivative is evaluated analytically and the $j$ derivative with finite differences, the Hessian matrix will have numerical noise in column indices but not in row indices, and thereby will be asymmetric, even though it is supposed to be symmetric. This is why it is a good idea to re-symmetrize the Hessian, $$ H_{ij} \to \frac {H_{ij} + H_{ji}} 2 $$.

Fully numerical scheme

What about the fully numerical Hessian? It is nothing more complicated than the above; you just use the finite difference equation again to expand $g_i$ to get

$$ H_{ij} \approx \frac {E({\bf x}+h\hat{e}_i+h\hat{e}_j)-E({\bf x}-h\hat{e}_i+h\hat{e}_j)-E({\bf x}+h\hat{e}_i-h\hat{e}_j)+E({\bf x}-h\hat{e}_i-h\hat{e}_j)} {4h^2}. $$

Even though this example only used two-point rules, you can use higher-order ones as well, such as the ones presented in Wikipedia.

  • $\begingroup$ Can I contact you in chat? I have few more questions. $\endgroup$ Feb 5, 2023 at 9:54
  • $\begingroup$ @HemanthHaridas I've copied over the discussion to a chat room. $\endgroup$
    – Tyberius
    Feb 5, 2023 at 13:50

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .