運動エネルギー積分を詳細に計算します

2023年12月30日2024年11月18日

このページでは量子化学計算を行うにあたり、必要な基底関数（原子軌道とイメージしてください）の運動エネルギー積分を計算する詳細について解説します。

(s|-1/2∇²|s)について

まずはs軌道同士の運動量積分を計算します。位置$\boldsymbol{A}$にある原子Aのs軌道（$G(a,\boldsymbol{A})$）と位置$\boldsymbol{B}$にある原子Bのs軌道（$G(b,\boldsymbol{B})$）の運動量積分$(s|-\frac{1}{2}\nabla^2|s)$を考えます。

運動積分ではガウス関数を二階微分してもガウス関数の形（に係数をかけたもの）になることを利用します。

この関係を利用してs軌道同士の運動エネルギー積分を行ってみましょう。

$\nabla^2=\frac{\partial^2}{\partial{x^2}}+\frac{\partial^2}{\partial{y^2}}+\frac{\partial^2}{\partial{z^2}}$なので、$(s|-\frac{1}{2}\nabla^2|s)$のまえに$(s|-\frac{1}{2}\frac{\partial^2}{\partial{x^2}}|s)$を考えます。

ここで、$\xi=\frac{ab}{a+b}$です。

結局、$(s|-\frac{1}{2}\nabla^2|s)=(s|-\frac{1}{2}\frac{\partial^2}{\partial{x^2}}|s)+(s|-\frac{1}{2}\frac{\partial^2}{\partial{y^2}}|s)+(s|-\frac{1}{2}\frac{\partial^2}{\partial{z^2}}|s)$であるので、$(s|-\frac{1}{2}\nabla^2|s)$はつぎのように$(s|s)$を利用して簡単な形にまとめることができます。

$(s|-\frac{1}{2}\nabla^2|s)$は以下のサブルーチン「CALC_TSS」内のコードで表現しています。「VAL」が$(s|-\frac{1}{2}\nabla^2|s)$を表しており、サブルーチンの返り値です。「SSS」は$(s|s)$のことで、サブルーチンの引数として用いています。

      SUBROUTINE CALC_TSS(NDimI,NDimJ,PA,PB,PP,ZETA,GZAI,d2,SSS,VAL)

      IMPLICIT NONE
      INTEGER  NDimI,NDimJ
      DOUBLE PRECISION  ZETA,GZAI,d2,SSS,VAL
      DOUBLE PRECISION  PA(3),PB(3),PC(3),PD(3) 
      DOUBLE PRECISION  PP(3),PQ(3),PR(3) 

      VAL = GZAI* (3.0d0 - 2.0d0*GZAI*d2)*SSS

      ENDSUBROUTINE

(p|-1/2∇²|s)について

同様にp軌道とs軌道の運動エネルギー積分$(p_x|-\frac{1}{2}\nabla^2|s)$を考えていきましょう。

$(p_x|-\frac{1}{2}\nabla^2|s)$では$(p|s)$と同様に$(x-A_x)G(a,\boldsymbol{A})=\frac{1}{2a}\frac{\partial}{\partial{A_x}}\exp(-a\boldsymbol{r}_A)$を利用します。

さらに、$\nabla$が$\frac{\partial^2}{\partial{x^2}}+\frac{\partial^2}{\partial{y^2}}+\frac{\partial^2}{\partial{z^2}}$であることから丁寧に展開し、$(p_x|-\frac{1}{2}\nabla^2|s)$では$\frac{\partial^2}{\partial{y^2}}$と$\frac{\partial^2}{\partial{z^2}}$は微分結果が奇関数となり積分値に影響しないことに注意します。

さらに、$G(a,\boldsymbol{A})$と$G(a,\boldsymbol{B})$の積が$K_{AB}G(p,\boldsymbol{P})$となることから、$(x-P_x)$のべきで整理することで、$(x-P_x)$の1次と3次のべきは奇関数となることからこちらも積分値に影響しないこともうまく利用します。

最終的に計算結果は↓のようになります。

↑の計算式で(省略)としているのは$(P_x-A_x)(s|-\frac{1}{2}\nabla^2|s)$です。一行が長くなりすぎるので省略しています。

$(p|-\frac{1}{2}\nabla^2|s)$は以下のサブルーチン「CALC_TPS」内のコードで表現しています。$(p_x|s)$が$(P_x-A_x)\times(s|s)$で表せることを利用しています。

      SUBROUTINE CALC_TPS(NDimI,NDimJ,PA,PB,PP,
     1                              ZETA,GZAI,d2,SSS,TSS,VAL)

      IMPLICIT NONE
      INTEGER  NDimI,NDimJ
      DOUBLE PRECISION  ZETA,GZAI,d2,SSS,TSS,VAL
      DOUBLE PRECISION  PA(3),PB(3),PC(3),PD(3) 
      DOUBLE PRECISION  PP(3),PQ(3),PR(3) 

      VAL=
     1  (PP(NDimI)-PA(NDimI))*TSS
     1      +2.0d0*GZAI*(PP(NDimI)-PA(NDimI))*SSS

      ENDSUBROUTINE

(p|-1/2∇²|p)について

$(p_x|-\frac{1}{2}\nabla^2|p_y)$も$(p_x|-\frac{1}{2}\nabla^2|s)$と同じように奇関数の積分値が$0$になることを利用することで計算をおこなうことができます。

ただし、展開する項数がやはり多くなるので、プラスマイナスを含めてひとつひとつ丁寧に計算しないと間違えます。（私は計算に4時間かかりました。）計算の羅列が並びますので、載せるかすこし悩みましたが、計算の詳細を伝えるページですので、しっかり載せます。

皆さんも時間があればチャレンジしてみてください。特に、理論系研究を志す方は計算を実際にしてみてください。

計算量がやや多いので、見たい人だけ見てください。

$(p|-\frac{1}{2}\nabla^2|p)$はサブルーチン「CALC_PP」と「CALC_TPS」を利用した↓のコードで表現しています。コードにするとあっけないですが、結構計算はめんどくさいです。

            DO NDimI = 1,3
            DO NDimJ = 1,3
            IF((N BAS_TYPE(I_BAS).EQ.NDimI)
     1                           .AND.(NBAS_TYPE(J_BAS).EQ.NDimJ)) THEN
            CALL CALC_PP(NDimI,NDimJ,PA,PB,PP,ZETA,GZAI,d2,SSS,VAL2)
            CALL CALC_TPS(NDimI,0,PA,PB,PP,ZETA,GZAI,d2,SSS,TSS,VAL1)

            t_t = t_t + CXIJ * (
     1                     (PP(NdimJ)-PB(NDimJ))*VAL1+2.0d0*GZAI*VAL2)


            IF(NDimI.EQ.NDimJ) THEN
            t_t = t_t + CXIJ *0.50d0/ZETA*TSS
            ENDIF
            ENDIF


            ENDDO
            ENDDO