Analytic Hessians for meta-GGA functionals

Question

In many of the free/open-source QM programs like GAMESS, ORCA, NWChem, the calculation of analytic hessians is not possible when a meta-GGA density functional is used, like M06-2X. My first question is — why is this? Why do meta-GGA functionals not have analytic hessians?

However, when I used Gaussian 16, I have found that it can calculate analytic hessians for M06-2X as well. My first guess was that Gaussian was automatically shifting to numerical calculation, but after reading the output file, it seems that Gaussian is performing CPHF on M06-2X (or some variation of CPHF, I am new to all of this). This should not be possible. So how does Gaussian do it?

Does it have something to with the coding of the programs or the theory of the m-GGA method itself?

Edit: At the time I wrote the question, Orca 4.2.1 did not have analytic hessians for m-GGA functionals. However, the most recent update, Orca 5.0, made m-GGA analytic hessians available.

Answer 1

To clarify the statements at the beginning of the question, at variance to NWChem which truly is free software,

• although the source code for GAMESS is available if you sign the license agreement, it is not open-source or free software, while
• ORCA is neither open source nor free software (it is available at no cost for academic users)

As to the claim about the lack of analytic Hessians for meta-GGAs in free software, this is true: in addition to NWChem, analytic Hessians for meta-GGAs are also missing in PySCF and Psi4.

There's no a priori reason for this: there's nothing that prevents analytic Hessians in Gaussian-basis codes, since it is perfectly possible to compute gradients and Hessians from the variational energy functional expressed in the atomic-orbital basis, see e.g. Pople et al in Chem. Phys. Lett. 199, 557 (1992). The key point here is that since the orbitals are expanded in terms of analytical functions, we can compute whatever derivatives we want!

If you do the calculation, it turns out that gradients can be computed just from the knowledge of the self-consistent field wave function just like in Hartree-Fock, whereas to compute Hessians you also need to solve for the perturbation in the wave function with the coupled-perturbed approach (this is what you see in Gaussian).

Now, the reason why meta-GGA Hessians are not available in most codes is that the implementation gets quite hairy. If you have an LDA, the inputs are the spin densities, $n_\uparrow$ and $n_\downarrow$. For GGAs, you add the reduced gradients $\gamma_{\sigma \sigma'} = \nabla n_\sigma \cdot \nabla n_{\sigma'}$; that is, $\gamma_{\uparrow \uparrow}$, $\gamma_{\uparrow \downarrow} = \gamma_{\downarrow \uparrow}$ and $\gamma_{\downarrow \downarrow}$, and for meta-GGAs also $\tau_\uparrow$ and $\tau_\downarrow$. This means that you have 2 input variables for LDAs, 5 for GGAs and 7 for meta-GGAs, meaning you have 2, 5, and 7 first derivatives for LDAs, GGAs, and meta-GGAs, respectively, which will be contracted with something that looks like the electron density to get the nuclear gradient.

For the Hessian, you get more terms. For meta-GGAs you can pick the first index in 7 ways, and the second index in 6 ways; eliminating the permutation of the two gives 21 second derivatives of the exchange-correlation function. (Compare this to GGAs that only have 10 terms, or LDAs that only have 1 term!) But, this is not the whole story: the above pertains only to the change in the exchange-correlation energy; you also get the coupling term of the density response (i.e. first derivative of electron density with respect to perturbation) with the gradient of the exchange-correlation energy; this is not symmetric so it would appear to add 7*7=49 more terms, bringing the total up to 70 terms.

The implementation of analytical Hessians is painstaking work, and since there's no new science in it, it just hasn't been a priority in free software programs. Analytical Hessians are available in commercial codes like Gaussian and Q-Chem, because they are needed for many routine applications in industry.

Things might change in the future, though: although laborious, the implementation of the analytic Hessians is straightforward. It's just a question of someone dedicating a few weeks of their time to write down the equations, implement them, and check that they are correct by finite differences...

Answer 2

Meta-GGA density functionals are dependent on the kinetic energy density $$ \tau_\sigma (\mathbf{r}) = \sum_i^\text{occ} \frac{1}{2} \lvert \nabla \phi_{\sigma,i} (\mathbf{r}) \rvert^2 $$ where $\phi_{\sigma,i} (\mathbf{r})$ are Kohn-Sham orbitals.$^1$ Certain spatial derivatives of this quantity will be necessary for the calculation of the Hessian and I have been told that achieving the necessary accuracy in the grid-based calculation is non-trivial.

Personally, I believe that this is down to most people assuming that little value is gained by mGGA frequencies over GGA frequencies. Of course, one should ideally use the Hessian to verify that one has reached a local minimum in one's geometry optimization, but beyond this, the usefulness may be limited. Therefore, development resources are directed elsewhere.

---

$^1$ Piotr de Silva, Clémence Corminboeuf, "Communication: A new class of non-empirical explicit density functionals on the third rung of Jacob’s ladder," J. Chem. Phys. 2015, 143, 111105 (PDF), doi:10.1063/1.4931628. Note that some authors use a different prefactor.