Restricted Hartree–Fock for Beryllium

Source code

Helium /  Lithium /  Beryllium

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
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

Define STO function

$STO$$ \; function \; has \; format \; of: \; Nr^{n-1}e^{-r\zeta}$

1
2
3
4
5
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)

Overlap Integrate

$S = \int_{0}^\infty f_1^* f_2 \; r^2dr$

1
2
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^2dr$

$H = \int_{0}^\infty f_1 ((-\dfrac{1}{2}) \nabla^2 - \dfrac{Z_{\alpha}}{r})f_2 \; r^2 dr$

$Laplace \; 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

1
2
3
# 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))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# 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$

1
2
3
4
5
6
7
8
9
10
11
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

1
2
3
4
5
6
7
8
9
10
# 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$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
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$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# 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} $

1
2
3
4
5
6
7
8
# 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

Part 2 Hartree Fork Iteration

  • 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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
# 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#################################################
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
# output
------------------------------ Initialization ------------------------------
------------------------- Ignore repulsion integral ------------------------
𝜁1=5.591𝜁2=3.355𝜁3=1.011𝜁4=0.610
π‘‚π‘Ÿπ‘π‘–π‘‘π‘Žπ‘™π‘ :
πœ™1=𝑐11πœ’1+𝑐21πœ’2+𝑐31πœ’3+𝑐41πœ’4
πœ™2=𝑐12πœ’1+𝑐22πœ’2+𝑐32πœ’3+𝑐42πœ’4
𝑐11=βˆ’0.275𝑐21=βˆ’0.751𝑐31=0.052𝑐41=βˆ’0.031
𝑐12=0.127𝑐22=0.113𝑐32=βˆ’1.542𝑐42=0.683

πœ–1π‘“π‘œπ‘Ÿπœ™1=βˆ’7.985
πœ–2π‘“π‘œπ‘Ÿπœ™1=βˆ’1.774
π»π‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’ πΉπ‘œπ‘Ÿπ‘˜π‘Žπ‘‘π‘œπ‘šπ‘’π‘›π‘’π‘Ÿπ‘”π‘¦=βˆ’19.51846 β„Žπ‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’=βˆ’531.11686 𝑒𝑉
------------------------------ Iteration 1 ------------------------------
------- Iteration 1 needs more time to caculate Repulsion Integral ------
𝜁1=5.591𝜁2=3.355𝜁3=1.011𝜁4=0.610
π‘‚π‘Ÿπ‘π‘–π‘‘π‘Žπ‘™π‘ :
πœ™1=𝑐11πœ’1+𝑐21πœ’2+𝑐31πœ’3+𝑐41πœ’4
πœ™2=𝑐12πœ’1+𝑐22πœ’2+𝑐32πœ’3+𝑐42πœ’4
𝑐11=0.156𝑐21=0.856𝑐31=0.003𝑐41=βˆ’0.001
𝑐12=βˆ’0.004𝑐22=0.252𝑐32=βˆ’1.173𝑐42=0.166

πœ–1π‘“π‘œπ‘Ÿπœ™1=βˆ’4.422
πœ–2π‘“π‘œπ‘Ÿπœ™1=βˆ’0.245
π»π‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’ πΉπ‘œπ‘Ÿπ‘˜π‘Žπ‘‘π‘œπ‘šπ‘’π‘›π‘’π‘Ÿπ‘”π‘¦=βˆ’14.28744 β„Žπ‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’=βˆ’388.77544 𝑒𝑉
Time used: 65.3 s
------------------------------ Iteration 2 ------------------------------
𝜁1=5.591𝜁2=3.355𝜁3=1.011𝜁4=0.610
π‘‚π‘Ÿπ‘π‘–π‘‘π‘Žπ‘™π‘ :
πœ™1=𝑐11πœ’1+𝑐21πœ’2+𝑐31πœ’3+𝑐41πœ’4
πœ™2=𝑐12πœ’1+𝑐22πœ’2+𝑐32πœ’3+𝑐42πœ’4
𝑐11=0.155𝑐21=0.858𝑐31=βˆ’0.002𝑐41=0.001
𝑐12=βˆ’0.005𝑐22=0.230𝑐32=βˆ’1.013𝑐42=βˆ’0.020

πœ–1π‘“π‘œπ‘Ÿπœ™1=βˆ’4.654
πœ–2π‘“π‘œπ‘Ÿπœ™1=βˆ’0.290
π»π‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’ πΉπ‘œπ‘Ÿπ‘˜π‘Žπ‘‘π‘œπ‘šπ‘’π‘›π‘’π‘Ÿπ‘”π‘¦=βˆ’14.50425 β„Žπ‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’=βˆ’394.67518 𝑒𝑉
Time used: 1.0 s
------------------------------ Iteration 3 ------------------------------
𝜁1=5.591𝜁2=3.355𝜁3=1.011𝜁4=0.610
π‘‚π‘Ÿπ‘π‘–π‘‘π‘Žπ‘™π‘ :
πœ™1=𝑐11πœ’1+𝑐21πœ’2+𝑐31πœ’3+𝑐41πœ’4
πœ™2=𝑐12πœ’1+𝑐22πœ’2+𝑐32πœ’3+𝑐42πœ’4
𝑐11=0.155𝑐21=0.858𝑐31=βˆ’0.003𝑐41=0.002
𝑐12=βˆ’0.004𝑐22=0.223𝑐32=βˆ’0.963𝑐42=βˆ’0.076

πœ–1π‘“π‘œπ‘Ÿπœ™1=βˆ’4.708
πœ–2π‘“π‘œπ‘Ÿπœ™1=βˆ’0.304
π»π‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’ πΉπ‘œπ‘Ÿπ‘˜π‘Žπ‘‘π‘œπ‘šπ‘’π‘›π‘’π‘Ÿπ‘”π‘¦=βˆ’14.55105 β„Žπ‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’=βˆ’395.94867 𝑒𝑉
Time used: 1.0 s
------------------------------ Iteration 4 ------------------------------
𝜁1=5.591𝜁2=3.355𝜁3=1.011𝜁4=0.610
π‘‚π‘Ÿπ‘π‘–π‘‘π‘Žπ‘™π‘ :
πœ™1=𝑐11πœ’1+𝑐21πœ’2+𝑐31πœ’3+𝑐41πœ’4
πœ™2=𝑐12πœ’1+𝑐22πœ’2+𝑐32πœ’3+𝑐42πœ’4
𝑐11=0.155𝑐21=0.858𝑐31=βˆ’0.003𝑐41=0.002
𝑐12=βˆ’0.003𝑐22=0.221𝑐32=βˆ’0.948𝑐42=βˆ’0.093

πœ–1π‘“π‘œπ‘Ÿπœ™1=βˆ’4.726
πœ–2π‘“π‘œπ‘Ÿπœ™1=βˆ’0.308
π»π‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’ πΉπ‘œπ‘Ÿπ‘˜π‘Žπ‘‘π‘œπ‘šπ‘’π‘›π‘’π‘Ÿπ‘”π‘¦=βˆ’14.56616 β„Žπ‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’=βˆ’396.35976 𝑒𝑉
Time used: 1.0 s
------------------------------ Iteration 5 ------------------------------
𝜁1=5.591𝜁2=3.355𝜁3=1.011𝜁4=0.610
π‘‚π‘Ÿπ‘π‘–π‘‘π‘Žπ‘™π‘ :
πœ™1=𝑐11πœ’1+𝑐21πœ’2+𝑐31πœ’3+𝑐41πœ’4
πœ™2=𝑐12πœ’1+𝑐22πœ’2+𝑐32πœ’3+𝑐42πœ’4
𝑐11=0.155𝑐21=0.858𝑐31=βˆ’0.003𝑐41=0.002
𝑐12=βˆ’0.003𝑐22=0.220𝑐32=βˆ’0.944𝑐42=βˆ’0.097

πœ–1π‘“π‘œπ‘Ÿπœ™1=βˆ’4.731
πœ–2π‘“π‘œπ‘Ÿπœ™1=βˆ’0.309
π»π‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’ πΉπ‘œπ‘Ÿπ‘˜π‘Žπ‘‘π‘œπ‘šπ‘’π‘›π‘’π‘Ÿπ‘”π‘¦=βˆ’14.57061 β„Žπ‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’=βˆ’396.48075 𝑒𝑉
Time used: 1.0 s
------------------------------ Iteration 6 ------------------------------
𝜁1=5.591𝜁2=3.355𝜁3=1.011𝜁4=0.610
π‘‚π‘Ÿπ‘π‘–π‘‘π‘Žπ‘™π‘ :
πœ™1=𝑐11πœ’1+𝑐21πœ’2+𝑐31πœ’3+𝑐41πœ’4
πœ™2=𝑐12πœ’1+𝑐22πœ’2+𝑐32πœ’3+𝑐42πœ’4
𝑐11=0.155𝑐21=0.858𝑐31=βˆ’0.003𝑐41=0.002
𝑐12=βˆ’0.003𝑐22=0.220𝑐32=βˆ’0.942𝑐42=βˆ’0.099

πœ–1π‘“π‘œπ‘Ÿπœ™1=βˆ’4.732
πœ–2π‘“π‘œπ‘Ÿπœ™1=βˆ’0.309
π»π‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’ πΉπ‘œπ‘Ÿπ‘˜π‘Žπ‘‘π‘œπ‘šπ‘’π‘›π‘’π‘Ÿπ‘”π‘¦=βˆ’14.57187 β„Žπ‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’=βˆ’396.51521 𝑒𝑉
Time used: 1.1 s
------------------------------ Iteration 7 ------------------------------
𝜁1=5.591𝜁2=3.355𝜁3=1.011𝜁4=0.610
π‘‚π‘Ÿπ‘π‘–π‘‘π‘Žπ‘™π‘ :
πœ™1=𝑐11πœ’1+𝑐21πœ’2+𝑐31πœ’3+𝑐41πœ’4
πœ™2=𝑐12πœ’1+𝑐22πœ’2+𝑐32πœ’3+𝑐42πœ’4
𝑐11=0.155𝑐21=0.858𝑐31=βˆ’0.003𝑐41=0.002
𝑐12=βˆ’0.003𝑐22=0.220𝑐32=βˆ’0.942𝑐42=βˆ’0.099

πœ–1π‘“π‘œπ‘Ÿπœ™1=βˆ’4.733
πœ–2π‘“π‘œπ‘Ÿπœ™1=βˆ’0.309
π»π‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’ πΉπ‘œπ‘Ÿπ‘˜π‘Žπ‘‘π‘œπ‘šπ‘’π‘›π‘’π‘Ÿπ‘”π‘¦=βˆ’14.57223 β„Žπ‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’=βˆ’396.52493 𝑒𝑉
Time used: 1.2 s
------------------------------ Iteration 8 ------------------------------
𝜁1=5.591𝜁2=3.355𝜁3=1.011𝜁4=0.610
π‘‚π‘Ÿπ‘π‘–π‘‘π‘Žπ‘™π‘ :
πœ™1=𝑐11πœ’1+𝑐21πœ’2+𝑐31πœ’3+𝑐41πœ’4
πœ™2=𝑐12πœ’1+𝑐22πœ’2+𝑐32πœ’3+𝑐42πœ’4
𝑐11=0.155𝑐21=0.858𝑐31=βˆ’0.003𝑐41=0.002
𝑐12=βˆ’0.003𝑐22=0.220𝑐32=βˆ’0.942𝑐42=βˆ’0.099

πœ–1π‘“π‘œπ‘Ÿπœ™1=βˆ’4.733
πœ–2π‘“π‘œπ‘Ÿπœ™1=βˆ’0.309
π»π‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’ πΉπ‘œπ‘Ÿπ‘˜π‘Žπ‘‘π‘œπ‘šπ‘’π‘›π‘’π‘Ÿπ‘”π‘¦=βˆ’14.57233 β„Žπ‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’=βˆ’396.52766 𝑒𝑉
Time used: 1.1 s
------------------------------ Iteration 9 ------------------------------
𝜁1=5.591𝜁2=3.355𝜁3=1.011𝜁4=0.610
π‘‚π‘Ÿπ‘π‘–π‘‘π‘Žπ‘™π‘ :
πœ™1=𝑐11πœ’1+𝑐21πœ’2+𝑐31πœ’3+𝑐41πœ’4
πœ™2=𝑐12πœ’1+𝑐22πœ’2+𝑐32πœ’3+𝑐42πœ’4
𝑐11=0.155𝑐21=0.858𝑐31=βˆ’0.003𝑐41=0.002
𝑐12=βˆ’0.003𝑐22=0.220𝑐32=βˆ’0.942𝑐42=βˆ’0.099

πœ–1π‘“π‘œπ‘Ÿπœ™1=βˆ’4.733
πœ–2π‘“π‘œπ‘Ÿπœ™1=βˆ’0.309
π»π‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’ πΉπ‘œπ‘Ÿπ‘˜π‘Žπ‘‘π‘œπ‘šπ‘’π‘›π‘’π‘Ÿπ‘”π‘¦=βˆ’14.57236 β„Žπ‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’=βˆ’396.52843 𝑒𝑉
Time used: 1.2 s
------------------------------ Iteration 10 ------------------------------
𝜁1=5.591𝜁2=3.355𝜁3=1.011𝜁4=0.610
π‘‚π‘Ÿπ‘π‘–π‘‘π‘Žπ‘™π‘ :
πœ™1=𝑐11πœ’1+𝑐21πœ’2+𝑐31πœ’3+𝑐41πœ’4
πœ™2=𝑐12πœ’1+𝑐22πœ’2+𝑐32πœ’3+𝑐42πœ’4
𝑐11=0.155𝑐21=0.858𝑐31=βˆ’0.003𝑐41=0.002
𝑐12=βˆ’0.003𝑐22=0.220𝑐32=βˆ’0.942𝑐42=βˆ’0.099

πœ–1π‘“π‘œπ‘Ÿπœ™1=βˆ’4.733
πœ–2π‘“π‘œπ‘Ÿπœ™1=βˆ’0.309
π»π‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’ πΉπ‘œπ‘Ÿπ‘˜π‘Žπ‘‘π‘œπ‘šπ‘’π‘›π‘’π‘Ÿπ‘”π‘¦=βˆ’14.57237 β„Žπ‘Žπ‘Ÿπ‘‘π‘Ÿπ‘’π‘’=βˆ’396.52865 𝑒𝑉
Time used: 1.4 s