今回は、女性の方、就学していない方、さらにこれから就学する小学生の方を対象に、量子化学入門を企画しました。プログラムは Python で作成しています。
量子化学は Linus Pauling の混成軌道 (Hybrid Orbital) に始まります。1個の炭素原子 C が、4個の水素原子 H と結合しているメタン分子 CH4 が最もシンプルな例です。H2 水素分子が最もシンプルな分子ですが、水素分子は量子化学が分岐する以前、初期量子力学の段階から研究されていて、水素分子はどちらかというと量子力学の方に研究の大部分があります。まず量子化学計算の最初期は Hückel MO 法です。しかし今回は Hückel MO 法から少し進化した電子スピンおよび 2 電子間反発も考慮した ppp 近似の SCF 量子化学計算です。1930～1960 年代に Roothaan や Pople らによってなされた最新版です。しかし現在では、核スピン、磁場、および電場の影響を考慮した量子計算の研究が進行中であり、研究としては少し古い感があります。しかし少し古いと言いながら、この SCF 量子化学計算から突然少し難解になりハードルが上がります。入門者の方に、これを乗り越えてもらおうという意味合いもあります。まずプログラムを見て、どういう計算をやっているのかを把握していただくだけでも、読んでいただく価値はあると思います。
最初に Python プログラムを掲載しています。1980年代に出版された「量子化学入門（上）」の p292-p299, p333-p334 から計算式とパラメーターを引用しています。
その後で同じ参考書から引用した理論を少しだけ書いています。理論をより深く追求したい方は参考書をご覧ください。（丁寧に説明が付記された本でしたが少し古いので書籍自体の入手がむつかしい場合は、多くの図書館に所蔵がありますので図書館のコピーサービスが便利です。）
Roothaan の SCF 計算プログラム
閉殻系（closed shell）のアリルカチオン、そして開殻系（opened shell）のアリルラジカルの RHF、この 2 種類の SCF 計算をするプログラムです。共通パラメーターはプログラム上部に列挙しています。F11, F12, F13, F22 ... だけが個別のパラメーターです。原子積分の値は半経験的手法で得た既知の値を使用します。一方、原子積分の値を非経験的手法で得る方法は後述の「参考：非経験的な解法」をご覧ください。以下、Python で作成したプログラムです。Python IDLE 3.7 (64bit) 版で動作確認をしています。また有効数字の取り方が手計算と少し異なるので値がズレるかもしれません。
方程式部分は、整数型で高次方程式を解く... を参考にさせていただきました。楽しみな試みだと思いました。
from multiprocessing import Manager, Process
from decimal import Decimal, ROUND_HALF_UP
## - - - - - - - - - - - - - - - - - - - -
kaisu = 10 # repeat MAX
ji_suu = 3
r_f = 10 # resolution
y_zoom = 0
## - - - - - - - - - - - - - - - - p295 -
W2p = -11.42
I11 = W2p - 13.15
I22 = W2p - 15.04
I33 = I11
I12 = - 2.85
I13 = - 0.45
I23 = I12
r11 = 10.84 # (11|11)
r22 = r11 # (22|22)
r33 = r11 # (33|33)
r12 = 7.52 # (11|22)
r23 = r12 # (22|33)
r13 = 5.63 # (11|33)
r21 = r12
r31 = r13
r32 = r12
## - - - - - - - - - - - - - - - - p296 -
Cpre = [ 0.5000, 0.7071, 0.5000,\
0.7071, 0, -0.7071,\
0.5000, -0.7071, 0.5000]
## - - - - - - - - - - - - - - - - - - - -
def f0(List0, R_f, Ji_su, Y_zoom, X_init, X_end, Ksu, Keta, Fugo):
Ypre =0
for X in range( int(X_init), int(X_end)+1 ):
X *= Fugo
Y=0
for j in range(0, Ji_su+1):
Y += Ksu[j]*(X**j)
if X == -int(X_init):
Ypre =Y
List0.pop(0)
if Ypre * Y > 0:
Ypre=Y
else:
if Y == 0:
Y=Y
else:
X_tmp = X-1
Y_tmp = 0
for k in range(0, Ji_su+1):
Y_tmp += Ksu[k]*(X_tmp**k)
List0.append(( X_tmp/(R_f*Keta), X_tmp, Y_tmp/((R_f*Keta)**(Ji_su-Y_zoom)), Y_tmp ))
Ypre=Y
if Fugo < 0:
List0.reverse()
if __name__=="__main__":
CpreDecimal=[]
for i in range(len(Cpre)):
CpreDecimal.append( Decimal(str(Cpre[i])).quantize(Decimal('0.0001'), rounding=ROUND_HALF_UP) )
print(f"Original φi :")
for i in range( int(len(CpreDecimal)/3) ):
print('{: }'.format(CpreDecimal[3*i+0]), f"χ1 ", '{:+}'.format(CpreDecimal[3*i+1]), f"χ2 ", '{:+}'.format(CpreDecimal[3*i+2]), f"χ3")
with Manager() as manager:
for x in range(2):
if x == 0:
Calc_type = 'Aryl Cation'
if x == 1:
Calc_type = 'Aryl Radical RHF'
for xx in range(kaisu):
if xx == 0:
print(f" ")
print(f"－－－－－－－－－－－－－－－－－－－")
print(Calc_type, f"SCF Calculation Start.")
print(f"－－－－－－－－－－－－－－－－－－－")
## Ptu Frs
if Calc_type == 'Aryl Cation':
keta = 10000
p11 = 2*Cpre[0]*Cpre[0]
p22 = 2*Cpre[1]*Cpre[1]
p33 = p11
p12 = 2*Cpre[0]*Cpre[1]
p23 = p12
p13 = p11
F11 = I11 + p11*r11*0.5 + p22*r12 + p33*r13
F22 = I22 + p22*r22*0.5 + p11*r12 + p33*r23
F33 = I33 + p33*r33*0.5 + p22*r23 + p11*r13
F12 = I12 - p12*r12*0.5
F13 = I13 - p13*r13*0.5
F23 = I23 - p23*r23*0.5
F21 = F12
F31 = F13
F32 = F23
elif Calc_type == 'Aryl Radical RHF':
keta = 10000
#F11 = I11 + ( (Cpre[0]**2) - (1/2)*(Cpre[3]**2) + 2*( (Cpre[0]*Cpre[3])**2 ) + (Cpre[3]**4) ) *r11 +\
# ( 2*(Cpre[1]**2) + (Cpre[4]**2) + Cpre[3]*Cpre[4]*( 2*Cpre[0]*Cpre[1] + Cpre[3]*Cpre[4] ) ) *r12 +\
# ( 2*(Cpre[2]**2) + (Cpre[5]**2) + Cpre[3]*Cpre[5]*( 2*Cpre[0]*Cpre[2] + Cpre[3]*Cpre[5] ) ) *r13
#F11 = I11 + 0.5*r11 + r12 + r13
F11 = I11 + ( Cpre[0]**2 - 0.5*(Cpre[3]**2) + 2*(Cpre[0]*Cpre[3])**2 + Cpre[3]**4 ) *r11 +\
2*Cpre[1]**2*r12 +\
( 2*Cpre[2]**2 + Cpre[5]**2 + Cpre[3]*Cpre[5]*( 2*Cpre[0]*Cpre[2] + Cpre[3]*Cpre[5] ) ) *r13
#F22 = I22 + ( (Cpre[1]**2) - (1/2)*(Cpre[4]**2) + 2*( (Cpre[1]*Cpre[4])**2 ) + (Cpre[4]**4) ) *r22 +\
# ( 2*(Cpre[0]**2) + (Cpre[3]**2) + Cpre[4]*Cpre[3]*( 2*Cpre[1]*Cpre[0] + Cpre[4]*Cpre[3] ) ) *r21 +\
# ( 2*(Cpre[2]**2) + (Cpre[5]**2) + Cpre[4]*Cpre[5]*( 2*Cpre[1]*Cpre[2] + Cpre[4]*Cpre[5] ) ) *r23
#F22 = I22 + 0.5*r22 + 2*r12
F22 = I22 + Cpre[1]**2*r22 + ( 2*Cpre[0]**2 + Cpre[3]**2 + 2*Cpre[2]**2 + Cpre[5]**2 ) *r12
F33 = F11
#F12 = I12 + ( Cpre[0]*Cpre[1] + (1/2)*Cpre[3]*Cpre[4] ) * ( Cpre[3]**2*r11 + Cpre[4]**2*r22 ) +\
# ( -Cpre[0]*Cpre[1] - (3/2)*Cpre[3]*Cpre[4] + (1/2)*Cpre[3]*Cpre[4]*( 2*(Cpre[0]**2) + (Cpre[3]**2) + 2*(Cpre[1]**2) + (Cpre[4]**2) ) ) *r12 +\
# Cpre[5]*Cpre[3] * ( Cpre[2]*Cpre[1] + (1/2)*Cpre[5]*Cpre[4] ) *r13 +\
# Cpre[5]*Cpre[4] * ( Cpre[2]*Cpre[0] + (1/2)*Cpre[5]*Cpre[3] ) *r23
#F12 = I12 + (32**(-0.5))*( r11 - 2*r12 -r13 )
F12 = I12 + Cpre[0]*Cpre[1] * (Cpre[3]**2) *r11 +\
(-Cpre[0]*Cpre[1]) *r12 +\
Cpre[5]*Cpre[3] * (Cpre[2]*Cpre[1]) *r13
#F13 = I13 + ( Cpre[0]*Cpre[2] + (1/2)*Cpre[3]*Cpre[5] ) * ( Cpre[3]**2*r11 + Cpre[5]**2*r22 ) +\
# ( -Cpre[0]*Cpre[2] - (3/2)*Cpre[3]*Cpre[5] + (1/2)*Cpre[3]*Cpre[5]*( 2*(Cpre[0]**2) + (Cpre[3]**2) + 2*(Cpre[2]**2) + (Cpre[5]**2) ) ) *r12 +\
# Cpre[4]*Cpre[3] * ( Cpre[1]*Cpre[2] + (1/2)*Cpre[4]*Cpre[5] ) *r13 +\
# Cpre[4]*Cpre[5] * ( Cpre[1]*Cpre[0] + (1/2)*Cpre[4]*Cpre[3] ) *r23
#F13 = I13
F13 = I13 + ( Cpre[0]*Cpre[2] + (1/2)*Cpre[3]*Cpre[5] ) * ( Cpre[3]**2*r11 + Cpre[5]**2*r22 ) +\
( -Cpre[0]*Cpre[2] - (3/2)*Cpre[3]*Cpre[5] + (1/2)*Cpre[3]*Cpre[5]*( 2*(Cpre[0]**2) + (Cpre[3]**2) + 2*(Cpre[2]**2) + (Cpre[5]**2) ) ) *r12
F23 = F12
F21 = F12
F31 = F13
F32 = F23
## Eq Calculation, 4 core
enen = [ -F11, -F12, -F13,\
-F21, -F22, -F23,\
-F31, -F32, -F33]
enenDecimal=[]
for i in range(len(enen)):
enenDecimal.append( Decimal(str(enen[i])).quantize(Decimal('0.0001'), rounding=ROUND_HALF_UP) )
enenReFloat=[]
for i in range(len(enenDecimal)):
enenReFloat.append( float(enenDecimal[i]) )
enen.clear()
enen = enenReFloat.copy()
Enen = []
for i in range(len(enen)):
Enen.append( int(enen[i]*keta) )
ksu_0 = []
tmp_0 = Enen[0]*Enen[4]*Enen[8] + Enen[3]*Enen[7]*Enen[2] + Enen[6]*Enen[1]*Enen[5] - Enen[6]*Enen[4]*Enen[2] - Enen[0]*Enen[7]*Enen[5] - Enen[3]*Enen[1]*Enen[8]
ksu_0.append(tmp_0)
tmp_1 = Enen[0]*Enen[4] + Enen[0]*Enen[8] + Enen[4]*Enen[8] - Enen[6]*Enen[2] - Enen[7]*Enen[5] - Enen[1]*Enen[3]
ksu_0.append(tmp_1)
tmp_2 = Enen[0] + Enen[4] + Enen[8]
ksu_0.append(tmp_2)
tmp_3 = 1
ksu_0.append(tmp_3)
ji_su =int(ji_suu)
ksu=[]
for i in range(len(ksu_0)):
ksu.append( int(ksu_0[i]*r_f**(ji_su-i)) )
listMinus = manager.list(range(1))
fugo = -1
x_from = 16
x_to = 25
Xfrom = int(keta*r_f*x_from)
Xto = int(keta*r_f*x_to)
Xrange = Xto - Xfrom
p1 = Process(target=f0, args=(listMinus,r_f,ji_su,y_zoom,Xfrom,Xto,ksu,keta,fugo))
p1.start()
listMinus2 = manager.list(range(1))
fugo = -1
x_from = 8
x_to = 16
Xfrom = int(keta*r_f*x_from)
Xto = int(keta*r_f*x_to)
Xrange = Xto - Xfrom
p2 = Process(target=f0, args=(listMinus2,r_f,ji_su,y_zoom,Xfrom,Xto,ksu,keta,fugo))
p2.start()
listMinus3 = manager.list(range(1))
fugo = -1
x_from = 0
x_to = 8
Xfrom = int(keta*r_f*x_from)
Xto = int(keta*r_f*x_to)
Xrange = Xto - Xfrom
p3 = Process(target=f0, args=(listMinus3,r_f,ji_su,y_zoom,Xfrom,Xto,ksu,keta,fugo))
p3.start()
listPlus = manager.list(range(1))
fugo = 1
x_from = 0
x_to = 8
Xfrom = int(keta*r_f*x_from)
Xto = int(keta*r_f*x_to)
Xrange = Xto - Xfrom
p4 = Process(target=f0, args=(listPlus,r_f,ji_su,y_zoom,Xfrom,Xto,ksu,keta,fugo))
p4.start()
p1.join()
p2.join()
p3.join()
p4.join()
## after calc
listAll=[]
if len(listMinus) > 0:
for i in range(len(listMinus)):
listAll.append(listMinus[i])
if len(listMinus2) > 0:
for i in range(len(listMinus2)):
listAll.append(listMinus2[i])
if len(listMinus3) > 0:
for i in range(len(listMinus3)):
listAll.append(listMinus3[i])
if len(listPlus) > 0:
for i in range(len(listPlus)):
listAll.append(listPlus[i])
listXall=[]
for i in range(len(listAll)):
listXall.append( listAll[i][0] )
listXallW2p=[]
for i in range(len(listXall)):
listXallW2p.append( listXall[i] - W2p )
listXallW2pDecimal=[]
for i in range(len(listXallW2p)):
listXallW2pDecimal.append( Decimal(str(listXallW2p[i])).quantize(Decimal('0.00001'), rounding=ROUND_HALF_UP) )
print(f"New Energy :")
for i in range(len(listXallW2pDecimal)):
print( f"ε-W2p =", '{: }'.format(listXallW2pDecimal[i]) )
## New Cir
e1 = listXall[0]*keta
e2 = listXall[1]*keta
e3 = listXall[2]*keta
C11N = ( Enen[1]**2 / ( 2*(Enen[1]**2) + (Enen[0] + e1 + Enen[2])**2 ) )**0.5
C12N = ( (Enen[0] + e1 + Enen[2])**2 / ( 2*(Enen[1]**2) + (Enen[0] + e1 + Enen[2])**2 ) )**0.5
C13N = C11N
C21N = ( Enen[2]**2 / ( (Enen[2]**2) + (Enen[0] + e2)**2 ) )**0.5
C22N = 0
C23N = -C21N
C31N = ( Enen[1]**2 / ( 2*(Enen[1]**2) + (Enen[0] + e3 + Enen[2])**2 ) )**0.5
C32N = -( (Enen[0] + e3 + Enen[2])**2 / ( 2*(Enen[1]**2) + (Enen[0] + e3 + Enen[2])**2 ) )**0.5
C33N = C31N
listCnew=[C11N,C12N,C13N,C21N,C22N,C23N,C31N,C32N,C33N]
listCnewDecimal=[]
for i in range(len(listCnew)):
listCnewDecimal.append( Decimal(str(listCnew[i])).quantize(Decimal('0.0001'), rounding=ROUND_HALF_UP) )
listCnewReFloat=[]
for i in range(len(listCnewDecimal)):
listCnewReFloat.append( float(listCnewDecimal[i]) )
print(f"New φi :")
for i in range( int(len(listCnewDecimal)/3) ):
print('{: }'.format(listCnewDecimal[3*i+0]), f"χ1 ", '{:+}'.format(listCnewDecimal[3*i+1]), f"χ2 ", '{:+}'.format(listCnewDecimal[3*i+2]), f"χ3")
## Decision
f_all=0
for i in range(len(listCnewReFloat)):
if listCnewReFloat[i]-Cpre[i] > 0.0001:
f_all=1
if f_all == 0:
print(f"Result", xx+1, f": Under 0.0001 all")
print(f"－－－－－－－－－－－－－－－－－－－")
print(Calc_type, f"SCF Calculation Complete.")
print(f" Ｏ Ｏ Ｏ Ｏ Ｏ Ｏ Ｏ Ｏ Ｏ Ｏ Ｏ Ｏ")
break
else :
print(f"Result", xx+1, f": Over 0.0001")
print(f"－－－－－－－－－－－－－－－－－－－")
if xx == kaisu-1:
print(f"Repeat Count Over !")
print(Calc_type, f"SCF Calculation Incomplete.")
print(f" Ｘ Ｘ Ｘ Ｘ Ｘ Ｘ Ｘ Ｘ Ｘ Ｘ Ｘ Ｘ")
break
else:
print(f"SCF Calculation is repeating ...")
Cpre.clear()
Cpre = listCnewReFloat.copy()
#New Energy :
#ε-W2p = -9.81440
#ε-W2p = 1.75868
#ε-W2p = 5.76657
#New φi :
# 0.4966 χ1 +0.7119 χ2 +0.4966 χ3
# 0.7071 χ1 +0.0000 χ2 -0.7071 χ3
# 0.5034 χ1 -0.7023 χ2 +0.5034 χ3
#Result 5 : Under 0.0001 all
#－－－－－－－－－－－－－－－－－－－
#Aryl Cation SCF Calculation Complete.
#
#
#New Energy :
#ε-W2p = -1.28964
#ε-W2p = 5.88678
#ε-W2p = 11.69602
#New φi :
# 0.5121 χ1 +0.6896 χ2 +0.5121 χ3
# 0.7071 χ1 +0.0000 χ2 -0.7071 χ3
# 0.4876 χ1 -0.7242 χ2 +0.4876 χ3
#Result 5 : Under 0.0001 all
#－－－－－－－－－－－－－－－－－－－
#Aryl Radical RHF SCF Calculation Complete.
アリルカチオン (Aryl Cation) および アリルラジカル (Aryl Radical UHF) とも、SCF 計算は、収束しました。アリルカチオンは、参考書１の pp296-297 と全く同じ値になります。アリルラジカルUHF は、参考書１の pp333-334 と同じ値になりませんでした。
※ プログラム下方のNew Cirについて： 0.5乗の符号の取り方は、下図に示した初期Huckel法の軌道の節（符号の反転する点）に沿っています。π電子近似計算ですから、各原子上のπ電子の符号を見ます。C11N,C12N,C13Nは符号＋。C21Nは符号＋、C22Nは零、C23Nは符号－。C31Nは符号＋、C32Nは符号－、C33Nは符号＋。これはノーベル化学賞を受賞された福井謙一先生のフロンティア軌道理論、HOMO、LUMOといった電子軌道の別分野の理論に発展していきます。
以上、SCF 計算プログラムを Python で作成しました。一部参考書と同じ値にはなりませんでしたが、SCF 計算は収束しました。
参考： SCF 計算の理論
F φ_i = ε_i φ_i \quad or \quad (F-ε_i)φ_i = 0
これが Roothaan の SCF 方程式 です。F は Fock の演算子です。$ε_i$ は分子軌道 $φ_i$ のエネルギーです。今回は分子軌道 $φ_i$ に LCAO（Liner Combination Atomic Orbital）を使います。次の式です。
φ_i = \sum_{r=1}^mC_{ir}χ_r
このLCAO型の分子軌道 $φ_i$ を上述の Roothaan の SCF 方程式 の右側の式に代入すると次のような連立方程式が導かれます。
\left.
\begin{array}{ll}
(F_{11}-ε_i)C_{i1} &+\,(F_{12}-ε_iS_{12})C_{i2} &+\quad... &+\,(F_{1m}-ε_iS_{1m})C_{im} = 0\\
(F_{12}-ε_iS_{12})C_{i1} &+\,(F_{22}-ε_i)C_{i2} &+\quad... &+\,(F_{2m}-ε_iS_{2m})C_{im} = 0\\
............. &............. &\:\,\,\quad... &............... \\
(F_{1m}-ε_iS_{1m})C_{i1} &+\,(F_{2m}-ε_iS_{2m})C_{i2} &+\quad ... &+\,(F_{mm}-ε_iS_{1m})C_{im} = 0\\
\end{array}
\right\}
次の行列式に書き換えることもできます。
\begin{vmatrix}
F_{11}-ε_i & F_{12}-ε_iS_{12} & ..... & F_{1m}-ε_iS_{1m} \\
F_{12}-ε_iS_{12} & F_{12}-ε_i & ..... & F_{2m}-ε_iS_{2m} \\
..... & ..... & ..... & ..... \\
F_{1m}-ε_iS_{1m} & F_{1m}-ε_iS_{2m} & ..... & F_{mm}-ε_i
\end{vmatrix}
\times
\begin{pmatrix}
C_{i1} \\
C_{i2} \\
... \\
C_{im}
\end{pmatrix}
= 0
ここで $S_{rs}$ は、ある原子軌道 $χ_r$ とある原子軌道 $χ_s$ 間の重なり積分です。今回はppp近似を使うので0です。
$F_{rs}$ は、Fock の行列要素といわれます。後述の全エネルギーの中の原子積分 $I_{rs} と (rs|tu)$ 、全電子密度 $P_{tt}$ 、結合次数 $P_{tu}$ を使うと次のように表せます。今回のアリルカチオンで使いました。
F_{rs} = I_{rs} + \sum_{t=1}^{m}\sum_{u=1}^{m}p_{tu}\left\{\left(rs|tu\right)-\left(rt|su\right)\right\}
一方、アリルラジカル RHF では、 全電子密度 $P_{tt}$ 、結合次数 $P_{tu}$ を使わず、$F_{rr}^{RHF}$ と $F_{rs}^{RHF}$ を求めます。閉殻系のアリルカチオンと同じπ電子近似を使います。ただし開殻系で2重項を持つラジカルでは、zero-differentialのπ電子近似と呼ばれるものを使います。$F_{rr}^{RHF}$ と $F_{rs}^{RHF}$ の数値は参考書 p333 の式 (6.118) である下式で求めます。
F_{rr}^{RHF} = \quad I_{rr} + \left[\sum_i^cC_{ir}^2
-\frac{1}{2}C_{mr}^2+2\sum_i^c(C_{ir}C_{mr})^2+C_{mr}^4\right]γ_{rr} \\
+ \sum_{s≠r} \left[2\sum_i^cC_{is}^2+C_{ms}^2 \\
+ C_{mr}C_{ms}\left\{2\sum_i^cC_{ir}C_{is}+C_{mr}C_{ms}\right\}\right]γ_{rs}
F_{rs}^{RHF} = \quad I_{rs} + \left[\sum_i^cC_{ir}C_{is}
+ \frac{1}{2}C_{mr}C_{ms}\right](C_{mr}^2γ_{rr}+C_{ms}^2γ_{ss}) \\
+ \left[-\sum_i^cC_{ir}C_{is}-\frac{3}{2}C_{mr}C_{ms} \\
+ \frac{1}{2}C_{mr}C_{ms}\left\{2\sum_i^cC_{ir}^2+C_{mr}^2+2\sum_i^cC_{is}^2+C_{ms}^2\right\}\right]γ_{rs} \\
+ \sum_{t≠r,\,s}C_{mt}C_{mr} \left[\sum_i^cC_{it}C_{is}+\frac{1}{2}C_{mt}C_{ms}\right]γ_{rt} \\
+ \sum_{t≠r,\,s}C_{mt}C_{ms} \left[\sum_i^cC_{it}C_{ir}+\frac{1}{2}C_{mt}C_{mr}\right]γ_{st}
上式で、$m$ は開殻系の分子軌道（ラジカルでは 1 個）。ここでのアリルラジカルでは分子軌道 $φ_2$ にだけ 1 個のラジカル電子が入っているので $m=2$ だけ。
$i$ は閉殻系の分子軌道。ここでのアリルラジカルでは分子軌道 $φ_1$ にだけ 2 個の電子が入っているので $i=1$ だけ。
$\sum_i^c$ は 2 電子被占のエネルギー準位について総和をとることを意味する。
つまりここでは 2 電子被占のエネルギー準位は ε1 だけであり、これは分子軌道 $φ_1$ のエネルギー準位である。
ゆえに、ここでのアリルラジカルでは 2 電子被占準位は分子軌道 $φ_1$ だけなので $i=1$ だけ。
これらの値を代入すると $F_{rr}^{RHF}$ と $F_{rs}^{RHF}$ は下式のようになります。
F_{rr}^{RHF} = \quad I_{rr} + \left[C_{1r}^2-\frac{1}{2}C_{2r}^2+2(C_{1r}C_{2r})^2+C_{2r}^4\right]γ_{rr} \\
+ \sum_{s≠r} \left[2C_{1s}^2+C_{2s}^2+ C_{2r}C_{2s}\left\{2C_{1r}C_{1s}+C_{2r}C_{2s}\right\}\right]γ_{rs}
F_{rs}^{RHF} = \quad I_{rs} + \left[C_{1r}C_{1s}
+ \frac{1}{2}C_{2r}C_{2s}\right](C_{2r}^2γ_{rr}+C_{2s}^2γ_{ss}) \\
+ \left[C_{1r}C_{1s}-\frac{3}{2}C_{2r}C_{2s} + \frac{1}{2}C_{2r}C_{2s}\left\{2C_{1r}^2+C_{2r}^2+2C_{1s}^2+C_{2s}^2\right\}\right]γ_{rs} \\
+ \sum_{t≠r,\,s}C_{2t}C_{2r} \left[C_{1t}C_{1s}+\frac{1}{2}C_{2t}C_{2s}\right]γ_{rt} \\
+ \sum_{t≠r,\,s}C_{2t}C_{2s} \left[C_{1t}C_{1r}+\frac{1}{2}C_{2t}C_{2r}\right]γ_{st}
これより $F_{11}^R$ ，$F_{11}^R$ ，$F_{12}^R$ ，$F_{13}^R$ を求めます。また 分子軌道 $φ_i$ はそれぞれ 3 項しかありません。
$F_{11}^R$ は r=1 で s=2, 3 です。$F_{22}^R$ は r=2 で s=1, 3 です。
$F_{12}^R$ は r=1, s=2 で t=3 だけです。$F_{13}^R$ は r=1, s=3 で t=2 だけです。
F_{11}^{R} = \quad I_{11} + \left[C_{11}^2-\frac{1}{2}C_{21}^2+2(C_{11}C_{21})^2+C_{21}^4\right]γ_{11} \\
+ \left[2C_{12}^2+C_{22}^2+ C_{21}C_{22}\left\{2C_{11}C_{12}+C_{21}C_{22}\right\}\right]γ_{12} \\
+ \left[2C_{13}^2+C_{23}^2+ C_{21}C_{23}\left\{2C_{11}C_{13}+C_{21}C_{23}\right\}\right]γ_{13}
F_{22}^{R} = \quad I_{22} + \left[C_{12}^2-\frac{1}{2}C_{22}^2+2(C_{12}C_{22})^2+C_{22}^4\right]γ_{22} \\
+ \left[2C_{11}^2+C_{21}^2+ C_{22}C_{21}\left\{2C_{12}C_{11}+C_{22}C_{21}\right\}\right]γ_{21} \\
+ \left[2C_{13}^2+C_{23}^2+ C_{22}C_{23}\left\{2C_{12}C_{13}+C_{22}C_{23}\right\}\right]γ_{23} \\
∵γ_{21}=γ_{23}=γ_{12}
F_{12}^{R} = \quad I_{12} + \left[C_{11}C_{12}
+ \frac{1}{2}C_{21}C_{22}\right](C_{21}^2γ_{11}+C_{22}^2γ_{22}) \\
+ \left[C_{11}C_{12}-\frac{3}{2}C_{21}C_{22} + \frac{1}{2}C_{21}C_{22}\left\{2C_{11}^2+C_{21}^2+2C_{12}^2+C_{22}^2\right\}\right]γ_{12} \\
+ C_{23}C_{21} \left[C_{13}C_{12}+\frac{1}{2}C_{23}C_{22}\right]γ_{13} \\
+ C_{23}C_{22} \left[C_{13}C_{11}+\frac{1}{2}C_{23}C_{21}\right]γ_{23} \\
∵γ_{23}=γ_{12}
F_{13}^{R} = \quad I_{13} + \left[C_{11}C_{13}
+ \frac{1}{2}C_{21}C_{23}\right](C_{21}^2γ_{11}+C_{23}^2γ_{33}) \\
+ \left[C_{11}C_{13}-\frac{3}{2}C_{21}C_{23} + \frac{1}{2}C_{21}C_{23}\left\{2C_{11}^2+C_{21}^2+2C_{13}^2+C_{23}^2\right\}\right]γ_{13} \\
+ C_{22}C_{21} \left[C_{12}C_{13}+\frac{1}{2}C_{22}C_{23}\right]γ_{12} \\
+ C_{22}C_{23} \left[C_{12}C_{11}+\frac{1}{2}C_{22}C_{21}\right]γ_{32} \\
∵γ_{32}=γ_{23}=γ_{12}
これらの $F_{11}^R$ ，$F_{11}^R$ ，$F_{12}^R$ ，$F_{13}^R$ にプログラム上部の原子積分の数値および係数 Cir を代入して使います。
参考書と同じ p296 の Hückel 法の係数から始めています。ただ参考書１ p334 の式 (6.119a) は、ある近似法で求めた表現になっていますが、数値自体は同じです。また今回は有効数字の取り方が異なるので、計算は収束しますが結果が少し違います。有効数字の取り方で係数が 0.0005 ズレると、途中計算で 2 乗して足して平方根をとり、この計算を 5 回繰り返すと、係数の誤差が 0.003 出ています。エネルギー準位 ε では 0.05 (eV) ほどの誤差がでます。全体値 1.234 (eV) に対して 0.05 (eV) ですから、4% 以上の誤差が出ています。計算結果に有効数字の取り方は大きな意味があります。ただ今回のような近似で小数点以下 4 桁目にどれほどの意味があるかについて疑問もあります。しかしプログラムでは Frs の式にある近似法を使わず、そのまま原子積分および係数の数値を小数点以下 4 桁で代入して計算しています。
SCF 計算のアルゴリズム
アリルカチオンの例で示します。
(1) 一組の分子軌道係数 $C_{ir}$ を仮定します。この $C_{ir}$ で全電子密度 $P_{tt}$ 、結合次数 $P_{tu}$ を計算します。ここでは p299 に示されるように φ1 の係数だけで計算します。
P_{tu} = 2C_{1t}C_{1u}
(2) $P_{tt}$ と $P_{tu}$ を代入して $F_{rs}$ を求めます。今回は PPP 近似で簡略化した次式を使います。
F_{rr} = I_{rr} + \frac{1}{2}p_{rr}(rr|rr) + \sum_{r≠s}P_{ss}(rr|ss)
F_{rs} = I_{rs} + \frac{1}{2}p_{rs}(rr|ss)
(3) $F_{rs}$ を連立方程式に代入し、連立方程式で $C_{ir}(r=1,2,3)$ のすべてがゼロにならない条件から次の永年方程式が得られます。
\begin{vmatrix}
F_{11}-ε_i & F_{12}-ε_iS_{12} & F_{13}-ε_iS_{13} \\
F_{21}-ε_iS_{21} & F_{22}-ε_i & F_{23}-ε_iS_{23} \\
F_{31}-ε_iS_{31} & F_{32}-ε_iS_{32} & F_{33}-ε_i
\end{vmatrix}
= 0
さらに PPP 近似で $S=0$ ですから、さらに次のように簡略化されます。
\begin{vmatrix}
F_{11}-ε_i & F_{12} & F_{13} \\
F_{21} & F_{22}-ε_i & F_{23} \\
F_{31} & F_{32} & F_{33}-ε_i
\end{vmatrix}
= 0
この永年方程式を解いて $ε_i$ を求めます。
さらに求めた εi を代入した連立方程式を解き、係数の規格化条件
C_{i1}^2 + C_{i2}^2 + C_{i3}^2 = 1
を使い、新しい $φi$ の新しい $C_{ir}$ を求めます。
判定：求めた新しい $C_{ir}$ と (1) で仮定した古い $C_{ir}$ とを比較します。一致していれば新しい $C_{ir}$ が確定で、SCF計算完了です。
収束した新しい分子軌道 φi が求まりましたので進みます。一致していなければ、(1)に戻り、新しい $C_{ir}$ で、古い $C_{ir}$ を置き換えて、(1)～(3)を繰り返します。
一致していれば、一応 SCF 分子軌道 φi のエネルギーと名づけられている $ε_i$ が確定します。
ここからはプログラム中には無いですが $ε_i$ は次式で表せます。
ε_i = H_i + \sum_{j=1}^{n}\left(2J_{ij} - K_{ij}\right)
ただし、この εi で分子全体の2n個の電子で合計すると、後述する全電子エネルギー E0 と一致しません。しかし全電子エネルギー $E_0$ と SCF 分子軌道 $φ_i$ のエネルギー $ε_i$ とは次の関係式が成立しています。第2項は電子間相互作用に関する項です。
E_0 = 2\sum_{i=1}^{n}ε_i - \sum_{i,j=1}^{n}\left(2J_{ij} - K_{ij}\right)
一応と前置きされているようにSCF 分子軌道 φi のエネルギー $ε_i$ は、全電子エネルギー $E_0$ と比較して、電子間相互作用を 2 倍見積もっています。Roothaan の SCF 計算をする場合は、この違いに注意が必要です。
参考：全電子エネルギー
測定される分子の電子スペクトルは、基底状態の電子配置から電子が1個励起され、励起状態の電子配置へ、この1個の電子が遷移した時、2つの状態のエネルギー準位の差
⊿E = hν, h:プランク定数 ν:振動数
が、振動数 $ν$ として観測されます。ゆえに理論および計算式が合っているのかどうか、まず基底状態および励起状態のエネルギー準位を計算することは非常に重要です。ここでは分子の全電子エネルギーの式を示します。
電子スピンおよび2電子間反発も考慮した分子の全エネルギー $E_0$ です。ただし核スピン、磁場、および電場の影響は考慮していません。
E_0 = 2\sum_{i=1}^{n}H_i + \sum_{i,j=1}^{n}\left(J_{ij} - K_{ij}\right) ：全電子エネルギー
参考書1, pp282-283 式(6.11) ~ 式(6.16)より
分子積分
H_i = \intφ_i\left(1\right)\left[-\frac{ℏ^2}{2m}⊿_1+V\left(1\right)\right]φ_i\left(1\right)dτ_1 ：分子コア積分
J_{ij} = \iintφ_i\left(1\right)φ_i\left(1\right)\frac{e^2}{r_{12}}φ_j\left(2\right)φ_j\left(2\right)dτ_1dτ_2 ：分子クーロン積分
K_{ij} = \iintφ_i\left(1\right)φ_j\left(1\right)\frac{e^2}{r_{12}}φ_i\left(2\right)φ_j\left(2\right)dτ_1dτ_2 ：分子交換積分
式中の$φ_i$は分子軌道です。
原子積分で表した積分
LCAO 法での分子軌道 $φ_i$ は、下式のような原子軌道 $χ_r$ の線型結合で表せます。これを上述の分子積分の式に代入します。
φ_i = \sum_{r=1}^{m}C_{ir}χ_r
さらに途中で下式の原子積分 $I_{rr}$、$I_{rs}$、$\left(rs|tu\right)$ を略記します。
I_{rr}= \intχ_r\left(1\right)\left[-\frac{ℏ^2}{2m}⊿_1+V\left(1\right)\right]χ_r\left(1\right)dτ_1
I_{rs}= \intχ_r\left(1\right)\left[-\frac{ℏ^2}{2m}⊿_1+V\left(1\right)\right]χ_s\left(1\right)dτ_1
\left(rs|tu\right)= \iintχ_r\left(1\right)χ_s\left(1\right)\frac{e^2}{r_{12}}χ_t\left(2\right)χ_u\left(2\right)dτ_1dτ_2
分子積分は、次のように原子積分で表せます。
H_i = \sum_{r=1}^{m}\left(C_{ir}\right)^2I_{rr} + 2\sum_{r>s}^{m}C_{ir}C_{is}I_{rs} ：コア積分
J_{ij} = \sum_{r=1}^{m}\sum_{s=1}^{m}\sum_{t=1}^{m}\sum_{u=1}^{m}C_{ir}C_{is}C_{jt}C_{ju}\left(rs|tu\right) ：クーロン積分
K_{ij} = \sum_{r=1}^{m}\sum_{s=1}^{m}\sum_{t=1}^{m}\sum_{u=1}^{m}C_{ir}C_{js}C_{it}C_{ju}\left(rs|tu\right) ：交換積分
さらに電子密度 $q_r=p_{rr}$ および結合次数 $p_{rs}$ 次のようにおき
q_r = p_{rr} = \sum_i^{occ}(C_{ir})^2, p_{rs} = \sum_i^{occ}C_{ir}C_{is}
全電子エネルギー $E_0$ の式に代入すると次のように表せます。
E_0 = \quad \sum_{r=1}^{m}p_{rr}I_{rr} + 2\sum_{r>s}^{m}p_{rs}I_{rs} \\
+ \sum_{r=1}^{m}\sum_{s=1}^{m}\sum_{t=1}^{m}\sum_{u=1}^{m}(p_{rs}p_{tu} - \frac{1}{2}p_{rt}p_{su})\left(rs|tu\right)
今回は、この式と似た手順で求めたFockの行列要素 $F_{rs}$ に、半経験的な手法で得た電子密度 $p_{rr}$、結合次数 $p_{rs}$、原子積分 $I_{rr}$、$I_{rs}$、$\left(rs|tu\right)$ の値を代入してSCF計算をしました。
参考：非経験的な解法
今回はRoothaanおよびPopleらのSCF計算アルゴリズムのプログラムの流れを作成しました。そして計算途中では半経験的な手法で得た電子密度 $p_{rr}$、結合次数 $p_{rs}$、原子積分 $I_{rr}$、$I_{rs}$、$\left(rs|tu\right)$ の定数値を使用しました。π 電子近似の半経験的な方法は参考書１の三訂量子化学入門(上)の p312-p317 をご覧ください。さらにCNDO 法および INDO 法などの σ 電子を含めた全価電子を考慮した半経験的 SCF 法については同 p317-p328 をご覧ください。
一方、非経験的な手法で原子積分の定数値を得ることもできます。これも SCF 計算アルゴリズムとは別個にあらかじめ求めておくものです。SCF 計算途中で使用する原子積分の定数値が、半経験的な手法でも、非経験的な手法でも、PPP 近似であれば重なり積分 S を０とする今回の SCF 計算アルゴリズムのプログラムが使えます。ただ非経験的な手法では重なり積分 S を０と近似しませんので少し複雑になります。非経験的な手法は、参考書２の三訂量子化学入門(下)の最初 p373 から始まります。
ここでは、今回は実施しなかった非経験的な手法で原子積分の定数値を得る方法を、代表的な STO-3G を例に少し述べます。現在の量子化学計算プログラムはほとんどガウス型でこの非経験的な手法で計算済の原子積分の定数値をあらかじめ既に持っています。変数がズレているのでまず 4 中心を 2 中心への変換が少し大変ですが、今回の積分にわざわざ導入したガウス型の特徴で項数が多いだけで積分自体はあっさり解け、解もシンプルです。一方、以下では分子軌道に使う基底関数がガウス型の例を示します。
STO-NG の一般式
STO \,\, 1s \,\, : \,\, χ_{1s}^{STO}(r,ζ_1) \,\, = \,\, \left(ζ_1^3/π\right)^{1/2}exp(-ζ_1r)
STO-NG \, 1s \,\,\,\, : \quad χ_{1s}^{STO-NG}(r,ζ_1) \quad = \quad\quad\quad \\
\sum_{k=1}^{N}d_{1sk}(2α_{1k}ζ_1^2/π)^{3/4} × exp [-α_{1k}(ζ_1r)^2]
参考書2, p382 式(7.10)より
次では具体例として、スレーター型である STO 1s と、ガウス型である STO-NG の N=3 の時つまり STO-3G 1s の例を見てみます。
STO-3G の具体例
STO \,\, 1s \,\, : \,\, χ_{1s}^{STO}(r, 1) \, = \, \left(1/π\right)^{1/2}e^{-r}
STO-3G \, 1s \,\,\,\, : \,\,\, χ_{1s}^{STO-3G}(r, 1) \,\, = \quad \\
0.154329(2×2.22766/π)^{3/4}e^{-2.22766r^2} \\
+ 0.535328(2×0.405771/π)^{3/4}e^{-0.405771r^2} \\
+ 0.444635(2×0.109818/π)^{3/4}e^{-0.109818r^2}
参考書2, p383 表7-2の値を、p382 式(7.10)に代入
スレーター型である STO 1s では水素原子の 1s 軌道を使います。項数は 1 個です。これに対してガウス型の STO-3G は項数が 3 個もあります。何故わざわざガウス型を使うのでしょうか。以下の引用を見てみます．
－－－ 参考書 2, pp379-380 より引用 －－－
ＧＴＯ の積分は ＳＴＯ に比べ，著しく簡単である．たとえば，Ａにある 1s 関数と Ｂにある 1s 関数との交換積分 $\left(ab|ab\right)$ の結果だけ見ると，ＳＴＯ では
\left(ab|ab\right) \quad = \\
\frac{1}{R}\left[\left(\frac{5}{8}a - \frac{23}{20}a^2 - \frac{3}{5}a^3 - \frac{13}{15}a^4\right)e^{-2a} \\
+ \frac{5}{6}\left\{\left(0.5772+log\,a\right)S^2(a) \\
-2E(-2a)\,S(a)\,S(-a) \\
+E(-4a)\,S^2(-a)\right\}\right]
ただし，a = ζR, R はＡＢ間の距離（bohr 単位)
S(u) = \left(1 + u + \frac{1}{3}u^2 \right)e^{-u}
E(-v) = -\int_{v}^{∞}\frac{e^{-t}}{t}dt
という複雑な形をしているが，ＧＴＯだと
\left(ab|ab\right) = 2\sqrt{\frac{α}{π}}e^{-αR^2}
と簡単になる。ＳＴＯでは，他の多中心積分の計算はますますたいへんになるのに比べ，ＧＴＯではｐやｄ関数が入っても，多中心積分でも大して複雑にならない．そのため，ごく大ざっぱに、言って，一つの二電子積分を計算する時間は，ＧＴＯ では ＳＴＯ の １/１０００ 以下ですんでしまう．
ＧＴＯ で項数を ＳＴＯ の 3 倍使うと仮定すると，必要な積分の数は $3^4 = 81$ 倍になる．しかし，それぞれが １/１０００ ですむので，積分の計算時間は，ＧＴＯ では ＳＴＯ の ８１/１０００ ～ １/１０ ですむことになる．というわけで，実際には ＧＴＯ が ＳＴＯ を駆逐してしまったのである．
－－－ 以上、参考書 2, pp379-380 より引用 －－－
ガウス型が導入された理由は、シンプルに積分時間が節約できるからです。
ここで引用文中に出てきた原子積分 (ab|ab) について、もう少しだけ詳しく見てみます。今回はアリルカチオンの 2 電子の π 電子近似でしたが、以下では水素分子の 2 電子を LCAO-MO 法で見てみます。
電子１ 電子２
◎－－－－◎
│ ∖
│ ∖
│ ∖
●－－－－－－－－●
χＡ χＢ
φ_+ = C_{A_+}χ_A + C_{B_+}χ_B : 被占分子軌道\\
φ_- = C_{A_-}χ_A - C_{B_-}χ_B : 空分子軌道
分子軌道は電子 1 と電子 2 が逆スピンで入った被占分子軌道 $φ_+$ と空軌道で表せます。ここでは積分部分だけに注目していますので、仮に各原子軌道の係数は１にします。また以下でスピン項は途中で消えるので省略しています。分子の波動関数は以下になります。
Φ=φ_{+(1)}φ_{+(2)}
水素分子のエネルギー $E$ は分子軌道および原子軌道で以下のように表せます。$Φ$ が実数関数であることより、その複素共役関数 $Φ^*$ も同じ実数関数になります。
E=\intΦ^*HΦdτ \\
= \iintφ_{+(1)}φ_{+(2)}Hφ_{+(1)}φ_{+(2)}dτ_1dτ_2 \\
= \iint\left(χ_{A(1)}+χ_{B(1)}\right)\left(χ_{A(2)}+χ_{B(2)}\right)H \\
\quad \quad \left(χ_{A(1)}+χ_{B(1)}\right)\left(χ_{A(2)}+χ_{B(2)}\right)dτ_1dτ_2 \\
∵H=-\frac{ℏ^2}{2m}(⊿_1+⊿_2) \\
+\frac{e^2}{4πε_0}(-\frac{1}{r_{A1}}-\frac{1}{r_{A2}}-\frac{1}{r_{B1}}-\frac{1}{r_{B2}}+\frac{1}{r_{AB}}+\frac{1}{r_{12}})
H はハミルトニアンです。ハミルトニアンの第 1 項は電子 1 と電子 2 の運動エネルギー項です。第 2 項が距離 r 離れた 2 粒子間の相互作用ポテンシャル項です。この中で問題となるのが第 2 項中の最後の $1/r_{12}$ です。エネルギー E の積分を個々に展開した時、以下の積分がここでの話題の (ab|ab) の原子積分です。本によっては <AB|BA> とすることもありますが同じものです。
\iintχ_{A(1)}χ_{B(2)}\frac{1}{r_{12}}χ_{B(1)}χ_{A(2)}dτ_1dτ_2 : <AB|BA> \\
\iintχ_{A(1)}χ_{B(1)}χ_{A(2)}χ_{B(2)}\frac{1}{r_{12}}dτ_1dτ_2 : (ab|ab)
χA に電子 1、χA に電子 2、χB に電子 1、χB に電子 2 である状態が混在しています。4 中心の 2 電子積分といわれ、引用文中にあった複雑な積分 (ab|ab) です。スレーター型の複雑な積分計算の場合、永遠に計算が終わらないかもしれません。ゆえに 2 原子分子まではスレーター型、3 原子分子以上はガウス型と使い分けているようです。
そこで導入されたのがガウス型です。個々の原子軌道 χ を基底関数と呼びかえて、基底関数 χA と χB 、それぞれに更にガウス型 LCAO を適用して次のようにおきます。ここでは積分部分だけに注目していますので、仮に各ガウス型の係数は１にしています。
χ_A = χ_{A_{Gα}} + χ_{A_{Gβ}} + χ_{A_{Gγ}} \\
χ_B = χ_{B_{Gα}} + χ_{B_{Gβ}} + χ_{B_{Gγ}}
これを(ab|ab)の積分に代入すると次のようになります。
(ab|ab) \quad
= \iintχ_{A(1)}χ_{B(2)}\frac{1}{r_{12}}χ_{B(1)}χ_{A(2)}dτ_1dτ_2 \\
= \iint
\left\{{χ_{A_{Gα}(1)}+χ_{A_{Gβ}(1)}+χ_{A_{Gγ}(1)}}\right\}
\left\{{χ_{B_{Gα}(2)}+χ_{B_{Gβ}(2)}+χ_{B_{Gγ}(2)}}\right\} \\
× \frac{1}{r_{12}} × \\
\left\{{χ_{B_{Gα}(1)}+χ_{B_{Gβ}(1)}+χ_{B_{Gγ}(1)}}\right\}
\left\{{χ_{A_{Gα}(2)}+χ_{A_{Gβ}(2)}+χ_{A_{Gγ}(2)}}\right\}dτ_1dτ_2
積分の数は増えますが、ガウス型の特徴で 4 中心を 2 中心にすることが可能です。展開すると 81 個の積分がありますが、その中の 1 つを次に示します。
\iintχ_{A_{Gα}(1)}χ_{B_{Gβ}(1)}χ_{A_{Gβ}(2)}χ_{B_{Gγ}(2)}\frac{1}{r_{12}}dτ_1dτ_2
ここからが、Samuel Francis Boys が始めた 4 中心積分を 2 中心積分に変換する公式です。1 例として、まずそれぞれのガウス型関数を、仮に次のように具体的に決めます。
χ_{A_{Gα}(1)}=e^{-α|r_1-A|^2} \quad ∵ |r_1-A|=r_{A1} \\
χ_{B_{Gβ}(1)}=e^{-β|r_1-B|^2} \quad ∵ |r_1-B|=r_{B1} \\
χ_{A_{Gβ}(2)}=e^{-β|r_2-A|^2} \quad ∵ |r_2-A|=r_{A2} \\
χ_{B_{Gγ}(2)}=e^{-γ|r_2-B|^2} \quad ∵ |r_2-B|=r_{B2} \\
そして電子 1 同士のガウス型関数の積をとり、公式を適用します（公式はうろ覚えなので間違っているかもしれません。あったはずの元文を捜索中です）。P というのは点 A と点 B の α：β の内分点（逆内分点 β：α かも！）。$r_{P1}$ は点 P と電子 1 の距離でこれが電子 1 の変数です。
χ_{A_{Gα}(1)}χ_{B_{Gβ}(1)} \\
=e^{-α|r_1-A|^2}e^{-β|r_1-B|^2}
=e^{-\frac{αβ}{α+β}|A-B|^2}×e^{-(α+β)|r_1-P|^2} \\
∵ |A-B|=r_{AB} \, , \quad |r_1-P|=r_{P1} \\
∵ P=(βA+αB)/(α+β) \,:\, interior \, division \, point \\
電子 2 でも同様に積をとり、公式を適用します。Q というのは点 A と点 B の β：γ の内分点です。$r_{Q2}$ は点 Q と電子 2 の距離でこれが電子 2 の変数です。
χ_{A_{Gβ}(2)}χ_{B_{Gγ}(2)} \\
=e^{-β|r_2-A|^2}e^{-γ|r_2-B|^2}
=e^{-\frac{βγ}{β+γ}|A-B|^2}×e^{-(β+γ)|r_2-Q|^2} \\
∵ |A-B|=r_{AB} \, , \quad |r_2-Q|=r_{Q2} \\
∵ Q=(γA+βB)/(β+γ) \,:\, interior \, division \, point \\
前述の積分に代入します。
\iintχ_{A_{Gα}(1)}χ_{B_{Gβ}(1)}χ_{A_{Gβ}(2)}χ_{B_{Gγ}(2)}\frac{1}{r_{12}}dτ_1dτ_2 \\
=\iint{e^{-α|r_1-A|^2}}e^{-β|r_1-B|^2}×e^{-β|r_2-A|^2}e^{-γ|r_2-B|^2}×\frac{1}{r_{12}}dτ_1dτ_2 \\
=e^{-\frac{αβ}{α+β}|A-B|^2}e^{-\frac{βγ}{β+γ}|A-B|^2}
\iint{e^{-(α+β)|r_1-P|^2}}e^{-(β+γ)|r_2-Q|^2}\frac{1}{r_{12}}dτ_1dτ_2
元の上側の２行は、核 A と核 B に、電子 1 と電子 2 が混在した 4 中心の 2 電子積分 (ab|ab) です。一方、最後の行は、核 A と核 B を α：β に内分した点 P に数学的に合成した擬似核 P と電子 1 だけの関数、核 A と核 B を β：γ に内分した点 Q に数学的に合成した擬似核 Q と電子 2 だけの関数、そして 1/ r12 の積の積分です。これは (PP|QQ) 、より (P|Q) と表現した方がいいかもしれませんが、2 中心電子間反発積分に変わっています。ガウス型関数の性質を応用した変換公式です。2 中心電子間反発積分であれば、81 個の積分があっても、個々の積分を解くことはそれほど困難ではないです（スレーター型に比べて）。ただ今回は指数関数の符号と係数に注意する事とあります。解けるはずですが、まだ解いていないので解けるかどうかはわかりません。
コロナ禍で、次々と官庁閉鎖の恐れのある状態にあがなえないながら、フラットな視点から小さな突破口をこじ開け、散った桜の涙を注入し、新緑あふれる並木道を歩けることを願って。
お忙しい中、最後までお読みいただき、ありがとうございました。
参考書
- 三訂 量子化学入門（上）, 米澤貞次郎ら，pp275-285 p292-299 pp310-317 pp328-337, 化学同人，1983
- 三訂 量子化学入門（下）, 米澤貞次郎ら，pp373-383, 化学同人，1983