量子化学において量子効果を生む源泉を知っていますか？答えは二電子積分にあります。二電子積分があることでクーロン力からは発生しえない量子効果が生まれます。
しかし、量子化学計算ソフトの中で二電子積分がどのようにして計算されているか知っている人は少ないのではないでしょうか？
このページでは量子化学の二電子積分を計算方法を解説します。二電子積分は計算量を爆発的に増大させている根本原因であると同時に、量子効果を生むゆえんでもあります。
- 二電子積分が重要な理由
- ガウス関数の二電子積分を計算する方法
二電子積分の大変なところ
量子化学計算で二電子積分がネックになる原因は、「基底関数が$10$倍になると、必要な電子積分は$10^4=10,000$倍に増える。」というものです
対称性を利用すると、この$1/8$程度に収まりますが、いずれにしても膨大です。
例えば、ベンゼンの量子化学計算でメモリが$8$GB必要だったとします。メチル基を増やしてトルエンの量子化学計算を行う際に、もし基底関数が$1.3$倍になった場合、メモリは$8\times1.3^4=22.8$GBまで必要になってしまいます。
この時点で一般的なノートパソコンでは計算ができないね。
一方で、二電子積分を形成する交換積分は量子的効果から出てくるもので、量子化学の量子効果を生む源泉です。
つまり、二電子積分は量子化学計算を行う上で避けては通れないのか。
では、二電子積分をプログラミングに乗せるための計算表式を得る手続きを見ていきましょう。
s軌道の二電子積分について
早速、$s$型原子軌道の二電子積分$(ss|ss)$を考えましょう。二電子積分の計算は次の3つのポイントで解決することができます。
- ガウス関数の積もまたガウス関数
- $1/r$はラプラス変換する
- 誤差関数を用いる
結論として、指数係数が$a,b,c,d$、中心座標が$r_A,r_B,r_C,r_D$の4つの$s$軌道$G(a,A)～G(d,D)$による二電子成分$(ss|ss)$の結果がこちらです。
いきなり結論を言われても、訳が分からないよ。
順番に説明していきますね。まずは「ガウス関数の積もまたガウス関数」を使います。
ガウス関数の積もまたガウス関数
指数のべき数が$a$、中心座標が$r_A$で与えられる$s$軌道は$G(a)=\exp(-a({r-r_A})^2)$で与えることができます。このガウス関数については、指数の性質から、ガウス関数の積もまたガウス関数です。
そのため、$(ss|ss)=\int{G(a)G(b)}\frac{1}{r}G(c)G(d)\text{d}r$は新しい指数と中心座標を用いて$(ss|ss)=\int{G(p)}\frac{1}{r}G(q)\text{d}r$と簡略化できます。
一気に式の長さが半分くらいなったけど、まだ$\frac{1}{r}$の処理が残っているよ。。
そこで、つぎは「$1/r$はラプラス変換する」を使います。
$1/r$はラプラス変換しよう
一見すると、$\frac{1}{r}$はやっかいな形ですが、ラプラス変換すると$\frac{1}{r}=\int^\infty_0e^{-r^2u^2}\text{d}u$と変換できます。
よくみると、$e^{-r^2u^2}$もガウス関数だね。
その通りです。二電子積分も結局ガウス関数の積分に落ち着きます。
$(ss|ss)$における$\frac{1}{r}$の部分をラプラス変換させ、式をまとめると、最終的には$(ss|ss)=K\int^\infty_0(\frac{\rho}{\rho+u^2})^{3/2}e^{-\frac{\rho{u^2}}{\rho+u^2}}\text{d}u$（$K$は定数）と指数関数の積分の形に落とし込むことができます。
こんな複雑な積分、実際にできるの？
じつは、$u^2=\frac{\eta{u^2}}{1-t^2}$と置換積分を行うとすごいきれいに積分できるよ。
誤差関数$F(\alpha)$を用いる
置換積分でうまく積分を行うと$(ss|ss)=K\int^1_0\exp(-\alpha{t^2})\text{d}t$の形まですっきりと変形できます。
この$\int^1_0\exp(-\alpha{t^2})\text{d}t$部分はガウス積分の一部を切り取ったものです。
ガウス積分は区間$[0～\infty]$だけど、この積分は$[0～1]$なんだね。
この中途半端な積分は誤差関数やBoys関数と呼ばれています。「$\int^1_0\exp(-\alpha{t^2})\text{d}t=F(\alpha)$」と表すことで、$(ss|ss)=K\times{F}(\alpha)$（$K$は定数）と二電子積分の最終的な標識を得ることができます。
実際の計算式をごらんになりたい方は下のアコーディオンを展開してください。
二電子積分(ss|ss)の詳細な証明はこちら
計算がやや細かいですが、重要な結論として二電子積分もガウス関数の性質をうまく使うことで計算することができます。
二電子積分でもエラー関数$F_0(x)=\int^1_0\exp(-x^2)\text{d}x$の積分が必要です。この部分はBOYS関数とも呼ばれ、繰り返し計算を行うことで精度よく計算できます。
BOYS関数のプログラミング例はこちら
BOYS関数は次のようにプログラミングすることで得ることができます。（$x$の部分を$ddd$としており、0次の結果$\int^1_0\exp(-x^2)\text{d}x$をF0であたえ、1次の積分結果$\int^1_0(\text{d}/\text{d}(x^2))\exp(-x^2)\text{d}x$をF1であたえます。）
SUBROUTINE BOYS(ddd,F0,F1,F2,F3,F4)
IMPLICIT NONE
INTEGER I,J
DOUBLE PRECISION ddd,PI,F0,F1,F2,F3,F4
DOUBLE PRECISION ZP,ZQ,ZW,ZX
F0 = 0.0d0
F1 = 0.0d0
F2 = 0.0d0
F3 = 0.0d0
F4 = 0.0d0
I = 0
IF(ddd.LE.50.0d0) THEN
ZP = 1.0d0
ZQ = (2.0d0*REAL(0)+1.0d0)
ZW = 0.0d0
ZX = 1.0d0
DO I= 1, 100
ZP = ZP * 2.0d0 *ddd
ZQ = ZQ * (2.0d0*REAL(I)+1.0d0)
ZW = ZP / ZQ
ZX = ZX + ZW
IF(ZW.LE.(1.0d0/(10.0d0**7.0d0))) THEN
EXIT
ENDIF
ENDDO
F0 = exp(-ddd)*ZX
ZP = 1.0d0
ZQ = (2.0d0*REAL(1)+1.0d0)
ZW = 0.0d0
ZX = 1.0d0/(2.0d0*REAL(1)+1.0d0)
DO I= 1, 100
ZP = ZP * 2.0d0 *ddd
ZQ = ZQ * (2.0d0*REAL(I)+2.0d0*REAL(1)+1.0d0)
ZW = ZP / ZQ
ZX = ZX + ZW
IF(ZW.LE.(1.0d0/(10.0d0**7.0d0))) THEN
EXIT
ENDIF
ENDDO
F1 = exp(-ddd)*ZX
ZP = 1.0d0
ZQ = (2.0d0*REAL(2)+1.0d0)
ZW = 0.0d0
ZX = 1.0d0/(2.0d0*REAL(2)+1.0d0)
DO I= 1, 100
ZP = ZP * 2.0d0 *ddd
ZQ = ZQ * (2.0d0*REAL(I)+2.0d0*REAL(2)+1.0d0)
ZW = ZP / ZQ
ZX = ZX + ZW
IF(ZW.LE.(1.0d0/(10.0d0**7.0d0))) THEN
EXIT
ENDIF
ENDDO
F2 = exp(-ddd)*ZX
ZP = 1.0d0
ZQ = (2.0d0*REAL(3)+1.0d0)
ZW = 0.0d0
ZX = 1.0d0/(2.0d0*REAL(3)+1.0d0)
DO I= 1, 100
ZP = ZP * 2.0d0 *ddd
ZQ = ZQ * (2.0d0*REAL(I)+2.0d0*REAL(3)+1.0d0)
ZW = ZP / ZQ
ZX = ZX + ZW
IF(ZW.LE.(1.0d0/(10.0d0**7.0d0))) THEN
EXIT
ENDIF
ENDDO
F3 = exp(-ddd)*ZX
ZP = 1.0d0
ZQ = (2.0d0*REAL(4)+1.0d0)
ZW = 0.0d0
ZX = 1.0d0/(2.0d0*REAL(4)+1.0d0)
DO I= 1, 100
ZP = ZP * 2.0d0 *ddd
ZQ = ZQ * (2.0d0*REAL(I)+2.0d0*REAL(4)+1.0d0)
ZW = ZP / ZQ
ZX = ZX + ZW
IF(ZW.LE.(1.0d0/(10.0d0**7.0d0))) THEN
EXIT
ENDIF
ENDDO
F4 = exp(-ddd)*ZX
ELSE
PI = dacos(-1.0d0)
F0 = 0.50d0*dsqrt(PI/ddd)
F1 = 1.0d0/( 4.0d0*ddd )*dsqrt(PI/ddd)
F2 = 3.0d0/( 8.0d0*ddd**2)*dsqrt(PI/ddd)
F3 = 15.0d0/(16.0d0*ddd**3)*dsqrt(PI/ddd)
F4 =105.0d0/(32.0d0*ddd**4)*dsqrt(PI/ddd)
ENDIF
END(ss|ss)についてプログラミングするとこちらのようになります（プログラムから該当部分を抜粋しています）。
CXIJKL= CX(I_BAS,I)*CX(J_BAS,J)*CX(K_BAS,K)*CX(L_BAS,L)
d2_pq = (PP(1)-PQ(1))**2.0d0
1 + (PP(2)-PQ(2))**2.0d0
1 + (PP(3)-PQ(3))**2.0d0
Z_ab = dsqrt(2.0d0)*PI**(1.25d0)/ZETA1*EXP(-GZAI1*d2_ab)
Z_cd = dsqrt(2.0d0)*PI**(1.25d0)/ZETA2*EXP(-GZAI2*d2_cd)
CALL BOYS(RHO*d2_pq,F0,F1,F2,F3,F4)
SSSS(0) = 1.0d0/dsqrt(ZETA1+ZETA2)*Z_ab*Z_cd*F0
IF((NBAS_TYPE(I_BAS).EQ.0).AND.
1 (NBAS_TYPE(J_BAS).EQ.0).AND.
1 (NBAS_TYPE(K_BAS).EQ.0).AND.
1 (NBAS_TYPE(L_BAS).EQ.0)) THEN
C < s1s1 | s1s1 >
e_t = e_t + CXIJKL*SSSS(0)
ENDIF
p軌道を含む二電子積分について
p軌道までを量子化学計算で扱う場合、必要となる二電子積分の種類は$(ss|ss)$、$(p_is|ss)$、$(p_is|p_js)$、$(p_ip_j|ss)$、$(p_ip_j|p_ks)$、$(p_ip_j|p_kp_l)$の6種類（$i,j,k,l$は$x,y,z$のいずれか）です。
$p$軌道を含む二電子積分の注意点としては、これまでは$0$次の誤差関数$F_0(\alpha)=\int^1_0\exp(-\alpha{t^2})\text{d}t$しか登場してきませんでしたが、p軌道が関与すると高次の誤差関数が必要になってきます。
例えば1次の誤差関数$F_1(\alpha)$は$\int^1_0t^2\exp(-\alpha{t^2})\text{d}t$だよ。
$(p_xp_x|p_yp_y)$、$(p_xp_y|p_zp_z)$など4つのp軌道が関与する二電子積分を計算するためには最終的には$4$次の誤差関数$F_4(x)$まで必要になってきます。
まとめ
二電子積分の標識を得る途中の計算式は多少（かなり？）複雑でしたが、プログラミングに用いる最終的な標識はすっきりした形で与えられることがわかったのではないでしょうか？
業務として量子化学計算を行っていて、プログラミングまでは必要ない方も、「量子化学計算では原子軌道（基底関数）としてガウス関数を用いているので、二電子積分なども簡単にできる。」と理解しておいていただければと思います。
