Restricted Hartree-Fock theory

The RHF (restricted Hartree-Fock) wave function model is given by \[\begin{equation} \require{physics} \ket{\text{RHF}(\boldsymbol{\kappa})} = \exp(-\hat{\kappa}) \ket{\text{RHF}} \thinspace , \end{equation}\] where \(\ket{\text{RHF}}\) is a closed-shell determinant: \[\begin{align} \require{physics} \ket{\text{RHF}} & = \qty( \prod_i^{N_P} \hat{a}^\dagger_{i \alpha} \hat{a}^\dagger_{i \beta} ) \ket{\text{vac}} \\ &= \hat{A}_{\text{core}} \ket{\text{vac}} \thinspace , \end{align}\] i.e. a single Slater determinant built from real orthonormal spatial orbitals. Written with spin-adapted creation operators, the single Slater determinant may also be written as \[\begin{equation} \ket{\text{RHF}} = \frac{1}{ \sqrt{2^{N_P}} } \qty( \prod_i \hat{Q}^{0,0}_{ii} ) \ket{\text{vac}} \thinspace . \end{equation}\] For RHF theory, \(\hat{\kappa}\) is the generator of singlet occupied-virtual orbital rotations: \[\begin{equation} \hat{\kappa} = \sum_{ai} {}^{\text{R}} \kappa_{ai} \hat{E}^-_{ai} + i \sum_{ai} {}^{\text{I}} \kappa_{ai} \hat{E}^+_{ai} \thinspace , \end{equation}\] where the only non-redundant parameters are the occupied-virtual orbital generators. When working only with real orbitals, we may set \({}^{\text{I}} \kappa_{ai}\) to zero.

The expressions for the RHF density matrices are simple. The 1-DM is diagonal and only non-zero for occupied orbitals \(ij\): \[\begin{equation} D_{ij} = \ev{\hat{E}_{ij}}{\text{RHF}} = 2 \delta_{ij} \thinspace , \end{equation}\] while the 2-DM can be written as: \[\begin{equation} d_{ijkl} = \ev{\hat{e}_{ijkl}}{\text{RHF}} = 4 \delta_{ij} \delta_{kl} - 2 \delta_{il} \delta_{kj} \thinspace , \end{equation}\] such that we can write the RHF electronic energy as \[\begin{equation} E_{\text{RHF}} = 2 \sum_i h_{ii} + \sum_{ij} (2 g_{iijj} - g_{ijji}) \thinspace , \end{equation}\] in which we remember that \(i\) and \(j\) are indices of the doubly occupied orbitals.

For the RHF wave function model, the Fockian and super-Fockian matrix elements simplify considerably. For the Fockian matrix, we find: \[\begin{equation} \mathscr{F}_{ip} = 2 {}^{i} F_{pi} = 2 \qty[ h_{pi} + \sum_j ( 2 g_{pijj} - g_{pjji} ) ] \thinspace , \end{equation}\] where its first index must be an occupied index \(i\) and \({}^{i} F_{pi}\) is the inactive Fock matrix. The \(\mathscr{Z}_{pqrs}\)-tensor must sum only over occupied indices, and the first index must be occupied: \[\begin{equation} \mathscr{Z}_{ipqr}( \hat{g} ) = \sum_{jk} ( g_{pjrk} d_{ijqk} - g_{pqjk} d_{irjk} - g_{pjkq} d_{ijkr} ) \thinspace . \end{equation}\] In the RHF orbital Hessian elements, we will encounter these non-zero objects in two appearances. The first one is when the third index is occupied: \[\begin{equation} \mathscr{Z}_{iajb}( \hat{g} ) = 4 g_{aibj} - 2 g_{ajbi} \thinspace , \end{equation}\] and the second type is when the fourth index is occupied: \[\begin{equation} \mathscr{Z}_{iabj}( \hat{g} ) = -\delta_{ij} \sum_k ( 4 g_{abkk} - 2 g_{akkb} ) - ( 4 g_{aijb} - 2 g_{abji} ) \thinspace . \end{equation}\]

The parts of the RHF orbital gradient can be calculated as: \[\begin{align} {}^{\text{R}} \vb{E}^{(1)}_{ai} &= \ev{ \comm{\hat{E}^-_{ai}}{\hat{\mathcal{H}}} }{\text{core}} \\ &= - \qty[ \mathscr{F}_{ia}(\hat{\mathcal{H}}) + \mathscr{F}^*_{ia}(\hat{\mathcal{H}}) ] \\ &= -2 ( {}^{i} F_{ai} + {}^{i} F_{ia} ) \end{align}\] and \[\begin{align} {}^{\text{I}} \vb{E}^{(1)}_{ai} &= i \ev{ \comm{\hat{E}^+_{ai}}{\hat{\mathcal{H}}} }{\text{core}} \\ &= i \qty[ \mathscr{F}_{ia}(\hat{\mathcal{H}}) - \mathscr{F}^*_{ia}(\hat{\mathcal{H}}) ] \\ &= 2 i ( {}^{i} F_{ai} - {}^{i} F_{ia} ) \thinspace . \end{align}\] From the real part of the orbital gradient follows that, when the gradient is made to vanish, the inactive Fock matrix should be diagonal: \[\begin{equation} {}^{i} F_{pq} = \delta_{pq} \epsilon_p \thinspace , \end{equation}\] where \(\epsilon_p\) are the diagonal elements of the inactive Fock matrix, i.e. the orbital energies.

Furthermore, we find for the parts of the RHF orbital Hessian: \[\begin{equation} {}^{\text{R} \text{R}} \vb{E}^{(2)}_{ai,bj} = 2 \qty[ \delta_{ij} ( {}^{i} F_{ab} + {}^{i} F_{ba} ) - \delta_{ab} ( {}^{i} F_{ij} + {}^{i} F_{ji} ) + (1 + P_{ai})(1 + P_{bj}) ( 2 g_{aibj} - g_{ajbi} ) ] \thinspace , \end{equation}\] \[\begin{equation} {}^{\text{R} \text{I}} \vb{E}^{(2)}_{ai,bj} = -2i (1 - P_{ij,ab}) (2 g_{aibj} - g_{ajbi}) \end{equation}\] and \[\begin{equation} {}^{\text{I} \text{I}} \vb{E}^{(2)}_{ai,bj} = 2 \qty[ \delta_{ij} ( {}^{i} F_{ab} + {}^{i} F_{ba} ) - \delta_{ab} ( {}^{i} F_{ij} + {}^{i} F_{ji} ) - (1 - P_{ai})(1 - P_{bj}) ( 2 g_{aibj} - g_{ajbi} ) ] \thinspace . \end{equation}\]

The partial derivatives of the RHF \(1\)-DM are: \[\begin{equation} \eval{ \pdv{ D_{pq}(\boldsymbol{\kappa}) }{ {}^{\text{R}} \kappa_{ai} } }_{\boldsymbol{\kappa}_0} = -2 ( \delta_{aq} \delta_{ip} + \delta_{ap} \delta_{iq} ) \end{equation}\] and \[\begin{equation} \eval{ \pdv{ D_{pq}(\boldsymbol{\kappa}) }{ {}^{\text{I}} \kappa_{ai} } }_{\boldsymbol{\kappa}_0} = 2i ( \delta_{ap} \delta_{iq} - \delta_{aq} \delta_{ip} ) \thinspace . \end{equation}\]

The density can be written as: \[\begin{equation} \rho(\vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0})) = 2 \sum_i \phi^*_i(\vb{r}) \phi_i(\vb{r}) \thinspace . \end{equation}\]