Pythonで量子化学計算（適宜更新）
量子化学計算をやりたい
Gaussian とか VASP とか使いたいところですが、個人ではとても買えません。100 万円位かかります。まじで。 （大学生なら大学が契約してる場合もあるので無料で使えたりする）
そんなときに Python でフリーで使えるライブラリがあることを知ったので試してようかなと。
ライブラリについて
適当に検索したところ、Psi4 と PySCFが有名みたいです。それぞれ Ψ と自己無撞着場から来ているであろうネーミングで、いかにも量子化学計算用ライブラリな印象を受けます。
Star 数の推移だけ見ると若干 PySCF のほうが勢いがありそうなので PySCF を使います。Psi4 は日本語記事も充実してるみたいですし（「その他参考」にリンク貼った）。
とりあえずインストール
$ pip install pyscf でいけます（ちなみに psi4 は conda でしかだめらしい）
分子構造の定義
いくつかの方法があるみたいですが、一番手軽な一つを使うことに。
gto（おそらく Gauss 型軌道の略）モジュールを使用して、原子の座標と基底関数を指定できます。
座標の単位はデフォルトでは Å で、ボーア半径への変更も可能です（が使うかね？）。
xyz 座標系を用いていますが、Z-matrix も対応可能らしいです（これは嬉しいかも。molden 使うと Z-matrix 形式になるので）。
from pyscf import gto
mol = gto.M(
atom="""
O 0 0 0;
H 0 1 0;
H 0 0 1;
""",
basis="sto-3g"
)
座標指定についての補足
print(mol._atom) で各種設定が見れるのですが、なぜか Boar 単位ですね。謎。
[('O', [0.0, 0.0, 0.0]), ('H', [0.0, 1.8897261245650618, 0.0]), ('H', [0.0, 0.0, 1.8897261245650618])]
print(mol.atom_coords(unit="ANG"))とすれば Å で見れます。
[[0. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
ちなみに以下のように書く原子に番号をつけて区別することができたり、
mol = gto.M(
atom="""
O 0 0 0;
H1 0 1 0;
H2 0 0 1;
""",
basis="sto-3g"
)
print(mol._atom)
-> [('O', [0.0, 0.0, 0.0]), ('H1', [0.0, 1.8897261245650618, 0.0]), ('H2', [0.0, 0.0, 1.8897261245650618])]
外部ファイルから座標を読むのも簡単にできるみたいです（が上手くいかなかったです）。
mol = gto.M(atom="my_molecule.xyz")
他にもこんな指定方法もあるみたいです。スクリプト組んで色々したいときに便利そうで良きですね。
mol = gto.M(atom=[['O',(0, 0, 0)], ['H',(0, 1, 0)], ['H',(0, 0, 1)]])
基底関数についての補足
錯体の計算とかするときに、金属原子だけは基底関数変えるとかはよくある話ですが、それも対応してるみたいです。
書き方はこんな感じ。
basis = {'O': 'sto-3g', 'H': '6-31g'}
使える基底関数はおそらくここに一覧があります。当然 Pople 系の 6-311G**とかはあった（リンク）ので、試して遊ぶ分には困らなさそうです。
ダニングのなんとかなんとかみたいなのもパット見た感じありました。
エネルギーを計算
とりあえずハートリー・フォック, DFT のほか、MP2 とか CCSD(T)などのポスト HF 系もいけそうです。
とりあえずサンプルコードそのまま書いて実行してみます。実行してみると converged SCF energy = -74.9796525178623 と言われました。計算できてるみたいですね。
scf.xxx() で計算方法の指定、kernel()で実行のようです。
RHF はさすがに Restricted HF でしょうけど、なんで、mfとかdmという変数名にしているのかは謎です。
mf = scf.RHF(mol)
mf.kernel()
dm1 = mf.make_rdm1()
mf.mulliken_pop() とすると、以下のように Muliken 電荷が計算できるみたいです。
ちゃんと酸素が負電荷を、水素が正の電荷を帯びていて良さそうです。
** Mulliken pop **
pop of 0 O 1s 1.99810
pop of 0 O 2s 1.90489
pop of 0 O 2px 2.00000
pop of 0 O 2py 1.20533
pop of 0 O 2pz 1.20533
pop of 1 H 1s 0.53642
pop of 1 H 2s 0.30676
pop of 2 H 1s 0.53642
pop of 2 H 2s 0.30676
** Mulliken atomic charges **
charge of 0O = -0.31364
charge of 1H = 0.15682
charge of 2H = 0.15682
DFT は dft をインポートし、交換相関項の設定すれば、HF と同様の方法で計算可能です
とりあえず皆大好き B3LYP も使えます。
mf = dft.RKS(mol)
mf.xc = "b3lyp"
mf.kernel()
使用可能な汎関数はこちらにあるみたいです（カスタムも可らしい。2 つのライブラリに依存しているらしいので、それぞれある）
libxc ↓
xcfun ↓
Discussion