2.53: An Interactive SCF Calculation for the Helium Atom

This is an outline of the Mathcad implementation of a two‐electron SCF calculation published in JCE by Snow and Bills. [Snow, R. L.; Bills, J. L. J. Chem. Educ. 1975, 52, 506.]

Under the orbital approximation, Ψ(1,2) = Φ(1)Φ(2), the two‐electron Schrödinger equation for helium can be decoupled into two, one‐electron Hartree differential equations.The individual helium atom electrons are assumed to occupy an orbital which is a linear combination of two Slater 1s orbitals,

\[ \Phi = C_1 f_1 + C_2 f_2 \nonumber \]

\[ \begin{matrix} f_1 = \sqrt{ \frac{ \alpha^3}{ \pi}} \text{exp} ( - \alpha r) & \text{and} & f_2 = \sqrt{ \frac{ \beta^3}{ \pi}} \text{exp} ( - \beta r) \end{matrix} \nonumber \]

The effective Hamiltonian for the i^th electron is,

\[H_i = - \frac{1}{2r^2} \frac{d}{dr} \left( r^2 \frac{d}{dr} \blacksquare \right) - \frac{2}{r_i} + \int_0^{ \infty} \left( C_{j1} f_1 + C_{j2} f_2 \right) \frac{1}{r_{ij} } \left( C_{j1} f_1 + C_{j2} f_2 \right) d \tau \nonumber \]

where the Cjs are the coefficients of the j^th electron.

Assuming that the i^th electron is also in an orbital of the form given in equation (1), the variational method yields the following expression for the energy of the i^th electron.

\[ \varepsilon_i = \frac{ \int_0^{ \infty} \left( C_{j1} f_1 + C_{j2} f_2 \right) H_i \left( C_{j1} f_1 + C_{j2} f_2 \right) d \tau}{ \int_0^{ \infty} \left( C_{j1} f_1 + C_{j2} f_2 \right)^2 d \tau} = \frac{Ci_1^2 + H_{11} 2 Ci_i Ci_2 H_{12} + Ci_2^2 H_{22}}{Ci_1^2 + 2Ci_1 Ci_2 S_{12} + Ci_2^2} \nonumber \]

The optimum values for the orbital scale factors are given below. The user can change these values to demonstrate that they do indeed yield a minimum energy for the trial wavefunction chosen.

\[ \begin{matrix} \alpha = 1.45 & \beta = 2.90 \end{matrix} \nonumber \]

Evaluation of integrals given values for the scale factors, α and β, of the Slater type orbitals follows (See Snow and Bills for details):

\[ \begin{matrix} \tau = \frac{ \alpha - \beta}{ \alpha + \beta} & \tau 1 = \frac{ \alpha - \beta}{3 \alpha + \beta} & \tau 2 = \frac{ \alpha - \beta}{ \alpha + 3 \beta} \end{matrix} \nonumber \]

Kinetic energy integrals:

\[ \begin{matrix} T_{11} = \frac{ \alpha^2}{2} & T_{22} = \frac{ \beta^2}{2} & T_{12} = \frac{1}{8} \left( \alpha + \beta \right)^2 \left( 1 - \tau^2 \right)^{2.5} \end{matrix} \nonumber \]

Electron‐nucleus potential energy integrals:

\[ \begin{matrix} V_{11} = -2 \alpha & V_{22} = -2 \beta & V_{12} = - ( \alpha + \beta ) \left( 1 - \tau \right)^{1.5} \end{matrix} \nonumber \]

Electron‐electron repulsion integrals:

\[ \begin{matrix} V_{1111} = \frac{5}{8} \alpha & V_{2222} = \frac{5}{8} \alpha & V_{1212} = \frac{5}{16} \left( 1 - \tau^2 \right)^3 ( \alpha + \beta ) \end{matrix} \nonumber \]

\[ \begin{matrix} V_{1122} = \frac{1}{16} \left( 1 - \tau^2 \right) \left( 5 - \tau^2 \right) ( \alpha + \beta ) & V_{1222} = \frac{1}{32} \left( 1 - \tau^2 \right)^{1.5} \left( 1 - tau2^2 \right) \left( 5 - \tau 2^2 \right) ( \alpha + 3 \beta ) \\ V_{1112} = \frac{1}{32} \left( 1 - \tau^2 \right)^{1.5} \left( 1 - \tau 1^2 \right) \left( 5 - \tau 1^2 \right) ( 3 \alpha + \beta ) & \text{Overlap integral: } S_{12} = \left( 1 - \tau^2\ \right)^{1.5} \end{matrix} \nonumber \]

Having evaluated the integrals, the next step is the calculation of the matrix elements that appear in equation (4).

\[ \begin{matrix} H_{11} = T_{11} + V_{11} + Cj_1^2 V_{1111} + 2Cj_1 Cj_2 V_{1112} + Cj_2^2 V_{1122} \\ H_{12} = T_{12} + V_{12} + Cj_1^2 V_{1112} + 2Cj_1 Cj_2 V_{1212} + Cj_2^2 V_{1222} \\ H_{22} = T_{22} + V_{22} + Cj_1^2 V_{1122} + 2Cj_1 Cj_2 V_{1222} + Cj_2^2 V_{2222} \end{matrix} \nonumber \]

Given initial values for the coefficients of the jth electron (these are declared below in a more convient location using Mathcadʹs global equal sign), minimization of the orbital energy simultaneously with respect to the coefficients of the ith electron yields the orbital energy of the ith electron and its coefficients. These output coefficients become the input coefficients of the next iteration when the orbital energy of the j th electron is calculated. The procedure is repeated until self‐consistency is achieved. This occurs when the output coefficients are equal to the input coefficients, or to put it another way, when the coefficients of the two electrons are identical.

In the numeric mode Mathcad requires seed values for all variable which appear in the expression to be evaluated in a Given/Find solve block. The seed values for the coefficients shown below are arbitrarily set at 0.5.

\[ \begin{matrix} Ci_1 = .5 & Ci_2 = .5 \end{matrix} \nonumber \]

The variational integral for the electron orbital energy:

\[ \varepsilon \left( Ci_1,~ Ci_2 \right) = \frac{Ci_1^2 H_{11} + 2Ci_1 H_{12} + Ci_2^2 H_{22}}{Ci_1^2 + 2Ci_1 Ci_2 S_{12} + Ci_2^2} \nonumber \]

Minimization of the variational integral simultaneously with respect to the coefficients plus the normalization condition

\[ \begin{matrix} \text{Given} & \frac{d}{dCi_1} \varepsilon \left( Ci_1,~Ci_2 \right) = 0 & \frac{d}{dCi_2} \varepsilon \left( Ci_1,~Ci_2 \right) = 0 & Ci_1^2 + 2 Ci_1 Ci_2 S_{12} + Ci_2^2 = 1 \end{matrix} \nonumber \]

yields the output coefficients and the orbital energy:

\[ \begin{matrix} \begin{pmatrix} Ci_1 \\ Ci_2 \end{pmatrix} = \text{Find} \left( Ci_1,~Ci_2 \right) & \begin{pmatrix} Ci_1 \\ Ci_2 \end{pmatrix} = \begin{pmatrix} 0.809249 \\ 0.219060 \end{pmatrix} \end{matrix} \nonumber \]

In the SCF method the output coefficients become the input coefficients in the next iteration.

Display the orbital energy:

\[ \varepsilon \left( Ci_1,~Ci_2 \right) = -0.984326 \nonumber \]

The energy of the atom at this point in the calculation is the orbital energy of the ith electron plus the kinetic and nuclear potential energy of the jth electron.

\[ \begin{matrix} E_{atom} = \varepsilon \left( Ci_1,~Ci_2 \right) + Cj_1^2 \left( T_{11} + V_{11} \right) + 2C_1 Cj_2 \left( T_{12} + V_{12} \right) + Cj_2^2 \left( T_{22} + V_{22} \right) & E_{atom} = -2.833076 \end{matrix} \nonumber \]

Input coefficients:

\[ \begin{pmatrix} Cj_1 \\ Cj_2 \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \end{pmatrix} \nonumber \]

Paste new input coefficients after each iteration.

Summary of the SCF results for the following initial input coefficients:

\[ \begin{pmatrix} Cj_1 \\ Cj_2 \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \end{pmatrix} \nonumber \]

\[ \begin{pmatrix} \text{Iteration} & Ci_1 & Ci_2 & \varepsilon \left( Ci_1,~Ci_2 \right) & E_{atom} \\ 1 & 0.809249 & 0.21906 & -0.984326 & -2.833076 \\ 2 & 0.847034 & 0.176951 & -0.90556 & -2.860616 \\ 3 & 0.839638 & 0.185241 & -0.920653 & -2.86163 \\ 4 & 0.841091 & 0.183615 & -0.917676 & -2.86167 \\ 5 & 0.840806 & 0.183934 & -0.918258 & -2.861672 \\ 6 & 0.840862 & 0.183871 & -0.918144 & -2.861672 \\ 7 & 0.840851 & 0.183884 & -0.918167 & -2.861672 \\ 8 & 0.840852 & 0.183881 & -0.918161 & -2.861673 \\ 9 & 0.840852 & 0.183882 & -0.918164 & -2.861672 \\ 10 & 0.840852 & 0.183882 & -0.918164 & -2.861672 \end{pmatrix} \nonumber \]