The Hartree self-consistent field method

The Hartree self-consistent field method is an iterative procedure that optimises an approximate wavefunction using the variational principle. Its goal is to estimate the eigenvalues of a modified form of the non-relativistic multi-electron Hamiltonian. This method was developed by the English scientist Douglas Hartree in 1927.

In an earlier article, we showed that the three dimensional Hamiltonian operator $\hat{H}_T$ for an n-electron atom (excluding spin-orbit and other interactions) is:

$\hat{H}_T=-\frac{\hbar^2}{2m_e}\sum_{i=1}^{n}\nabla_i^2-\sum_{i=1}^{n}\frac{Ze^2}{4\pi\varepsilon_0r_i}+\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}\frac{e^2}{4\pi\varepsilon_0r_{ij}}\;\;\;\;\;\;\;\;(1)$

where $\nabla_i^2=\frac{\partial^2}{\partial x_i^2}+\frac{\partial^2}{\partial y_i^2}+\frac{\partial^2}{\partial z_i^2}$ , $r_i=\sqrt{x_i^2+y_i^2+z_i^2}$ and $r_{ij}=\vert\mathbf{r}_i-\mathbf{r}_j\vert=\sqrt{\left(x_i-x_j\right)^2-\left(y_i-y_j\right)^2-\left(z_i-z_j\right)^2}$ .

Our objective now is to solve eigenvalue equations for $\hat{H}_T$ . The eigenvalue equation of the hydrogen atom has an exact solution (i.e. a solution represented by a formula). However, for an atom with more than one electron, the inter-electronic repulsion term $\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}\frac{e^2}{4\pi\varepsilon_0r_{ij}}$ poses a problem, which is a three-body problem that has no exact solution. Therefore, approximation methods are needed.

Firstly, we simplify the Hamiltonian by rewriting it in a system devised by Hartree called “atomic units” (abbreviated to a.u.), where $\hbar=m_e=\vert e\vert=a_0=1$ ( $a_0$ is the Bohr radius, where $a_0=\frac{4\pi\varepsilon_0\hbar^2}{m_ee^2}$ ). For example, if the charge of a particle is 3.4 a.u., it is $3.4\times1.602\times10^{-19}C$ in S.I. units. Eq1 therefore becomes

$\hat{H}_T=-\frac{1}{2}\sum_{i=1}^{n}\nabla_i^2-\sum_{i=1}^{n}\frac{Z}{r_i}+\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}\frac{1}{r_{ij}}\;\;\;\;\;\;\;\;(2)$

with the eigenvalues of $\widehat{H}_T$ also in atomic units.

Secondly, we approximate the form a multi-electron wavefunction takes. As hydrogenic wavefunctions have exact solutions, it is logical to express the multi-electron wavefunction as products of orthonormal one-electron spatial wavefunctions $\phi_i(\mathbf r_i)$ :

$\psi(\mathbf r_1, \mathbf r_2,\cdots,\mathbf r_n)=\phi_1(\mathbf r_1)\phi_2(\mathbf r_2)\cdots\phi_n(\mathbf r_n)\;\;\;\;\;\;\;\;(3)$

Thirdly, we utilise the variational principle, which states that the calculated energy $E=\int\psi^*\widehat{H}\psi d\tau$ is always greater than or equal to the ground state $E_0$ of any system. This implies that we can calculate an upper bound to $E_0$ of any system by using any trial function $\psi$ . The closer the trial function is to the actual or best wavefunction of the system $\psi_0$ , the closer E will be to $E_0$ .

Using eq3 and eq2, we have

$E=\int\cdots\int\phi_1^*(\mathbf r_1)\cdots\phi_n^*(\mathbf r_n)\left(-\frac{1}{2}\sum_{i=1}^n\nabla_i^2-\sum_{i=1}^n\frac{Z}{r_i}+\sum_{i=1}^{n-1}\sum_{j=i+1}^n\frac{1}{r_{ij}}\right)\phi_1(\mathbf r_1)\cdots\phi_n(\mathbf r_n)d\tau_1\cdots d\tau_n\;\;\;\;\;\;\;\;(4)$

Simplifying the above equation and minimising E (see this article for proof), we obtain the following set of n ‘one-particle’ equations called the Hartree equations:

$\left(-\frac{1}{2}\nabla_i^2-\frac{Z}{r_i}+\sum_{j\neq i}\int\frac{\vert\phi_j\vert^2}{r_{ij}}d\mathbf r_j\right)\phi_i(\mathbf r_i)=\varepsilon_i\phi_i(\mathbf r_i)\;\;\;\;\;\;\;\;(5)$

where $i=1\cdots n$ .

Question

How do we interpret eq5?

Answer

Eq5 can be perceived as a set of n Schrodinger equations, each consisting of a modified one-electron non-relativistic Hamiltonian. Solving each of the n equations gives a minimised value of $\varepsilon_i$ corresponding to the i-th electron.

To further explain the modified inter-electronic repulsion term, we refer to eq2. Focussing on electron 1, we have

$\hat{H}_T=-\frac{1}{2}\nabla_1^2-\frac{Z}{r_1}+\sum_{j\neq 1}\frac{1}{r_{1j}}$

We regard electrons 2, 3, …, n as being smeared out into a continuous and spherically symmetrical charge distribution through which electron 1 moves. For example, the interaction $V_{12}$ between electron 1 and electron 2 is

$V_{12}=\frac{q_1}{4\pi\varepsilon_0}\int\frac{\rho_2}{r_{12}}d\mathbf r_2\;\;\;\;\;\;(in\; S.I.\; units)$

where $\rho_2$ is the charge density of (i.e. charge per unit volume) of electron 2.

Since $q_1=-e$ , $\int\rho_2 d\mathbf r_2=-e$ and $\int\vert\phi_2\vert^2d\mathbf r_2=1$ , we have $\int\left(-e\vert\phi_2\vert^2\right)d\mathbf r_2=-e$ and hence $\rho_2=-e\vert\phi_2\vert^2$ . So,

$V_{12}=\int\frac{\vert\phi_2\vert^2}{r_{12}}d\mathbf r_2\;\;\;\;\;\;(in \;atomic\; units)$

The total potential energy between electron 1 and all other electrons is therefore

$V_1\left(r_1,\theta_1,\phi_1\right)=V_{12}+V_{13}+\cdots+V_{1n}=\sum_{j=2}^n\int\frac{\vert\phi_j\vert^2}{r_{1j}}d\mathbf r_j=\sum_{j\neq 1}\int\frac{\vert\phi_j\vert^2}{r_{1j}}d\mathbf r_j$

Since eq5 is the outcome of minimising eq4, subject to the constraints of the orthonormality of one-electron wavefunctions that make up the overall wavefunction of a system, any solution set of $\phi_i(\mathbf r_i)$ from eq5 is associated with a set of $\varepsilon_i$ , each of which is also obtained from multiplying eq5 by $\phi_i(\mathbf r_i)$ and integrating over all space:

$\varepsilon_i=H_i+\sum_{j\neq i}J_{ij}\;\;\;\;\;\;\;\;(6)$

where $H_i=\int\phi_i^*(\mathbf r_i)\left(-\frac{1}{2}\nabla_i^2-\frac{Z}{r_i}\right)\phi_i(\mathbf r_i)d\mathbf r_i$ and $J_{ij}=\int\int\frac{\vert\phi_i\vert^2\vert\phi_j\vert^2}{r_{ij}}d\mathbf r_id\mathbf r_j$ .

Summing all $\varepsilon_i$ of eq6, we have

$\sum_{i=1}^n\varepsilon_i=\sum_{i=1}^nH_i+\sum_{i=1}^n\sum_{j\neq i}J_{ij}$

Since $\frac{1}{2}\sum_{j\neq i}\sum_{i=1}^n=\sum_{j>i}\sum_{i=1}^n$ for the Coulomb term (see this article for explanation),

$\sum_{i=1}^n\varepsilon_i=\sum_{i=1}^nH_i+2\sum_{J>i}\sum_{i=1}^nJ_{ij}\;\;\;\;\;\;\;\;(7)$

To establish the relation between E and $\varepsilon_i$ , we rewrite eq4 as

$E=\sum_{i=1}^nH_i+\sum_{J>i}\sum_{i=1}^nJ_{ij}\;\;\;\;\;\;\;\;(8)$

where we have used the fact that $\sum_{i=1}^{n-1}\sum_{j=i+1}^n$ is equivalent to $\sum_{j>i}\sum_{i}^n$ (see this article for explanation).

Combining eq7 and eq8,

$E=\sum_{i=1}^n\varepsilon_i-\sum_{J>i}\sum_{i=1}^nJ_{ij}\;\;\;\;\;\;\;\;(9)$

E is greater than the ground state $E_0$ of the system, except when the solution set $\varepsilon_i$ is associated with the set of one-electron wavefunctions $\phi_i(\mathbf r_i)$ that best describes the ground state of the system, in which case E is the best estimate of $E_0$ . This best set of $\phi_i(\mathbf r_i)$ and the corresponding set of $\varepsilon_i$ are found using the Hartree self-consistent field method, which involves the following steps:

Make guesses for all single-electron wavefunctions except for one, e.g. use trial functions for $\phi_2(\mathbf r_2)\phi_3(\mathbf r_3)\cdots \phi_n(\mathbf r_n)$ except $\phi_1(\mathbf r_1)$ . The corresponding equation to solve is:

$\left[-\frac{1}{2}\nabla_1^2-\frac{Z}{r_1}+\sum_{j\neq 1}\int\frac{\vert\phi_j\vert^2}{r_{1j}}d\mathbf r_j\right]\phi_1(\mathbf r_1)=\varepsilon_1\phi_1(\mathbf r_1)$

Solving the above equation means finding $\varepsilon_1$ and $\phi_1(\mathbf r_1)$ , which can only be done numerically because of the presence of unknown functions and variables (see this article for an example). The solutions are then the first approximations of $\varepsilon_1$ and $\phi_1(\mathbf r_1)$ , which leads to a first improvement to eq3.

To further improve the wavefunction of eq3, we solve for $\varepsilon_2$ and $\phi_2(\mathbf r_2)$ using

$\left[-\frac{1}{2}\nabla_2^2-\frac{Z}{r_2}+\sum_{j\neq 2}\int\frac{\vert\phi_j\vert^2}{r_{2j}}d\mathbf r_j\right]\phi_2(\mathbf r_2)=\varepsilon_2\phi_2(\mathbf r_2)$

As per the previous step, the same trial functions for $\phi_3(\mathbf r_3)\phi_4(\mathbf r_4)\cdots \phi_n(\mathbf r_n)$ are employed together with the computed $\phi_1(\mathbf r_1)$ . The solutions are then the 1^st approximation of $\varepsilon_2$ and $\phi_2(\mathbf r_2)$ , and hence a further improvement to eq3.

Next, $\phi_1(\mathbf r_1)$ and $\phi_2(\mathbf r_2)$ derived from the previous two steps are used along with trial functions for $\phi_4(\mathbf r_4)\phi_5(\mathbf r_5)\cdots \phi_n(\mathbf r_n)$ to solve for $\varepsilon_3$ and $\phi_3(\mathbf r_3)$ . The process is repeated until $\varepsilon_n$ and $\phi_n(\mathbf r_n)$ are solved. The entire cycle is then repeated, i.e. step 1 is repeated using improved wavefunctions of $\phi_2(\mathbf r_2)\phi_3(\mathbf r_3)\cdots \phi_n(\mathbf r_n)$ to solve for an even better approximation of $\phi_1(\mathbf r_1)$ , and then followed by step 2 and step 3. This cyclical process stops when all $\varepsilon_i$ and $\phi_i(\mathbf r_i)$ become invariant, which implies that the best approximation of eq3 for the particular system has been derived.

Finally, E is computed with eq9 using the derived $\varepsilon_i$ and $\phi_i(\mathbf r_i)$ .

As the numerical method of calculating E is usually quite complex, an analytical method is often used. Such a method involves finding analytical expressions for all terms in eq8 using trial one-electron wavefunctions. As per the numerical method, guess values are used for the variables in all trial one-electron wavefunctions except one, e.g. $\phi_1(\mathbf r_1)$ . Eq8 is then minimised to obtain the solutions for E and $\phi_1(\mathbf r_1)$ . The process is repeated until E and all $\phi_i(\mathbf r_i)$ become invariant.

In conclusion, the Hartree self-consistent field method produces relatively accurate results when compared to experiment data for the ground state energy of a small atom like helium. However, it is less accurate for finding E of larger atoms, as it does not take into consideration exchange forces due to spin exchange interactions. A more accurate method, call the Hartree-Fock method is needed for such atoms.

The Hartree self-consistent field method

Question

Answer

Next article: Lagrange method of undetermined multipliers

Content page of computational chemistry

Content page of advanced chemistry

Main content page

Leave a Reply Cancel reply