はじめに

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

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

今回は、教科書2章,図2.1で紹介されている、Cuが単純立方格子の場合に、格子定数を手動で変えた計算により、格子定数を最適化、また、Birch Murnagham状態方程式を使って平衡体積弾性率を見積もってみます。

前準備

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

# import libs
import os
import subprocess
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# improt ase
from ase.build import bulk
from ase.visualize import view
from ase.calculators.espresso import Espresso

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

関数の定義

モデルや計算のための関数を定義します。
今回のモデルでは、単純立方格子(sc)を使います。

*.ipynb

def make_crystal_lattice(a=None):
    # Build model for DFT calculations
    model = bulk(name='Cu', crystalstructure='sc', a=a, orthorhombic=False, cubic=False, basis=None)
    return model

def calc_dft_cu():
    # ASEを使ってQuantum EspressoでDFTするためのセッティング
    # 並列化する場合はcmdの"1"を適宜書き換え
    pseudopotentials = {'Cu':pp_file}
    #実行環境に合わせて、必要に応じて下記を変えてください。
    cmd = 'pw.x -in espresso.pwi > espresso.pwo'
    input_data = {
        'control':{'pseudo_dir':'./'},
        'system': {
            'occupations' : 'smearing',
            'smearing' : 'mp',
            'degauss' : 0.001,
            'ecutwfc': 21.461,     #21.461[Ry] ≒ 292[eV]
            'ecutrho': 85.0},      #ecutwfc * 4
        'disk_io': 'low'}
    calc = Espresso(command=cmd,
                    pseudopotentials=pseudopotentials,
                    kpts=(12, 12, 12),
                    input_data=input_data)
    return calc

def move_results_file(dir_for_result):
    # QEの計算ファイルを別のフォルダにまとめて格納する
    os.makedirs(dir_for_result, exist_ok=True)
    subprocess.call('mv espresso.* pwscf.save pwscf.xml '+dir_for_result+'/.', shell=True)
    subprocess.call('cp *.UPF '+dir_for_result+'/.', shell=True)

###DFT計算を行う
格子定数aを2.2～2.9Åで変量して、エネルギーを計算してみます。

*.ipynb

lattice_list = np.arange(2.2,2.95,0.05)
lattice_list = np.delete(lattice_list, [-2,-4,-6]) #教科書と揃えるため少し削除

r_ene=[]
for lattice in lattice_list:
    print(f"Lattice constant :{lattice} [Angs.]")
    Cu = make_crystal_lattice(lattice)
    calc = calc_dft_cu()
    Cu.set_calculator(calc)
    r_ene.append(Cu.get_potential_energy())
    print(f"   {r_ene[-1]} [eV]")
    fname="lattice_"+'{:.2f}'.format(lattice)
    move_results_file(fname) #各格子定数毎の結果をフォルダに移動

r_ene = np.array(r_ene)

###出力結果### 
    Lattice constant :2.2 [Angs.]
       -1654.2331551425614 [eV]
    Lattice constant :2.25 [Angs.]
       -1654.434029306143 [eV]
    ......

結果の解析

得られたエネルギーを使って解析します。
教科書に従って、格子定数の変化に対するエネルギーの変化を２次関数やBirch Murnagham状態方程式に従ってFittingして、最安定の格子定数や体積弾性率を求めてみます。

*.ipynb

#Difine Fitting functions
def poly2(lattice_a, a0, beta, E0):
    E = beta*(lattice_a**2) -2*a0*beta*lattice_a + E0 - beta*a0**2
    #E = a0*(lattice_a**2) -2*beta*lattice_a + E0    
    return E

def Birch_Murnagham_EoS(a, a0, B0, B02, E0):
    #a:lattice_a
    E = 9*(a0**3)*B0/16 * ( B02*( (a0/a)**2-1)**3 + ( ((a0/a)**2-1)**2 * (6-4*(a0/a)**2) ) ) +E0
    return E

まずは２次関数を使ってみます。やってみると分かりますが、すべてのデータを使うとあまりよくフィッティングできません。平衡位置周辺のみのデータを使ってフィットしてみます。

*.ipynb

# Fit the poly2 function and estimate energy
print('Fit Quadratic Polynomial')
fit_params, std = curve_fit(poly2, lattice_list[1:-4], r_ene[1:-4])
est_poly2 = poly2(lattice_list, *fit_params)
print(f' a0 = {fit_params[0]}  [Angstrom]')
print(f' B0 = {2/9*(1/fit_params[0])*fit_params[1]*160.21766208}  [GPa]')
#B0 [eV/Angstrom^3] -> [GPa]

###出力結果### 
    Fit Quadratic Polynomial
     a0 = 2.42  [Angstrom]
     B0 = 108.8  [GPa]

格子定数は教科書とほぼ一致、平衡体積弾性率は若干高め？でしょうか。

次は Birch Murnagham状態方程式を使ってみます。こちらは全データでFit。ただし、初期値を与えないと全然違うところでFitしてしまうので、初期値を与えます。

*.ipynb

# fit the Birch_Murnagham_EoS and estimate energy by it
print('Fit Birch Murnagham EoS')
init_param = np.array([2.2, 120, 0, 0]) #curve_fitの初期値与えないとうまく収束しない
fit_params, std = curve_fit(Birch_Murnagham_EoS, lattice_list, r_ene, p0=init_param, maxfev=10000)
est_BMeos = Birch_Murnagham_EoS(lattice_list, *fit_params)
print(f' a0 = {fit_params[0]}  [Angstrom]')
print(f' B0 = {fit_params[1]*160.21766208}  [GPa]')

###出力結果### 
    Fit Birch Murnagham EoS
     a0 = 2.41  [Angstrom]
     B0 = 100.0  [GPa]

a0もB0もそれっぽい値です。とはいえ教科書とずれているのは、アプリケーションや擬ポテンシャルが違うことが原因でしょうか。

最後にDFTの計算結果や関数によるFitting結果を図示してみましょう。

*.ipynb

# make a fig.
plt.rcParams["font.size"] = 12
fig = plt.figure(figsize=(5,4))
plt.xlabel("a [Angstrom]")
plt.ylabel("Etot/atom [eV]")

plt.scatter(lattice_list,r_ene,label="DFT")
plt.plot(lattice_list,est_poly2,label="fit_poly2",c='gray',linestyle="dashed")
plt.plot(lattice_list,est_BMeos,label="fit_BM-EoS",c='black')

plt.legend()
plt.show()
#fig.savefig("result.png")
#print("図2.1 単純立方格子Cuの格子定数aに対する全エネルギーEtot(a)")

　二次関数では平衡位置からずれると乖離が見られますが、状態方程式の方は、平衡位置からずれた位置でもエネルギーと格子定数の関係をうまく表現できていることが確認できます。

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

ASEとQuantum Espressoを使って、

他