Solution of the Hartree-Fock equations

We already took a look at the Hartree-Fock equations and introduced the Fock operator. This post can be seen as a continuation. We will now take a closer look at how actually to solve the Hartree-Fock equations starting from an Ansatz in which we expand our molecular orbitals (MOs) in atomic orbitals (AOs).

We expand our molecular orbitals (MOs) in a basis of atomic orbitals (AOs)

\begin{equation} \psi_i = \sum_{\nu} |{\nu}\rangle c_{\nu i} \quad\quad\text{with}\quad |{\nu}\rangle = \chi_{\nu}(\vec{r}) \cdot \Big\{ \begin{array}{c} \alpha \\ \beta \\ \end{array} \end{equation}
\begin{equation} \delta \psi_i = d\psi_i = \sum_{\nu} |{\nu}\rangle dc_{\nu i} \end{equation}

For the Hartree Fock equations in a basis set we get:

\begin{equation} \delta L = \big[ \sum_i \langle{\delta\psi_i}|\hat{F}|{\psi_i}\rangle - \sum_{ij} \varepsilon_{ij} \langle{\delta\psi_i|\psi_j}\rangle \big] + c.c. = 0 \end{equation}
\begin{equation} \frac{dL}{dc_{\mu i}} = \langle{\mu}| \hat{F} |{\nu}\rangle - \sum_j \varepsilon_{ij} \langle{\mu|\psi_j}\rangle = 0 \quad \forall_{\mu,i} \end{equation}

Introducing the matrices

\begin{equation} F_{\mu\nu} = \langle{\mu}|\hat{F} |{\nu}\rangle \end{equation}
\begin{equation} S_{\mu\nu} = \langle{\mu|\nu}\rangle \end{equation}

the Hartree Fock equations can be written in matrix form as

\begin{equation} \mathbf{F} \vec{c}_i = \sum_j \varepsilon_{ij} \mathbf{S} \vec{c}_j \end{equation}

In case of the canonical Hartree Fock equations

\begin{equation} \mathbf{F} \vec{c}_i = \varepsilon_{i} \mathbf{S} \vec{c}_i \end{equation}

The two most time consuming computational steps are

  • Calculation of the two electron AO integrals \(\big(\mu\nu|\kappa\lambda\big)\).
  • The construction of the Fock matrix using the two electron AO integrals \(\big(\mu\nu|\kappa\lambda\big)\) and the MO coefficients \(c_{\mu i}\)
    \begin{equation} F_{\mu\nu} = h_{\mu\nu} + J_{\mu\nu} - K_{\mu\nu} \end{equation}
    \begin{equation} J_{\mu\nu} = \langle{\mu}| \hat{J} |{\nu}\rangle = \sum_j \langle{\mu \psi_j|\nu\psi_j}\rangle = \sum_j \big(\mu\nu|\psi_j\psi_j \big) = \sum_{\kappa\lambda j} \big(\mu\nu|\kappa\lambda\big) c_{\lambda j} c_{\kappa j} \end{equation}
    \begin{equation} K_{\mu\nu} = \langle{\mu}| \hat{K} |{\nu}\rangle = \sum_j \langle{\mu \psi_j|\psi_j\nu}\rangle = \sum_j \big(\mu\psi_j|\psi_j\nu \big) = \sum_{\kappa\lambda j} \big(\mu\kappa|\lambda\nu\big) c_{\lambda j} c_{\kappa j} \end{equation}

Since the Fock matrix depends on the orbital coefficients \(c_{\mu i}\), which should be determined by solving the Hartree Fock equations, we obtain a non-linear problem, which needs to be solved iteratively. The standard approach is the self-consistent-field (SCF) method.

  • Guess start orbitals \(\psi_i^{(0)}\)
  • Calculate for \(\psi_i^{(m)}\) the matrices \(\mathbf{J}\) and \(\mathbf{K}\) and build \(\mathbf{F}\)
  • Solve \(\mathbf{F} \psi_i^{(m+1)}=\varepsilon_i^{(m+1)} \psi_i^{(m+1)}\) (diagonalize Fock matrix) to get a new set of orbitals.
  • Calculate the energy
    \begin{equation*} E_{\rm{HF}}^{(m+1)} = \sum_i \frac{1}{2} \left( h_{ii} + F_{ii} \right) \end{equation*}
  • Is \(|E_{\rm{HF}}^{(m+1)}-E_{\rm{HF}}^{(m)}| < T_1\) and \(||\psi^{(m+1)}-\psi^{(m)}||< T_2\), where \(T_1\) and \(T_2\) are convergence thresholds stop, else continue with 2.

Unrestricted and restricted Hartree-Fock

The Ansatz we used so far

  • Single Slater determinant
  • For all spin orbitals the spatial parts are optimized independent of each other. The only requirement is the ortohonormality of the spatial part and that the spin parts are eigenfunctions of \(\hat{S}_z\).
  • Thus the spatial part for \(\alpha\) and \(\beta\) can be different.

is denoted as unrestricted Hartree Fock (UHF). An other Ansatz is restricted Hartree Fock (RHF). It varies from UHF in two points.

  • For the orbitals the additional constraint is introduced that the spatial part for \(\alpha\) and \(\beta\) is identical (thus restricted).
    \begin{equation} \phi_i = \sum_{\nu} \chi_{\nu} c_{\nu i} \quad \phi^{\alpha}_i = \phi_i \alpha \quad \phi^{\beta}_i = \phi_i \beta \end{equation}
    Due to this restriction the parameters, which have to be optimized are reduced by one half, which reduces the computational costs.
  • The Ansatz for the wave function is one configuration state function (CSF). Only in the closed-shell and high-spin open-shell case this is possible with a single Slater determinant.

Using RHF the Fock matrix looks a bit different

\begin{equation} F_{\mu\nu} = h_{\mu\nu} + \sum_i \left[ 2\big(\mu\nu|ii\big) - \big(\mu i|i \nu\big) \right] \end{equation}

since due to spin integration we get different factors in front of the Coulomb and Exchange matrices. RHF is the preferred choice due to the reduced computational costs. Often it is an advantage to express the Fock matrix in terms of the density matrix

\begin{equation} D_{\mu\nu} = \sum_i c_{\mu i} c_{\nu i} \quad . \end{equation}

For the Fock matrix we then obtain

\begin{equation} \begin{split} F_{\mu\nu} & = h_{\mu\nu} + \sum_{\kappa\lambda} \left[ 2\big(\mu\nu|\kappa\lambda\big) D_{\kappa\lambda} - \big(\mu \kappa|\lambda \nu\big) D_{\kappa\lambda} \right] \\ & = h_{\mu\nu} + \sum_{\kappa\lambda} D_{\kappa\lambda} \left[ 2\big(\mu\nu|\kappa\lambda\big) - \big(\mu \kappa|\lambda \nu\big) \right] \end{split} \end{equation}

In a later article we will see how this formulation helps us to incrementally build up the Fock matrix using the change in the density matrix from iteration to iteration, which paves the way for efficient integral screening and integral direct implementations.

social