Advanced Chemistry - Mono Mole

December 7, 2022September 26, 2024

Numerical solution of the Hartree self-consistent field method for He

In this article, we shall calculate the ground state energy of a helium atom in one-dimension using an Excel spreadsheet. From eq5, the equations to be solved iteratively are:

$\left[-\frac{1}{2}\frac{d^2}{dx_1^2}-\frac{2}{x_1}+\int_0^{\infty}\frac{\vert\phi_2(x_2)\vert^2}{\vert x_1-x_2\vert}dx_2\right]\phi_1(x_1)=\varepsilon_1\phi_1(x_1)$

$\left[-\frac{1}{2}\frac{d^2}{dx_2^2}-\frac{2}{x_2}+\int_0^{\infty}\frac{\vert\phi_1(x_1)\vert^2}{\vert x_2-x_1\vert}dx_1\right]\phi_2(x_2)=\varepsilon_2\phi_2(x_2)$

The presence of multiple functions and variables in each of these equations requires us to devise a computation strategy, the first of which is to derive a trial wavefunction.

Question

Show that $y=4\sqrt {2}xe^{-2x}$ is a reasonable wavefunction describing the helium ion.

Answer

The one-dimensional Schrodinger equation for the helium ion is $\left(-\frac{1}{2}\frac{d^2}{dx^2}-\frac{2}{x}\right)Y(x)=EY(x)$ , which can be written as

$P(x)y^{''}+Q(x)y=0\;\;\;\;\;\;\;\;(27)$
where $P(x)=-\frac{1}{2}$ , $Q(x)=\left(-\frac{2}{x}+E\right)$ , $y=Y(x)$ and $y^{''}=\frac{d^2}{dx^2}Y(x)$ .

The general solutions of eq27 are of the form $y=e^{bx}$ or $y=xe^{bx}$ . Let’s use $y=xe^{bx}$ , which is a modified form of a Slater-type orbital, where $b=2$ for the helium ion. Since $\int_0^{\infty}x^ne^{-mx}dx=\frac{n!}{m^{n+1}}$ , the normalisation constant N for y is $N=\sqrt\frac{m^{n+1}}{n!}=4\sqrt2$ . So, $y=4\sqrt {2}xe^{-2x}$ .

As the helium atom experiences an effective nuclear charge that is between that of a hydrogen atom ( $b=1$ ) and a helium ion ( $b=2$ ), it can expressed by the trial wavefunction $y=xe^{-3x/2}$ , which after normalisation is

$y=3\sqrt {\frac{3}{2}}xe^{-3x/2}\;\;\;\;\;\;\;\;(28)$
Next, the pair of equations mentioned at the top of the article have two distinct wavefunctions, each representing one of the two electrons. To simplify our computation, we assume that the two electrons can be expressed by the same spatial wavefunction. This is known as the restricted case. To iteratively solve for the common wavefunction, we need to transform the two equations into a new pair of equations. Rearranging eq5 for helium in one-dimension, we have

$\phi_i^{''}(x_i)=2\left[\int_0^{\infty}\frac{\phi_j^2(x_j)dx_j}{\vert x_i-x_j\vert}-\frac{2}{x_i}-\varepsilon_i\right]\phi_i(x_i)\;\;\;\;\;\;\;\;(29)$

where $x_i=x_0+dx_0$ and $\phi_i^{''}(x_i)=\frac{d^2}{dx^2}\phi_i(x_i)$ .

Substitute the definition of the 2^nd derivative $\phi_i^{''}(x_i)=\frac{\phi_i^{'}(x_i)-\phi_i^{'}(x_0)}{x_i-x_0}$ in eq29 and rearranging, we have

$\phi_i^{'}(x_i)=\phi_i^{'}(x_0)+2\left[\int_0^{\infty}\frac{\phi_j^2(x_j)dx_j}{\vert x_i-x_j\vert}-\frac{2}{x_i}-\varepsilon_i\right]\phi_i(x_i)(x_i-x_0)$

To prevent $\phi_i^{'}(x_i)$ from being undefined when $x_i=x_j$ , we modify the Coulomb integral by adding a small constant C to the denominator. It is found that $C=0.5$ gives reasonable computation results. So,

$\phi_i^{'}(x_i)=\phi_i^{'}(x_0)+2\left[\int_0^{\infty}\frac{\phi_j^2(x_j)dx_j}{\vert x_i-x_j\vert+0.5}-\frac{2}{x_i}-\varepsilon_i\right]$ $\times\phi_i(x_i)(x_i-x_0)\;\;\;\;\;\;\;\;(30)$

The definition of the 1^st derivative is $\phi_i^{'}(x_i)=\frac{\phi_i(x_i)-\phi_i(x_0)}{x_i-x_0}$
or:

$\phi_i(x_i)=\phi_i(x_0)+ \phi_i^{'}(x_i)(x_i-x_0)\;\;\;\;\;\;\;\;(31)$

Eq30 and eq31 are used in an iterative algorithm to find $\varepsilon_i$ and $\phi_i(x_i)$ . The procedure is as follows:

1) Set the boundary condition $\phi_i(x_{i=0})=0$ and $\phi_i(x_{i=\infty})=0$ . With reference to eq30, use eq28 as the initial guess wavefunction $\phi_j(x_j)$ for the Coulomb integral. Let’s also assume an initial guess gradient of $\phi_i^{'}(x_0)=3.7$ .

2) Let $0\leq x_i\leq 7$ with intervals of $dx_0=0.02$ (350 intervals in total). This implies that $x_i\rightarrow\infty$ at $x_i=7$ .

3) Use Simpson’s rule to integrate the Coulomb integral $V_{12}=\int_0^{\infty}\frac{\phi_j^2(x_j)dx_j}{\vert x_i-x_j\vert +0.5}$ for all $\boldsymbol{\mathit{x_i}}$ . An example of the excel table is shown below.

With reference to eq24, the formulae for cells C4, D4 and E4 are =(C2^2)*EXP(-2*1.5*C2)/(ABS(B4-C2)+0.5), =4*(D2^2)*EXP(-2*1.5*D2)/(ABS(B4-D2)+0.5) and =2*(E2^2)*EXP(-2*1.5*E2)/(ABS(B4-E2)+0.5) respectively. The formula for cell MP4 is =SUM(B4:MN4). Finally, the formula for cell MQ is 13.5*7/350/3*MP4, where 13.5 is the square of the normalisation constant of eq28. The formulae for the remaining cells within C4:MQ354 are populated accordingly.

4) The iterative table is formulated as follows:

The second row is for guess values and the error function. Rows 4 to 354 are for the algorithm. Cells MS2 and MV2 are the initial guess values of $\phi_i^{'}(x_0)$ and b respectively. The formula for the normalisation constant in cell MU2 is =SQRT(((-MV2*2)^3)/2). Cell MS4 is obviously zero. The formula for cell MT4 is =MS2. With reference to eq29 and eq28, the formulae for cells MS5 and MT5 are, =MS4+MT4*(B5-B4) and =MT4+(MQ5-(2/B5)-MT2)*2*MS5*(B5-B4) respectively. The formulae for the remaining cells within MS4:MT354 are populated accordingly.

Next, we plot a graph using values in cells MS4:MS354 as vertical coordinates and B4:B354 as horizontal coordinates. By visually inspecting the graph, the value in cell MT2 is manually adjusted so that $\phi_i(x_i)\rightarrow 0$ when $x_i\rightarrow\infty$ . To optimise the values in cells MS2 and MV2, we include a trendline in the plot, where the formulae for cells MU4 and MU5 is zero and =MU2*B5*EXP(MV2*B5) respectively. The formulae for the remaining cells within MU4:MU354 are populated accordingly.

The error function formula in cell MW2 is {=SQRT(SUM((MS5:MS365-MU5:MU354)^2))}. The curly brackets {} are obtained by first typing the formula within {} in the cell and pressing CTRL+SHIFT+ENTER. Lastly, the Excel Solver application is use to optimise the values in cells MS2 and MV2, with the following settings:

Set objective: MW2
To: Min
By changing variable cells: MS2,MV2
Selecting a solving method: GRG Nonlinear

5) For the next set of iterations, repeat steps 3 and 4, where the guess values of MS2 and MV2 are the optimised values of the previous set of iterations. The value for MU is no longer the formula =SQRT(((-MV2*2)^3)/2) but the optimised normalisation constant of the previous set of iterations. Similarly, the values of and used in the Coulomb integral in step 3 are no longer $N=\sqrt {13.5}$ and $b=1.5$ but the optimised values of MU2 and MV2 of the previous set of iterations.

6) The cycle continues until the values of $\varepsilon_i$ , N and b are invariant up to five decimal points.

The results are summarized below:

To find $E_{He}$ we use eq9, where

$E_{He}=2\varepsilon_i-\int_0^{\infty}\phi_i(x_i)\left[\int_0^{\infty}\frac{\phi_j^2(x_j)dx_j}{\vert x_i-x_j\vert+0.5}\right]\phi_i(x_i)dx_i\;\;\;\;\;\;\;\;(32)$

$\varepsilon_i$ is taken from the 12^th iteration set in the table above, while the double integral in eq32 is evaluated using form 2 of eq25. We being by computing the inner integral $\int_0^{\infty}\frac{\phi_j^2(x_j)dx_j}{\vert x_i-x_j\vert+0.5}$ as per step 3, using $N=4.838127$ and $b=-1.802042$ . Next, we apply the Simpson rule a second time for the outer integral. The Excel table is as follows:

The formula for cell NA4 is =(4.838127^2)*(B4^2)*EXP(-1.802042*2*B4))*MQ4. With reference to form 2 of eq25, the formula in NC4 is =NA4*NB4*7/350/3, representing the double integral for $x_i=0$ . The formulae for the remaining cells within NA4:NB354 are populated accordingly. Summing NC4:NC354, we obtain the total value of the double integral of 1.1386946, which when substituted in eq32 together with $\varepsilon_i=-0.854799$ gives $E_{He}=-2.848293$ . This theoretical value has a deviation of about 1.91% versus the experiment data of the ground state of helium.

Next article: The angle between two vectors in spherical coordinates

Previous article: Slater-type orbitals

Content page of computational chemistry

Content page of advanced chemistry

Main content page

December 7, 2022April 18, 2025

Slater-type orbitals

Slater-type orbitals (STO) are mathematical functions that resemble hydrogenic wavefunctions. Introduced by John Slater in 1930, they are used as trial wavefunctions in computational chemistry to approximate the energies of molecular systems. The general form of an STO is:

$\psi(r,\theta,\phi)=Nr^{n-1}e^{-\zeta r}Y_{l,ml}(\theta,\phi)\;\;\;\;\;\;\;\;(26)$

where

N is a normalisation constant

n, $l$ and $m_l$ are quantum numbers

$r$ is the distance between an electron and the nucleus of an atom.

$\zeta$ is the effective charge of the nucleus

$Y_{l,ml}(\theta,\phi)$ represents the spherical harmonics

STOs with different n but the same $l$ and $m_l$ are not orthogonal to one another. When working with such STOs, the Gram-Schmidt process is used to convert the non-orthogonal orbitals to orthogonal ones.

Question

Show that the normalisation constant of a 1s orbital is $\sqrt\frac{\zeta^3}{\pi}$ .

Answer

$N^2\int_0^{2\pi}d\phi\int_0^{\pi}sin\theta d\theta\int_0^{\infty}e^{-2\zeta r}r^2dr=1$
Using the identity $\int_0^{\infty}x^ne^{-\alpha x}dx=\frac{n!}{\alpha^{n+1}}$ , we have $\sqrt\frac{\zeta^3}{\pi}$ .

Question

Using the identity $\int_0^{\infty}x^ne^{-\alpha x}dx=\frac{n!}{\alpha^{n+1}}$ (see this article for proof), show that the 1s orbitals are orthogonal to each other but not orthogonal to the 2s orbital.

Answer

$\frac{\zeta_1^3}{\pi}\int e^{-\zeta_1r}e^{-\zeta_1r}d\tau=1$
$\sqrt\frac{\zeta_1^3\zeta_3^5}{3\pi^2}\int e^{-\zeta_1r}re^{-\zeta_3r}d\tau=\frac{\sqrt{192\zeta_1^3\zeta_3^5}}{(\zeta_1+\zeta_3)^4}$

For the orbitals to be orthogonal to each other, either $\zeta_1$ or $\zeta_3$ needs to be zero, which is forbidden because the effective nuclear charge is not zero. Hence, the orbitals are not orthogonal to each other.

Next article: Numerical solution of the Hartree self-consistent method for He

Previous article: Simpson’s rule