## Roothaan-Hall equations¶

The Hartree-Fock equations need to be converted into a tabular form (matrix form) using the basis set of atomic orbitals to allow for a solution to be determined using the computer. We insert our basis set of orbitals $\sum_i c_i \phi^{CGF}_i(\boldsymbol{r}_i)$ (we will drop the superscript and just have the orbitals denoted as $\phi_i$) We are also going to be building up a table by combining different basis set indices so will use the subscripts $\nu$ and $\mu$ to denotes the indices as we build up the table/matrix. We will also replace $(\boldsymbol{r}_\nu)$ with $(1)$ for ease of reading.

$$ \mathscr{f}_i\psi_i(\boldsymbol{r}_i) = \epsilon_{i}\psi_i(\boldsymbol{r}_i) $$$$ \mathscr{f}_i\sum^K_{\nu =1} c_{\nu 1} \phi_{\nu}(1) = \epsilon_{i}\sum^K_{\nu =1} c_{\nu i} \phi_{\nu}(1) $$

We left multiply the other basis set index of our table and integrate.

$$ \sum^K_{\nu =1}c_{\nu 1} \int d\nu_1 \phi_{\mu}(1) \mathscr{f}_i \phi_{\nu}(1) = \epsilon_{i}\sum^K_{\nu =1} c_{\nu i}\int d\nu_1 \phi_{\mu}(1)\phi_{\nu}(1) $$

This last term $\int d\nu_1 \phi_{\mu}(1)\phi_{\nu}(1)$ is called the overlap integral written as $S_{\mu \nu}$. This term will help to illustrate the population of a table if we want to know how much the 1s orbital on the He interacts with the 1s orbital on the H atom we would calculate the values on the diagonal of the overlap matrix (for our calculation below this is 0.45077041 in the code the variable is called S12). However, the overlap of the basis function with itself is one by definition.

$$S_{\mu \nu}= \left[ \begin{matrix} 1.0 & 0.45077041 \\
0.45077041 & 1.0
\end{matrix}\right]
$$

We have another matrix called the Fock matrix $F_{\mu \nu}$ which is made up of the core contributions $\mathscr{H}^{CORE}$ and the effective potential $v_{HF}$ lets see how the matrices are constructed from these.

$$
\int d\nu_1 \phi_{\mu}(1)\mathscr{H}^{CORE}\phi_{\nu}(1) = \int d\nu_1 \phi_{\mu}(1)\left[ -\sum_i \frac{\nabla_i^{2}}{2}- \sum_{A,i}\frac{Z_{A}}{r_{Ai}} \right]\phi_{\nu}(1) = H^{CORE}_{\mu \nu}
$$

Where $H^{CORE}_{\mu \nu}$ is another matrix that has to be calculated containing a kinetic energy and potential energy contribution. These are called one electron integrals (as they are the interactions of a single electron in a molecular orbital). In the code the core contributions are broken up into kinetic terms (T11, T12 and T22, note we do not need T21 as the matrix is symmetrical) and potential terms (V11A, V12A, V22A for atom one and V11B, V12B, V22B for the second atom, note the sum is also over the nuclei for this term)

The effective Hartree potential is a bit tricker to calculate.

$$
\int d\nu_1 \phi_{\mu}(1)v_{HF}\phi_{\nu}(1) = \sum_{j=1}^{N/2}\int d\nu_1 \phi_{\mu}(1)\left[ 2\mathscr{J}_j - \mathscr{K}_j \right]\phi_{\nu}(1) = G^{2e}_{\mu \nu \lambda \sigma}
$$

Where $N$ is the total number of electrons and for closed shell molecules we want half of these electrons, due to two electrons residing in each occupied molecular orbital. This integrate is tricky as it involves the impact of exchange which requires that electrons are indistingishable and therefore instead of a matrix we have 3D matrix or a cube of values (often called a tensor in mathematics). We therefore have to introduce two new indices $\lambda$ and $\sigma$ to denote these indistinguishable versions of the electrons. Which creates quite a construction.

$$
G^{2e}_{\mu \nu \lambda \sigma} = \sum^{N/2}_{j=1}\sum^K_{\lambda=1}\sum^K_{\sigma=1}c_{\lambda j}c_{\sigma j}\left[ \begin{matrix}
2\int d\nu_1d\nu_2\phi_\mu(1)\phi_\nu(1)\frac{1}{r_{12}}\phi_\lambda(2)\phi_\sigma(2)\\
-\int d\nu_1\nu_2 \phi_\mu(1)\phi_\lambda(1)\frac{1}{r_{12}}\phi_{\nu}(2)\phi_\sigma(2)
\end{matrix}
\right]
$$

Or in short hand

$$
G^{2e}_{\mu \nu \lambda \sigma} = \sum^{N/2}_{j=1}\sum^K_{\lambda=1}\sum^K_{\sigma=1}c_{\lambda j}c_{\sigma j}\left[2(\mu\nu|\lambda\sigma)-(\mu\lambda|\nu\sigma) \right]
$$

The code uses the variable names V1111, V2111, V2121, V2211, V2221 and V2222 and then constructs $G^{2e}_{\mu \nu \lambda \sigma}$ with the variable name G.

Another shorthand that will be helpful is to consider the charge density matrix $P$

$$
P_{\mu \nu}=2\sum^{N/2}_{i=1}c_{\mu i}c_{\nu i}\; and\; P_{\lambda \sigma}=2\sum^{N/2}_{i=1}c_{\lambda i}c_{\sigma i}
$$

The electron density can be determined from the density matrix as

$$
\rho(\boldsymbol{r})=\sum^K_{\mu=1}\sum^K_{\nu=1}c_{\mu}(\boldsymbol{r})c_{\nu}(\boldsymbol{r})
$$

The Fock matrix $F_{\mu \nu}$ is then calculated as for closed shell molecules

$$
F_{\mu \nu} = H^{CORE}_{\mu \nu} + \sum^K_{\lambda=1}\sum^K_{\sigma=1}P_{\lambda \sigma}\left[2(\mu\nu|\lambda\sigma)-(\mu\lambda|\nu\sigma) \right] = H^{CORE}_{\mu \nu} + G^{2e}_{\mu \nu \lambda \sigma}
$$

In the code G is used when calculating this term

Putting this back into the Hartree-Fock equation we have now transformed into five matrices/tables or the Roothaan-Hall equations

$$ \boldsymbol{F} \boldsymbol{C} = \boldsymbol{S} \boldsymbol{C} \boldsymbol{E} $$

The coefficient matrix $\boldsymbol{C}$ is a table with $K\times K$ entries with coefficients where $K$ are the number of electrons.

$$
\boldsymbol{C} = \left( \begin{matrix} c_{1,1} & c_{1,2} & ... & c_{1,K} \\
c_{2,1} & c_{2,2} & ... & c_{2,K} \\
\vdots & \vdots & & \vdots \\
c_{K,1} & c_{K,2} & ... & c_{K,K}
\end{matrix} \right)
$$

$\boldsymbol{E}$ is a diagonal matrix of the energies of the molecular orbitals

$$
\boldsymbol{E} = \left( \begin{matrix} \epsilon_1 & 0 & ... & 0 \\
0 & \epsilon_2 & ... & 0 \\
\vdots & \vdots & & \vdots \\
0 & 0 & ... & \epsilon_K
\end{matrix} \right)
$$

We are trying to find out the values of $\boldsymbol{C}$ which minimise the energy however $\boldsymbol{C}$ is on both sides of the equation so we have to iterative solve the equations until the $\boldsymbol{C}$ on both sides are equal and the energy is the lowest possible this is called the self consistent field method or SCF for short.

Computers are best at solving equations in the form of $\boldsymbol{FC}=\boldsymbol{CE}$, however, that overlap integral is a bit of a problem ($\boldsymbol{F} \boldsymbol{C} = \boldsymbol{S} \boldsymbol{C} \boldsymbol{E}$) so we have to do some matrix algebra to rearrange the equation into an equivalent equation $\boldsymbol{F'C'}=\boldsymbol{C'E}$ (where $\boldsymbol{F'}=\boldsymbol{S}^{-1/2}\boldsymbol{F}\boldsymbol{S}^{-1/2}$ and $\boldsymbol{C'}=\boldsymbol{S}^{-1/2}\boldsymbol{C}$ (In the code $\boldsymbol{S}^{-1/2}$ and $\boldsymbol{S}^{1/2}$ are denoted as X and X.T )