Siを例に、フォノンの簡単な計算を実行してみます。

===== SCF計算 =====

事前にpw.xによるSCF計算を行います。

<file - Si.scf.in>
&control
   calculation='scf',
   prefix='Si'
   pseudo_dir = './pseudo/',
   outdir='./tmp/'
/
&system
   ibrav=2, celldm(1)=10.26, nat=2, ntyp=1,
   ecutwfc=18.0
/
&electrons
/
ATOMIC_SPECIES
 Si  28.086  Si.pz-vbc.UPF
ATOMIC_POSITIONS (alat)
 Si 0.00 0.00 0.00
 Si 0.25 0.25 0.25
K_POINTS automatic
 8 8 8 1 1 1
</file>

  * 金属の場合はsmearingを指定する必要がありますが、テトラヘドロン法はフォノン計算に使えません（以下のph.x実行時にエラーが出ます）。

<code>
$ pw.x < Si.scf.in > Si.scf.out
</code>


===== Γ点のフォノン =====

電子状態の計算が終われば、フォノンの振動数を計算できるようになる。 まず、k=(0,0,0)のΓ点におけるフォノンを計算してみる。 以下のファイルを用意する。

<file - Si.phG.in>
phonons of Si at Gamma
 &inputph
  tr2_ph=1.0d-14,
  prefix='Si',
  epsil=.true.,
  outdir='./tmp/',
  fildyn='Si.dynG',
 /
0.0 0.0 0.0
</file>

第1行はコメント行とみなされ無視されます。
また、最後の行は計算する波数（デカルト座標で、単位は格子定数aを使って2π/a）を表します。

各引数の意味は次のとおりです。

^変数^初期値^説明^
|tr2_ph|1.0d-12|フォノンの計算の収束条件。経験的に1.0d-14を使ったほうがよい|
|prefix|pwscf|SCF計算で使用したプレフィックスと同じにする|
|epsil|.false.|誘電率を計算する。Γ点の計算に必要。ただし.true.にする場合は、半導体かつΓ点の計算でないとエラーが出る|
|outdir|./|出力ファイルの場所。SCF計算のものと同じにする|
|fildyn|matdyn|結果であるdynamical matrixを出力するファイル名|

  * 金属の場合は、epsil=.false.として、計算するk点を 0.01 0.0 0.0 のようにΓ点から少しずらしてやるとうまくいきます。

実行は次のようにします。

<code>
$ ph.x < Si.phG.in > Si.phG.out
</code>

結果はfildynで指定したもの（今の場合はSi.dynG）に出力されます。

<file - Si.dynG>
（中略）
 **************************************************************************
     freq (    1) =      -0.170404 [THz] =      -5.684069 [cm-1]
 (  0.274931  0.000000 -0.504256  0.000000 -0.412479  0.000000 ) 
 (  0.274931  0.000000 -0.504256  0.000000 -0.412479  0.000000 ) 
     freq (    2) =      -0.170404 [THz] =      -5.684069 [cm-1]
 ( -0.161796  0.000000  0.380828  0.000000 -0.573404  0.000000 ) 
 ( -0.161796  0.000000  0.380828  0.000000 -0.573404  0.000000 ) 
     freq (    3) =      -0.170404 [THz] =      -5.684069 [cm-1]
 (  0.631059  0.000000  0.317327  0.000000  0.032690  0.000000 ) 
 (  0.631059  0.000000  0.317327  0.000000  0.032690  0.000000 ) 
     freq (    4) =      15.296893 [THz] =     510.249417 [cm-1]
 ( -0.463984  0.000000  0.270574  0.000000  0.459900  0.000000 ) 
 (  0.463984  0.000000 -0.270574  0.000000 -0.459900  0.000000 ) 
     freq (    5) =      15.296893 [THz] =     510.249417 [cm-1]
 (  0.527053  0.000000  0.327503  0.000000  0.339053  0.000000 ) 
 ( -0.527053  0.000000 -0.327503  0.000000 -0.339053  0.000000 ) 
     freq (    6) =      15.296893 [THz] =     510.249417 [cm-1]
 ( -0.083269  0.000000  0.565271  0.000000 -0.416575  0.000000 ) 
 (  0.083269  0.000000 -0.565271  0.000000  0.416575  0.000000 ) 
 **************************************************************************
</file>

  * 6つの振動モードがあって、それぞれのΓ点における振動数がわかる。
  * 音響モード（-0.170404 THz とあるが、これは0と見てよい）と光学モード（15.296893 THz）がそれぞれ三重に縮退している。
  * 音響モードの振動数を0にしたければ構造緩和を行うとよい。
  * 原子の変位ベクトルは、例えばfreq 4であれば1番目のSiが(-0.463984, 0.270574, 0.459900)方向、2番目のSiが(0.463984, -0.270574, -0.459900)である（デカルト座標）。
  * 振動モードの既約表現はSi.phG.out内に出力されている。


フォノンの振動に対応する変位ベクトルはXcrysDenを使って見ることができます。
次のファイルを用意します。

<file - Si.dynmatG.in>
&inputfil
  dyn ='Si.dynG',
  filxsf = 'dynmatG.axsf'
  asr ='simple',
  lperm = .true.,
  q(1)=1.0,
  q(2)=0.0,
  q(3)=0.0
/
</file>

  * Γ点のフォノンは q->0 の極限で定義されます。lpermは q=(q(1), q(2), q(3)) に沿って q->0 の極限をとることを指定します。
    * Si.phG.inで波数を 0.01 0.0 0.0 とした場合はasrおよびlpermの行を削除して q(1) = 0.01 としてください。
  * asrはacoustic sum ruleの略で、'simple'にするとΓ点で音響フォノンの振動数がゼロになるようにします。デフォルトは'no'（何もしない）です。このほか、次のものが指定できます
    * zero-dim ... 分子の計算で使います
    * one-dim ... 1次元系の計算（カーボンナノチューブなど）に使います
    * crystal ... simpleよりも計算精度があがりますが、その分コストも増えます。simpleでうまく行かない場合はこちらにします。

実行は次のようにします。
<code>
$ dynmat.x < Si.dynmatG.in > Si.dynmatG.out
</code>

すると dynmatG.axsf というファイルができるので、以下の手順で変位ベクトルを表示します。

  - XCrysDenを起動して、[File] -> [Open Structure] -> [Open AXSF] よりdynmatG.axsfを開く
  - ウィンドウが開くので、矢印ボタンを押して変位ベクトルを表示したいフォノンを選ぶ（エネルギーの低い順です）
  - いったん [Hide] ボタンを押してウィンドウを最小化する
  - [Display] から [Forces] を選ぶ
  - このままだと矢印が長すぎるので、[Modify] -> [Force Settings] でLength Factorの大きさを調節する（デフォルトは200だが30くらいがちょうどよい）

===== X点のフォノン =====

今度はX点(2π/a,0,0)のフォノンを調べます。
手順はさきほどと同じですが、誘電率の計算が不要です。

<file - Si.phX.in>
phonons of Si at X
 &inputph
  tr2_ph=1.0d-14,
  prefix='Si',
  outdir='./tmp/',
  fildyn='Si.dynX',
 /
1.0 0.0 0.0
</file>

  $ ph.x < Si.phX.in > Si.phX.out

結果は Si.dynX に出力されます。

<file -Si.dynX>
（中略）
 **************************************************************************
     freq (    1) =       4.207901 [THz] =     140.360482 [cm-1]
 ( -0.000000  0.000000 -0.023972  0.000000  0.706700  0.000000 ) 
 ( -0.000000  0.000000  0.706700  0.000000 -0.023972  0.000000 ) 
     freq (    2) =       4.207901 [THz] =     140.360482 [cm-1]
 ( -0.000000  0.000000 -0.706700  0.000000 -0.023972  0.000000 ) 
 ( -0.000000  0.000000 -0.023972  0.000000 -0.706700  0.000000 ) 
     freq (    3) =      12.236850 [THz] =     408.177391 [cm-1]
 ( -0.999882  0.000000  0.000000  0.000000  0.000000  0.000000 ) 
 ( -0.015351  0.000000 -0.000000  0.000000  0.000000  0.000000 ) 
     freq (    4) =      12.236850 [THz] =     408.177391 [cm-1]
 ( -0.015351  0.000000  0.000000  0.000000  0.000000  0.000000 ) 
 (  0.999882  0.000000 -0.000000  0.000000  0.000000  0.000000 ) 
     freq (    5) =      13.746712 [THz] =     458.540957 [cm-1]
 (  0.000000  0.000000 -0.000000  0.000000  0.707107  0.000000 ) 
 ( -0.000000  0.000000 -0.707107  0.000000  0.000000  0.000000 ) 
     freq (    6) =      13.746712 [THz] =     458.540957 [cm-1]
 (  0.000000  0.000000  0.707107  0.000000  0.000000  0.000000 ) 
 (  0.000000  0.000000  0.000000  0.000000 -0.707107  0.000000 ) 
 **************************************************************************
</file>

2重に縮退したものが3つあることがわかります。

変位ベクトルを計算する場合は次のファイルを用意して、Γ点と同じように計算してください。


<file - Si.dynmatX.in>
&inputfil
  dyn ='Si.dynX',
  filxsf = 'dynmatX.axsf'
/
</file>

<code>
$ dynmat.x < Si.dynmatX.in > Si.dynmatX.out
</code>