ASEとQuantum Espressoを使って、手動でCuの単位格子最適化して、平衡体積弾性率を見積る

### はじめに
　Jupyter上で、Pythonライブラリである[ASE(Atomic Simulation Environment)](https://wiki.fysik.dtu.dk/ase/)を介して、第一原理計算ソフト[Quantum Espresso](https://www.quantum-espresso.org/)で、第一原理計算を実施してみます。  
　本記事は、書籍『[密度半関数理論入門(D.S.ショール，J.A.ステッケル 著、佐々木泰造，末原茂 共訳](http://www.yoshiokasyoten.sakura.ne.jp/single_titles/ISBN978-4-8427-0365-7.html)』(以下、教科書)の例題や課題を参考に、プログラムを作って計算します。

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

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


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

```python
# 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)を使います。

```python:*.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Åで変量して、エネルギーを計算してみます。

```python:*.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して、最安定の格子定数や体積弾性率を求めてみます。

```python:*.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
```

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

```python:*.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してしまうので、初期値を与えます。

```python:*.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結果を図示してみましょう。

```python:*.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)")
```

![output_15_0.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/302038/312b376d-4914-2b14-5a42-3527620b6f80.png)


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

-----

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

**ASEとQuantum Espressoを使って、**
- [手動でCuの単位格子最適化して、平衡体積弾性率を見積る](https://qiita.com/egoo/items/3b68017fedd3d160014a)
- [Ptが単純立方構造、fcc、hpcのいずれの構造をとるかを調べる](https://qiita.com/egoo/items/8fa70a073df10782caaa)
- [六方晶Hfの格子定数aとcを最適化する](https://qiita.com/egoo/items/003669a4bc5da4f503ca)
- [H2Oなどの３原子分子を構造を最適化して、安定構造を見つける](https://qiita.com/egoo/items/b6c4b8a1f5280fbfd718)

**他**
- [ASEでモデルの図を保存する,描写する,モデルを出力する(png,cif,traj,pov)](https://qiita.com/egoo/items/55be3a27cd8a83c26ef2)
- [Quantum Espressoを色々なライブラリでコンパイルして、自分のPCに最速な状態にする](https://qiita.com/egoo/items/7fc55b4f7fa9c4b2f87d)