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

#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のバンド計算では、以下の手順にて計算を実行します。

SCFまたは構造最適化計算を実施
Brillouin Zoneを網羅するq点群にてDFPTの計算を実施
q点のFFT
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サイトをご覧ください。

Quantum ESPRESSO入力ファイル作成手順シリーズ
１．結晶構造の作成
 ２．SCF計算の設定
 ３．擬ポテンシャルファイルの選択方法
 ４．原子座標のみの最適化
 ５．格子ベクトル及び格子内部の原子座標の最適化
 ６．状態密度の計算
 ７．局所状態密度の計算
 ８．バンド構造の計算
 ９．バンド数の設定
 １０．van der Waals相互作用
 １１．LDA+U法
 １２．NEB法
 １３．Phonon計算（特定q点）
１４．Phonon計算（バンド構造）
１５．擬ポテンシャルの作成
 １６．擬ポテンシャルのテスト
 １７．SCF計算の入力ファイルの実例