はじめに

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

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

今回は、教科書3章,練習6で紹介されている、
『六方晶Hfの格子定数aとcを最適化する』
問題に取り組んでみます。

前準備

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

*.ipynb

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

from ase import Atoms
from ase.build import bulk
from ase.optimize import BFGS, GPMin
from ase.constraints import UnitCellFilter
from ase.visualize import view
from ase.io.trajectory import Trajectory
from ase.calculators.espresso import Espresso

#今回はPAW、PBEを使用
#!wget http://pseudopotentials.quantum-espresso.org/upf_files/Hf.pbe-spn-kjpaw_psl.1.0.0.UPF
hf_pp = "Hf.pbe-spn-kjpaw_psl.1.0.0.UPF"

関数定義

*.ipynb

def make_crystal_lattice(type=None, a=None,c=None):
    #結晶構造のモデルを作る関数
    model = bulk(name='Hf', crystalstructure=type, a=a, c=c,orthorhombic=False, cubic=False, basis=None)
    return model

def calc_dft(Mk):
　　#QEで計算を行う関数
    pseudopotentials = {'Hf':hf_pp}
    cmd = 'pw.x -in espresso.pwi > espresso.pwo'
    input_data = {
         'control':{'pseudo_dir':'./'},
         'system':{
            'occupations' : 'smearing',
            'smearing' : 'mp',
            'degauss' : 0.001, 
             'ecutwfc': 46.0,  #Suggested minimum value by the potential
             'ecutrho': 264.0}, #Suggested minimum value by the potential
         'disk_io':'low'}
    calc = Espresso(command=cmd,
                    pseudopotentials=pseudopotentials,
                    kpts=(Mk, Mk, Mk),
                    tprnfor=True,#ASEを使った構造最適化時に必要
                    tstress=True,
                    input_data=input_data)
    return calc

六方晶Hfの構造を作る

*.ipynb

model = make_crystal_lattice(type='hcp', a=3.0, c=3.5)

#できたモデルの格子ベクトルをチェック
print(model.cell)
###出力結果###
Cell([[3.0, 0.0, 0.0], [-1.5, 2.598076211353316, 0.0], [0.0, 0.0, 3.5]])

#必要に応じて可視化
#view(model,viewer='ngl')

ここで、六方晶の格子ベクトルは下記で、

\boldsymbol{a} = (a ,0 ,0)　　\boldsymbol{b} = (-\frac{1}{2}a ,\frac{\sqrt{3}}{2}a ,0)　　\boldsymbol{c} = (0,0,c)

で表されます。上記で出力された格子ベクトルと比較すると、作ったモデルが、指示通りとa=3.0, c=3.5で作成されていることが確認できます。

格子定数最適化する

では、上で作ったモデルを最適化してみます。
今回は、aseの格子定数最適化機能、 UnitCellFilterを用います。

*.ipynb

#最適化の実行
calc = calc_dft(Mk=12)
model.set_calculator(calc)
opt_unitcell = UnitCellFilter(model)
opt = BFGS(opt_unitcell, trajectory='qn.traj',logfile='qn.log')
opt.run(fmax=0.01)
#少し時間がかかります。

#最適化された構造の格子定数cとaを得る
opt_a = model.cell[0][0]
opt_c = model.cell[2][2]
print(f'a={opt_a:.2f}, c={opt_c:.2f}, c/a={opt_c/opt_a:.2f}')

### 実行結果 ###
a=3.18, c=5.04, c/a=1.58

c/aは1.58という値が得られました。

実験値[Ref]では、a = 0.3188nm, c = 0.5042 nm, c/a = 1.58
とのことですので、よく一致した結果が得られています。
[Ref] Markin V.Y. et al, Inorg. Mater. v.19,p.434, (1983) by NIMS AtomWorks

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

ASEとQuantum Espressoを使って、

他