Why is CPHF/CPKS necessary for calculating second derivatives?

This question is coming from an answer to one of my previous questions. During optimizations, QM programs usually compute the gradient(first derivative) analytically, and take a guess of the hessian (second derivatives). If the hessian is needed a coupled-perturbed hartree fock (CPHF) or coupled-perturbed Kohn-Sham (CPKS) is usually required which is very computationally expensive. From what I have been told, the gradient does not need CPHF, it can be calculated directly from the SCF.

My question is why the second derivative is so much more expensive than the first derivative? For a compound I was working on, the first derivatives took about 2 minutes, while the second derivatives took almost 15 minutes to run. That's more than 7 times! What I don't understand is that if the SCF solution can be differentiated once w.r.t the coordinates, then why can't it be differentiated twice?

It comes down to the fact that HF and KS both are variational methods. This short article by Julien Toulouse gives a great description of ways to compute static/dynamic response properties. I'll just summarize the relevant portion.

We can compute derivatives of the energy with respect to any variable $$x$$ as: $$\frac{dE}{dx}=\frac{\partial E}{\partial x}+\sum_i \frac{\partial E}{\partial p_i}|_{\mathbf{p}=\mathbf{p}^0} \frac{\partial p_i^0}{\partial x}$$ Here we are writing the derivative in two terms. The first is due to the explicit dependence of the energy on the variable $$x$$. The latter term is due to implicit dependence, with the energy depending on particular wavefunction parameters $$\mathbf{p}$$, which in turn may depend on $$x$$. For SCF methods, these parameters are just the MO coefficients $$C$$.

For a general method, this would require some type of response calculation to solve, as we typically don't have an explicit form for $$\frac{\partial p_i^0}{\partial x}$$. However, since the energy for HF/KS is variational $$\frac{\partial E}{\partial p_i}|_{\mathbf{p}=\mathbf{p}^0}=0$$, which zeros out this term.

So to compute the forces with HF/KS, we only need to consider the explicit dependence of the energy on the nuclear positions. However, once we want to compute the Hessian, we can no longer ignore this implicit term. If we write the Hessian as a derivative of the force in the same form as above, the force is not variationally optimal and so this terms doesn't cancel. Thus we need to determine the dependence of the MO coefficients on the nuclear positions, which is typically done using CPHF/CPKS.

You can get a rough order of magnitude estimate of a 2nd order property from just the explicit term. For example, when computing the polarizability (2nd derivative of the energy with respect to an applied electric field) Gaussian will print out an approximate polarizability, which it calculates by just contracting the dipole with itself rather than the perturbed electric density. While this can sometimes be close to the final result, I'm not aware of any formal bound on the size of the implicit contribution, so in general it would be a major approximation to neglect the perturbation of the density.

• What if we ignore the term coming from the MO coefficients? I guess what I am trying to understand is which term makes the biggest contribution to the second derivative? Nov 26 '20 at 21:51
• @ShoubhikRMaiti I added a bit to address your comment. My experience with this is mainly with computing polarizabilities, so it may differ slightly for geometric derivatives.
– Tyberius
Nov 27 '20 at 0:24