9. フォノンの計算

フォノンの計算方法について紹介します。 quantum ESPRESSO では密度汎関数摂動論(Density Functional Perturbation Theory (DFPT))を用いたフォノンの計算が可能です。 単純な Frozen phonon 法では、フォノンの波数にあわせて大きなユニットセルをとらなければなりませんが、DFPT では 常に最も小さなユニットセル内の計算だけでフォノンが計算できます。(ただし、それでもフォノン計算はかなり重くなります。)

フォノンのバンド構造を計算する流れとしては、

1 pw.x を使って scf 計算

2 ph.x を使って粗いｑ点メッシュ上の dynamical matrix を計算

3 q2r.x を使って 2 の結果を実空間にフーリエ変換し、force constantを計算

4 3 の結果を matdyn.x を使って任意の波数にフーリエ変換

となります。3,4 の過程ではacoustic sum rule を用いることで、計算誤差による imaginary phonon mode がでないようにすることができます。

注釈

imaginary phonon mode とは振動数ωの2乗が計算の結果負になってしまったモードのことを指します。(出力では - \sqrt{\lvert \omega \rvert} の値が出力されます) 可能性としては計算がちゃんと収束していない場合と、構造がそのモードに対して不安定な場合の2通りが考えられます。

以下では、diamond のフォノン計算を例に 1-4 のステップの説明をしていきます。 (graphene でのフォノン計算は収束させるための条件がきついので、チュートリアル向きではありません。興味があれば以下を参考に試してみて下さい。)

9.1. scf 計算

フォノン計算に向けて特に気をつけるべきことはありません。 diamond の scf インプットファイル(diamond.scf.in)は以下の通りです。

&control
    calculation = 'scf'
    restart_mode='from_scratch',
    prefix='diamond',
    tstress = .true.
    tprnfor = .true.
    pseudo_dir = './',
    outdir='./work/'
 /
 &system
    ibrav=  2,
    celldm(1) = 6.6559,
    nat=  2,
    ntyp= 1,
    ecutwfc = 30.0,
    ecutrho = 150.0,
 /
 &electrons
    mixing_mode = 'plain'
    mixing_beta = 0.7
    conv_thr =  1.0d-8
 /
ATOMIC_SPECIES
 C  12.0107  C.pz-van_ak.UPF
ATOMIC_POSITIONS
C        0.00   0.00   0.00
C        0.25   0.25   0.25
K_POINTS {automatic}
 8 8 8 1 1 1

9.2. phonon 計算

次にフォノンの計算に移ります。 計算にはこれまでの pw.x ではなく ph.x を用いるため、pw.x のものとは全く異なります。 オプションの詳細は、パッケージに含まれる INPUT_PH.html を参考にして下さい。

具体的には以下のようなインプットファイル(diamond.ph.in)を用意します。

diamond phonon
&inputph
	tr2_ph=1.0d-14,
	prefix='diamond',
	outdir='./work/'
	fildyn='diamond.dyn',
	amass(1)=12
	ldisp=.true.
	nq1 = 2
	nq2 = 2
	nq3 = 2
/

フォノンの計算では一行目はコメントになります。2行目以降のそれぞれのオプションの意味は以下の通りです。

	tr2_ph=1.0d-14,

収束条件です。default は tr2_ph=1.d-12 ですが、私の経験上 tr2_ph=1.0d-14 はあった方がよいと思います。

	prefix='diamond',
	outdir='./work/'

これらは、scf 計算と同じものを指定します。

	fildyn='diamond.dyn',

dynamical matrix の計算結果を記録するファイル名です。

	amass(1)=12

炭素原子の質量です。 何も指定しなければ、scf 計算のときに指定した質量が計算に使われます。

	ldisp=.true.
	nq1 = 2
	nq2 = 2
	nq3 = 2

フォノンのq点として 2x2x2 のメッシュ状の点をとることを意味しています。 scf計算に用いるk点と違い、必ずΓ点を含むようなメッシュをとります。

このインプットファイルを以下のコマンドで実行します。 (それなりに時間がかかります。)

% ph.x < diamond.ph.in > diamond.ph.out

これで、diamond.ph.out には2x2x2 のメッシュ状のq点でのフォノンの情報が出力されます。 このメッシュ状の点の情報だけでよければ、ここで計算は終わりになります。 また、あるq点だけの計算がしたければ、ldisp=.false.にして、インプットファイルの最後に計算したいq点のリストを書くということも可能です。例えばΓ点だけの計算でよければ、以下のようなインプットで計算できます。

diamond phonon
&inputph
	tr2_ph=1.0d-14,
	prefix='diamond',
	outdir='./work/'
	fildyn='diamond.dyn',
	amass(1)=12
	ldisp=.false.
/
0.0  0.0  0.0

9.3. dynamical matrixのフーリエ変換(q->r)

次に、ph.x で計算した dynamical matrix を実空間にフーリエ変換します。 以下のインプットファイル(q2r.in)を用意します。

&input
  fildyn = 'diamond.dyn'
  flfrc = 'diamond.fc'
/

fildyn は diamond.ph.in で指定した dynamical matrix のファイル名、flfrc はそれをフーリエ変換した出力を保存するファイル名です。 このインプットを用いて以下のコマンドを実行します。

% q2r.x < q2r.in > q2r.out

q2r.x へのインプットは、 INPUT_Q2R.html を参考にして下さい。

9.4. dynamical matrixのフーリエ変換(r->q)

最後に任意の波数のフォノンの情報を matdyn.x を用いて計算します。 matdyn.x の入力ファイルの詳細は INPUT_MATDYN.html を参照してください。

これによりフォノンのDOSやバンドが計算できます。

DOSを計算する場合は以下のインプット(matdyn.dos.in)を使います。

&input
  flfrc = 'diamond.fc'
  asr = 'crystal'
  dos = .true.
  nk1 = 16,
  nk2 = 16,
  nk3 = 16,
/

まず q2r.x で作成したファイル名を flfrc で指定します。 asr は acoustic sum rule の条件を課すことを明示しており、これにより、音響モードが q=0 でomega=0 にいくようになります。（計算の精度の問題で、asrを課さないとq=0 でも完全には omega=0 になりません。） なお、asr としてはナノチューブのような一次元系の場合の回転に対するモードも考慮することができます。(asr='one-dim')

dos=.true. でDOSの計算をすることを指定し、nk1,nk2,nk3でDOS計算のためのメッシュを指定します。

このファイルを使って

% matdyn.x < matdyn.dos.in

とやれば、phonon の DOSが matdyn.dos というファイルに出力されます。 DOS 計算におけるエネルギーステップはdeltaE というオプションで指定でき、default は 1cm^-1 になっています。

band を計算する場合は、以下のようなインプット(matdyn.freq.in)を使います。

&input
  asr = 'crystal'
  flfrc = 'diamond.fc'
  dos=.false.
  q_in_band_form = .true.
/
3
0.50    0.50     0.50   25
0.00    0.00     0.00   25
0.00    0.00     1.00   25

DOS計算との違いは、dos=.false.として、計算するq点の情報を入力する点です。 q_in_band_form = .true.とすると、band計算のときと同様に高対称点の情報を入れていけばよいことになります。 座標の単位は、q_in_cryst_coord オプションで指定できます。（デフォルトは q_in_cryst_coord = .false.）

% matdyn.x < matdyn.freq.in

とすれば、matdyn.freq というファイルにバンド計算の結果が出力されます。 また、matdyn.freq.gp にはプロット用の出力も得られます。 これをプロットするとフォノンのバンドは以下のようになります。

_images/band3.png