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

#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計算は、以下の手順で実施される。

SCFまたは構造最適化計算を実施
DFPTの計算を実施
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サイトをご覧ください。

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