ASEとGPAWでDFT計算を行う

10942 ワード
JupyterLab Python ASE GPAW DFT Python テキストリンク

ASEとGPAWを用いてDFT計算を行ってみます．簡単な例として水素分子の原子化エネルギーを求めましょう．

  • ASEは原子シミュレーション用のPythonモジュールです．多くの機能が実装されています．Calculatorクラスをインターフェースとして各種コードを呼び出すことができます．
  • GPAWはPAW法とASEに基づいたDFT計算コードです．C言語とPythonにより記述されています．広く使用されています．

以下の例はJupyterLabを用いて実行しました．

インストール

Anacondaで構築したPython環境に，pipでASEとGPAWをインストールします．

$ pip install ase --user
$ pip install gpaw --user

GPAWのPAWのデータセットもダウンロードしてください．

$ gpaw install-data $HOME/gpaw-data

インストールは以上で終わりです．
この段階で私の環境は次のようになっています．

$ python --version
Python 3.7.6
$ pip list installed | grep -e ase -e gpaw
ase                                3.19.1
gpaw                               20.1.0

GPAWの計算条件

次の条件でGPAWクラスをインスタンス化します．

  • カットオフエネルギー: 300eV
  • 汎関数: PBE
from gpaw import GPAW, PW

gpaw = GPAW(mode=PW(300), xc='PBE')

水素分子のモデルを作成

ase.build モジュールには小分子のプリセットが含まれます．
molecule関数に分子式を渡すだけで，3次元の分子モデルであるAtomsオブジェクトを作成できます．
今回は孤立系の計算を行いたいので，作成したAtomsオブジェクトのcenter メソッドを呼び出して，モデルに3.0Aの真空層を追加しています．

from ase.build import molecule

atoms_h2 = molecule('H2')
atoms_h2.center(vacuum=3.)

view関数で可視化することができるので，分子の構造を確認しておきましょう．

from ase.visualize import view

view(atoms_h2)

計算を実行する

AtomsオブジェクトにGPAWのCalculatorを関連付けて，get_potential_energy関数を呼び出すと，GPAWによる計算が実行されます．

atoms_h2.set_calculator(gpaw)
e_h2 = atoms_h2.get_potential_energy() # takes several seconds

水素原子の計算

水素原子についても同様に計算を行ってやります．

atoms_h = molecule('H')
atoms_h.center(vacuum=3.)
atoms_h.set_calculator(GPAW(mode=PW(300), xc='PBE', hund=True))
e_h = atoms_h.get_potential_energy()

結果の表示

最後に結果をまとめてみましょう．
H2の原子化エネルギー$\Delta E$は次式で得られます．

$\Delta E = 2 E_{\rm H} - E_{\rm H_2}$

同時に分子構造もNotebook上に表示してみます．
matplotlibのラッパであるase.visualize.plotモジュールを使用します．

import matplotlib.pyplot as plt
from ase.visualize.plot import plot_atoms

delta_e = 2 * e_h - e_h2
print(f'Atomization Energy: {delta_e: 5.2f} eV')

fig, ax = plt.subplots(1, 2, figsize=(8,6))

title = f'Calculation Result for {atoms_h2.get_chemical_formula()}\n' + \
        f'Total Energy : {atoms_h2.get_potential_energy(): 5.2f} eV'
ax[0].set_title(title)
ax[0].axis('off')
plot_atoms(atoms_h2, ax=ax[0], rotation='90x')

title = f'Calculation Result for {atoms_h.get_chemical_formula()}\n' + \
        f'Total Energy : {atoms_h.get_potential_energy(): 5.2f} eV'
ax[1].set_title(title)
ax[1].axis('off')
plot_atoms(atoms_h, ax=ax[1], rotation='90x')

plt.show()

参考