Coupled-perturbed Hartree-Fock equations
The coupled-perturbed Hartree-Fock (CPHF) equations are used to solve the response of molecular orbitals to each external perturbation (the \(U^x\) matrices below). Before working on the derivations, it is worth noting that the CPHF equations are not explicitly solved for first-order energy derivatives (i.e., gradients) in quantum chemistry. Please see the comments at the end of this note for details.
The Hartree-Fock (or density functional theory) eigenvalue-equations remain true in the presence of external perturbations \(x\), i.e.,
\[F^{AO}(x)C(x)=S^{AO}(x)C(x)E(x)\]The unperturbed coefficients \(C(0)\) represent the unperturbed molecular orbitals \(\phi_p(0)=\sum_\mu C_{\mu p}(0)f_\mu\), where \(f_\mu\) are the atomic orbitals. With perturbation, the molecular orbitals can be expanded in the basis of unperturbed molecular orbitals
\[\phi_p(x)=\sum_{q}U_{qp}(x)\phi_q(0)=\sum_{\mu}\sum_{q}U_{qp}(x)C_{\mu q}(0)f_\mu\]It is obvious that the perturbed coefficients \(C(x)\) can be written as
\[C(x)=C(0)U(x)\]With defining \(F(x)\equiv C^\dagger(0)F^{AO}C(0)\) and \(S(x)\equiv C^\dagger(0)S^{AO}C(0)\), we have \(F(0)=E(0)\), \(U(0)=S(0)=I\), and
\[F(x)U(x)=S(x)U(x)E(x)\]as well as the orthonormality condition
\[U^\dagger(x)S(x)U(x)=I\]Evaluate the derivatives of the above equations with respect to \(x\) at \(x=0\), use the notation \(O^x\equiv\frac{\partial O}{\partial x}\), and we have
\[F^x(0)+E(0)U^x(0)=S^x(0)E(0) + U^x(0)E(0) + E^x(0)\] \[U^{\dagger x}(0) + U^x(0) = - S^x(0)\]Ommiting the \(0\), we have the following equations
\[F^x_{pq} +E_pU^x_{pq} = S^x_{pq}E_q + U^x_{pq}E_q + E^x_{p}\delta_{pq}\] \[U^{x \ast}_{qp} + U^x_{pq} = - S^x_{pq}\]for \(p=q\)
\[Re[U^x_{pp}] = -\frac{1}{2}S^x_{pp}\]for \(p\neq q\)
\[U^x_{pq} = \frac{F^x_{pq}-S^x_{pq}E_q}{E_q-E_p}\]To solve \(U^x\), we have to know \(F^x\) and \(S^x\). The latter is easier and can be obtained as \(S^x(0)=C^\dagger(0) S^{AO,x} C(0)\). For Fock operator, we have
\[\begin{gathered} F_{pq}(x)=\sum_{\mu\nu}C_{\mu p}^{*}(0)C_{\nu q}(0)\left[h_{\mu\nu}(x)+\sum_{i\rho\sigma}C_{\rho i}^{*}(x)C_{\sigma i}(x)\langle\mu\rho||\nu\sigma\rangle(x)\right] \\ =\sum_{\mu\nu}C_{\mu p}^{*}(0)C_{\nu q}(0)\left[h_{\mu\nu}(x)+\sum_{i\rho\sigma}\sum_{ts}C_{\rho t}^{*}(0)U_{ti}^{*}(x) C_{\sigma s}(0)U_{si}(x)\langle\mu\rho||\nu\sigma\rangle(x)\right] \\ =h_{pq}(x)+\sum_{i}\sum_{ts}U_{ti}^{*}(x)U_{si}(x)\langle pt||qs\rangle(x) \end{gathered}\] \[\begin{gathered} F_{pq}^{x}=h_{pq}^{x}(0)+\sum_{i}\left[\langle pi||qi\rangle^{x}(0)+\sum_{t}U_{ti}^{x*}(0)\langle pt||qi\rangle(0)+\sum_{s}U_{si}^{x}(0)\langle pi||qs\rangle(0)\right] \\ =f_{pq}^{x}+\sum_{it}[U_{ti}^{x*}\langle pt||qi\rangle+U_{ti}^{x}\langle pi||qt\rangle] \\ =f_{pq}^{x}+\sum_{ij}\bigl[U_{ji}^{x*}\langle pj||qi\rangle+U_{ji}^{x}\langle pi||qj\rangle\bigr]+\sum_{ia}\bigl[U_{ai}^{x*}\langle pa||qi\rangle+U_{ai}^{x}\langle pi||qa\rangle\bigr] \\ =f_{pq}^{x}-\sum_{ij}S_{ij}^{x}\langle pj||qi\rangle+\sum_{ia}[U_{ai}^{x*}\langle pa||qi\rangle+U_{ai}^{x}\langle pi||qa\rangle] \end{gathered}\]As a result, we have
\[U_{pq}^x=\frac{G_{pq}+\sum_{ia}[U_{ai}^{x*}\langle pa||qi\rangle+U_{ai}^x\langle pi||qa\rangle]}{E_q-E_p}\]where
\[G_{pq}=f_{pq}^x-S_{pq}^xE_q-\sum_{ij}S_{ij}^x\langle pj||qi\rangle\]We note that these equations are coupled only for \(U^x_{ai}\) block. After solving \(U^x_{ai}\), all the other blocks can be calculated without any iteration. The CPHF equations for \(U^x_{ai}\) can be reformulated as
\[\begin{bmatrix}\Delta+A&B\\B^*&\Delta+A^*\end{bmatrix}\begin{bmatrix}U^x\\U^{x*}\end{bmatrix}=-\begin{bmatrix}G^x\\G^{x*}\end{bmatrix}\]where
\[\Delta_{bj,ai}=\bigl(E_{b}-E_{j}\bigr)\delta_{ba}\delta_{ji}\] \[A_{bj,ai}=\langle bi||ja\rangle, B_{bj,ai}=\langle ba||ji\rangle\]One can see that the CPHF equation has to be solved for each perturbation (or each \(x\)) independently and this would be computationally expensive.
But as we mentioned before, the CPHF equations are not explicitly solved for first-order energy derivatives. This is consistent with the famous “2n+1” rule that to calculate the 2n+1 order property, one only needs to solve the n order response of the wave function. In practice, one can show that to calculate the first-order energy derivatives, one only needs to solve the contributions with the following form
\[\sum_{ia}\left(U_{ai}^{x}I_{ai}^\ast+U_{ai}^{x\ast}I_{ai}\right)\]It can be reformulated as
\[\left[\begin{array}{cc}I^*&I\end{array}\right]\left[\begin{array}{c}U^x\\U^{x*}\end{array}\right]=-\left[\begin{array}{cc}I^*&I\end{array}\right]\left[\begin{array}{cc}\Delta+A&B\\B^*&\Delta+A^*\end{array}\right]^{-1}\left[\begin{array}{c}G^x\\G^{x*}\end{array}\right]\]If we can solve the following equation to obtain \(D\) intermediates
\[\left[\begin{array}{cc}D^* & D\end{array}\right]\left[\begin{array}{c}\Delta+A&B\\B^*&\Delta+A^*\end{array}\right]=-\left[\begin{array}{cc}I^* & I\end{array}\right]\]the original contributions can be written as
\[\sum_{ia}\left(U_{ai}^{x}I_{ai}^\ast+U_{ai}^{x\ast}I_{ai}\right) = \sum_{ia}\left(G_{ai}^{x}D_{ai}^\ast+G_{ai}^{x\ast}D_{ai}\right)\]This is so-called “Z-vector” method in quantum chemistry. We emphasize that the Z-vector equations to solve \(D\) intermediates are perturbation-independent. It only needs to be solved once and can be used to calculate any first-order energy derivatives.