量子化学計算のプログラムを自作してみた

普段、量子化学の計算を行う際は、GaussianやGAMESS他のソフトが用いられることが殆どで、自分でコードを書くのは計算手法の研究をされている一部の方に限られるのではないかと思います。

近年では、B3LYP汎関数の使用を始めとする、密度汎関数理論(Density Functional Theory; DFT)が用いられるケースが殆どですが、それに先立って確立された手法として、Hartree Fock法があります。 Hartree Fock (HF)法は、計算結果のエネルギーの見積もりの不正確さ*1から、今ではこの手法が単体で使用される場面は多くはありません。しかし、近似の仕方が明確なため*2、量子化学の基礎として多くの方が一度は勉強される理論かと思います。

私も、数年前に教科書でHFの理論を勉強してから、一度は自分でコードを書いておきたいという気持ちがずっとあったのですが、毎回途中やめになっており、このたび漸くそのminimalな実装をC++で書ききりました。実際に書いてみると、式を読んで分かった気になっていたものの、実際の理解は曖昧だったりといったことがよくわかりました。

コードは、下記のリポジトリにおいていますので、気になる方、これから同じことを勉強しようと言う方はどうぞ。

github.com

感想としては、結構大変でした。Gaussianなどの開発チームだったり、もしくは日々新たな計算手法を開発されている方たちが神に見えてきます。以下に気づき、感想、苦労した点を列挙します。

1. Gauss型の短縮基底関数の正規化

これは、自分で書いてみるまで気がつかなかったことです。（といっても書き始めたのが２年くらいは前なので、その頃に気がついたのですが。。。）例えば、STO-3Gという最も基本的な基底関数系で、水素の1S軌道は、軌道指数(Orbital exponent, よくαと書かれます)と計数を強調して、

$H_{1s} (r) = 0.44635 g(-0.168856) + 0.535328 g(-0.623913) + 0.154329 g(-3.42525)$

という感じでよく表されます。ただ、これを $g(\alpha) = \exp(-\alpha r^{2})$ と考えて素直に計算すると、一番簡単なはずの重なり積分から数字が合いません。

実際には、（後で気づけば当たり前のことではあるのですが...）それぞれのプリミティブガウス関数に、それぞれを全空間で積分した際の値が１になるように、規格化因子がかけられています。また、さらに、それらを重ね合わせた関数（短縮ガウス関数）にも、規格化因子が掛けられます。したがって、より正しい表記は、

$H_{1s} (r) = \frac{1}{0.999999 } (0.44635 \frac{\exp(-0.168856 r^{2})}{0.187736} + 0.535328 \frac{\exp(-0.623913 r^{2})}{0.500326} + 0.154329\frac{\exp(-3.42525 r^{2})}{1.79444} )$

となります。これらの規格化因子は、プリミティブガウス関数のそれぞれについてまず求め、その後、短縮ガウス関数について求める必要があります。なお、p軌道以上の角運動量を持つ軌道には、 $x^{l}y^{m}z^{n}$ という因子がかかるので、これも考慮に入れる必要があります。

2.積分の実装

積分は、1966年に竹田、藤永、大旗の３名が書いた論文を参考にしました。これの実装がとにかく難しかったです。慣れれば、重なり積分と、運動エネルギー積分はまだなんとかなります。難しいのは、核引力積分(Nuclear Attraction Integral)と電子間反発（４中心）積分の部分です。

これらは、とにかく展開方法もなかなか理解できませんし、非常に入れ子の深いループになってしまいます。

現状は上記の論文をそのまま参考にしており、STO-3G基底での水素分子やHeH+イオン程度は一瞬で計算は終わりますが、試しにp軌道を分極関数として追加してみると、一気に遅くなりました。

近年のGaussianを始めとする量子化学計算プログラムは、PRISMアルゴリズムなど、積分ごとに20種類くらいのパスを使い分けるような非常に複雑な実装*3がなされているそうですが、その重要性がよーく分かりました。今後は、できればDKRの積分方法や、Obara-Saikaの展開公式なども実装してみたいと思います（できれば、です。２電子積分の展開など、まだ全然理解できていません。）

3. 何を使って行列演算を行うか

王道はBLASなどのライブラリを使うことなのだろうと思いますが、恥ずかしながら私は他のプログラムからリンクさせる為にコンパイルした程度の経験しかありませんでした。

以前にnumpyを使って今回と同じようなものを書こうとしたことがありました*4ので、今回は比較的それと似たような感覚で使えそうなEigenというテンプレートライブラリを使用しました。

おそらくBLASなどのような高度なチューニングは難しいでしょうが、今回のような目的には非常に使い勝手が良かったと思います。

今後

現状はHartree FockのRHF計算のみ実装しましたので、今後はUHF計算をまずは実装したいと思います。あとは、スパコンが使える今のうちに、少しくらいMPIチューニングなどで遊んでみるのも良いかもしれません。その次は、DFTか微分計算ではありますが、おそらくUHFを書いてみるころには学生生活が終了している気がします...苦笑

2017.07.15 追記

UHFの計算はSzaboの本を見ながらすぐにかけてしまいました。今はDFTを実装するべくいろいろ調べていますが、交換・相関ポテンシャルの空間積分のためのグリッド生成がなかなか難しいです。またここにメモをする予定です。

追記その２ Qiitaにメモしました。まだ積分は書いていなくて、ハードコーディングしてます。 qiita.com