Electronic problem: Hartree-Fock

In this subsection we will try to show the main topics about the task that has occupied theoretical chemists during more than seventy years, the solution of electronic Schrödinger equation 1.6 for molecules. Obviously no analytical solution for this equation exists. Actually the many-body interacting problem is a very common problem that exists in many fields of physics, from the subparticles world to the astronomy.

The usual strategy is to consider a non-interacting system whose solution is usually known and then try to solve the real problem by means of perturbation theory or variational theory.

Independent particle model: Hydrogen atom:
The hydrogen atom was solved analytically by Schrödinger in 1926. Its solution serves as a basis to consider the molecular case[8].

$$\hat h\Phi=\epsilon\Phi \quad with \quad \Phi(r\theta\phi)=R_{nl}(r)Y_{lm}(\theta \phi )$$ (2.16)

Where $\Phi$ is the atomic orbital (radial and angular) and $\hat h$ is the mono-electronic Hamiltonian

$$\hat h = -\frac{1}{2}\Delta - \frac{1}{r}$$ (2.17)

Variational Principle:
The variational principle states that any wavefunction of the space accomplishes that its expectation value of the energy is an upper bound to the ground state energy of the system $E_0$.

$$E_0 \leq \frac{\langle\Psi\vert\hat H \vert\Psi\rangle}{\langle\Psi\vert\Psi\rangle}$$ (2.18)

For the sake of clarity the superscript $electronic$ will be avoided from now on.

A test function: Slater determinant
An intuitive initial wavefunction would be considering the independent particle model. That is, a molecule constituted by non-interacting electrons so that every electron will be in an atomic spin-orbital (from here on $\phi_i$ for atoms) or in a molecular spin-orbital ($\Phi$ for molecules). This choice characterizes the molecular orbital theory. An alternative to this choice are the Valence Bond Methods where the test function is a combination of configurations of atomic orbitals. This strategy was named the Heitler-London-Pauling-Slater [8] and during the thirties it was a good alternative to the molecular orbital theory. However, despite valence bond methods give a chemical vision of the wavefunction easier to interpret, the numerical convergence is low and quantitative results are expensive.

The operator of non-interacting electrons is a sum of independent one-particle operators and its eigenfunction is the product of the one-particle wavefunctions, this is the so-called Hartree product

$$\vert\Psi^{hartree}\rangle= \prod_i^n \Phi(i)$$ (2.19)

Since the electrons are fermions they must be antisymmetric when interchanged. To antisymmetrize the Hartree product we will introduce the Slater determinant

$$\vert\Psi^{slater}\rangle= \frac{1}{\sqrt {n!}} \begin{vmatrix}\P... ...ddots & \vdots \\ \Phi_1(n) & \Phi_2(n) & \ldots & \Phi_n(n) \\ \end{vmatrix}$$ (2.20)

This can be a good test function, even though it is important to note that we are dealing with only one function of the electronic space (equation 1.16 has infinite solutions), that is, the exact wavefunction for the real system is found in the basis that considers all electrons in all orbitals. It is defined a Slater determinant as a unique configuration. Then the exact solution cannot come from here but when we consider all possible configurations. In section 1.2.1.7 this option is commented as Configuration Interaction.

In any case we will see that a Slater determinant is not that bad. So, let us put this Slater determinant in equation 1.18. Considering the orthonormality of the molecular orbitals we obtain the following expressions known as Slater rules[28].

$$\langle\Psi^{slater}\vert\hat H \vert\Psi^{slater}\rangle = \sum_i^n h_{ii} + \sum_i^n\sum_{j>i}^n(J_{ij}-K_{ij})$$ (2.21)

where

$$h_{ii} = \langle \Phi_i(1) \vert-\frac{1}{2}\Delta_i - \sum_K^N \frac{Z_K}{R_{iK}}\vert\Phi_i(1)\rangle$$

$$= \langle \Phi_i(1) \vert\hat h_1 \vert\Phi_i(1) \rangle$$ (2.22)

$$J_{ij} = \langle \Phi_i(1)\Phi_j(2) \vert\frac{1}{r_{12}}\vert\Phi_i(1)\Phi_j(2)\rangle$$ (2.23)

$$K_{ij} = \langle \Phi_i(1)\Phi_j(2) \vert\frac{1}{r_{12}}\vert\Phi_j(1)\Phi_i(2)\rangle$$ (2.24)

The coulombic integrals $J_{ij}$ and the exchange integrals $K_{ij}$ can be rewritten in an operator form that will be useful later.

$$\hat J_{i}\vert\Phi_j(2)\rangle = \langle \Phi_i(1)\vert\frac{1}{r_{12}}\vert\Phi_i(1)\rangle \vert\Phi_j(2)\rangle$$ (2.25)

$$\hat K_{i}\vert\Phi_j(2)\rangle = \langle \Phi_i(1)\vert\frac{1}{r_{12}}\vert\Phi_j(1)\rangle \vert\Phi_i(2)\rangle$$ (2.26)

In this way, it can be seen the coulombic $J_{ij}$ and exchange integrals $K_{ij}$ as the matrix elements representation of operators $\hat J_i$ and $\hat K_i$ in the molecular orbitals basis. Obviously, the $\hat J_i$ and $\hat K_i$ depend explicitly on the molecular orbitals themselves.

Minimization of the energy:
In order to obtain the best molecular orbitals that give the lowest energy we have to minimize the energy with respect to the orbitals, provided that these orbitals cannot be zero or linear dependent, that is, constraining the orbitals to be orthogonal and to have a norm of one.

As in many other situations, the optimization of a functional with constraints is carried out by the Lagrange Multipliers method. Building the Lagrangian $\cal L$, differentiating with respect to the molecular orbitals and imposing the stationary conditions ( $\delta\cal L$= 0) we could see that the molecular orbitals can be obtained by solving the final Hartree-Fock (HF) equations:

$$(\hat h_i + \sum_j^n (\hat J_j - \hat K_j))\vert\Phi_i\rangle=\epsilon_i\vert\Phi_i\rangle$$ (2.27)

or defining the Fock operator as

$$\hat F=\hat h_i + \sum_j^n (\hat J_j - \hat K_j)$$ (2.28)

we have the condensed form

$$\hat F\vert\Phi_i\rangle=\epsilon_i\vert\Phi_i\rangle \quad ; \quad \epsilon_i=\langle \Phi_i\vert\hat F\vert\Phi_i\rangle$$ (2.29)

where $\epsilon_i$ are the Lagrange multipliers and can be interpreted somewhat as the energy of the $ith$ molecular orbital (Koopmans theorem gives an alternative and sometimes useful definition[9]). However we must recall that the Hartree-Fock equations are just an strategy to obtain the optimized molecular orbitals and that the energy comes from the equation 1.21.

Equation 1.29 is a pseudo-eigenvalue equation because the Fock operator cannot be known unless we know the molecular orbitals, and molecular orbitals are obtained by the Fock operator. This dependence (exemplified in equations 1.25 and 1.26) forces that the Hartree-Fock equations must be solved iteratively until self-consistency. This procedure is called Self-Consistent-Field (SCF).

Introduction of a basis: Roothan-Hall equations:
Hartree-Fock equations can only be solved numerically (mapping the orbitals on a set of grid points) for highly symmetric systems (mostly atoms). In molecular systems we must introduce a known basis set { $\vert\phi_j\rangle$} to span our molecular orbitals as a combination of these $m$ basis functions.

$$\vert\Phi_i\rangle=\sum_r^m c_{ri} \vert\phi_{r}\rangle$$ (2.30)

This is the Linear Combination of Atomic Orbitals method (LCAO). The term atomic is used because it has been seen that taking atomic-like orbitals centered on every nucleus as a basis set the numerical results are optimal. The implementation of the LCAO method (equation 1.30) to the HF equations (equation 1.29) leads the so-called Roothan-Hall equations[10]

$$\hat F\vert\sum_r^m c_{ri}\vert\phi_{r}\rangle=\epsilon_i\sum_r^m c_{ri}\vert\phi_{r}\rangle$$ (2.31)

Multiplying equation 1.31 by $\langle \phi_{s}\vert$ with $s=1,\ldots,m$ and writing the equation in matrix form

$${\bf FC=SC\epsilon}$$ (2.32)

where S is the overlap matrix and is not diagonal since the basis set are not orthogonal. An orthogonalization of the basis can lead to an eigenvalue equation ${\bf F'C'=C'\epsilon}$. Remember that the fock operator $\hat F$ depends on the molecular orbitals, so if we substitute equation 1.30 in 1.28 we have an expression for the elements of fock matrix $F_{pq}$ in the basis set representation

$$F_{pq} = \langle \phi_p\vert\hat F \vert\phi_q\rangle$$

$$= \langle \phi_p\vert\hat h \vert\phi_q\rangle + \sum_i^n \langle \phi_p\vert\hat J_i - \hat K_i \vert\phi_q\rangle$$ (2.33)

$$= \langle \phi_p\vert\hat h \vert\phi_q\rangle$$

$$+ \sum_i^n\sum_r^m\sum_s^m c_{ri}c_{si}( \langle \phi_p\phi_r \vert... ...\phi_s\rangle - \langle \phi_p\phi_r \vert\frac{1}{r}\vert\phi_s\phi_q\rangle )$$ (2.34)

$$F_{pq} = h_{pq} + \sum_{r,s}P_{rs}G_{pqrs}$$ (2.35)

The term $h_{pq}$ are the core integrals and $G_{pqrs}$ the two-electron integrals. The two electron integrals nomenclature is usually simplified as $\langle pr\vert sq\rangle$. Depending on the atom where the four different functions are centered it will be a two-electron integral of one, two, three or four centers. The three and four centers are the most common integrals and then the most expensive to evaluate. The matrix $P_{rs} \equiv\sum_i^n c_{ri}c_{si}$ is usually called the density matrix. Its elements will be the coefficients to optimize in order to have a diagonal representation of the Fock matrix and therefore the pursued coefficients to obtain the molecular orbitals.

The method described above is the Restricted Hartree-Fock (RHF), every spin-orbital is a molecular orbital occupied by two electrons with spin function $\alpha$ and $\beta$ so that the whole expectation value of spin operator $\hat S^2$ is zero. For open shell systems there are the Restricted Open Shell HF (ROHF) or the more adequate Unrestricted HF method (UHF). The corresponding Roothan-Hall equations for the UHF case are the Pople-Nesbet equations.

Few comments about the basis set:
The analytical solution of Hydrogen atom of equation 1.16 has a radial part represented by the Laguerre polynomials and an angular part represented by the Spherical Harmonics which are basically the associate Legendre Polynomials [8,9]. The basis that we usually introduce in the Roothan-Hall equations will have a radial and an angular part. The angular is almost never commented because it is always the same (s,p,d,f...). It is the radial part which decides how far or how close is the electron with respect to the nucleus. The radial function can be represented by the Slater Type Orbitals (STO). However, despite of the adequacy of the STO basis, the integrals are not analytical using STO. Then to avoid expensive numerical integrals it is more common to use Gaussian functions (GTO) which integrals are analytical or easier to evaluate. In any case some modern methods still employ the STO functions. Furthermore, the semiempirical methods used in this thesis are carried out by STO as well. Many textbooks can be found covering the basis set, the different types, their efficiency and the computational requirements[11].

The SCF process:
We can summarize all the HF process through the Roothan-Hall equations as follows:

1. Input data: composition ($Z_k$) and geometry ($X_K$) of the system and functions basis set $\{\vert\phi_r\rangle\}$
2. Integral computing: One electron ($h_{pq}$ and $S_{pq}$) and two electron ($G_{pqrs}$)
3. Orthogonalize the basis set: ${\bf C}={\bf LC'}$ so that ${\bf S}\to{\bf 1}$
4. Give an initial guess to obtain a tentative molecular orbitals
5. Build the Fock matrix as core integrals $+$ density matrix $\times$ two-integrals
6. Transform the Fock matrix to the orthogonal basis set ${\bf F'=L^TFL}$
7. Diagonalize the Fock matrix and obtain the coefficients ${\bf C'}\to{\bf C}$ and the new molecular orbitals
8. Compute the HF energy and check for convergence.
9. If it is not converged compute the new density matrix and go back to point 5.

This procedure has some bottlenecks that could make the whole process not feasible. The computation and the storage of integrals, several diagonalization and matrix transformations combined with a slow convergence process. A discussion on how expensive is every step and the different possibilities to overcome these problems may be found in the literature. Particularly the linear scaling techniques[29] deal with the problem of converting the whole computational requirements into a process that scales linearly with the size and/or with the basis set.

It is important to note that the HF method is the best method for a one-configuration wavefunction because we have used the exact electronic Hamiltonian to develop the equations. However there are many improvements to the HF method. Some of them will be outlined very briefly in section 1.2.1.7