【物理】クーロン積分

原子の周りの電子の密度を計算しています。これが計算できれば有効ポテンシャル
（クーロンポテンシャルが電子間の相互作用の影響でポテンシャルの形が変わる）が計算できます。
電磁波と原子とのトムソン散乱断面積は電子密度のフーリエ変換から計算できるらしい。
つまり、可視光の範囲で計算すれば各波長毎の散乱断面積、吸収・反射率といった
ＢＲＤＦに相当するものが原子番号から計算できるはずです。分子・結晶とかになるとまた難しいですが。
とりあえず大きな目標は、金属に光をあててレンダリングして、Auだったら金色に光って
Agだったら銀色に光るといった計算ができれば万々歳なのですが・・・。
リアルフォトンマップになります。無理かもしれません。
　
原子1個の場合のハートリーフォック方程式（ＨＦ方程式）は次となります。
H \phi(x)=E\phi(x), H:=-\frac{1}{2}\nabla^2 - \frac{Z}{|{\bf r}|}+\int d^3r |\phi({\bf r})|^2/|{\bf r}|
ハミルトニアンHの第一項は運動ポテンシャル、第二項は原子-電子間のクーロンポテンシャル、
第三項は電子-電子間のクーロン相互作用です。
　
状態関数を\phi({\bf r})=C e^{-\alpha r^2},{\bf r}=(x,y,z)^T(1s軌道のみ)と
ガウス関数で近似した場合を考える。この場合、ＨＦ方程式に\phiを
掛けて積分を取れば、係数Cに関する行列固有方程式の形になるので解けます。
　
重なり項、運動エネルギー項は上式を代入してガウス積分すれば行列要素が計算できますが、
クーロンポテンシャルのところは分母に\bf rが出てくるため、r=0の時
発散して積分できないのでは？と困ってましたが、うまく解析的に計算できるようです。
　
今回、1sの場合のクーロンポテンシャルの部分の積分（クーロン積分というらしい）を
計算してみます。結果は分かっていて\alpha=0.5のときに2\piになれば
いいです。
前回の日記のリンクの本の方はフーリエ変換とデルタ関数使って計算してるようでしたが、
グリーンの定理を使って計算してみました。
　
F(r)=\int e^{-r'^2} \frac{1}{|r-r'|} d^3r'

グリーンの定理により実はこの積分は次のポアソン方程式を解けばいいらしい。
\nabla^2 F(r)=-4\pi e^{-r^2}
極座標系のラプラシアンは
\nabla^2 = \frac{\partial^2}{\partial r^2}+\frac{1}{r}\frac{\partial}{\partial r}+\frac{1}{r^2}\Lambda
となる。F(r)は球対称なので最後の作用は無視できて、
\frac{1}{r}\frac{\partial^2(rF(r))}{\partial r^2}=-4\pi e^{-r^2}
に帰着する。
ここで、U(r)=rF(r)と置いてU(r)について解いて最終的に
F(r)=(erf(r)\pi^{3/2})/r
を得る。特に、\lim_{r\rightarrow 0}erf(r)/r=2/\pi^{1/2}より、
F(0)=2\pi
こういう計算は普段やらないので大変でした。数式tex打つのに時間かかりすぎました。
ＨＦ方程式の第3項の電子-電子間の積分はF(r)をさらに積分しないといけないですし、
p軌道,d軌道の計算もしないといけません。p軌道も同じような感じで計算できそうです。
　

■GTOのs,p軌道系の計算まとめ(11/2)
maximaに助けてもらいながら地味にp軌道系までの各種積分計算してました。
意外とこういう軌道ごとの直接の計算式は本には載ってないようです。
間違ってたらすいません。d軌道となると計算が複雑すぎて手動では無理そうです。
p,qはパラメータで、a=p+q,b=r+s,1s,1pは軌道を示してます。

【重なり積分】
<1s,p|1s,q>=\int d^3r\left{ e^{-p r^2}e^{-q r^2}\right}=\frac{\pi^{3/2}}{(p+q)^{3/2}}

<1p,p|1p,q>=\int d^3r\left{ x^2e^{-p r^2}e^{-q r^2}\right}=\frac{\pi^{3/2}}{2(p+q)^{5/2}}

【運動エネルギー積分】
<1s,p|\nabla^2|1s,q>=\int d^3r\left{ \nabla e^{-p r^2} \nabla e^{-q r^2}\right}=\frac{6\pi^{3/2}pq}{(p+q)^{5/2}}

<1p,p|\nabla^2|1p,q>=\frac{5\pi^{3/2}pq}{(p+q)^{7/2}}

【クーロン積分】
<1s,p|\frac{1}{|r-r'|}|1s,q>=\frac{erf(a^{1/2}r)\pi^{3/2}}{r a^{3/2}}

<1p,p|\frac{1}{|r-r'|}|1p,q>=\frac{1}{3}\left(-\frac{\pi e^{-ar^2}}{a^2}+\frac{3\pi^{3/2}}{2a^{5/2}}\frac{erf(a^{1/2}r)}{r} \right)

r=0を代入するとクーロン項の行列要素が計算できます。
【２電子積分】
クーロン積分結果に、基底関数との内積をとります。
<1s,1s|\frac{1}{|r-r'|}|1s,1s>=2\pi^{5/2}\frac{1}{a b (a+b)^{1/2}}

<1p,1s|\frac{1}{|r-r'|}|1s,1p>=\frac{\pi^{5/2}}{3a b^2}\left(\frac{b}{(a+b)^{3/2}}+\frac{2}{(a+b)^{1/2}}\right)

<1s,1s|\frac{1}{|r-r'|}|1p,1p>=\frac{\pi^{5/2}}{3}\left(-\frac{1}{a^2(a+b)^{3/2}}+\frac{1}{a b^2 (a+b)^{1/2}}\right)

<1p,1p|\frac{1}{|r-r'|}|1p,1p>=\frac{\pi^{5/2}}{3a^2}\left(-\frac{1}{2(a+b)^{5/2}}+\frac{1}{b^2(a+b)^{1/2}}\right)