ASEとQuantum Espressoを使って、H2Oなどの３原子分子を構造を最適化して、安定構造を見つける。

はじめに

　Jupyter上で、PythonライブラリであるASE(Atomic Simulation Environment)を介して、第一原理計算ソフトQuantum Espressoで、第一原理計算を実施してみます。
　本記事は、書籍『密度半関数理論入門(D.S.ショール，J.A.ステッケル 著、佐々木泰造，末原茂 共訳』(以下、教科書)の例題や課題を参考に、プログラムを作って計算します。

　なお、固体DFTの専門家ではありませんので、計算条件、結果等はご自身でご確認/判断ください。

今回は、教科書の3章の練習問題5、
『H2Oなどの３原子分子を構造最適化して、安定構造を見つける』
に取り組んでみたいと思います。

前準備

計算に必要になるライブラリや、擬ポテンシャルを用意します。

*.jpynb

# import libs
import os
import subprocess
import numpy as np
import matplotlib.pyplot as plt

# import ase
from ase import Atoms
from ase.optimize import BFGS
from ase.visualize import view
from ase.io.trajectory import Trajectory
from ase.calculators.espresso import Espresso

#教科書の例にならい、pw91の擬ポテンシャルを使います。
#!wget http://nninc.cnf.cornell.edu/psp_files/O.pw91-van_ak.UPF
#!wget http://nninc.cnf.cornell.edu/psp_files/H.pw91-van_ak.UPF
o_pp = "O.pw91-van_ak.UPF"
h_pp = "H.pw91-van_ak.UPF"

関数定義

繰り返し使う部分を関数定義しておきます。

*.jpynb

def calc_dft(cutoff=300.0,Mk=3):
    # QEでDFT計算をする. cutoff:[eV] 
    eV_to_Ry =  1/13.605698#[eV/Ry] 
    pseudopotentials = {'O':o_pp,'H':h_pp}
    #実行環境に合わせて、必要に応じて下記を変えてください。
    cmd = 'pw.x -in espresso.pwi > espresso.pwo'

    input_data = {
        'control':{'pseudo_dir':'./'},
        'system': {
            'ecutwfc': cutoff*eV_to_Ry,      #[eV]->[Ry]
            'ecutrho': cutoff*4*eV_to_Ry},   #ecutwfc * 4
        'disk_io': 'low'}

    calc = Espresso(command=cmd,
                    pseudopotentials=pseudopotentials,
                    kpts=(Mk, Mk, Mk),
                    tprnfor=True,#ASEを使った構造最適化時に必要
                    input_data=input_data)
    return calc

def make_model(molc,L,a,b=0):
    # 座標を指定して分子モデルを作成する
    # 原点に１つ、対称位置に原子を２つ配置して作成
    model = Atoms(molc, 
                positions=[(0,0,0),(a,b,0),(-a,b,0)],
                cell=[L,L,L],
                pbc=[1,1,1])
    return model

def calc_length(a,b):
    #結合長を計算する
    return np.linalg.norm((a-b))

def calc_angle(center,a0,b0):
    #結合角を計算する。centerに中心の原子の座標を指定。
    a = a0 - center
    b = b0 - center
    cos = np.inner(a,b)/(np.linalg.norm(a) * np.linalg.norm(b))
    deg = np.rad2deg(np.arccos(np.clip(cos, -1.0, 1.0)))
    return deg

水分子を最適化してみる

まずは初期の分子モデルを作成して、moleculeに格納します。格納された構造や中身も確認します。

*.ipynb

#分子モデルを作成
molecule = make_model(molec='OH2', L=10.0, a=1.2, b=0.5)

#原子と座標を確認
for i in range(len(molecule)):
    print(f'elem : {molecule.symbols[i]}, position : {molecule.positions[i]}')
###出力結果###  
  #原点がO(酸素)なのが確認できます
  elem : O, position :[0. 0. 0.]
  elem : H, position :[1.2 0.5 0. ]
  elem : H, position :[-1.2  0.5  0. ]

#初期構造の結合長、結合角も確認
pos=molecule.get_positions()
print(f'Length : {calc_length(pos[0],pos[1])}, {calc_length(pos[0],pos[2])} ')
print(f'Angle  : {calc_angle(pos[0],pos[1],pos[2])}' )
###出力結果###
  Length : 1.3, 1.3 
  Angle  : 134.76027010391914

#必要に応じて可視化して確認
#view(molecule ,viewer='ngl')

最適化

では、上記moleculeとして作ったH2Oを、最適化してみます。

*.ipynb

calc = calc_dft(cutoff=300,Mk=3) #
molecule.set_calculator(calc)
dyn = BFGS(molecule, trajectory='qn.traj',logfile='qn.log')
dyn.run(fmax=0.01) #fmax eV/A, max force for optimizations

しばらくすると、計算が終わります。（おそらくTrueの文字が返ってくるかと）

結果の処理

初期構造を設定したmoleculeの構造が最適化されたものに変わっています。確認してみましょう。

*.ipynb

# 原子の座標をposに格納する
pos=molecule.get_positions()
# 分子のエネルギーをenergyに格納する
energy = molecule.get_potential_energy()
print(f'Length : {calc_length(pos[0],pos[1])}, {calc_length(pos[0],pos[2])} ')
print(f'Angle  : {calc_angle(pos[0],pos[1],pos[2])}' )
print(f'Energy : {energy}')

###出力結果###
  Length : 0.972895505008442, 0.9728955050084419 
  Angle  : 104.61973639267624
  Energy : -468.72102581446654

#必要に応じて可視化して確認
#view(molecule ,viewer='ngl')

最適化によって、

• 結合長は1.3 →　0.97 Angstrom
• 結合角は134 → 104.6°

と最適化されていることが確認できます。

qn.logにエネルギーの変化が保存されています。確認すると、エネルギーが順調に下がっているのが分かります。
qn.trajに最適化の各ステップでの構造が保存されています。

初期構造を変えて最適化

今回使った最適化手法は、ASEに実装されているBFGS法です。教科書にもありますが、これらの最適化（QE本体の最適化法も）は基本、初期構造の近傍の局所安定構造を求めるにすぎません。そこで、初期構造を変えて計算し、結果を確認してみましょう。
モデル作成部の
　　molc = make_model(molec='OH2', L=10.0, a=1.2, b=0.5)
のa,bを適宜変更して最適化してみます。

case	初期 O-H結合	初期 H-O-H角	最適化 O-H結合	最適化 H-O-H角	Energy
1	1.3	134.8	0.973	104.62	-468.721
2	0.95	179.9	0.973	104.63	-468.721
3	1.3	180.0	0.944	180.00	-467.475

ケース１、２はほぼ同じ構造になってます。３は、b＝0として作った初期構造ですが、おそらく最適化時にデフォルトで対称性を維持する設定になっていることで、結合角が変わらない最適化がなされていると思われます。
Energyを見ると、H-O-Hが曲がっている１,2が安定であることがわかります。

ケース２の最適化の前後を図示しておきます

実験値との比較

念のため実験値と比較してみましょう。CCCBDBのH2Oのデータを参照します。

• O-H結合長　：0.9578 Å
• H-O-H結合角：104.478 degree

ケース１，２では角度はよく再現できていますが、結合長が若干長いですね。

---

ASEとQuantum Espresso関連の投稿記事リンク

ASEとQuantum Espressoを使って、

• 手動でCuの単位格子最適化して、平衡体積弾性率を見積る
• Ptが単純立方構造、fcc、hpcのいずれの構造をとるかを調べる
• 六方晶Hfの格子定数aとcを最適化する
• H2Oなどの３原子分子を構造を最適化して、安定構造を見つける

他

• ASEでモデルの図を保存する,描写する,モデルを出力する(png,cif,traj,pov)
• Quantum Espressoを色々なライブラリでコンパイルして、自分のPCに最速な状態にする