Task 2: Restricted Hartree-Fock ground state solver

We start from the time-independent Schrödinger equation

\[\begin{align*} \hat{H} | \Psi \rangle = E | \Psi \rangle, \end{align*}\]

where \(| \Psi \rangle\) is a many-body wave function of \(N\) particles. We know that in the non-interacting case the exact ground state of the Hamiltonian will be a single Slater determinant (which we denote \(| \Phi \rangle\)) of the \(N\) lowest single-particle eigenstates. This motivates the ansatz of using a single Slater determinant for the full many-body problem

\[\begin{align*} | \Psi \rangle \approx | \Phi \rangle = \sqrt{N!} \hat{A} | \phi_1 \phi_2 \dots \phi_N \rangle, \end{align*}\]

where the \(N\) lowest molecular orbitals (MO) (also known as the Hartree-Fock orbitals) from \(\{|\phi_p\rangle \mid p = 1, \dots, L \}\) (\(N \leq L\)) are the primary unknowns, \(| \phi_1 \phi_2 \dots \phi_N \rangle = | \phi_1 \rangle \otimes | \phi_2 \rangle \otimes \cdots \otimes | \phi_N \rangle\) is a product state of the \(N\) first molecular orbitals, and \(\hat{A}\) is the anti symmetrizer making the product state into an anti symmetrized wave function.

Note

We use a convention where we split up the set of molecular orbitals into a set of occupied and virtual orbitals. That is

\[\begin{align*} \{ \phi_p \}_{p = 1}^{L} = \{ \phi_i \}_{i = 1}^{N} \cup \{ \phi_a \}_{a = N + 1}^{L}, \end{align*}\]

where we let the indices \(i, j, k, l, \dots \in \{1, \dots, N\}\) denote the occupied orbitals, i.e., the orbitals that are contained in the Hartree-Fock Slater determinant, \(a, b, c, d, \dots \in \{N + 1, \dots, L \}\) are the indices of the virtual orbitals that are not contained in the reference determinant, and \(p, q, r, s, \dots \in \{1, \dots, L\}\) the full set of orbitals.

The Hartree-Fock method method finds the Slater determinant that minimizes the energy, i.e., gets closest to the “true” many-body ground state energy \(E_{gs}\) of the problem using the variational principle

\[\begin{align*} E_{gs} \leq E[\Phi, \Phi^{*}] = E[\phi_1, \dots, \phi_N] = \langle \Phi | \hat{H} | \Phi \rangle, \end{align*}\]

subject to the constraint that the molecular orbitals are orthonormal, \(\langle \phi_i | \phi_j \rangle = \delta_{ij}\). This in turn implies that the Slater determinant is normalized, i.e., \(\langle \Phi | \Phi \rangle = 1\). To include this requirement in the variational principle we use Lagrange’s method of undetermined multipliers

(5)\[ L[\phi_1, \dots, \phi_N, \lambda] = E[\phi_1, \dots, \phi_N] - \sum_{i = 1}^{N} \sum_{j = 1}^{N} \lambda_{ji} \left( \langle \phi_i | \phi_j \rangle - \delta_{ij} \right), \]

where \(\lambda_{ji}\) are the Lagrange multipliers, one for each constraint. We require the Lagrangian, \(L\), to be real. As the constraints we use in general are hermitian, we require that that the Lagrange multipliers be hermitian as well, that is, \(\lambda_{ji} = \lambda^{*}_{ij}\). Expressed in terms of the molecular orbitals the Hartree-Fock energy is

(6)\[ E \equiv E[\phi_1, \dots, \phi_N] = \langle \Phi | \hat{H} | \Phi \rangle = \sum_{i = 1}^{N} \langle \phi_i | \hat{h} | \phi_i \rangle + \frac{1}{2} \sum_{i = 1}^{N} \sum_{j = 1}^{N} \left( \langle \phi_i \phi_j | \hat{u} | \phi_i \phi_j \rangle - \langle \phi_i \phi_j | \hat{u} | \phi_j \phi_i \rangle \right). \]

The notation here is general, but we will for the ground state solver use the one-body Hamiltonian from equation (2) , and the two-body Hamiltonian from equation (3).

Types of orbital systems

Before proceeding, we need to decide on which type of single-particle functions, and hence which type of Hartree-Fock method, we intend to use. In particular, we need to decide on how we wish to treat spin. We will be looking at two electrons with spin \(1/2\), and as the Hamiltonian in equation (1) is independent of spin, we know that the “true” ground state will be a singlet state. However, as we are working with an approximation to the full ground state wave function it is not necessarily guaranteed that our solution will be an eigenstate of the \(\hat{S}^2\)-operator. For spin \(s = 1/2\) we have the two spin-\(z\) projections \(m_s = \pm 1/2\). We will for brevity denote the two orthonormal spin-states as \(| \sigma \rangle = | s, m_s \rangle\) with \(| 0 \rangle = | 1/2, +1/2 \rangle\) for spin up and \(| 1 \rangle = | 1/2, -1/2 \rangle\) for spin down. To be able to control the spin symmetry of the problem we typically look at three different kinds of orbitals.

1. The general spin-orbital is on the form

   \[\begin{align*} | \phi_p \rangle = \sum_{\sigma = 0}^{1} a_{p\sigma} | \varphi_{p\sigma} \rangle | \sigma \rangle = a_{p0} | \varphi_{p0} \rangle | 0 \rangle + a_{p1} | \varphi_{p1} \rangle | 1 \rangle, \end{align*}\]

   where \(| \varphi_{p\sigma} \rangle\) are normalized spatial orbitals for a given spin-direction \(\sigma\). The coefficients \(\alpha_{p\sigma}\) are included to ensure an overall normalization, i.e., \(|\alpha_{p0}|^2 + |\alpha_{p1}|^2 = 1\).

2. The unrestricted spin-orbital is on the form

   \[\begin{align*} | \phi_p \rangle = | \varphi_{p} \rangle | \sigma \rangle = | \varphi_{\lfloor p / 2 \rfloor, \sigma} \rangle | \sigma \rangle, \end{align*}\]

   where an unrestricted spin-orbital has a definite spin-direction, but allows for the spatial part to be different between spin-directions. Here \(\lfloor p / 2 \rfloor\) is used to denote integer division.

3. The restricted spin-orbital is on the form

   (7)\[ | \phi_p \rangle = | \varphi_{\lfloor p / 2 \rfloor} \rangle | \sigma \rangle, \]

   where the spin-orbitals have a definite spin-direction with the spatial part equal between the two spin-directions. That is, each spatial orbital is doubly occupied.

The naming can be a bit misleading, but it is related to the three most common types of Hartree-Fock approximations. For the general Hartree-Fock method (GHF) using general spin-orbitals we only care about the energy. This means that the Hartree-Fock ground state \(|\Phi\rangle\) does not need to be an eigenstate of \(\hat{S}_z\) nor \(\hat{S}^2\). The unrestricted Hartree-Fock method (UHF) using unrestricted spin-orbitals will give a ground state that is an eigenstate of \(\hat{S}_z\), but not \(\hat{S}^2\) in general. Finally, the restricted Hartree-Fock method (RHF) using restricted spin-orbitals ensures that the wave function is an eigenstate of both \(\hat{S}_z\) and \(\hat{S}^2\). For two particles the RHF method will give a singlet ground state.

In this project we will opt for the RHF-method with restricted spin-orbitals.

Problem a)

Assuming that \(N\) is an even number (which it will be for us), show that the restricted Hartree-Fock energy will be given by

\[\begin{align*} E = 2 \sum_{i = 1}^{N / 2} \langle \varphi_i | \hat{h} | \varphi_i \rangle + \sum_{i = 1}^{N / 2} \sum_{j = 1}^{N / 2} \left( 2 \langle \varphi_i \varphi_j | \hat{u} | \varphi_i \varphi_j \rangle - \langle \varphi_i \varphi_j | \hat{u} | \varphi_j \varphi_i \rangle \right), \end{align*}\]

when inserting the restricted spin-orbitals from equation (7) into equation (6).

In coordinate representation we write the state as

\[\begin{align*} \Phi(x_1, \dots, x_N) = \frac{1}{\sqrt{N!}} \begin{vmatrix} \phi_1(x_1) & \dots & \phi_N(x_1) \\ \vdots & \ddots & \vdots \\ \phi_1(x_N) & \dots & \phi_N(x_N) \end{vmatrix}. \end{align*}\]

Problem a)

Add derivation of the (canonical) RHF-equations.

Minimization of the Lagrangian in equation (5) yields the Hartree-Fock eigenvalue equation

(8)\[ \hat{f} | \phi_p \rangle = \varepsilon_p | \phi_p \rangle, \]

where \(\hat{f}\) is the single-particle Fock-operator. Note that this equation is not limited to the occupied states, but to all the molecular orbitals. The matrix elements of the single-particle Fock operator are given by

(9)\[\begin{align} \langle \psi | \hat{f} | \xi \rangle = \langle \psi | \hat{h} | \xi \rangle + \sum_{i = 1}^{n} \langle \psi \phi_i | \hat{u} | \xi \phi_i \rangle - \sum_{i = 1}^{n} \langle \psi \phi_i | \hat{u} | \phi_i \xi \rangle = \langle \psi | \hat{h} | \xi \rangle + \sum_{i = 1}^{n} \langle \psi \phi_i | \hat{u} | \xi \phi_i \rangle_{AS}, \end{align}\]

where \(|\psi\rangle\) and \(|\xi \rangle\) are two arbitrary single-particle states.

The Roothan-Hall equations

The way we solve the Hartree-Fock equation in (8) is by expanding the Hartree-Fock orbitals \(\{ \phi_p \}_{p = 1}^k\) in an atomic orbital basis \(\{ \psi_{\mu} \}_{\mu = 1}^{l}\) (which we here take to be the orthonormal harmonic oscillator basis states with spin-doubling), and then projecting the equations onto the known basis set.

Note

In general \(k \leq l\) for the molecular orbitals. The consequence of this is that Hartree-Fock can be used as a method to truncate the basis set. This is common when using computationally more complex post Hartree-Fock methods such as e.g., coupled-cluster, to lower simulation times. However, in this project we limit our attention to \(k = l\).

We express a basis transformation from the AO basis to the MO basis by

(10)\[ |\phi_p\rangle = \sum_{\mu = 1}^{l} C_{\mu p} | \psi_{\mu} \rangle, \]

where the coefficients \(C_{\mu p}\) now become our primary unknowns in the Hartree-Fock method. Projecting the AO basis onto equation (8) we have

(11)\[\begin{gather} \langle \psi_{\mu} | \hat{f} | \phi_p \rangle = \varepsilon_{p} \langle \psi_{\mu} | \phi_p \rangle \implies \sum_{\nu = 1}^{l} f_{\mu \nu} C_{\nu p} = \varepsilon_{p} \sum_{\nu = 1}^{l} \delta_{\mu \nu} C_{\nu p} = C_{\mu p} \varepsilon_{p} \\ \implies \mathbf{F} \mathbf{C} = \mathbf{C} \boldsymbol{\varepsilon}, \end{gather}\]

where \([\mathbf{F}]_{\mu \nu} = f_{\mu \nu} = \langle \psi_{\mu} | \hat{f} | \psi_{\nu} \rangle\) are the Fock matrix elements in the AO basis, \([\mathbf{C}]_{\mu p} = C_{\mu p}\) the basis transformation coefficient matrix, and \(\boldsymbol{\varepsilon} = [\varepsilon_1, \dots, \varepsilon_k]^T\) the vector of Hartree-Fock single-particle eigenenergies (not to be confused with the full many-body energy). The last eigenvalue equation is known as the Roothan-Hall equations. Here both \(\mathbf{C}\) and \(\boldsymbol{\varepsilon}\) are the primary unknowns. We now need to construct the Fock matrix in the atomic orbital basis, that is

(12)\[\begin{align} f_{\mu \nu} &= \langle \psi_{\mu} | \hat{f} | \psi_{\nu} \rangle = \langle \psi_{\mu} | \hat{h} | \psi_{\nu} \rangle + \sum_{i = 1}^{n} \langle \psi_{\mu} \phi_i | \hat{u} | \psi_{\nu} \phi_i \rangle_{AS}. \end{align}\]

Replacing the sum over the occupied molecular orbitals in the two-body elements by equation (10) we get

(13)\[\begin{align} f_{\mu \nu} &= \langle \psi_{\mu} | \hat{h} | \psi_{\nu} \rangle + \sum_{i = 1}^{n} C^{*}_{\kappa i} C_{\lambda i} \langle \psi_{\mu} \psi_{\kappa} | \hat{u} | \psi_{\nu} \psi_{\lambda} \rangle_{AS} = h_{\mu \nu} + D_{\lambda \kappa} u^{\mu \kappa}_{\nu \lambda}, \end{align}\]

where we have defined the density matrix \(D_{\lambda \kappa} = \sum_{i = 1}^{n} C^{*}_{\kappa i} C_{\lambda i}\).

A complicating factor with the Roothan-Hall equations is that the Fock matrix \(\mathbf{F}\) depends on the coefficient matrix \(\mathbf{C}\). To get around this we use a technique called self-consistent field iterations.