Why is CPHF/CPKS necessary for calculating second derivatives?

Question

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?

Answer 1

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.

Answer 2

In order to clear up the confusion, there are in fact two questions to be answered:

1. What are actually computed in HF energy gradients?
2. What does CPHF theory solve?

Let us begin with direct differentiation of HF energy with respect to an arbitrary nuclear coordinate, $\mathbf{R}$, \begin{equation} \label{direct differentiation of HF energy} \frac{\partial E^{\mathrm{HF}}}{\partial \mathbf{R}} = \sum_{\mu \nu} \frac{\partial P_{\nu \mu}}{\partial \mathbf{R}} F_{\mu \nu} + P_{\nu \mu} \frac{\partial F_{\mu \nu}}{\partial \mathbf{R}} + \frac{\partial V_{\mathrm{nuc}}}{\partial \mathbf{R}}, \end{equation} where $F_{\mu \nu}$ and $P_{\nu \mu}$ are Fock matrix and spin-orbital density matrix, respectively. Generally, explicit evaluations of $\partial P_{\nu \mu}/ \partial \mathbf{R}$ is not easy and can be avoided if HF energy gradient is all one desires. What do I mean by that? If we recall the minimisation of HF energy with respect to coefficients, $C$ (subject to the orthonomality conditions, $\sum_{\mu \nu} C^{*}_{\mu p} S_{\mu \nu} C_{\nu q} = \delta_{pq}$) leads to the Fock-type equations \begin{equation} \sum_{\nu} (F_{\mu \nu} - \epsilon_{p} S_{\mu \nu})C_{\nu p} = 0. \end{equation} Then the first term of \ref{direct differentiation of HF energy} can be expanded as \begin{equation} \label{rewrite dPdR} \begin{aligned} \sum_{\mu \nu} \frac{\partial P_{\nu \mu}}{\partial \mathbf{R}} F_{\mu \nu} &= \sum_{\mu \nu} \sum_{i} \frac{\partial C^{*}_{\mu i}}{\partial \mathbf{R}} F_{\mu \nu} C_{\nu i} + \mathrm{CC} \\ &= \sum_{\mu \nu} \sum_{i} \frac{\partial C^{*}_{\mu i}}{\partial \mathbf{R}} \epsilon_{i} S_{\mu \nu} C_{\nu i} + \mathrm{CC}, \end{aligned} \end{equation} where $\mathrm{CC}$ stands for the complex conjugate form of the previous term. Furthermore, differentiation of orthonomality condition (with $p = q = i$) results in \begin{equation}\label{differentiation of orthonomality condition} \sum_{\mu \nu} \left[ \frac{\partial C^{*}_{\mu i}}{\partial \mathbf{R}} S_{\mu \nu} C_{\nu i} + C^{*}_{\mu i} \frac{\partial S_{\mu \nu}}{\partial \mathbf{R}} C_{\nu i} + C^{*}_{\mu i} S_{\mu \nu} \frac{\partial C_{\nu i}}{\partial \mathbf{R}} \right] = 0. \end{equation} We may now use above equation to eliminate the coefficient derivatives and therefore we arrive at the final expression, \begin{equation} \frac{\partial E^{\mathrm{HF}}}{\partial \mathbf{R}} = \sum_{\mu \nu} P_{\nu \mu} \frac{\partial F_{\mu \nu}}{\partial \mathbf{R}} + \frac{\partial V_{\mathrm{nuc}}}{\partial \mathbf{R}} - \sum_{\mu \nu} W_{\nu \mu} \frac{\partial S_{\mu \nu}}{\partial \mathbf{R}}, \end{equation} where $W_{\nu \mu} = \sum_{i} \epsilon_i C^{*}_{\mu i} C_{\nu i}$ is an energy-weighted density matrix. Above equation is the answer to the first question, it is much cheaper than directly computing $\partial P /\partial \mathbf{R}$.

We may obtain the second order derivative of the HF energy as \begin{equation} \begin{aligned} \frac{\partial^2 E^{\mathrm{HF}}}{\partial \mathbf{R}_1 \partial \mathbf{R}_2} =& \sum_{\mu \nu} P_{\nu \mu} \frac{\partial^2 F_{\mu \nu}}{\partial \mathbf{R}_1 \partial \mathbf{R}_2} + \frac{\partial^2 V^{\mathrm{nuc}}}{\partial \mathbf{R}_1 \partial \mathbf{R}_2} - \sum_{\mu \nu} W_{\nu \mu} \frac{\partial^2 S_{\mu \nu}}{\partial \mathbf{R}_1 \partial \mathbf{R}_2} \\ &+ \sum_{\mu \nu} \frac{\partial P_{\nu \mu}}{\partial \mathbf{R}_2} \frac{\partial F_{\mu \nu}}{\partial \mathbf{R}_1} - \sum_{\mu \nu} \frac{ \partial W_{\nu \mu}}{\partial \mathbf{R}_2} \frac{\partial S_{\mu \nu}}{\partial \mathbf{R}_1}. \end{aligned} \end{equation} Here the terms involving computations of $\partial P_{\nu \mu} / \partial \mathbf{R}$ and $\partial W_{\nu \mu} / \partial \mathbf{R}$ can no longer be avoided, and therefore it is where CPHF theory kicks in.

I will finish up the second question in detail when I have spare time, but here's the general idea of CPHF theory for you to chew on: In CPHF theory, the SCF equation is perturbed via a unitary transformation where its unitary matrix is unknown and has to be solved, then the perturbed equation is approximated via Taylor expansion, the perturbed Fock matrix is then "somehow" a linear function of the derivative of perturbed unitary matrix, $U'(\mathbf{R})$, and $\partial P_{\nu \mu}/ \partial \mathbf{R} = U^{'}_{qp}(\mathbf{R}) + U^{'*}_{qp}(\mathbf{R})$.