Quantum ESPRESSO入力ファイル作成手順１４．Phonon計算（バンド構造）

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

#0. はじめに
第一原理計算ソフトQuantunm ESPRESSO <https://www.quantum-espresso.org> の入力ファイルの書き方について、複数回にわたって解説していきます。今回は第１４回目です。Phononのバンド計算について解説します。

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


#1. 計算手順
前回の記事では、ph.xを用いて単一q点でのPhonon計算について解説しました(<https://qiita.com/xa_member/items/77e464c01c5db1e847ea>)。今回のPhononのバンド計算では、以下の手順にて計算を実行します。

  1. SCFまたは構造最適化計算を実施
  2. Brillouin Zoneを網羅するq点群にてDFPTの計算を実施
  3. q点のFFT
  4. FFTメッシュによる内挿法にてPhononのバンド構造作成

先ずは、SCF計算を実施する。ただし、Phonon計算では最適化された構造が必要であるため、手元に最適化構造がない場合には SCF計算の代わりに構造最適化計算を実施する。その後、Brillouin Zoneを網羅する一様なq点メッシュを用いて、DFPT計算を実施する。このDFPT計算には、ph.xという実行モジュールを使用する。さらに、q2r.xという実行モジュールを用いて、各q点上のエネルギーをFFTにて実空間に変換する。FFTメッシュによる内挿法を適用すると、任意のq点にてPhononのエネルギーを計算可能となる。Brillouin Zone内のパスを指定して matdyn.xという実行モジュールを使用すると、任意のパスに沿ったPhononのバンド構造が計算される。

#2. Brillouin Zoneを網羅するq点群でのDFPT計算
DFPT計算(ph.x)の入力ファイルは、以下の通りである。

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

基本的な内容は前回(単一q点でのPhonon計算)と同様である。。ただし、ldispが.TRUE.となっている点が異なる。当該設定により、一様メッシュによるBrillouin Zoneを網羅するq点群でのPhonon計算が実行される。q点のメッシュ数は、nq1, nq2, nq3 にて指定する。多数のq点を含むため、計算には比較的に多くの時間が掛かる。


#3. q点のFFT
一様メッシュのq点群にてDFPT計算が完了した後、q2r.xにて各q点におけるエネルギーを実空間表示に変換する。q2r.xの入力ファイルの以下の通りである。

```
&INPUT
    fildyn = "pwscf.dyn"
    flfrc  = "pwscf.ifc"
    zasr   = "crystal"
/
```

fildynは、ポスト処理などで使用するデータを保存するファイルの名前であり、ph.x実行時と同様の名前を指定する。flfrcは、力の定数などを保存するファイル名である。zasrは、並進運動の自由度(音響フォノン)に対する Acoustic Sum Ruleの適用有無です。zasr = "crystal" を設定すると、並進運動を射影分離して計算精度が若干向上します。


#4. FFTメッシュによる内挿法にてPhononのバンド構造作成
実空間表示したPhononエネルギーを再びFFTすることで、任意のq点に対してエネルギー内挿することができる。当該処理には、matdyn.xを用いる。matdyn.xの入力ファイルは以下の通りである。

```
&INPUT
    amass(1)       =  28.0855
    asr            = "crystal"
    flfrc          = "pwscf.ifc"
    flfrq          = "pwscf.freq"
    q_in_band_form = .TRUE.
/
12
gG     20
X      20
W      20
K      20
gG     20
L      20
U      20
W      20
L      20
K      0
U      20
X      0
```

amassは各元素の質量であり、ph.xの実行時に指定したものと同様の値を設定する。asrは Acoustic Sum Ruleである。適用する場合には、"crystal"を指定する。flfrcは、q2r.xにて指定したものと同様の名前を設定する。flfrqは、振動数に関するデータが出力されるファイル名である。q_in_band_formは、Phononのバンドを出力するというフラグである。当該フラグには必ず.TRUE.を設定します。&INPUTのブロックより以降の行には、Brillouin Zone内のパスを記述します。パスの指定の仕方は、電子のバンド構造計算におけるK_POINTSと同様である。

上述の一連の計算を実行すると、Phononのバンド構造がファイルに出力される。出力されるファイル名は、flfrqにて指定したファイル名＋”.gp”という名前である。この例では、”pwscf.freq.gp”というファイル名である。


#5. Acoustic Sum Rule
q2r.x および matdyn.x の実行時に、Acoustic Sum Rule（ASR)の設定がありました。このでは、ASRの適用要否について、考察します。そもそも、ASRとは 系の並進運動の自由度を事前に取り除き、振動モードのみを抽出するという手法です。しかしながら、DFPT法にて計算されるHessian行列を対角化すると、その固有ベクトルは自ずと並進運動の自由度と振動モードに分離されます。つまり、本来的にはASRによる補正は不要なわけです。ではなぜ、Quantum ESPRESSOにはASRが実装されているのでしょうか？

Phonon計算においては、事前に構造最適化計算を実施します。この構造最適化計算において 原子に働く力が”完全に”ゼロに収束していれば、Phonon計算において 並進運動と振動モードは完全に分離されます。しかしながら、多くの場合、原子に働く力は完全にはゼロにならず、微小ながらも有限値を持ちます。また、Hessian行列の並進運動に関する固有値はゼロです。このため、並進運動の固有ベクトルは微小な力の収束誤差に影響を受けて、容易に並進運動と振動モードが混成されます。このコンタミネーションを回避するために、ASRが必要となります。構造最適化計算における力の収束閾値が十分に小さい場合には、並進運動は十分によく分離されるためASRは不要です。しかしながら、Quantum ESPRESSOのデフォルトの閾値程度で最適化している際には、ASRは適用した方がよさそうです。

---------------------------------------
なお、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)