【化学でPython】GPAW：実空間グリッド法を用いた本格的な第一原理計算を手軽に実行する
はじめに
この「化学でPython」シリーズでは、化学の分野で有用な Python ライブラリを紹介しています。
今回紹介するのは、GPAW です。
GPAW とは？
GPAW は、実空間グリッド（real-space grid）法 で密度汎関数理論（DFT）計算を実行するためのライブラリです。
Python の原子・分子シミュレーション環境である ASE (Atomic Simulation Environment) と統合されているので、Python スクリプトの中で直感的に計算を実行できます。
インストール
pip で簡単にインストールできます。
pip install gpaw
基本的な使い方
まずは、最も基本的な使い方として、水素分子のエネルギー計算 を行ってみましょう。
GPAW は ASE の Calculator （計算エンジン）として使えます。ASE で原子（Atoms）を作り、それに GPAW 計算エンジン（GPAW）をセットして計算を実行します。
from ase import Atoms
from gpaw import GPAW
# 1. 分子の定義 (水素分子 H2, 結合長 0.74 Å)
# pbc=False は周期境界条件なし（孤立分子）を意味します
# 孤立分子の場合でも、シミュレーションセル（真空領域）を定義する必要があります。
# ここでは、分子から 2.5 Å の真空領域を確保して、一辺が 5 Å 程度のセルを定義します。
h2 = Atoms('H2', positions=[[0, 0, 0], [0, 0, 0.74]], pbc=False)
h2.center(vacuum=2.5)
# 2. 計算手法の設定
# mode='lcao': 原子軌道線形結合法（高速・定性的）
# basis='dzp': Double Zeta Polarized 基底
# xc='LDA': 局所密度近似（交換相関汎関数）
# txt='h2.txt': 計算結果を保存するファイル名
h2.calc = GPAW(mode='lcao', basis='dzp', xc='LDA', txt='h2.txt')
# 3. ポテンシャルエネルギーの計算
energy = h2.get_potential_energy()
print(f"H2 molecule energy: {energy:.3f} eV")
H2 molecule energy: -6.412 eV
このように、ASE のインターフェースを介して、簡単に DFT 計算を実行できます。
実践例: 一酸化炭素(CO)の構造最適化
実践的な例として、「一酸化炭素（CO）分子の結合長を第一原理計算で決定する」 タスクをやってみましょう。
1. 初期構造の作成と計算設定
まず、適当な結合長（ここでは 1.0 Å）で CO 分子を作成します。
今回は精度を重視して、GPAW の特徴である 実空間グリッドモード (mode='fd') を使ってみます。
from ase import Atoms
from ase.optimize import BFGS
from gpaw import GPAW
# 1. CO分子の初期構造 (結合長 1.0 Å)
co = Atoms('CO', positions=[[0, 0, 0], [0, 0, 1.0]], pbc=False)
co.center(vacuum=2.5)
# 2. 計算機の設定 (実空間グリッド法)
# h=0.2: グリッド間隔 (Å)
# xc='PBE': 一般的な GGA 汎関数
calc = GPAW(mode='fd', h=0.2, xc='PBE', txt='co.txt')
co.calc = calc
print(f"Initial bond length: {co.get_distance(0, 1):.3f} Å")
Initial bond length: 1.000 Å
2. 構造最適化の実行
ASE の BFGS アルゴリズムを使って、力がゼロになる（安定な）構造まで原子を動かします。
# 3. 構造最適化
# trajectory='co_opt.traj': 最適化の履歴を保存
opt = BFGS(co, trajectory='co_opt.traj')
# fmax=0.05: 原子にかかる力が 0.05 eV/Å 以下になるまで計算
opt.run(fmax=0.05)
# 4. 結果の確認
final_distance = co.get_distance(0, 1)
print(f"Optimized bond length: {final_distance:.3f} Å")
# 実験値との比較
exp_value = 1.128
diff = final_distance - exp_value
print(f"Difference from Exp: {diff:.3f} Å")
Step Time Energy fmax
BFGS: 0 23:59:45 -13.602541 25.980744
BFGS: 1 23:59:49 -12.781053 12.872487
BFGS: 2 23:59:54 -14.323885 9.830642
BFGS: 3 23:59:57 -6.823306 86.294735
BFGS: 4 00:00:01 -14.682890 7.872000
BFGS: 5 00:00:05 -14.884819 5.924364
BFGS: 6 00:00:09 -15.042435 3.426015
BFGS: 7 00:00:12 -15.080746 0.824305
BFGS: 8 00:00:14 -15.083460 0.224957
BFGS: 9 00:00:17 -15.084544 0.320481
BFGS: 10 00:00:21 -15.089016 0.984184
BFGS: 11 00:00:24 -15.094244 1.136725
BFGS: 12 00:00:27 -15.100616 0.847572
BFGS: 13 00:00:29 -15.103207 0.244326
BFGS: 14 00:00:33 -15.103371 0.043250
Optimized bond length: 1.133 Å
Difference from Exp: 0.005 Å
計算された結合長は約 1.133 Å となりました。実験値（1.128 Å）との誤差はわずか 0.005 Å 程度です。
PBE 汎関数は結合長を少し過大評価する傾向がありますが、非常に良い精度で構造を予測できていることがわかります。
3. 最適化過程の可視化
最適化の履歴 (co_opt.traj) を読み込んで、結合長とエネルギーの変化をグラフにしてみましょう。
import matplotlib.pyplot as plt
from ase.io import read
import numpy as np
# 履歴ファイルの読み込み
traj = read('co_opt.traj', index=':')
# 結合長とエネルギーの抽出
distances = [atoms.get_distance(0, 1) for atoms in traj]
energies = [atoms.get_potential_energy() for atoms in traj]
steps = range(len(traj))
# グラフの描画
fig, ax1 = plt.subplots()
color = 'tab:red'
ax1.set_xlabel('Step')
ax1.set_ylabel('Bond Length (Å)', color=color)
ax1.plot(steps, distances, marker='o', color=color)
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx() # 2軸目（右側）を作成
color = 'tab:blue'
ax2.set_ylabel('Potential Energy (eV)', color=color)
ax2.plot(steps, energies, marker='x', color=color, linestyle='--')
ax2.tick_params(axis='y', labelcolor=color)
plt.title('CO Geometry Optimization')
fig.tight_layout()
plt.show()
出力結果
ステップが進むにつれて結合長が 1.133 Å に収束し、ポテンシャルエネルギーが低下していく様子が確認できます。
まとめ
今回は GPAW を紹介しました。
- Point 1: 実空間グリッド法を用いた DFT コードで、基底関数の選択に悩まされにくいです。
- Point 2: ASE と完全に統合されており、Python で柔軟に計算フローを構築できます。
- Point 3: 分子だけでなく、金属表面や半導体などの固体物理の計算にも活用できます。
GPAW は「手軽に第一原理計算をやってみたい」という場合に最適なツールの一つ。
ぜひ試してみてください。
Discussion