Quantum ESPRESSO入力ファイル作成手順１３．Phonon計算（特定q点）

クロスアビリティ Winmostarサポートチームです。
WinmostarはQuautum ESPRESSOを簡単かつ高度に利用するための統合GUI環境となっています。学生は無料で利用でき、学生以外も無料トライアルを入手することができますので、是非 [WinmostarのWebサイト](https://winmostar.com/jp/) をご覧ください。

#0. はじめに
第一原理計算ソフトQuantunm ESPRESSO <https://www.quantum-espresso.org> の入力ファイルの書き方について、複数回にわたって解説していきます。今回は第１３回目です。Phonon計算について解説します。Phonon計算には、特定のq点にて振動モードを計算するものと、Brillouinゾーン内のすべてのq点について計算を実行してPhononの分散関係を計算するものがあります。本記事では、前者が対象です。後者については、次回の記事で解説します。

> この記事では、初歩的な内容しか紹介しません。
> 入力ファイルの詳細について知りたい方は、マニュアル http://www.quantum-espresso.org/Doc/INPUT_PH.html をご参照下さい。


#1. Phonon
Phononは原子核の運動を量子化した擬似粒子です。このため、Phononを取り扱う際には、原子座標の変位に対する全エネルギーの変化を計算する必要があります。具体的には、全エネルギーを原子座標で2階微分する必要があります。1階微分である力は、Hellmann-Feynman定理に従って SCF計算完了後に容易に評価できます。しかしながら、2階微分の計算は容易ではありません。一般には 有限差分法による数値微分と、DFPTによる解析的微分の2種類の計算方法があるのですが、Quantum ESPRESSOでは後者が採用されています。DFPTでは、SCF計算完了後に、原子座標の変位を用いた摂動計算を実行します。


#2. q点
原子座標を変位させる場合、すぐに思いつくのが、全てのセルに含まれる原子を同一方向に動かすというやり方です。しかしながら、周期境界条件においては、これ以外にも変位の方法があります。例えば、隣接するセルの原子を互いに逆並行に変位させる方法です。この場合においても変位の周期性は担保されており、変位は \pi / a の波数を持つことになります。さらに一般化すれば、波数ベクトルqを有する変位を定義することが出来ます。ただし、波数ベクトルqの定義域は、電子波動関数におけるk点と同様に、Brillouinゾーンの内側に限定されます。k点と同様のアナロジーで、この波数ベクトルをq点と呼びます。Quantum ESPRESSOでは、任意のq点に対してDFPTが適用可能となっています(VASPなどの他のソフトではΓ点のみの場合がある。)。


#3. 計算手順
単一q点でのPhonon計算は、以下の手順で実施される。

  1. SCFまたは構造最適化計算を実施
  2. DFPTの計算を実施
  3. DFPTのポスト処理

先ずは、SCF計算を実施する。ただし、Phonon計算では最適化された構造が必要であるため、手元に最適化構造がない場合には SCF計算の代わりに構造最適化計算を実施する。その後、DFPT計算を実施する。当該計算には、実行モジュールをph.xを用いる。DFPT計算が完了した後、ポスト処理を実行すると振動モードに関するデータが生成される。ポスト処理には、dynmat.xという実行モジュールを使用する。


#4. DFPT計算
DFPT計算(ph.x)の入力ファイルは、以下の通りである。

```
This is DFPT calculation
&INPUTPH
    fildyn       = "pwscf.dyn"
    ldisp        = .FALSE.
    epsil        = .FALSE.
    amass(1)     =  28.0855
    alpha_mix(1) =  0.5
/
  0.000000   0.000000   0.000000
```

１行目はコメント行である。２～８行目がDFPTの計算条件の設定である。fildynは、ポスト処理などで使用するデータを保存するファイルの名前である。ldispは、Phononの分散関係を計算するためのフラグである。今回は単一q点での計算なので、FALSEとする。epsilがTRUEのとき、Phononと同時に、格子誘電率テンソルが計算される。ただし、格子誘電率テンソルは、系が絶縁体(occupations = "fixed")で 且つ q点がΓ点の場合にのみ計算可能である。amass(1, 2, ...)は各原子種の質量である。単位は、amuである。alpha_mix(1)は、DFPTの繰り返し計算における電荷密度混合率である。0.0 ~ 1.0 の値が設定可能である。収束性が悪い場合には、少し小さな値(0.2 ~ 0.3)を設定すると良い。入力ファイルの最終行は、q点ベクトルの値である。単位は、2pi/aである(aは格子定数)。これにより、指定された単一のq点におけるPhononが計算できる。

DFPT計算が完了すると、標準出力にPhononの振動数が表示される。Phononの振動数モードは、原子座標の自由度に相当する。一つの原子当たりX,Y,Zの3つの自由度があるため、ユニットセル内に2原子がある場合には振動モードの数は 2 x 3 = 6 である。

```
     q = (    0.000000000   0.000000000   0.000000000 ) 

 **************************************************************************
     freq (    1) =       0.332357 [THz] =      11.086246 [cm-1]
     freq (    2) =       0.332357 [THz] =      11.086246 [cm-1]
     freq (    3) =       0.332357 [THz] =      11.086246 [cm-1]
     freq (    4) =      15.282554 [THz] =     509.771114 [cm-1]
     freq (    5) =      15.282554 [THz] =     509.771114 [cm-1]
     freq (    6) =      15.282554 [THz] =     509.771114 [cm-1]
 **************************************************************************
```


#5. DFPTポスト処理
DFPT計算(ph.x)が完了した後、ポスト処理(dynmat.x)を実行すると、可視化可能な振動モードのファイルが出力される。dynmat.xの入力ファイルは、以下の通りである。

```
&INPUT
    asr    = "crystal"
    fildyn = "pwscf.dyn"
    filxsf = "pwscf.dyn.axsf"
/
```

asrは、Acoustic Sum Ruleを意味しており、q点がΓ点の場合にのみ利用可能である。計算対象が結晶系の場合には、asr = "crystal"と設定すると並進運動の自由度(音響フォノン)を振動モードから除去して、純粋な振動(光学フォノン)のみを抽出してくれる。また、直線形の孤立分子の場合には asr = "one-dim"、非直線形の孤立分子の場合には asr = "zero-dim" とすることで、並進運動および回転運動の自由度を除去できる。fildynは、ph.xにて指定したものと同様のファイル名を設定する。filxsfは、振動モードのデータを出力するファイルの名前である。当該ファイルは、XCrySDenにて可視化可能である。

---------------------------------------
なお、Winmostarを利用するとQuautum ESPRESSOをGUI上から簡単かつ高度に活用することができます。学生は無料で利用でき、学生以外も無料トライアルを入手することができますので、是非 [WinmostarのWebサイト](https://winmostar.com/jp/) をご覧ください。

Quantum ESPRESSO入力ファイル作成手順シリーズ
[１．結晶構造の作成](https://qiita.com/xa_member/items/727c1a62930611babaf7) 
[２．SCF計算の設定](https://qiita.com/xa_member/items/a098e94d633a1d29f1c0) 
[３．擬ポテンシャルファイルの選択方法](https://qiita.com/xa_member/items/006d9dd44c662a17903d) 
[４．原子座標のみの最適化](https://qiita.com/xa_member/items/dc860b1656f7dfe93c8a)
[５．格子ベクトル及び格子内部の原子座標の最適化](https://qiita.com/xa_member/items/4c667817e50facd8a770)
[６．状態密度の計算](https://qiita.com/xa_member/items/7aa75304b2ce7f90c8d4)
[７．局所状態密度の計算](https://qiita.com/xa_member/items/95c7bd64981f7329b94b)
[８．バンド構造の計算](https://qiita.com/xa_member/items/a519c53b63df75a68894)
[９．バンド数の設定](https://qiita.com/xa_member/items/d923ea537328cdba8b10)
[１０．van der Waals相互作用](https://qiita.com/xa_member/items/bd64362329c7235bf248)
[１１．LDA+U法](https://qiita.com/xa_member/items/ab3e1333e92bc424f5a7)
[１２．NEB法](https://qiita.com/xa_member/items/e3b63a1fdca12f95e8cb)
[１３．Phonon計算（特定q点）](https://qiita.com/xa_member/items/77e464c01c5db1e847ea)
[１４．Phonon計算（バンド構造）](https://qiita.com/xa_member/items/5c02fcae08431588659d)
[１５．擬ポテンシャルの作成](https://qiita.com/xa_member/items/5feb8f432170141c89cd)
[１６．擬ポテンシャルのテスト](https://qiita.com/xa_member/items/3902c1660fb7064a4c34)
[１７．SCF計算の入力ファイルの実例](https://qiita.com/xa_member/items/febde40c8a17b7b83cc5)