双中心用椭圆坐标
用线性组合原子轨道LCAO方法把电子波函数设为
假设这个轨道是两个轨道的线性组合,φA,φB用氢原子的1s轨道代入
轨道能量E
将E的表达式展开
因为cA,cB是实系数,因此可以进一步化简
设能量
重叠积分
原式可以化简为
由于对称性
得到
把cA看作变量,微分得到
再把cB看作变量,对等式两边微分得到
得到行列式
展开
由于对称性
解得
将HAA展开
电子对原子核A的动能和势能+电子对原子核B的势能+两个原子核的势能。
因为
E1s氢原子基态能级=-0.5,原式可以化简为
用同样的办法展开HAB
设
得到
S是重叠积分=0.586
J是库仑积分=0.473
K是交换积分=0.406
因此
因此得到轨道能量
实验值是-0.60264*Hartrees,因此计算值是实验值的91.9%。
Peter Atkins Ronald S. Friedman,Molecular Quantum Mechanics P264
Michael P. Mueller,Fundamentals of Quantum Chemistry Molecular Spectroscopy and Modern Electronic Structure Computations P225
Gupta, V. P,Principles and Applications of Quantum Chemistry P73
H Eyring, J Walter, G E Kimball,Quantum Chemistry P197
*Bates, D.R., Ledsham, K, Stewart, A, L. Wave Functions of the Hydrogen Molecular Ion. Philos. Trans. Roy. Soc. London, Ser. A. 246, 215 (1953)
import sympy
import math
from sympy import symbols, cancel
a = sympy.Symbol('a')
e = sympy.Symbol('e')
m = sympy.Symbol('m')
h = sympy.Symbol('h')
l = sympy.Symbol('l')
lp = sympy.Symbol('lp')
r = sympy.Symbol('r')
EE = sympy.Symbol('EE')
R = sympy.Symbol('R')
r1 = sympy.Symbol('r1')
r2 = sympy.Symbol('r2')
r3 = sympy.Symbol('r3')
μ = sympy.Symbol('μ')
v = sympy.Symbol('v')
α = sympy.Symbol('α')
β = sympy.Symbol('β')
x = sympy.Symbol('x')
y = sympy.Symbol('y')
z = sympy.Symbol('z')
θ1= sympy.Symbol('θ1')
θ2= sympy.Symbol('θ2')
Φ1= sympy.Symbol('Φ1')
Φ2= sympy.Symbol('Φ2')
θ= sympy.Symbol('θ')
Ψ= sympy.Symbol('Ψ')
Φ= sympy.Symbol('Φ')
pi=sympy.Symbol('pi')
E=sympy.Symbol('E')
I=sympy.Symbol('I')
sin=sympy.Symbol('sin')
cos=sympy.Symbol('cos')
diff=sympy.Symbol('diff')
integrate=sympy.Symbol('integrate')
pi=sympy.pi
E=sympy.E
sin=sympy.sin
cos=sympy.cos
diff=sympy.diff
integrate=sympy.integrate
z1=1
z2=1
#原子核距离
R=2.0
fx1=(z1)**(1.5)*sympy.exp(-z1*R*(μ+v)/2 )*pi**(-0.5)
fx2=(z2)**(1.5)*sympy.exp(-z2*R*(μ-v)/2 )*(pi)**(-0.5)
#重叠积分
sab=( integrate( ( integrate( integrate( fx1*fx2 *(1/8)*R*R*R *(μ*μ-v*v) , (μ , 1 , float('inf') ) ) , (v, -1 , 1 ) ) ) , (Φ,0,2*pi) ) )
print ( "sab", sab )
#sab 0.586452894025322
#库仑积分
j=( integrate( ( integrate( integrate( ( pi**(-1)*sympy.exp(- R*(μ+v) ) / (R*(μ-v)/2 ) )*(R*R*R/(8) )*(μ*μ-v*v), (μ , 1 , float('inf') ) ) , (v, -1 , 1 ) ) ) , (Φ,0,2*pi) ) )
print ( "j", j )
#交换积分
k=( integrate( ( integrate( integrate( ( pi**(-1)*sympy.exp(- R*(μ) ) / (R*(μ+v)/2 ) )*(R*R*R/(8) )*(μ*μ-v*v), (μ , 1 , float('inf') ) ) , (v, -1 , 1 ) ) ) , (Φ,0,2*pi) ) )
print ( "k", k )
Haa=-0.5-j+1/R
Hab=-0.5*sab-k+0.5*sab
print ( "Haa", -0.5-j+1/R )
print ( "Hab", -0.5*sab-k+0.5*sab )
print ( "E+", (Haa+Hab)/(1+sab) , (Haa+Hab)/(1+sab)*27.2 )
print ( "E-", (Haa-Hab)/(1-sab) , (Haa-Hab)/(1-sab)*27.2 )
#E+ -0.553771495318483 -15.0625846726627
#E- -0.160853965596688 -4.37522786422990
#Haa -0.472526541666899
#Hab -0.406005849709838
#######