【物理】ハートリーフォック法で電子軌道計算

物理の勉強がてら、ハートリーフォック法で電子軌道と基底状態のエネルギーを
計算しました。GTOという状態関数を16個のガウス系関数の和で近似する方法です。
s軌道系はχ(r)=exp(-αr*r)、
p軌道系はχ(r)=x*exp(-αr*r)という感じです。
何が大変かというと、交換項の積分計算で
∫∫∫dr1dr2 χ(r1)χ(r2)(1/|r1-r2|)χ*(r2)χ*(r1)
というやつで、真中に1/|r1-r2|が入っているのがいやらしすぎる。
ガウス系だと漸化式を駆使して解析的に解けるらしいですが
結局どうやって導出するのか分からず本に載っている計算式を
使いました。それでもバグチェックが大変でした。
基本的に行列要素を計算して、あとは固有方程式を解くだけです。
LAPACKを使えば一撃です。
 
計算結果です。
 
O原子の電子軌道

Ne原子の電子軌道

線の軌道に±のスピンをもつ電子が2個づつ入ります。
 
参照値と比較しました。微妙に値が違いますが
大体あってます。

     計算結果[a.u] Ref[a.u]

                                                                    • -

[Ne]1s-2s-2p(6) -125.05  -128.35
[O]1s-2s-2p(4)  -76.45  -74.77
[C]1s-2s-2p(2)  -40.73  -37.66
[Be]1s-2s(2)   -16.09  -14.27
[He]1s(2)     -2.85  -2.85

                                                                    • -

この調子で頑張って散乱断面積を計算していけば、原子レベルのレンダリングができるかと思います。

参考:
[1]計算物理学
http://books.google.co.jp/books?id=bFXfyOsd4wsC&pg=PP16&dq=%E8%A8%88%E7%AE%97%E7%89%A9%E7%90%86&as_brr=3#v=onepage&q=&f=false

[2]Modern techniques in computational chemistry
http://books.google.co.jp/books?id=vZlR9sXdGowC&dq=Modern+techniques+in+computational+chemistry&printsec=frontcover&source=bl&ots=dnaNadzHnT&sig=_kP04RJn0FBG4T2jtHDhMGeJQ2c&hl=ja&ei=WNbjSuSFPITq6gOn8dDtAQ&sa=X&oi=book_result&ct=result&resnum=1&ved=0CA0Q6AEwAA#v=onepage&q=&f=false