用类氢轨道计算交换积分和Li原子2S谱项能级

Li原子核外有3个电子，让他们的轨道分别是(1，0，0）和(2，0，0)，用类氢轨道代入，3个电子的轨道分别是

其中z是核电荷数等于3，a0是第一层轨道半径等于1，a1是第二层轨道半径等于4.

将能量的狄拉克符号展开

三个电子的动能+势能+库仑排斥能-交换能

其中E1=E2，J13=J23，K13=K23

先计算动能和势能

计算库仑积分，因为3个电子都是s轨道所以可以用公式

因此

 

同样的方法J13=J23=0.18

 

计算交换积分

设Ψ3是α轨道，则电子3的β轨道是

因为都是s轨道所以同样可以用

因此

 

K23=K13

 

因此能级

E=-4.5-4.5-0.4921+0.18+0.18+1.875-0.034-0.034=-7.32Hartrees

用Hartree-Fork方法计算得到能量是-7.4327*，计算值是HF方法的98.5%

*徐光宪-《量子化学》p873

python代码
 

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

#Li原子基态能级

def hin( fx1 ,z ):
    fx = fx1
    z=z
    # 拉普拉斯算符
    f1 = (1 / (r * r)) * diff((r * r * diff(fx, r)), r)

    f2 = (1 / (r * r * sin(θ))) * diff((sin(θ) * diff(fx, θ)), θ)

    f3 = (1 / (r * r * sin(θ) * sin(θ))) * diff(fx, Φ, Φ)

    f8 = (-1 / 2) * (f1 + f2 + f3) * fx

    # print   (   f1 )
    # print   (   f2 )
    # print   (   f3 )

    # print    ( f8 )

    # 球坐标积分  动能
    f9 = (integrate((integrate(integrate(f8 * r * r * sin(θ), (r, 0, float('inf'))), (θ, 0, pi))), (Φ, 0, 2 * pi)))

    # print(f9)

    f10 = fx * (-z / r) * fx
    # 势能
    f11 = (integrate((integrate(integrate(f10 * r * r * sin(θ), (r, 0, float('inf'))), (θ, 0, pi))), (Φ, 0, 2 * pi)))

    # print(f11)
    print("H", f9 + f11)

    return f9 + f11

def jin (fr1 ,fr2 ,z):

    z=z
    f21 = fr1 * fr2 * (1 / r1) * fr1 * fr2 * r1 * r1 * r2 * r2

    f22 = fr1 * fr2 * (1 / r2) * fr1 * fr2 * r1 * r1 * r2 * r2

    f23 = (integrate(f21, (r2, 0, r1)))

    f24 = (integrate(f22, (r2, r1, float('inf'))))

    f25 = (16 * pi ** (2) * integrate(f24 + f23, (r1, 0, float('inf'))))

    # print("f23",f23)
    # print("f24",f24)

    print("J", f25)

    return f25

def kin (fr1,fr2,fr3,fr4 ,z):

    z=z
    # 交换积分
    #fr1 = (z) ** (1.5) * sympy.exp(-z * r1) * pi ** (-0.5)
    #fr2 = (z) ** (1.5) * sympy.exp(-z * r2) * pi ** (-0.5)

    #fr3 = (z) ** (1.5) * sympy.exp(-z * r2) * pi ** (-0.5)
    #fr4 = (z) ** (1.5) * sympy.exp(-z * r1) * pi ** (-0.5)

    f21 = fr1 * fr2 * (1 / r1) * fr3 * fr4 * r1 * r1 * r2 * r2

    f22 = fr1 * fr2 * (1 / r2) * fr3 * fr4 * r1 * r1 * r2 * r2

    f23 = (integrate(f21, (r2, 0, r1)))

    f24 = (integrate(f22, (r2, r1, float('inf'))))

    f36 = (16 * pi ** (2) * integrate(f24 + f23, (r1, 0, float('inf'))))

    # print("f23",f23)
    # print("f24",f24)

    print("K", f36)

    return f36


z=3

a0=1
a1=4*a0

fx1=(z/a0)**(1.5)*2*sympy.exp(-z*r/a0  )*(4*pi)**(-0.5)
fx2=(z/( 2*a1))**(1.5)*(2-z*r/a1)*sympy.exp(-z*r/(2*a1)  )*(4*pi)**(-0.5)

fj1 = (z/a0)**(1.5)*2*sympy.exp(-z*r1/a0  )*(4*pi)**(-0.5)
fj2 = (z/a0)**(1.5)*2*sympy.exp(-z*r2/a0  )*(4*pi)**(-0.5)

fj3=(z/( 2*a1))**(1.5)*(2-z*r1/a1)*sympy.exp(-z*r1/(2*a1)  )*(4*pi)**(-0.5)
fj4=(z/( 2*a1))**(1.5)*(2-z*r2/a1)*sympy.exp(-z*r2/(2*a1)  )*(4*pi)**(-0.5)



fh=hin(fx1,z)*2+hin(fx2,z)

fj=jin(fj1,fj2,z)+jin(fj2,fj3,z)*2

fk=kin(fj1,fj3,fj2,fj4,z)*2

print(   fh+fj-fk , (fh+fj-fk)/7.4327   )