【物理】ハートリーフォック法・p軌道2電子積分

前回の計算は間違っていました。ティッセン本と同じやり方（フーリエ変換とデルタ関数をうまく使う方法）
でp軌道の2電子積分を計算しました。テクニックとしてはフーリエ変換、デルタ関数、部分積分、球面座標変換を
使いますが結構面倒です。

結果は次の通りです。
$\left\langle 1p_{x},\gamma;1p_{x},\alpha|\frac{1}{|r_{1}-r_{2}|}|1p_{x},\beta;1p_{x},\delta\right\rangle=2\left(\frac{pq}{p+q}\right)^{1/2}\pi^{1/2}\left(\frac{\pi}{p}\right)^{3/2}\left(\frac{\pi}{q}\right)^{3/2}\left(\frac{1}{4pq}+\frac{3}{20(p+q)^{2}}-\frac{1}{12}\left(\frac{1}{p(p+q)}-\frac{1}{q(p+q)}\right)\right)$
p軌道系はpx,py,pzの3種類あるので、それらの2電子積分は以下となります。
$\left\langle 1p_{x},\gamma;1p_{x},\alpha|\frac{1}{|r_{1}-r_{2}|}|1p_{y},\beta;1p_{x},\delta\right\rangle=2\left(\frac{pq}{p+q}\right)^{1/2}\pi^{1/2}\left(\frac{\pi}{p}\right)^{3/2}\left(\frac{\pi}{q}\right)^{3/2}\left(\frac{1}{4pq}+\frac{1}{20(p+q)^{2}}-\frac{1}{12}\left(\frac{1}{p(p+q)}-\frac{1}{q(p+q)}\right)\right)$
$\left\langle 1p_{y},\gamma;1p_{x},\alpha|\frac{1}{|r_{1}-r_{2}|}|1p_{x},\beta;1p_{y},\delta\right\rangle=2\left(\frac{pq}{p+q}\right)^{1/2}\pi^{1/2}\left(\frac{\pi}{p}\right)^{3/2}\left(\frac{\pi}{q}\right)^{3/2}\left(\frac{1}{20(p+q)^{2}}\right)$

プログラムも置きます。VC++2008(無料)で開けます。
Ne原子の基底状態のエネルギーとKoopmanの定理で第１イオン化エネルギーを計算します。後者は微妙です。
http://g0307.hp.infoseek.co.jp/hartree_fock.zip

実行するとこんな感じです。

■上式の導出

原子の場合のｐ軌道系の2電子積分は以下となる。
$\displaystyle \left\langle 1p,\gamma;1p,\alpha|\frac{1}{|r_{1}-r_{2}|}|1p,\beta;1p,\delta\right\rangle=\int {}^{3}dr_{1}{}^{3}dr_{2}\left(x_{1}e^{-\alpha r_{1}^{2}}\times x_{2}e^{-\gamma r_{2}^{2}}\right)\frac{1}{|r_{1}-r_{2}|}\left(x_{1}e^{-\beta r_{1}^{2}}\times x_{2}e^{-\delta r_{2}^{2}}\right)$
$ p=\alpha+\beta,q=\gamma+\delta$とすると、
$\displaystyle \left\langle 1p,\gamma;1p,\alpha|\frac{1}{|r_{1}-r_{2}|}|1p,\beta;1p,\delta\right\rangle=\int {}^{3}dr_{1}{}^{3}dr_{2}\left(x_{1}^{2}e^{-pr_{1}^{2}}\frac{1}{|r_{1}-r_{2}|}x_{2}^{2}e^{-qr_{2}^{2}}\right)$
となる。
$x_{1}^{2}e^{-pr_{1}^{2}}$のフーリエ逆変換は$e^{-\frac{k^{2}}{4p}}\left(\frac{\pi}{p}\right)^{3/2}\left(\frac{k_{x}^{2}}{4p^{2}}-\frac{1}{2p}\right)$で、$\displaystyle \frac{1}{|r_{1}-r_{2}|}$のフーリエ逆変換は$4\pi k^{-2}$より、
$\displaystyle \left\langle 1p,\gamma;1p,\alpha|\frac{1}{|r_{1}-r_{2}|}|1p,\beta;1p,\delta\right\rangle=\int {}^{3}dr_{1}{}^{3}dr_{2}\left(\int d^{3}k_{1}e^{-\frac{k_{1}^{2}}{4p}}\left(\frac{\pi}{p}\right)^{3/2}\left(\frac{k_{1,x}^{2}}{4p^{2}}-\frac{1}{2p}\right)e^{-ik_{1}r_{1}}\right)$
$\times\left(\int d^{3}k4\pi k^{-2}e^{-ik(r_{1}-r_{2})}\right)$
$\times\left(\int d^{3}k_{3}e^{-\frac{k_{3}^{2}}{4q}}\left(\frac{\pi}{q}\right)^{3/2}\left(\frac{k_{3,x}^{2}}{4q^{2}}-\frac{1}{2q}\right)e^{-ik_{3}r_{2}}\right)$
となる。
ここで、デルタ関数$\displaystyle \delta(r)=(2\pi)^{-3}\int d^{3}r e^{-ikr}$を利用する。
例えば、$\displaystyle \delta(r_{1})=(2\pi)^{-3}\int d^{3}r_{1} e^{-ik_{1}r_{1}}$を代入して$r_{1}$と$r_{2}$に関する積分項を減らせる。
また、デルタ関数の作用$\displaystyle \int d^{3}r f(r)\delta(r-a)=f(a)$を用いると、
$\displaystyle \left\langle 1p,\gamma;1p,\alpha|\frac{1}{|r_{1}-r_{2}|}|1p,\beta;1p,\delta\right\rangle=4\pi\left(\frac{\pi}{p}\right)^{3/2}\left(\frac{\pi}{q}\right)^{3/2}(2\pi)^{-6}\int d^{3}k\left(\frac{k_{x}^{2}}{4p^{2}}-\frac{1}{2p}\right)\left(\frac{k_{x}^{2}}{4q^{2}}-\frac{1}{2q}\right)k^{-2}e^{-k^{2}\left(\frac{p+q}{4pq}\right)}$
となる。
積分の中は以下のように展開できる(変数を$r=k,x=k_{x}$とした)。
$\displaystyle \int e^{-r^{2}\left(\frac{p+q}{4pq}\right)}\left(\frac{x^{4}}{16p^{2}q^{4}}-\frac{x^{2}}{8pq^{2}}-\frac{x^{2}}{8p^{2}q}+\frac{1}{4pq}\right)r^{-2}$
それぞれの項目に対して球面座標に展開して積分を行うと、
$x=r\sin(\theta)\cos(\phi)$
より、
第1項目は、
$\displaystyle \int e^{-\left(\frac{p+q}{4pq}\right)r^{2}}x^{4} d^{3}r$
$=\displaystyle \int e^{-\left(\frac{p+q}{4pq}\right)r^{2}}r^{4} dr\int\left(\sin(\theta)\cos(\phi)\right)^{4}\sin(\theta)d\theta d\phi$
$=12\displaystyle \pi^{1/2}\left(\frac{pq}{p+q}\right)^{5/2}\times\frac{4}{5}\pi$
第２，３項目は、
$\displaystyle \int e^{-\left(\frac{p+q}{4pq}\right)r^{2}}x^{2} d^{3}r$
$=\displaystyle \int e^{-\left(\frac{p+q}{4pq}\right)r^{2}}r^{2} dr\int\left(\sin(\theta)\cos(\phi)\right)^{2}\sin(\theta)d\theta d\phi$
$=2\displaystyle \pi^{1/2}\left(\frac{pq}{p+q}\right)^{3/2}\times\frac{4}{3}\pi$
第４項目は、
$\displaystyle \int e^{-\left(\frac{p+q}{4pq}\right)r^{2}}d^{3}r$
$=\displaystyle \int e^{-\left(\frac{p+q}{4pq}\right)r^{2}}dr\int\sin(\theta)d\theta d\phi$
$=\pi^{1/2}\left(\frac{pq}{p+q}\right)^{1/2}\times 4\pi$
となる。

以上をまとめると以下の結果が得られる。
$\left\langle 1p_{x},\gamma;1p_{x},\alpha|\frac{1}{|r_{1}-r_{2}|}|1p_{x},\beta;1p_{x},\delta\right\rangle=2\left(\frac{pq}{p+q}\right)^{1/2}\pi^{1/2}\left(\frac{\pi}{p}\right)^{3/2}\left(\frac{\pi}{q}\right)^{3/2}\left(\frac{1}{4pq}+\frac{3}{20(p+q)^{2}}-\frac{1}{12}\left(\frac{1}{p(p+q)}-\frac{1}{q(p+q)}\right)\right)$
また、p軌道系はpx,py,pzの3種類あるので、それらの2電子積分は以下となる。
$\left\langle 1p_{x},\gamma;1p_{x},\alpha|\frac{1}{|r_{1}-r_{2}|}|1p_{y},\beta;1p_{x},\delta\right\rangle=2\left(\frac{pq}{p+q}\right)^{1/2}\pi^{1/2}\left(\frac{\pi}{p}\right)^{3/2}\left(\frac{\pi}{q}\right)^{3/2}\left(\frac{1}{4pq}+\frac{1}{20(p+q)^{2}}-\frac{1}{12}\left(\frac{1}{p(p+q)}-\frac{1}{q(p+q)}\right)\right)$
$\left\langle 1p_{y},\gamma;1p_{x},\alpha|\frac{1}{|r_{1}-r_{2}|}|1p_{x},\beta;1p_{y},\delta\right\rangle=2\left(\frac{pq}{p+q}\right)^{1/2}\pi^{1/2}\left(\frac{\pi}{p}\right)^{3/2}\left(\frac{\pi}{q}\right)^{3/2}\left(\frac{1}{20(p+q)^{2}}\right)$

終わり。

【補足】ガウス関数形のフーリエ変換

ガウス関数形のフーリエ変換(1次元)：

■$e^{-px^{2}}$のフーリエ変換
$\displaystyle \int e^{-px^{2}}e^{-ixk}dx=\left(\frac{\pi}{p}\right)^{1/2}e^{-\frac{k^{2}}{4p}}$

■$xe^{-px^{2}}$のフーリエ変換
$ \int xe^{-px^{2}}e^{-ixk}dx=\int(-\frac{1}{2p}e^{-px^{2}})^{\prime}e^{-ixk}dx$
$=-\frac{ik}{2p}\int e^{-px^{2}}e^{-ixk}dx$
$=-\frac{ik}{2p}\left(\frac{\pi}{p}\right)^{1/2}e^{-\frac{k^{2}}{4p}}$

■$x^{2}e^{-px^{2}}$のフーリエ変換
$\int x^{2}e^{-px^{2}}e^{-ixk}dx=\int(-\frac{1}{2p}e^{-px^{2}})^{\prime}xe^{-ixk}dx$
$=-[\displaystyle \frac{1}{2p}e^{-px^{2}}xe^{-ixk}]_{-\infty}^{\infty}-\frac{1}{2p}\int e^{-px^{2}}(e^{-ixk}-ikxe^{-ixk})dx$
$=-\displaystyle \frac{1}{2p}\left(\frac{\pi}{p}\right)^{1/2}e^{-\frac{k^{2}}{4p}}+\frac{ik}{2p}\int e^{-px^{2}}xe^{-ixk}dx$
$=-\displaystyle \frac{1}{2p}\left(\frac{\pi}{p}\right)^{1/2}e^{-\frac{k^{2}}{4p}}-\frac{ik}{2p}\frac{ik}{2p}\left(\frac{\pi}{p}\right)^{1/2}e^{-\frac{k^{2}}{4p}}$
$=-\displaystyle \frac{1}{2p}\left(\frac{\pi}{p}\right)^{1/2}e^{-\frac{k^{2}}{4p}}+\frac{k^{2}}{(2p)^{2}}\left(\frac{\pi}{p}\right)^{1/2}e^{-\frac{k^{2}}{4p}}$
$=e^{-\frac{k^{2}}{4p}}\left(\frac{\pi}{p}\right)^{1/2}\left(\frac{k^{2}}{4p^{2}}-\frac{1}{2p}\right)$

これより、3次元ガウス関数系のフーリエ変換は一次元ガウス関数のフーリエ変換結果の積より以下となる。
$\displaystyle \int x^{2}e^{-pr^{2}}e^{-irk}d^{3}r=e^{-\frac{k^{2}}{4p}}\left(\frac{\pi}{p}\right)^{3/2}\left(\frac{k_{x}^{2}}{4p^{2}}-\frac{1}{2p}\right)$
$\displaystyle \int x^{2}e^{-qr^{2}}e^{-irk}d^{3}r=e^{-\frac{k^{2}}{4q}}\left(\frac{\pi}{q}\right)^{3/2}\left(\frac{k_{x}^{2}}{4q^{2}}-\frac{1}{2q}\right)$