【計算化学】【Python】対話型のエネルギー計算プログラムを作ってみた
以前、拙著エントリにてpsi4のことについてほんのちょっとだけ触れたのだが、その時は結局itertoolsというライブラリを紹介して終わりになった。
今回こそはpsi4をやろう。せっかくだから、難しいコマンドは抜きにして、計算条件を入力すればあとは結果を返してくれるようなものが欲しいので、それを作ってみることにした。ただし、今回扱うのは一点計算のみであり、構造最適化はまた別の機会に回すことにする。
psi4とは?
psi4と書いて「サイフォー」と読む。Python上で動かせる、フリーのab initio量子化学計算プログラムとして位置づけられている。
量子化学計算を行う場合、計算用ソフトウェアとしてのfirst choiceはやはりGaussianということになろう。機能性、計算速度、業界での知名度など、あらゆる面で優秀な存在だ（と私は思っている）。しかし、そんな甘い蜜をタダで吸えるわけもなく、それなりの投資が必要になる。加えて、どれだけGaussianが優秀と言えども、計算リソースが伴わなければ無用の長物だ。やはり究極的には、フリーのソフトウェアで計算して、浮いた予算をリソースに充当したくなるというものだ。
psi4は、それを可能にしてくれるかもしれない。なんといっても無料である。言うまでもないがPython本体も無料だ。これでGaussianの購入資金をまるまる節約できる。
実際、Gaussianは高機能性が却って災いして、ポテンシャルを十分に引き出せないというケースは多いのではないだろうか。QM/MMとかONIOMとかNBOとか言われても、私のような非専門家からしてみれば「？」の嵐だ。
そんな複雑な計算はいらん、単にopt.とfreq.がしたいだけなんだという時に、psi4は良い選択肢になるだろう。というか、psi4が他の有償ソフトウェアと比べて限定的な機能しか有していないのかと問われると、そんなこともなさそうだ。詳しく調査してはいないがIRC計算もできるようだし、非共有結合性相互作用の評価に優れているとも聞く。
そしてこれは、psi4がどうのというより私個人の癖なのだが、GUIではなくCUIでPCに指令を出して結果を得るというのはなんとも格好がつく。ハッタリがきいている。だから、やる。そういう不純な動機でもいいじゃないか。いざ計算。
インスト―ル
以下のコマンドでインストール可能。このコマンドではpythonのバージョンが3.6である場合を想定してこのような記述になっているが、ここはお使いの環境に合わせて適宜変更していただきたい。
conda install psi4 psi4-rt python=3.6 -c psi4
エネルギー計算をしてみる
途中の過程を全て放り出しているような気もするが、以下のようなコードを書いてみた。計算に使用するスレッド数、メモリ数などを対話的に入力して行けば、エネルギーが計算される。ひとまずコード全体を俯瞰した後に、ひとつひとつ分割して考えていく。私の拙い英語は黙殺してください。
import psi4 print("Enter the number of Threads to be used.") nt = input() psi4.set_num_threads(nthread = int(nt)) print("Enter the amount memory to be used.") mem = input() psi4.set_memory(mem) print("Enter the name of the output file.") filename = input() psi4.set_output_file("{}.log".format(filename)) print("Enter the level of theory (e.g. hf, b3lyp) of this calculation.") lot = input() print("Enter the basis set (e.g. 6-31G) of this calculation.") basis = input() level = lot + "/" + basis print("Enter the charge number of your molecule.") charge = int(input()) print("Enter the spin multiplicity of your molecule.") spm = int(input()) print("Paste the input of your molecule.") car = (input()) molec = psi4.geometry(''' {} {} #charge and spin multiplicity {} '''.format(charge, spm, car)) e = psi4.energy(level, molecule = molec) print("Calculated Energy: {} a.u.".format(e))
では、上から順番に見ていこう。
import psi4 print("Enter the number of Threads to be used.") nt = input() psi4.set_num_threads(nthread = int(nt)) #ここでスレッド数を指定する
まずもって行うべきはimportだ。下の三行はまぁ見ればわかるのだけど、Python側が「使うスレッド数教えてね」と尋ねてくるので、こちらはそれに対して適当な数字を返す。が、inputメソッドで取得した値はstr扱いになってしまうので、無理やりintに変換している。
次に、使用メモリ量を指定する。
print("Enter the amount memory to be used.") mem = input() psi4.set_memory(mem) #ここでメモリを指定する
スレッド数の場合と同じ構造なのだが、今回はinputメソッドの値(mem)をint化していない。というのも、メモリについては"4GB"のように、単位付きで指定しなければならないからだ。
次は計算結果のoutput fileの名前入力だ。
print("Enter the name of the output file.") filename = input() psi4.set_output_file("{}.log".format(filename))
formatメソッドは最近知ったのだが、コード中においてどの部分が変数なのかわかりやすくなるので、覚える価値ありだろう。そのうち短めのエントリを書く予定。
次は計算の核心部分なので2つまとめて紹介する。汎関数と基底関数を指定する。つい汎関数という言葉を使ってしまうのは私が日ごろDFTをしているからなのだが、Hartree-Fockも可能。
print("Enter the level of theory (e.g. hf, b3lyp) of this calculation.") lot = input() print("Enter the basis set (e.g. 6-31G) of this calculation.") basis = input() level = lot + "/" + basis
例えば汎関数がB3LYPで基底関数が6-31Gなら、末尾の部分で"B3LYP/6-31G"というLevel of Theoryが出来上がることになる。これをlevelという変数に代入している。なお、psi4で使用可能な汎関数および基底関数については以下参照。繰り返しになるが、psi4はDFT以外にもいろいろ計算可能。単に私がDFTしか知らないだけだ。
www.psicode.org
だいぶ道のりが長くなってきたが、次に電荷数、スピン多重度、およびinputする分子構造を指定する。formatメソッドを使っている関係上、まとめて解説した方が良いだろう。
print("Enter the charge number of your molecule.") charge = int(input()) #電荷数の指定 print("Enter the spin multiplicity of your molecule.") spm = int(input()) #スピン多重度の指定 print("Paste the input of your molecule.") car = (input()) #分子構造の指定 Z-matrixかXYZ形式の座標をペーストする molec = psi4.geometry(''' {} {} #charge and spin multiplicity {} '''.format(charge, spm, car))
電荷数、スピン多重度、および座標の3つをまとめてinput分子を指定する。これらの値は、三重引用符で囲んでいることに注意。
ここまで入力したら、あとはPCが仕事をする番だ。
e = psi4.energy(level, molecule = molec) print("Calculated Energy: {} a.u.".format(e))
画面上に、計算で求められたエネルギーが出力される。単位はa.u.であることに注意。
画面上にエネルギーが表示されて終わりというのはなんともショボくさく見えるかもしれないが、実際にはカレントディレクトリにoutput fileが出力されているので、詳細な計算結果についてはそちらを確認すればいい。
