Derivation of the GHF SCF equations through Lagrange multipliers

Using an expansion of the spinors in known scalar bases for the components, the energy can be seen as a function of the expansion coefficients: \[\begin{equation} \require{physics} E(\vb{C}_I) = \sum_I^N h_{II}(\vb{C}_I) + \frac{1}{2} \sum_{IJ}^N (g_{IIJJ}(\vb{C}_I) - g_{IJJI}(\vb{C})_I) \thinspace . \end{equation}\] Upon close inspection, we see that only the coefficients that correspond to the occupied spinors (i.e. with indices \(I\)) contribute to the GHF energy. The optimizable parameters are thus the coefficients \(\vb{C}_I\), which is short-hand for the rectangular matrix that collects the columns of the first \(I\) occupied spinor expansion coefficients.

If we were to minimize this energy with respect to the expansion coefficients, nothing restricts us from finding expansion coefficients that do not represent an orthonormal spinor basis. Therefore, we should seek to optimize the Lagrangian \[\begin{equation} \mathscr{L} ( \vb{C}_I, \boldsymbol{\epsilon} ) = E(\vb{C}_I) - \sum_{IJ}^N \epsilon_{JI} \qty[ \vb{C}^\dagger \vb{S} \vb{C} - \vb{I} ]_{IJ} \thinspace , \end{equation}\] whose goal is to minimize the energy subject to the constraint that the orbitals remain orthonormal. The reason for the switched indices for the Lagrangian multipliers \(\epsilon_{JI}\) is such that the resulting equations can be represented by matrix multiplications. We now expand the spinor matrix elements in those of the underlying scalar bases, and we require the derivatives \[\begin{equation} \eval{ \pdv{ \mathscr{L} ( \vb{C}_I, \boldsymbol{\epsilon} ) }{ C^{\varepsilon *}_{\gamma L} } }_{\vb{C}_I^\star} = 0 \end{equation}\] to vanish at the optimal coefficients \(\vb{C}_I^\star\), in order to find equations that optimize the Lagrangian. We should note that \(L\) is again an occupied index, so \(C^{\varepsilon}_{\gamma L}\) represents one of the optimizable parameters. We therefore call these equations the stationary conditions. (We do not bother with the other derivative, which is its complex conjugate.) As an itermediary result, we find: \[\begin{multline} \sum_\sigma \sum_\mu^{K_\sigma} h^{\varepsilon \sigma}_{\gamma \mu} C^{\sigma \star}_{\mu L} + \sum_\sigma \sum_{\mu \nu}^{K_\sigma} \sum_{\rho}^{K_\varepsilon} P^{\sigma \sigma}_{\mu \nu}(\vb{C}_I^\star) C^{\varepsilon \star}_{\rho L} (\gamma \varepsilon \rho \varepsilon | \mu \sigma \nu \sigma) \\ - \sum_\sigma \sum_{\mu \nu}^{K_\sigma} \sum_{\rho}^{K_\varepsilon} P^{\sigma \varepsilon}_{\mu \rho}(\vb{C}_I^\star) C^{\sigma \star}_{\nu L} (\gamma \varepsilon \rho \varepsilon | \mu \sigma \nu \sigma) = \sum_{I}^N \epsilon_{IL} \sum_{\mu}^{K_\varepsilon} S^{\varepsilon \varepsilon}_{\gamma \mu} C^{\varepsilon \star}_{\mu I} \thinspace , \end{multline}\] by introducing the scalar basis GHF density matrix \(\vb{P}\): \[\begin{equation} P^{\sigma \tau}_{\mu \nu}(\vb{C}_I) = \sum_{I}^N C^{\sigma *}_{\mu I} C^{\tau}_{\nu I} \thinspace . \end{equation}\] We have already mentioned that our penultimate goal would be to find a matrix equation. In order to do so, we define the so-called direct, or Coulomb matrix \(\vb{J}\) \[\begin{equation} \vb{J}(\vb{C}_I) = \begin{pmatrix} \vb{J}^{\alpha \alpha}(\vb{C}_I) & \vb{0} \\ \vb{0} & \vb{J}^{\beta \beta}(\vb{C}_I) \end{pmatrix} \thinspace , \end{equation}\] with elements \[\begin{equation} J^{\sigma \sigma}_{\mu \nu}(\vb{C}_I) = \sum_\tau \sum_{\rho \lambda}^{K_\tau} P^{\tau \tau}_{\rho \lambda}(\vb{C}_I) ( \mu \sigma \nu \sigma | \rho \tau \lambda \tau ) , \end{equation}\] and the exchange matrix \(\vb{K}\) \[\begin{equation} \vb{K}(\vb{C}_I) = \begin{pmatrix} \vb{K}^{\alpha \alpha}(\vb{C}_I) & \vb{K}^{\alpha \beta}(\vb{C}_I) \\ \vb{K}^{\beta \alpha}(\vb{C}_I) & \vb{K}^{\beta \beta}(\vb{C}_I) \end{pmatrix} \thinspace , \end{equation}\] with elements \[\begin{equation} K^{\sigma \tau}_{\mu \nu}(\vb{C}_I) = \sum_{\rho}^{K_\sigma} \sum_{\lambda}^{K_\tau} P^{\tau \sigma}_{\lambda \rho}(\vb{C}_I) ( \mu \sigma \rho \sigma | \lambda \tau \nu \tau ) \thinspace . \end{equation}\] Using these definitions, we can rewrite the stationary equations as: \[\begin{equation} \sum_\sigma \sum_\mu^{K_\sigma} h^{\varepsilon \sigma}_{\gamma \mu} C^{\sigma \star}_{\mu L} + \sum_{\rho}^{K_\varepsilon} J^{\varepsilon \varepsilon}_{\gamma \rho}(\vb{C}_I^\star) C^{\varepsilon \star}_{\rho L} - \sum_{\sigma} \sum_{\nu}^{K_\sigma} K^{\varepsilon \sigma}_{\gamma \nu}(\vb{C}_I^\star) C^{\sigma \star}_{\nu L} = \sum_I^N \epsilon_{IL} \sum_{\mu}^{K_\varepsilon} S^{\varepsilon \varepsilon}_{\gamma \mu} C^{\varepsilon \star}_{\mu I} \thinspace , \end{equation}\] which, upon closer inspection, reveals a matrix multiplication. If we now define the GHF Fock matrix as \[\begin{equation} \vb{F}(\vb{C}_I) = \vb{H} + \vb{J}(\vb{C}_I) - \vb{K}(\vb{C}_I) \thinspace , \end{equation}\] the stationary equations can be written as \[\begin{equation} \vb{F}(\vb{C}_I^\star) \vb{C}_I^\star % = \vb{S} \vb{C}_I^\star \boldsymbol{\epsilon} % \thinspace , \end{equation}\] which are the so-called (SCF) equations. Concerning the dimensions of this generalized eigenvalue problem, the Fock and overlap matrix are of dimension \((M \times M)\), and the occupied spinor coefficients \(\vb{C}_I\) has dimension \((M \times N)\). The optimized occupied spinor coefficients \(\vb{C}^\star_I\) thus are the lowest \(N\) generalized eigenvectors of the Fock matrix (expressed in the underlying scalar orbital bases).

We should note that, in general, the matrix of the Lagrangian multipliers \(\boldsymbol{\epsilon}\) is not required to be diagonal in order to find a Hartree-Fock solution. We can, however, diagonalize the matrix of Lagrangian multipliers, which ultimately would lead to the spinors, i.e. those for which the GHF Fock matrix (in the orthonormal spinor basis) is diagonal.