pythonの量子化学計算用モジュールで遊んでみた。最初は、量子化学計算みたいなコストのかかる計算をスクリプト言語で書いてどうすんの？みたいに思ったが、速くしなきゃならん部分だけCで書いてあとはPythonでもいいのかも？いや逆にトータルの生産性はこっちのほうがいいのかもしれないゾと思ったり。
あとは、モパクーもガウシーもゲムスーもそもそもインプットからしてややこしいから、よりシンプルな記述で計算できるとベター。だからOOPな感じで書けるのは有用かもん。
GUIでインプットのややこしさをラップすればいいじゃん的ナ発想もあるんだろうが、複雑なオプション与えたいときに、GUIの操作ほうが余計複雑になっちゃったりとか、結局自分でインプットいじらなあかんていうことが良くあるしな。なにより、バッチとか、スクリプトに組み込むときメンドイ。
特にUserGuideの例にあるように、同じ分子をアルゴリズム変えたり、基底関数変えて計算したいとかいうときには、スクリプトのほうが便利だ。
ユーザーガイドには水素分子の計算が例として挙げられているが、ちょっとだけいじった、水分子のHF計算をのせておく
from PyQuante.hartree_fock import rhf
from PyQuante.Molecule import Molecule
h2o = Molecule('H2O',
atomlist = [(8,(0,0,0)),
(1,(0.959,0,0)),
(1,(-.230,0.930,0))],
units = 'Angstrom')
en,orbe,orbs = rhf(h2o,verbose=True)
モジュールのインポートをしたら、moleculeオブジェクトで水分子を作ったら、rhfでh2oを計算。これだけ。
お手軽!!
$ python test.py
Nbf = 25
Nclosed = 5
Optimization of HF orbitals
0 -68.2710764842
1 -71.4107153837
2 -73.8915175997
3 -75.0373719156
4 -75.6459826198
5 -75.8850728221
6 -75.9751560983
7 -76.0064588121
8 -76.0173702077
9 -76.0210974045
10 -76.0223774136
11 -76.0228145355
12 -76.0229641837
13 -76.023015328
結果はこんな感じで。verbose=Trueのオプションつけたからずらずら出てくるが、普通はエネルギーを出力させるくらいでいいでしょう。
計算時間は、C3マシンで実行したので3分程度かかったが、今のスペックならほとんど一瞬で終わるのではないかな。
Chemruby,PerlMol,CDKなどで、ある程度まともな3次元構造立ち上げられるようになったら、
構造立ち上げ -> MMで構造最適化 -> 半経験的MO構造最適化 -> いい基底関数とMP2レベルでシングルポイントの電子状態計算
という、トラディッショナルな構造最適ルーチンをスクリプトで実行したくなるだろうなぁ。
Corinaレベルの構造立ち上げライブラリでオープンなものってあるのだろうか？職場だったら、corina用モジュール書けばいいのだろうが。MMはTinker使えるような気がするが使ったことないからわからん。mopacもラッパー書けばいいじゃんって言われそうな気もするが。