# 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? Commented Nov 26, 2020 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
Commented Nov 27, 2020 at 0:24

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}$$, $$$$\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}},$$$$ 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 $$$$\sum_{\nu} (F_{\mu \nu} - \epsilon_{p} S_{\mu \nu})C_{\nu p} = 0.$$$$ Then the first term of \ref{direct differentiation of HF energy} can be expanded as \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} where $$\mathrm{CC}$$ stands for the complex conjugate form of the previous term. Furthermore, differentiation of orthonomality condition (with $$p = q = i$$) results in $$$$\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.$$$$ We may now use above equation to eliminate the coefficient derivatives and therefore we arrive at the final expression, $$$$\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}},$$$$ 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{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} 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})$$.

• Welcome to our community! We hope to see much more of you in the future!!! Would you like us to convert this to a comment? If not, why don't you expand this into a full answer? Commented Apr 22 at 15:40
• The $O(N^5)$ integral transformation terms only show up in the MO algorithm. In the AO algorithm (which is used by most programs by default, at least when the molecule is large), the program never converts the two-electron integrals to the MO basis even once. Rather, the CPHF equations are solved iteratively by Krylov subspace methods without explicitly forming the LHS matrix. Each iteration still costs $O(N^5)$, but that's only because there are $O(N)$ right hand side vectors. If there are only $O(1)$ RHS vectors (as in polarizability/TDDFT gradient) the cost per iteration would be $O(N^4)$. Commented Apr 23 at 11:24