# Restricted Hartree–Fock for Beryllium

Last updated: May 16th, 2019

## Restricted Hartree–Fock for Beryllium

$2018/12/07$

In [1]:
import sympy as sp
from sympy import oo
import numpy as np
from itertools import product
from scipy.linalg import eig
from sympy import diff
import time as time
import matplotlib.pyplot as plt
from sympy.plotting import plot
# %matplotlib notebook
%matplotlib inline
from IPython.display import Math
sp.init_printing()

r, r1, r2, zeta, zeta1, zeta2 = sp.symbols("r, r1, r2, zeta, zeta1, zeta2")
n = sp.Symbol('n',integer=True)


## Part 1 Function Definition

$STO$$\; function \; has \; format \; of: \; Nr^{n-1}e^{-r\zeta} In [2]: def STO(zeta, n, r=r): return (2*zeta)**n*(2*zeta/sp.factorial(2*n))**(1/2)*r**(n-1)*sp.exp(-zeta*r) print("For example: STO for 1s") STO(zeta, 1)  For example: STO for 1s  Out[2]:$$2 \zeta^{1.5} e^{- r \zeta}$$### Overlap Integrate¶$S = \int_{0}^\infty f_1^* f_2 \; r^2dr$In [3]: # S Overlap Integrate def S_int(f1, f2): return sp.integrate(f1*f2*r*r ,(r, 0, +oo))  ### Hamiltonian core¶ H core = kinetics energy + electron and nuclear potential energy$H = \int_{0}^\infty f_1 \hat{H} f_2 \; r^2drH = \int_{0}^\infty f_1 ((-\dfrac{1}{2}) \nabla^2 - \dfrac{Z_{\alpha}}{r})f_2 \; r^2 drLaplace \; operator: \nabla^2\therefore H = \int_{0}^\infty f_1 ((-\dfrac{1}{2}) \dfrac{1}{r} \dfrac{\partial}{\partial r} \dfrac{\partial}{\partial r} r f_2 - \dfrac{Z_{\alpha}}{r}r^2 f_2 )dr$reminder: change the Z of nucleus if you want to run for other atom In [4]: # H core = kinetics energy + electron and nuclear potential energy def H_int(f1, f2, Z): return sp.integrate(f1*(-((1/2)*(1/r)*diff(diff(r*f2, r), r))-((Z/r)*f2))*r*r, (r,0,+oo))  In [5]: # Returns the core hamiltonian matrix def H_matrix(fs, Z): H = np.zeros((len(fs),len(fs))) for i in range(len(fs)): for j in range(len(fs)): H[i, j] = H_int(fs[i], fs[j], Z) return H # Returns the overlap matrix def S_matrix(fs): S = np.zeros((len(fs),len(fs))) for i in range(len(fs)): for j in range(len(fs)): S[i, j] = S_int(fs[i], fs[j]) return S  ### Two elecron repulsion integral¶$(rs|tu) = \int \int \dfrac{f_r^*(1) f_s(1) f_t^*(2) f_u(2)}{r_{12}} \; dv_1 dv_2 $For 1s 2s orbital$(rs|tu) = \int_{0}^\infty \int_{0}^\infty \dfrac{f_r^*(1) f_s(1) f_t^*(2) f_u(2)}{r_{12}} \; r_1^2dr_1\; r_2^2dr_2 (rs|tu) = \int_{0}^\infty f_r^*(1) f_s(1) \; r_1^2dr_1\int_{0}^\infty \frac{ f_t^*(2) f_u(2)}{r_{12}}\; r_2^2dr_2 $From problem 9.14 in quantum_chemistry by levine$(rs|tu) = \int_{0}^\infty f_r^*(1) f_s(1) \; r_1^2dr_1\int_{0}^\infty \frac{ f_t^*(2) f_u(2)}{r_{>}}\; r_2^2dr_2 (rs|tu) = \int_{0}^\infty f_r^*(1) f_s(1) \; r_1^2dr_1(\int_{0}^{r_1} \frac{ f_t^*(2) f_u(2)}{r_{1}}\; r_2^2dr_2 + \int_{r_1}^\infty \frac{ f_t^*(2) f_u(2)}{r_{2}}\; r_2^2dr_2) Let \; B= \int_{0}^{r_1} \frac{ f_t^*(2) f_u(2)}{r_{1}}\; r_2^2dr_2 + \int_{r_1}^\infty \frac{ f_t^*(2) f_u(2)}{r_{2}}\; r_2^2dr_2(rs|tu) = \int_{0}^\infty f_r^*(1) f_s(1) B \; r_1^2 dr_1 $In [6]: def Repulsion_electron(zetas): f1=STO(zetas[0][0], zetas[0][1], r1) f2=STO(zetas[1][0], zetas[1][1], r1) f3=STO(zetas[2][0], zetas[2][1], r2) f4=STO(zetas[3][0], zetas[3][1], r2) fs = [f1, f2, f3, f4] B = (1/r1)*sp.integrate(f3*f4*r2*r2 ,(r2, 0, r1)) + sp.integrate((1/r2)*f3*f4*r2*r2 ,(r2, r1, +oo)) A = sp.integrate(f1*f2*r1*r1*B ,(r1, 0, +oo)) return A  ### Density matrix¶$ P_{tu} =2 \sum_{j=1}^{n/2}c_{tj}^* c_{uj} $Reminder: P need to be changed if the atom have unpaired electron In [7]: # Calculates Density matrix def P_matrix(Co): P = np.zeros([Co.shape[0], Co.shape[0]]) for t in range(Co.shape[0]): for u in range(Co.shape[0]): for j in range(int(Co.shape[0]/2)): P[t][u] += 2* Co[t][j]*Co[u][j] return P  ### Fock matrix¶$ F_{rs} = H_{rs}^{core} + \sum_{t=1}^{b} \sum_{t=1}^{b}P_{tu}[(rs|tu)- \frac{1}2(ru|ts)] G = \sum_{t=1}^{b} \sum_{t=1}^{b}P_{tu}[(rs|tu)- \frac{1}2(ru|ts)] F_{rs} = H_{rs}^{core} + G In\;G \;one\;is\;coulombic\;repulsion,\;another\;is\;exchange\;energy$In [8]: def R_matrix(zetas): R = np.zeros((len(zetas), len(zetas), len(zetas), len(zetas))) rs = list(product(range(len(zetas)),repeat=2)) tu = list(product(range(len(zetas)),repeat=2)) for r, s in rs: for t, u in tu: R[r,s,t,u] = Repulsion_electron((zetas[r], zetas[s], zetas[t], zetas[u])) return R # Caculate G Matrix def G_matrix(zetas, Co, R): G = np.zeros((Co.shape[0], Co.shape[0])) P = P_matrix(Co) rs = list(product(range( Co.shape[0]),repeat=2)) tu = list(product(range( Co.shape[0]),repeat=2)) for r, s in rs: g = 0 for t, u in tu: int1 = R[r, s, t, u] int2 = R[r, u, t, s] # print('({0}{1}|{2}{3}): {4}'.format(r, s, t, u, int1)) g+= P[t, u] * (int1 - 0.5 * int2) G[r, s] = g return G # Returns the Fock matrix def F_matrix(fs, Z, zetas, Co, R): return H_matrix(fs, Z) + G_matrix(zetas, Co, R)  ### Solve Hartree-Fork equation¶$det(F_{rs}-\epsilon_i S_{rs} = 0)The\;energy\;returned\;is\;the\;orbital\;energy\;for\;1\;electron$In [ ]: # slove secular equation, return the energy and improved coeffients # the energy here is orbital energy for 1 electron def secular_eqn(F, S): ei, C = eig(F, S) # sort eigvalue and eigvector from lower to higher idx = ei.argsort()[::1] ei = ei[idx] C = C[:,idx] # eigvector from scipy.linalg.eig is not normalized, which is a bug # this is to fix it Co = np.zeros((C.shape[0],C.shape[0])) inte = np.matmul(np.matmul(C.T, S), C) for i in range(C.shape[0]): for j in range(C.shape[0]): Co[j][i]=C[j][i]/np.sqrt(inte[i][i]) return ei, Co  ### Return atom energy¶$ E_{HF} = \sum_{i=1}^{2/n}\epsilon +\frac{1}2 \sum_{r=1}^{b} \sum_{s=1}^{b}P_{rs}H_{rs}+V_{NN} $In [12]: # return energy of atom def get_E0(e, P, H): E0 = 0 for i in range(int(e.shape[0]/2)): E0 += e[i].real E0 = E0 + 0.5*(P*H).sum() return E0  In [13]: # input # input for zeta zetas = [[5.59108, 1], [3.35538, 1], [1.01122, 2], [0.61000, 2]] # input nuclear charge (element number) Z = 4 # build basis function f1=STO(zetas[0][0], zetas[0][1]) f2=STO(zetas[1][0], zetas[1][1]) f3=STO(zetas[2][0], zetas[2][1]) f4=STO(zetas[3][0], zetas[3][1]) fs = [f1, f2, f3, f4] S = S_matrix(fs) S  Out[13]: array([[1. , 0.90780473, 0.09913506, 0.03600389], [0.90780473, 1. , 0.24088162, 0.10010333], [0.09913506, 0.24088162, 1. , 0.85384478], [0.03600389, 0.10010333, 0.85384478, 1. ]]) ## Part 2 Hartree Fork Iteration (Main Process) • Initialization 1. initializing Co (coefficients) without considering electron repulsion 2. Solve Hartree-Fork equation with H_matrix and S_matrix to get initial Co • Iteration 1. Using Co, we can get P_matrix (electron density) 2. Using P_matrix, H_matrix, G_matrix => F_matrix 3. Solve Hartree-Fork equation with F_matrix and S_matrix to get improved orbital energy and Co, which also means improved orbital functions. 4. Using improved Co, return to step 1 In [14]: # input # input for zeta zetas = [[5.59108, 1], [3.35538, 1], [1.01122, 2], [0.61000, 2]] # input nuclear charge (element number) Z = 4 # build basis function f1=STO(zetas[0][0], zetas[0][1]) f2=STO(zetas[1][0], zetas[1][1]) f3=STO(zetas[2][0], zetas[2][1]) f4=STO(zetas[3][0], zetas[3][1]) fs = [f1, f2, f3, f4] # initialization R = np.zeros((len(zetas), len(zetas), len(zetas), len(zetas))) H = H_matrix(fs, Z) S = S_matrix(fs) e, Co = secular_eqn(H, S) P = P_matrix(Co) scf_H = get_E0(e, P, H) ##############################################print information below################################################# print('-'*30, "Initialization", '-'*30) print('-'*25, "Ignore repulsion integral", '-'*24) display(Math('\zeta_1 = {0} \quad \zeta_2 = {1} \quad \zeta_3 = {2} \quad \zeta_4 = {3}'.format(format(zetas[0][0], '0.3f'), format(zetas[1][0], '0.3f'), format(zetas[2][0], '0.3f'),format(zetas[3][0], '0.3f')))) display(Math('Orbitals:')) display(Math(' \phi_1 = c_{11} \chi_1 + c_{21} \chi_2 + c_{31} \chi_3 + c_{41} \chi_4')) display(Math(' \phi_2 = c_{12} \chi_1 + c_{22} \chi_2 + c_{32} \chi_3 + c_{42} \chi_4')) display(Math('c11 = {0} \quad c21 = {1} \quad c31 = {2} \quad c41 = {3}'.format(format(Co[0][0], '0.3f'), format(Co[1][0], '0.3f'), format(Co[2][0], '0.3f'), format(Co[3][0], '0.3f')))) display(Math('c12 = {0} \quad c22 = {1} \quad c32 = {2} \quad c42 = {3}'.format(format(Co[0][1], '0.3f'), format(Co[1][1], '0.3f'), format(Co[2][1], '0.3f'), format(Co[3][1], '0.3f')))) # plot density graph colorlist = ['red', 'orange', 'yellow', 'green', 'blue', 'purple', 'black', 'red', 'orange', 'yellow', 'green', 'blue', 'purple', 'black'] phi1 = Co[0,0]*f1+Co[1,0]*f2+Co[2,0]*f3+Co[3,0]*f4 phi2 = Co[0,1]*f1+Co[1,1]*f2+Co[2,1]*f3+Co[3,1]*f4 density_1 = phi1*phi1*r*r density_2 = phi2*phi2*r*r p = plot((density_1, (r, 0, 5)), (density_2, (r, 0, 5)), show = False, legend = True) p[0].label = 'electron density$r^2 \phi_1^2$' p[1].label = 'electron density$r^2 \phi_2^2$' p[0].line_color = colorlist[0] p[1].line_color = 'blue' p.show() # print energy result display(Math(' \epsilon_1 \; for \; \phi_1 = {0} '.format(format(e[0].real, '0.3f')))) display(Math(' \epsilon_2 \; for \; \phi_1 = {0} '.format(format(e[1].real, '0.3f')))) display(Math(' Hartree \ Fork \; atom \; energy = {0} \ hartree = {1} \ eV'.format(format(scf_H, '0.5f'), format(scf_H*27.211, '0.5f')))) ##############################################print information above################################################# for i in range(10): print('-'*30, "Iteration", i + 1, '-'*30) if(i==0): print('-'*7, "Iteration 1 needs more time to caculate Repulsion Integral", '-'*6) start = time.time() R = R_matrix(zetas) else: start = time.time() F = F_matrix(fs, Z, zetas, Co, R) S = S_matrix(fs) e, Co = secular_eqn(F, S) P = P_matrix(Co) scf_H = get_E0(e, P, H) ##########################################print information below################################################# # print information display(Math('\zeta_1 = {0} \quad \zeta_2 = {1} \quad \zeta_3 = {2} \quad \zeta_4 = {3}'.format(format(zetas[0][0], '0.3f'), format(zetas[1][0], '0.3f'), format(zetas[2][0], '0.3f'),format(zetas[3][0], '0.3f')))) display(Math('Orbitals:')) display(Math(' \phi_1 = c_{11} \chi_1 + c_{21} \chi_2 + c_{31} \chi_3 + c_{41} \chi_4')) display(Math(' \phi_2 = c_{12} \chi_1 + c_{22} \chi_2 + c_{32} \chi_3 + c_{42} \chi_4')) display(Math('c11 = {0} \quad c21 = {1} \quad c31 = {2} \quad c41 = {3}'.format(format(Co[0][0], '0.3f'), format(Co[1][0], '0.3f'), format(Co[2][0], '0.3f'), format(Co[3][0], '0.3f')))) display(Math('c12 = {0} \quad c22 = {1} \quad c32 = {2} \quad c42 = {3}'.format(format(Co[0][1], '0.3f'), format(Co[1][1], '0.3f'), format(Co[2][1], '0.3f'), format(Co[3][1], '0.3f')))) # plot density graph phi1 = Co[0,0]*f1+Co[1,0]*f2+Co[2,0]*f3+Co[3,0]*f4 phi2 = Co[0,1]*f1+Co[1,1]*f2+Co[2,1]*f3+Co[3,1]*f4 density_1 = phi1*phi1*r*r density_2 = phi2*phi2*r*r p1 = plot((density_1, (r, 0, 5)), (density_2, (r, 0, 5)), show = False, legend = True) p1[0].label = None p1[1].label = None p1[0].line_color = colorlist[i+1] p1[1].line_color = colorlist[i+1] p.extend(p1) p.show() # print energy result display(Math(' \epsilon_1 \; for \; \phi_1 = {0} '.format(format(e[0].real, '0.3f')))) display(Math(' \epsilon_2 \; for \; \phi_1 = {0} '.format(format(e[1].real, '0.3f')))) display(Math(' Hartree \ Fork \; atom \; energy = {0} \ hartree = {1} \ eV'.format(format(scf_H, '0.5f'), format(scf_H*27.211, '0.5f')))) stop = time.time() print('Time used:',format(stop-start, '0.1f'),'s') ##########################################print information above#################################################  ------------------------------ Initialization ------------------------------ ------------------------- Ignore repulsion integral ------------------------ $\displaystyle \zeta_1 = 5.591 \quad \zeta_2 = 3.355 \quad \zeta_3 = 1.011 \quad \zeta_4 = 0.610\displaystyle Orbitals:\displaystyle \phi_1 = c_{11} \chi_1 + c_{21} \chi_2 + c_{31} \chi_3 + c_{41} \chi_4\displaystyle \phi_2 = c_{12} \chi_1 + c_{22} \chi_2 + c_{32} \chi_3 + c_{42} \chi_4\displaystyle c11 = -0.275 \quad c21 = -0.751 \quad c31 = 0.052 \quad c41 = -0.031\displaystyle c12 = 0.127 \quad c22 = 0.113 \quad c32 = -1.542 \quad c42 = 0.683\displaystyle \epsilon_1 \; for \; \phi_1 = -7.985 \displaystyle \epsilon_2 \; for \; \phi_1 = -1.774 \displaystyle Hartree \ Fork \; atom \; energy = -19.51846 \ hartree = -531.11686 \ eV$------------------------------ Iteration 1 ------------------------------ ------- Iteration 1 needs more time to caculate Repulsion Integral ------ $\displaystyle \zeta_1 = 5.591 \quad \zeta_2 = 3.355 \quad \zeta_3 = 1.011 \quad \zeta_4 = 0.610\displaystyle Orbitals:\displaystyle \phi_1 = c_{11} \chi_1 + c_{21} \chi_2 + c_{31} \chi_3 + c_{41} \chi_4\displaystyle \phi_2 = c_{12} \chi_1 + c_{22} \chi_2 + c_{32} \chi_3 + c_{42} \chi_4\displaystyle c11 = 0.156 \quad c21 = 0.856 \quad c31 = 0.003 \quad c41 = -0.001\displaystyle c12 = -0.004 \quad c22 = 0.252 \quad c32 = -1.173 \quad c42 = 0.166\displaystyle \epsilon_1 \; for \; \phi_1 = -4.422 \displaystyle \epsilon_2 \; for \; \phi_1 = -0.245 \displaystyle Hartree \ Fork \; atom \; energy = -14.28744 \ hartree = -388.77544 \ eV$Time used: 57.4 s ------------------------------ Iteration 2 ------------------------------ $\displaystyle \zeta_1 = 5.591 \quad \zeta_2 = 3.355 \quad \zeta_3 = 1.011 \quad \zeta_4 = 0.610\displaystyle Orbitals:\displaystyle \phi_1 = c_{11} \chi_1 + c_{21} \chi_2 + c_{31} \chi_3 + c_{41} \chi_4\displaystyle \phi_2 = c_{12} \chi_1 + c_{22} \chi_2 + c_{32} \chi_3 + c_{42} \chi_4\displaystyle c11 = 0.155 \quad c21 = 0.858 \quad c31 = -0.002 \quad c41 = 0.001\displaystyle c12 = -0.005 \quad c22 = 0.230 \quad c32 = -1.013 \quad c42 = -0.020\displaystyle \epsilon_1 \; for \; \phi_1 = -4.654 \displaystyle \epsilon_2 \; for \; \phi_1 = -0.290 \displaystyle Hartree \ Fork \; atom \; energy = -14.50425 \ hartree = -394.67518 \ eV$Time used: 1.2 s ------------------------------ Iteration 3 ------------------------------ $\displaystyle \zeta_1 = 5.591 \quad \zeta_2 = 3.355 \quad \zeta_3 = 1.011 \quad \zeta_4 = 0.610\displaystyle Orbitals:\displaystyle \phi_1 = c_{11} \chi_1 + c_{21} \chi_2 + c_{31} \chi_3 + c_{41} \chi_4\displaystyle \phi_2 = c_{12} \chi_1 + c_{22} \chi_2 + c_{32} \chi_3 + c_{42} \chi_4\displaystyle c11 = 0.155 \quad c21 = 0.858 \quad c31 = -0.003 \quad c41 = 0.002\displaystyle c12 = -0.004 \quad c22 = 0.223 \quad c32 = -0.963 \quad c42 = -0.076\displaystyle \epsilon_1 \; for \; \phi_1 = -4.708 \displaystyle \epsilon_2 \; for \; \phi_1 = -0.304 \displaystyle Hartree \ Fork \; atom \; energy = -14.55105 \ hartree = -395.94867 \ eV$Time used: 1.1 s ------------------------------ Iteration 4 ------------------------------ $\displaystyle \zeta_1 = 5.591 \quad \zeta_2 = 3.355 \quad \zeta_3 = 1.011 \quad \zeta_4 = 0.610\displaystyle Orbitals:\displaystyle \phi_1 = c_{11} \chi_1 + c_{21} \chi_2 + c_{31} \chi_3 + c_{41} \chi_4\displaystyle \phi_2 = c_{12} \chi_1 + c_{22} \chi_2 + c_{32} \chi_3 + c_{42} \chi_4\displaystyle c11 = 0.155 \quad c21 = 0.858 \quad c31 = -0.003 \quad c41 = 0.002\displaystyle c12 = -0.003 \quad c22 = 0.221 \quad c32 = -0.948 \quad c42 = -0.093\displaystyle \epsilon_1 \; for \; \phi_1 = -4.726 \displaystyle \epsilon_2 \; for \; \phi_1 = -0.308 \displaystyle Hartree \ Fork \; atom \; energy = -14.56616 \ hartree = -396.35976 \ eV$Time used: 1.3 s ------------------------------ Iteration 5 ------------------------------ $\displaystyle \zeta_1 = 5.591 \quad \zeta_2 = 3.355 \quad \zeta_3 = 1.011 \quad \zeta_4 = 0.610\displaystyle Orbitals:\displaystyle \phi_1 = c_{11} \chi_1 + c_{21} \chi_2 + c_{31} \chi_3 + c_{41} \chi_4\displaystyle \phi_2 = c_{12} \chi_1 + c_{22} \chi_2 + c_{32} \chi_3 + c_{42} \chi_4\displaystyle c11 = 0.155 \quad c21 = 0.858 \quad c31 = -0.003 \quad c41 = 0.002\displaystyle c12 = -0.003 \quad c22 = 0.220 \quad c32 = -0.944 \quad c42 = -0.097\displaystyle \epsilon_1 \; for \; \phi_1 = -4.731 \displaystyle \epsilon_2 \; for \; \phi_1 = -0.309 \displaystyle Hartree \ Fork \; atom \; energy = -14.57061 \ hartree = -396.48075 \ eV$Time used: 1.2 s ------------------------------ Iteration 6 ------------------------------ $\displaystyle \zeta_1 = 5.591 \quad \zeta_2 = 3.355 \quad \zeta_3 = 1.011 \quad \zeta_4 = 0.610\displaystyle Orbitals:\displaystyle \phi_1 = c_{11} \chi_1 + c_{21} \chi_2 + c_{31} \chi_3 + c_{41} \chi_4\displaystyle \phi_2 = c_{12} \chi_1 + c_{22} \chi_2 + c_{32} \chi_3 + c_{42} \chi_4\displaystyle c11 = 0.155 \quad c21 = 0.858 \quad c31 = -0.003 \quad c41 = 0.002\displaystyle c12 = -0.003 \quad c22 = 0.220 \quad c32 = -0.942 \quad c42 = -0.099\displaystyle \epsilon_1 \; for \; \phi_1 = -4.732 \displaystyle \epsilon_2 \; for \; \phi_1 = -0.309 \displaystyle Hartree \ Fork \; atom \; energy = -14.57187 \ hartree = -396.51521 \ eV$Time used: 1.3 s ------------------------------ Iteration 7 ------------------------------ $\displaystyle \zeta_1 = 5.591 \quad \zeta_2 = 3.355 \quad \zeta_3 = 1.011 \quad \zeta_4 = 0.610\displaystyle Orbitals:\displaystyle \phi_1 = c_{11} \chi_1 + c_{21} \chi_2 + c_{31} \chi_3 + c_{41} \chi_4\displaystyle \phi_2 = c_{12} \chi_1 + c_{22} \chi_2 + c_{32} \chi_3 + c_{42} \chi_4\displaystyle c11 = 0.155 \quad c21 = 0.858 \quad c31 = -0.003 \quad c41 = 0.002\displaystyle c12 = -0.003 \quad c22 = 0.220 \quad c32 = -0.942 \quad c42 = -0.099\displaystyle \epsilon_1 \; for \; \phi_1 = -4.733 \displaystyle \epsilon_2 \; for \; \phi_1 = -0.309 \displaystyle Hartree \ Fork \; atom \; energy = -14.57223 \ hartree = -396.52493 \ eV$Time used: 1.4 s ------------------------------ Iteration 8 ------------------------------ $\displaystyle \zeta_1 = 5.591 \quad \zeta_2 = 3.355 \quad \zeta_3 = 1.011 \quad \zeta_4 = 0.610\displaystyle Orbitals:\displaystyle \phi_1 = c_{11} \chi_1 + c_{21} \chi_2 + c_{31} \chi_3 + c_{41} \chi_4\displaystyle \phi_2 = c_{12} \chi_1 + c_{22} \chi_2 + c_{32} \chi_3 + c_{42} \chi_4\displaystyle c11 = 0.155 \quad c21 = 0.858 \quad c31 = -0.003 \quad c41 = 0.002\displaystyle c12 = -0.003 \quad c22 = 0.220 \quad c32 = -0.942 \quad c42 = -0.099\displaystyle \epsilon_1 \; for \; \phi_1 = -4.733 \displaystyle \epsilon_2 \; for \; \phi_1 = -0.309 \displaystyle Hartree \ Fork \; atom \; energy = -14.57233 \ hartree = -396.52766 \ eV$Time used: 1.3 s ------------------------------ Iteration 9 ------------------------------ $\displaystyle \zeta_1 = 5.591 \quad \zeta_2 = 3.355 \quad \zeta_3 = 1.011 \quad \zeta_4 = 0.610\displaystyle Orbitals:\displaystyle \phi_1 = c_{11} \chi_1 + c_{21} \chi_2 + c_{31} \chi_3 + c_{41} \chi_4\displaystyle \phi_2 = c_{12} \chi_1 + c_{22} \chi_2 + c_{32} \chi_3 + c_{42} \chi_4\displaystyle c11 = 0.155 \quad c21 = 0.858 \quad c31 = -0.003 \quad c41 = 0.002\displaystyle c12 = -0.003 \quad c22 = 0.220 \quad c32 = -0.942 \quad c42 = -0.099\displaystyle \epsilon_1 \; for \; \phi_1 = -4.733 \displaystyle \epsilon_2 \; for \; \phi_1 = -0.309 \displaystyle Hartree \ Fork \; atom \; energy = -14.57236 \ hartree = -396.52843 \ eV$Time used: 1.3 s ------------------------------ Iteration 10 ------------------------------ $\displaystyle \zeta_1 = 5.591 \quad \zeta_2 = 3.355 \quad \zeta_3 = 1.011 \quad \zeta_4 = 0.610\displaystyle Orbitals:\displaystyle \phi_1 = c_{11} \chi_1 + c_{21} \chi_2 + c_{31} \chi_3 + c_{41} \chi_4\displaystyle \phi_2 = c_{12} \chi_1 + c_{22} \chi_2 + c_{32} \chi_3 + c_{42} \chi_4\displaystyle c11 = 0.155 \quad c21 = 0.858 \quad c31 = -0.003 \quad c41 = 0.002\displaystyle c12 = -0.003 \quad c22 = 0.220 \quad c32 = -0.942 \quad c42 = -0.099\displaystyle \epsilon_1 \; for \; \phi_1 = -4.733 \displaystyle \epsilon_2 \; for \; \phi_1 = -0.309 \displaystyle Hartree \ Fork \; atom \; energy = -14.57237 \ hartree = -396.52865 \ eV\$
Time used: 1.3 s

In [ ]:



Reference

[1] Levine, Quantum Chemistry, 7th Edition, chapter 14
[2] C. Roetti and E. Clementi, J. Chem. Phys., 60, 4725 (1974)
[3] Acosta C R. Restricted closed shell Hartree Fock Roothaan matrix method applied to Helium atom using Mathematica[J]. European Journal of Physics Education, 2017, 5(1): 1-14.
[4] Simple Quantum Chemistry: Hartree-Fock in Python, 2018
[5] Helium hartree fork, prtkm, 2015

In [15]:
Co

Out[15]:
array([[ 1.55022874e-01, -3.03031822e-03,  4.06097130e-02,
2.56921068e+00],
[ 8.57664787e-01,  2.19993903e-01, -2.54907937e-01,
-2.52397487e+00],
[-3.09301538e-03, -9.41967018e-01,  1.79177857e+00,
7.40166137e-01],
[ 1.91274461e-03, -9.90329330e-02, -1.96939890e+00,
-4.71016053e-01]])
In [16]:
S

Out[16]:
array([[1.        , 0.90780473, 0.09913506, 0.03600389],
[0.90780473, 1.        , 0.24088162, 0.10010333],
[0.09913506, 0.24088162, 1.        , 0.85384478],
[0.03600389, 0.10010333, 0.85384478, 1.        ]])
In [20]:
Co.transpose() @ S @ Co

Out[20]:
array([[ 1.00000000e+00,  1.03614213e-16, -5.99612282e-16,
-1.85369662e-15],
[ 7.51265203e-17,  1.00000000e+00, -5.93392144e-16,
-2.31082821e-16],
[-5.34461285e-16, -6.03515972e-16,  1.00000000e+00,
1.95670965e-16],
[-1.45676379e-15, -1.97311684e-16,  1.47015765e-16,
1.00000000e+00]])
In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]: