クロスアビリティ 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計算の入力ファイルの実例