Вычисление атомных орбиталей водорода¶
Уравнение Шрёдингера для атомов, содержащих более одного электрона, не может быть решено в аналитическом виде. В связи с этим рассматривают приближённые методы, наиболее существенным из которых является метод самосогласованного поля (Self consistent field(SCF) method). Идея метода заключается в том, что каждый электрон в атоме рассматривается как движущийся в самосогласованном поле, создаваемом ядром вместе со всеми остальными электронами.
Построение самосогласованного поля может осуществляться либо методом последовательных приближений (изначально предложенным Хартри) или прямым вариационным методом.
Существенно, что вычисления методом самосогласованного поля весьма громоздки, особенно для сложных атомов. Для них применяются другие методы — метод Томаса-Ферми, метод функционала плотности а также различные приближённые методы решения уравнений Хартри-Фока — например, метод Хартри-Фока-Слейтера или Хартри-Фока-Рутана.
Метод Хартри-Фока — в квантовой механике приближённый метод решения уравнения Шрёдингера путём сведения многочастичной задачи к одночастичной в предположении, что каждая частица двигается в некотором усреднённом самосогласованном поле, создаваемом всеми остальными частицами системы.
Метод состоит из нескольких стадий. На первом этапе решается задача о движении электрона в определённом модельном потенциале, который должен как можно лучше отображать взаимодействие выбранного электрона с ядрами атомов и другими электронами. Найденные волновые функции используются для того, чтобы определить взаимодействие электрона с другими электронами и ядрами, уточняя потенциал. В дальнейшем опять решается задача нахождения волновых функций электрона для нового потенциала и нахождения из него следующего, более точного. Процедура продолжается до достижения сходимости.
Волновая функция многоэлектронной системы выбирается в виде детерминанта Слэтера. Уравнения Хартри — Фока являют собой одноэлектронные уравнения типа уравнения Шрёдингера, которым соответствуют орбитали $ {\displaystyle \varphi _{j}} \phi _{j} $, отвечающие минимальным значениям энергии молекулярной системы. В простейшем случае уравнения Хартри — Фока имеют вид:
$ \displaystyle {\hat {F}}[\{\phi _{j}\}](1)={\hat {H}}^{\text{core}}(1)+\sum _{j=1}^{N/2}[2{\hat {J}}_{j}(1)-{\hat {K}}_{j}(1)] $,
где фокиан $ \displaystyle {\hat {F}}[\{\phi _{j}\}](1) $ является оператором Гамильтона для одного электрона, находящегося в самосогласованном поле. Фокиан состоит из суммы одноэлектронного оператора $ {\hat {H}}^{\text{core}}(1) $, равного сумме оператора кинетической энергии электрона (1) и оператора потенциальной энергии его взаимодействия со всеми ядрами:
$ \displaystyle {\hat {H}}^{\text{core}}(1)=-{\frac {1}{2}}\nabla _{1}^{2}-\sum _{\alpha }{\frac {Z_{\alpha }}{r_{1\alpha }}} $,
и суммы операторов $ {\displaystyle (2{\hat {J}}_{j}(1)-{\hat {K}}_{j}(1))} $, определяющих взаимодействие рассматриваемого электрона (1) с усреднённым полем остальных электронов. Действие двух последних операторов на орбиталь $ \phi _{j} $ определяется следующими соотношениями:
$ \displaystyle {\widehat {J}}_{j}(1)f(1)=f(1)\int {\left|\phi _{j}(2)\right|}^{2}{\frac {1}{r_{12}}}\,dr_{2} $ — оператор Кулона, учитывающий взаимодействие с орбиталью $ {\displaystyle j} $-го электрона,
$ \displaystyle {\hat {K}}_{j}(x_{1})f(x_{1})=\phi _{j}(x_{1})\int {{\frac {\phi _{j}^{*}(x_{2})f(x_{2})}{r_{12}}}\mathrm {d} x_{2}} $ — обменный оператор
Метод Хартри-Фока-Рутана (ХФР) — это алгебраический подход к решению уравнений Хартри-Фока, в котором неизвестные одноэлектронные функции-орбитали ищутся в виде линейных комбинаций функций заданного вида — атомных орбиталей (приближение ЛКАО).
import numpy as np
import scipy.special as sp
import scipy.misc
from IPython.display import display,Image
import os, sys
Импортируем функции для квантовых расчетов:
import hfscf
! cat hfscf.py
Метод Хартри-Фока-Рутана¶
Базисные наборы представляют собой компромисс между вычислительными возможностями и желаемой точностью. Существует два типа базисных функций (атомных орбиталей) – орбитали слэйтеровского типа (STO – Slater Type Orbitals) и гауссовского типа (GTO – Gaussian Type Orbitals).
Базисные функции гауссова и слэтеровского типа, локализованные на разных атомах, не ортогональны друг другу. Пусть в качестве базисного набора выбран ряд таких неортогональных функций.
Основная идея:уменьшение количества интегралов, за счётобнуления перекрывания некоторых орбиталей.Если два ядра далеко друг от друга то вероятно центрированныек ядрам функции не перекрываются. Основной результат с точки зрения уравнений – это матрица обмена, $ S $ = 1.Тогда уравнение Рутан-Хола будет выглядеть как:
$ {\mathbf {F}}{\mathbf {C}}={\mathbf {S}}{\mathbf {C}}{\mathbf {\epsilon }} $
где $F$ – матрица оператора Фока, $S$ – матрица перекрытия для базисных функций.
Алгоритм расчета с помощью метода Хартри Фока¶
Входные данные: В качестве входящих данных задаются координаты атомов, заряды ядер, полное число электронов и базисные функции.
Получение матриц: Получение матриц, не зависящих от собственных векторов (матрица перекрытия $S_{pq} $, матрица одноэлектронного гамильтониана $H_{pq} $, двухэлектронные интегралы.
Приведение матрицы перекрытия к единичной: Расчет матрицы $V_{pq} $.
Диагонализация матрицы Фока: Используя матрицу Vpq, решается обобщенное уравнение на собственные значения.
Построение новой матрицы плотности: Используя полученные собственные векторы, строится новая матрица плотности.
x = Image(filename='plot.jpg')
display(x)
# Import some libraries to have the maths functions and plotting
import matplotlib.pyplot as plt
x = np.linspace(-5,5,num=1000)
r = abs(x)
zeta = 1.0
psi_STO = (zeta**3/np.pi)**(0.5)*np.exp(-zeta*r)
plt.figure(figsize=(4,3))
plt.plot(x,psi_STO)
Функция расчёта орбиталей на основе подхода Хартри-Фока-Рутана¶
def SCF (r = 1.4632, Z=[1,1], b1 = hfscf.GTO["H"], b2 = hfscf.GTO["H"], b = 3, vbs=False):
R = [0, r] # r - electronic degrees of freedom,
# R - nuclear degrees of freedom.
# Generate orbital overlap matrix S
if vbs: print("*) Generando matriz de traslape S.")
s_scf = hfscf.S(R, b1, b2, b)
# Generate Hamiltonian H for Fock operator
if vbs: print("\n*) Generando hamiltoniano H.")
h_scf = hfscf.H(R, Z, b1, b2, b)
# Diagonalize matrix S and find matrix X
if vbs: print("\n*) Diagonalizando matriz S y hallando matriz diagonal X.")
X = hfscf.diagon(s_scf)
Xa = X.getH()
# Estimate density matrix P
if vbs: print("\n*) Creando matriz de densidad P.")
p_scf = np.matrix([[0,0],[0,0]], dtype=np.float64) # Referencia (7) p. 148
# Begin iterative process
if vbs: print("\n*) Comenzando con el SCF.")
for iteracion in range(50):
# Build Fock matrix F
# F = H + G
# G = p*(J-0.5*K), p - p_psf, J - Coulomb operator, K - exchange operator
if vbs: print("\n**) Generando la matriz de Fock: calculando \ integrales de dos electrones.")
g_scf = hfscf.G(r, p_scf, b1, b2, b)
f_scf = h_scf + g_scf # Referencia (7) p. 141 eq. (3.154)
# Hartree-Fock equations to Hartree-Fock-Roothaan ones
# Build matrix F '
# F' = X_adj * F * X
# X_adj - hermitian adjoint from the complex Hilbert space
if vbs: print("**) Cambiando la base de F.")
f_tra = Xa * f_scf * X
# Diagonalize matrix F and construct matrix C'
if vbs: print("**) Diagonalizando F' y generando C'.")
c_tra = hfscf.diagon2(f_tra)
# Build matrix C
# C = X * C'
if vbs: print("**) Construyendo matriz de coeficientes C.")
c_scf = X * c_tra
# Build matrix P from matrix C
if vbs: print("**) Recalculando matriz de densidad P.")
p_temp = hfscf.P(c_scf)
# Calculating of electron and total energies
electron_E = hfscf.ener_elec(p_temp, h_scf, f_scf)
total_E = hfscf.ener_tot(r, Z, 0)
# Plots
hfscf.orbital(c_tra, r, b1, b2, b) # 1D plot
hfscf.orbital2D(c_scf, X, f_scf, r, b1, b2, b, 0.02) # 2D plot
print("\nConcluida la " + str(iteracion + 1) + ". iteracion.\n")
# Review convergence
if np.linalg.norm(p_temp - p_scf) < 1E-4: # Referencia (7) p. 148
print("\n\n-->El campo autoconsistente SI ha convergido!") # The self-consistent
# field has converged!
return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}
else:
p_scf = p_temp
print("\n\n-->El campo autoconsistente NO ha convergido!\nRevisar supuestos.")
return {"S":s_scf,"H":h_scf,"X": X,"F":f_scf,"C":c_scf,"P":p_temp}
print(SCF())
Расчет самосогласованного поля¶
# Coeff is the d_n variable in the equation above
Coeff = np.array([[1.00000,0.0000000,0.000000],
[0.678914,0.430129,0.000000],
[0.444635,0.535328,0.154329]])
# Expon is the alpha variable in the equation above
Expon = np.array([[0.270950,0.000000,0.000000],
[0.151623,0.851819,0.000000],
[0.109818,0.405771,2.227660]])
psi_CGF_STO1G = Coeff[0,0]*(2*Expon[0,0]/np.pi)**(0.75)*np.exp(-Expon[0,0]*r**2)
psi_CGF_STO2G = Coeff[1,0]*(2*Expon[1,0]/np.pi)**(0.75)*np.exp(-Expon[1,0]*r**2) \
+ Coeff[1,1]*(2*Expon[1,1]/np.pi)**(0.75)*np.exp(-Expon[1,1]*r**2) \
+ Coeff[1,2]*(2*Expon[1,2]/np.pi)**(0.75)*np.exp(-Expon[1,2]*r**2)
psi_CGF_STO3G = Coeff[2,0]*(2*Expon[2,0]/np.pi)**(0.75)*np.exp(-Expon[2,0]*r**2) \
+ Coeff[2,1]*(2*Expon[2,1]/np.pi)**(0.75)*np.exp(-Expon[2,1]*r**2) \
+ Coeff[2,2]*(2*Expon[2,2]/np.pi)**(0.75)*np.exp(-Expon[2,2]*r**2)
# Plot the three functions
plt.figure(figsize=(5,3))
plt.title("Approximations to a STO with CGF")
plt.plot(x,psi_STO,label="STO")
plt.plot(x,psi_CGF_STO1G,label="STO-1G")
plt.legend()
plt.figure(figsize=(5,3))
plt.plot(x,psi_STO,label="STO")
plt.plot(x,psi_CGF_STO2G,label="STO-2G")
plt.legend()
plt.figure(figsize=(5,3))
plt.plot(x,psi_STO,label="STO")
plt.plot(x,psi_CGF_STO3G,label="STO-3G")
plt.legend()
Интегралы¶
def S_int(A,B,Rab2):
"""
Calculates the overlap between two gaussian functions
"""
return (np.pi/(A+B))**1.5*np.exp(-A*B*Rab2/(A+B))
def T_int(A,B,Rab2):
"""
Calculates the kinetic energy integrals for un-normalised primitives
"""
return A*B/(A+B)*(3.0-2.0*A*B*Rab2/(A+B))*(np.pi/(A+B))**1.5*np.exp(-A*B*Rab2/(A+B))
def V_int(A,B,Rab2,Rcp2,Zc):
"""
Calculates the un-normalised nuclear attraction integrals
"""
V = 2.0*np.pi/(A+B)*F0((A+B)*Rcp2)*np.exp(-A*B*Rab2/(A+B))
return -V*Zc
# Mathematical functions
def F0(t):
"""
F function for 1s orbital
"""
if (t<1e-6):
return 1.0-t/3.0
else:
return 0.5*(np.pi/t)**0.5*sp.erf(t**0.5)
def erf(t):
"""
Approximation for the error function
"""
P = 0.3275911
A = [0.254829592,-0.284496736,1.421413741,-1.453152027,1.061405429]
T = 1.0/(1+P*t)
Tn=T
Poly = A[0]*Tn
for i in range(1,5):
Tn=Tn*T
Poly=Poly*A[i]*Tn
return 1.0-Poly*np.exp(-t*t)
def TwoE(A,B,C,D,Rab2,Rcd2,Rpq2):
"""
Calculate two electron integrals
A,B,C,D are the exponents alpha, beta, etc.
Rab2 equals squared distance between centre A and centre B
"""
return 2.0*(np.pi**2.5)/((A+B)*(C+D)*np.sqrt(A+B+C+D))*F0((A+B)*(C+D)*Rpq2/(A+B+C+D))*np.exp(-A*B*Rab2/(A+B)-C*D*Rcd2/(C+D))
def Intgrl(N,R,Zeta1,Zeta2,Za,Zb):
"""
Declares the variables and compiles the integrals.
"""
global S12,T11,T12,T22,V11A,V12A,V22A,V11B,V12B,V22B,V1111,V2111,V2121,V2211,V2221,V2222
S12 = 0.0
T11 = 0.0
T12 = 0.0
T22 = 0.0
V11A = 0.0
V12A = 0.0
V22A = 0.0
V11B = 0.0
V12B = 0.0
V22B = 0.0
V1111 = 0.0
V2111 = 0.0
V2121 = 0.0
V2211 = 0.0
V2221 = 0.0
V2222 = 0.0
R2 = R*R
# The coefficients for the contracted Gaussian functions are below
Coeff = np.array([[1.00000,0.0000000,0.000000],
[0.678914,0.430129,0.000000],
[0.444635,0.535328,0.154329]])
Expon = np.array([[0.270950,0.000000,0.000000],
[0.151623,0.851819,0.000000],
[0.109818,0.405771,2.227660]])
D1 = np.zeros([3])
A1 = np.zeros([3])
D2 = np.zeros([3])
A2 = np.zeros([3])
# This loop constructs the contracted Gaussian functions
for i in range(N):
A1[i] = Expon[N-1,i]*(Zeta1**2)
D1[i] = Coeff[N-1,i]*((2.0*A1[i]/np.pi)**0.75)
A2[i] = Expon[N-1,i]*(Zeta2**2)
D2[i] = Coeff[N-1,i]*((2.0*A2[i]/np.pi)**0.75)
# Calculate one electron integrals
# Centre A is first atom centre B is second atom
# Origin is on second atom
# V12A - off diagonal nuclear attraction to centre A etc.
for i in range(N):
for j in range(N):
# Rap2 - squared distance between centre A and centre P
Rap = A2[j]*R/(A1[i]+A2[j])
Rap2 = Rap**2
Rbp2 = (R-Rap)**2
S12 = S12 + S_int(A1[i],A2[j],R2)*D1[i]*D2[j]
T11 = T11 + T_int(A1[i],A1[j],0.0)*D1[i]*D1[j]
T12 = T12 + T_int(A1[i],A2[j],R2)*D1[i]*D2[j]
T22 = T22 + T_int(A2[i],A2[j],0.0)*D2[i]*D2[j]
V11A = V11A + V_int(A1[i],A1[j],0.0,0.0,Za)*D1[i]*D1[j]
V12A = V12A + V_int(A1[i],A2[j],R2,Rap2,Za)*D1[i]*D2[j]
V22A = V22A + V_int(A2[i],A2[j],0.0,R2,Za)*D2[i]*D2[j]
V11B = V11B + V_int(A1[i],A1[j],0.0,R2,Zb)*D1[i]*D1[j]
V12B = V12B + V_int(A1[i],A2[j],R2,Rbp2,Zb)*D1[i]*D2[j]
V22B = V22B + V_int(A2[i],A2[j],0.0,0.0,Zb)*D2[i]*D2[j]
# Calculate two electron integrals
for i in range(N):
for j in range(N):
for k in range(N):
for l in range(N):
Rap = A2[i]*R/(A2[i]+A1[j])
Rbp = R - Rap
Raq = A2[k]*R/(A2[k]+A1[l])
Rbq = R - Raq
Rpq = Rap - Raq
Rap2 = Rap*Rap
Rbp2 = Rbp*Rbp
Raq2 = Raq*Raq
Rbq2 = Rbq*Rbq
Rpq2 = Rpq*Rpq
V1111 = V1111 + TwoE(A1[i],A1[j],A1[k],A1[l],0.0,0.0,0.0)*D1[i]*D1[j]*D1[k]*D1[l]
V2111 = V2111 + TwoE(A2[i],A1[j],A1[k],A1[l],R2,0.0,Rap2)*D2[i]*D1[j]*D1[k]*D1[l]
V2121 = V2121 + TwoE(A2[i],A1[j],A2[k],A1[l],R2,R2,Rpq2)*D2[i]*D1[j]*D2[k]*D1[l]
V2211 = V2211 + TwoE(A2[i],A2[j],A1[k],A1[l],0.0,0.0,R2)*D2[i]*D2[j]*D1[k]*D1[l]
V2221 = V2221 + TwoE(A2[i],A2[j],A2[k],A1[l],0.0,R2,Rbq2)*D2[i]*D2[j]*D2[k]*D1[l]
V2222 = V2222 + TwoE(A2[i],A2[j],A2[k],A2[l],0.0,0.0,0.0)*D2[i]*D2[j]*D2[k]*D2[l]
return
def Colect(N,R,Zeta1,Zeta2,Za,Zb):
"""
Takes the basic integrals and assembles the relevant matrices,
that are S,H,X,XT and Two electron integrals
"""
# Form core hamiltonian
H[0,0] = T11+V11A+V11B
H[0,1] = T12+V12A+V12B
H[1,0] = H[0,1]
H[1,1] = T22+V22A+V22B
# Form overlap matrix
S[0,0] = 1.0
S[0,1] = S12
S[1,0] = S12
S[1,1] = 1.0
# This is S^-1/2
X[0,0] = 1.0/np.sqrt(2.0*(1.0+S12))
X[1,0] = X[0,0]
X[0,1] = 1.0/np.sqrt(2.0*(1.0-S12))
X[1,1] = -X[0,1]
# This is the coulomb and exchange term (aa|bb) and (ab|ba)
TT[0,0,0,0] = V1111
TT[1,0,0,0] = V2111
TT[0,1,0,0] = V2111
TT[0,0,1,0] = V2111
TT[0,0,0,1] = V2111
TT[1,0,1,0] = V2121
TT[0,1,1,0] = V2121
TT[1,0,0,1] = V2121
TT[0,1,0,1] = V2121
TT[1,1,0,0] = V2211
TT[0,0,1,1] = V2211
TT[1,1,1,0] = V2221
TT[1,1,0,1] = V2221
TT[1,0,1,1] = V2221
TT[0,1,1,1] = V2221
TT[1,1,1,1] = V2222
# Self consistent field calculation (SCF):
# 1. Guess an initial density matrix $\boldsymbol{P}$ (for this example will use the initial guess $ P=0 $)
# 2. From $\boldsymbol{P}$ calculate the Fock matrix
# 3. Calculate $\boldsymbol{F'}=\boldsymbol{S}^{-1/2}\boldsymbol{F}\boldsymbol{S}^{1/2}$
# 4. Solve the eigenvalue problem using the secular equation $|\boldsymbol{F'}-\boldsymbol{E}\boldsymbol{I}|=0$ giving the eigenvalues $\boldsymbol{E}$ and eigenvectors $\boldsymbol{C'}$ which can be solved by diagonalising $\boldsymbol{F'}$
# 5. Calculate the molecular orbitals coefficients by reverting $\boldsymbol{C'}=\boldsymbol{S}^{-1/2}\boldsymbol{C}$
# 6. Calculate the new density matrix $\boldsymbol{P}$ from the matrix $\boldsymbol{C}$
# 7. Check to see if the energy has converged if not then head back to step 6 and repeat.
def SCF(N,R,Zeta1,Zeta2,Za,Zb,G):
"""
Performs the SCF iterations
"""
Crit = 1e-11 # Convergence critera
Maxit = 250 # Maximum number of iterations
Iter=0
######## STEP 1. Guess an initial density matrix ########
# Use core hamiltonian for initial guess of F, I.E. (P=0)
P = np.zeros([2,2])
Energy = 0.0
while (Iter<Maxit):
Iter += 1
print(Iter)
######## STEP 2. calculate the Fock matrix ########
# Form two electron part of Fock matrix from P
G = np.zeros([2,2]) # This is the two electron contribution in the equations above
for i in range(2):
for j in range(2):
for k in range(2):
for l in range(2):
G[i,j]=G[i,j]+P[k,l]*(TT[i,j,k,l]-0.5*TT[i,j,k,l])
# Add core hamiltonian H^CORE to get fock matrix
F = H+G
# Calculate the electronic energy
Energy = np.sum(0.5*P*(H+F))
print('Electronic energy = ',Energy)
######## STEP 3. Calculate F' (remember S^-1/2 is X and S^1/2 is X.T) ########
G = np.matmul(F,X)
Fprime = np.matmul(X.T,G)
######## STEP 4. Solve the eigenvalue problem ########
# Diagonalise transformed Fock matrix
Diag(Fprime,Cprime,E)
######## STEP 5. Calculate the molecular orbitals coefficients ########
# Transform eigen vectors to get matrix C
C = np.matmul(X,Cprime)
######## STEP 6. Calculate the new density matrix from the old P ########
Oldp = np.array(P)
P = np.zeros([2,2])
# Form new density matrix
for i in range(2):
for j in range(2):
#Save present density matrix before creating a new one
for k in range(1):
P[i,j] += 2.0*C[i,k]*C[j,k]
######## STEP 7. Check to see if the energy has converged ########
Delta = 0.0
# Calculate delta the difference between the old density matrix Old P and the new P
Delta = (P-Oldp)
Delta = np.sqrt(np.sum(Delta**2)/4.0)
print("Delta",Delta)
#Check for convergence
if (Delta<Crit):
# Add nuclear repulsion to get the total energy
Energytot = Energy+Za*Zb/R
print("Calculation converged with electronic energy:",Energy)
print("Calculation converged with total energy:",Energytot)
print("Density matrix", P)
print("Mulliken populations",np.matmul(P,S))
print("Coeffients",C)
break
def FormG():
"""
Calculate the G matrix from the density matrix and two electron integals
"""
for i in range(2):
for j in range(2):
G[i,j]=0.0
for k in range(2):
for l in range(2):
G[i,j]=G[i,j]+P[k,l]*(TT[i,j,k,l]-0.5*TT[i,j,k,l])
def Mult(A,B,C_,IM,M):
"""
Multiples two square matrices A and B to get C
"""
for i in range(M):
for j in range(M):
for k in range(M):
C_[i,j] = A[i,j]*B[i,j]
def Diag(Fprime,Cprime,E):
"""
Diagonalises F to give eigenvectors in C and eigen values in E, theta is the angle describing the solution
"""
#
import math
# Angle for heteronuclear diatonic
Theta = 0.5*math.atan(2.0*Fprime[0,1]/(Fprime[0,0]-Fprime[1,1]))
#print('Theta', Theta)
Cprime[0,0] = np.cos(Theta)
Cprime[1,0] = np.sin(Theta)
Cprime[0,1] = np.sin(Theta)
Cprime[1,1] = -np.cos(Theta)
E[0,0] = Fprime[0,0]*np.cos(Theta)**2+Fprime[1,1]*np.sin(Theta)**2+Fprime[0,1]*np.sin(2.0*Theta)
E[1,1] = Fprime[1,1]*np.cos(Theta)**2+Fprime[0,0]*np.sin(Theta)**2-Fprime[0,1]*np.sin(2.0*Theta)
if (E[1,1] <= E[0,0]):
Temp = E[1,1]
E[1,1] = E[0,0]
E[0,0] = Temp
Temp = Cprime[0,1]
Cprime[0,1] = Cprime[0,0]
Cprime[0,0] = Temp
Temp = Cprime[1,1]
Cprime[1,1]=Cprime[1,0]
Cprime[1,0]=Temp
return
def HFCALC(N,R,Zeta1,Zeta2,Za,Zb,G):
"""
Calculates the integrals constructs the matrices and then runs the SCF calculation
"""
# Calculate one and two electron integrals
Intgrl(N,R,Zeta1,Zeta2,Za,Zb)
# Put all integals into array
Colect(N,R,Zeta1,Zeta2,Za,Zb)
# Perform the SCF calculation
SCF(N,R,Zeta1,Zeta2,Za,Zb,G)
return
"""
Let's set up the variables and perform the calculations
"""
global H,S,X,XT,TT,G,C,P,Oldp,F,Fprime,Cprime,E,Zb
H = np.zeros([2,2])
S = np.zeros([2,2])
X = np.zeros([2,2])
XT = np.zeros([2,2])
TT = np.zeros([2,2,2,2])
G = np.zeros([2,2])
C = np.zeros([2,2])
P = np.zeros([2,2])
Oldp = np.zeros([2,2])
F = np.zeros([2,2])
Fprime = np.zeros([2,2])
Cprime = np.zeros([2,2])
E = np.zeros([2,2])
Energy = 0.0
Delta = 0.0
N = 3
R = 1.4632
Zeta1 = 2.0925
Zeta2 = 1.24
Za = 2.0
Zb = 1.0
HFCALC(N,R,Zeta1,Zeta2,Za,Zb,G)
Молекулярные орбитали¶
C = np.matmul(X,Cprime)
P = np.array([[ 1.28614168,0.54017322],
[ 0.54017322 ,0.22687011]])
x = np.linspace(-8,8,num=1000)
r1 = abs(x+R/2)
r2 = abs(x-R/2)
psi_CGF_STO3G_He = Coeff[2,0]*(2*Expon[2,0]/np.pi)**(0.75)*np.exp(-Expon[2,0]*r1**2) \
+ Coeff[2,1]*(2*Expon[2,1]/np.pi)**(0.75)*np.exp(-Expon[2,1]*r1**2) \
+ Coeff[2,2]*(2*Expon[2,2]/np.pi)**(0.75)*np.exp(-Expon[2,2]*r1**2)
psi_CGF_STO3G_H = Coeff[2,0]*(2*Expon[2,0]/np.pi)**(0.75)*np.exp(-Expon[2,0]*r2**2) \
+ Coeff[2,1]*(2*Expon[2,1]/np.pi)**(0.75)*np.exp(-Expon[2,1]*r2**2) \
+ Coeff[2,2]*(2*Expon[2,2]/np.pi)**(0.75)*np.exp(-Expon[2,2]*r2**2)
density = np.zeros(x.shape)
density = density + P[0,0]*psi_CGF_STO3G_He*psi_CGF_STO3G_He
density = density + P[1,1]*psi_CGF_STO3G_H*psi_CGF_STO3G_H
density = density + 2*P[0,1]*psi_CGF_STO3G_He*psi_CGF_STO3G_H
plt.figure(figsize=(7,4))
plt.title("Molecular orbitals")
plt.plot(x,C[0,0]*psi_CGF_STO3G_He+C[1,0]*psi_CGF_STO3G_H,label="Bonding")
plt.plot(x,C[0,1]*psi_CGF_STO3G_He+C[1,1]*psi_CGF_STO3G_H,label="Anti-bonding")
plt.legend(loc=4)
plt.figure(figsize=(7,4))
plt.title("Electron density $\Psi^{2}$")
plt.plot(x,density)
