仕事で化合物のHOMO LUMO GAPを計算したくなったんですが、計算ツールが手元になかったので、pythonで行おうと思いました。
量子化学計算を行えるツールとして、psi4(非経験的分子軌道法のツール)、mopac(半経験的分子軌道法)のパッケージを選びました。
今回はpsi4についてメモを残します。
１ psi4をインストールしました。
conda install -c psi4 psi4
２ 分子情報を入力しました。
geometry.py
import psi4
#分子構造(ベンゼン)の入力
#psi4.geometry()は引数としてcartesian座標形式、またはz座標形式を受け付けます。
#ここでは前者を選択しました。
#ベンゼンはCとHをそれぞれ6つずつ持っています。
#それぞれの原子の3次元位置情報(x,y,z)が記されています。(1行目の0、1はそれぞれチャージと多重度(一重項とか三重項とか))
#cartesian座標はavogadro(無料の分子描画ソフトウェア)にベンゼンの構造式を入れ、「extensions」から「PSI4」を選択し、そこからコピペしてきました。
molecule = psi4.geometry("""
0 1
C 0.16410 -1.37260 -0.00020
C -1.10660 -0.82840 0.00010
C -1.27070 0.54420 0.00000
C -0.16410 1.37260 0.00000
C 1.10660 0.82840 0.00060
C 1.27070 -0.54420 -0.00070
H 0.29230 -2.44490 0.00440
H -1.97120 -1.47560 0.00000
H -2.26350 0.96930 -0.00050
H -0.29230 2.44490 -0.00080
H 1.97120 1.47560 -0.00010
H 2.26350 -0.96930 0.00020
""")
３ 構造最適化計算します。((特に構造が複雑な場合に)conformerをいくつか作ってそちらも計算し、最もエネルギーが低い結果を採用すべし。)
optimize.py
#計算のメソッドにB6LYPを、基底関数に6-31G(d)を選択しました。
#psi4.optimize()を使いました。
energy, wfn = psi4.optimize("B6LYP/6-31G(d)", molecule, return_wfn = True)
#普通のPCで普通に実行すると恐らく数時間かかります(ベンゼンなら数分？)
#ハイスペックのPCで以下のメモリとスレッドの設定をした方がいいと思います。
#例 (会社のPC使いました。こんなに必要かは知りません笑)
psi4.set_memory(1000GB)
psi4.set_num_threads(80)
４ wfnからHOMO LUMOの情報を取り出します。
homolumogap.py
#HOMO軌道とLUMO軌道のエネルギー値(単位はhartree)を取り出し、27.2114をかけてeVの単位に変換しました。
HOMO = wfn.epsilon_a_subset('AO', 'ALL').np[wfn.nalpha()-1]*27.2114
LUMO = wfn.epsilon_a_subset('AO', 'ALL').np[wfn.nalpha()]*27.2114
#gapの計算です。
HLG = LUMO - HOMO
#表示
print(HOMO, LUMO, HLG)
安定化した構造情報を保存して、avogadro等で描画させることもできます。
以上