• Siを例にフォノンの分散と状態密度を計算します。
• 計算には超格子の作成が必要です。
• 化合物で見られる極性半導体ではLO-TO分裂というものを考えないといけませんが、それは別のページで解説します。
• 事前に cif2cell をインストールしてください。
  • CIF2Cell

SiのCIFファイルと擬ポテンシャルを入手します。 今回は結晶構造は Crystallography Open Database から、擬ポテンシャルはQEのウェブサイトから入手しました。

• CIFファイル … http://www.crystallography.net/cod/9008565.html
• 擬ポテンシャル … https://www.quantum-espresso.org/pseudopotentials/ps-library/si

擬ポテンシャルは Si.pbe-n-kjpaw_psl.1.0.0.UPF を使用します。作業ディレクトリ内にpseudoという名前のディレクトリを作り、そこに入れておきます。

CIFファイル名を9008565.cifとします。 cif2structを使って、普通の入力ファイルと、超格子の入力ファイルを作成します。

$ cif2cell -p pwscf -f 9008565.cif -o Si.pw.in
$ cif2cell --supercell=[3,3,3] -p pwscf -f 9008565.cif -o Si_333.pw.in

• 超格子のサイズは最初は2×2×2など小さいものから、結果が収束するまで徐々に大きくしていきます

超格子のインプットファイルに足りない部分を書き加えます。 cif2cellのデフォルトで格子定数はA（オングストローム単位）で指定されますが、celldm(1)（ボーア半径単位）に直しておきます。 また、原子に加わる力の大きさを出力するために tprnfor=.true. を忘れないようにします。

Si_333.pw.in

&control
   calculation='scf',
   prefix='Si_333',
   pseudo_dir = './pseudo/',
   outdir='./tmp/'
   tprnfor=.true.
/

&SYSTEM
  ibrav = 0
  celldm(1) = 10.2625
  nat = 54
  ntyp = 1
  ecutwfc = 40
  ecutrho = 320
/

&electrons
  conv_thr = 1.0d-8
/

K_POINTS automatic
 4 4 4 1 1 1

CELL_PARAMETERS {alat}
  1.500000000000000   1.500000000000000   0.000000000000000 
  1.500000000000000   0.000000000000000   1.500000000000000 
  0.000000000000000   1.500000000000000   1.500000000000000 
ATOMIC_SPECIES
  Si   28.08500  Si.pbe-n-kjpaw_psl.1.0.0.UPF
ATOMIC_POSITIONS {crystal}
Si   0.000000000000000   0.000000000000000   0.000000000000000 
Si   0.083333333333333   0.083333333333333   0.083333333333333 
（略）
Si   0.750000000000000   0.416666666666667   0.416666666666667

フォノンの計算に必要な原子間力定数(Interatomic force constant)の計算を行います

超格子の原子位置を少しだけずらしたQEの入力ファイルを作成します。 &cellフィールドと&positionフィールドには、さきほどの超格子のインプットファイルで得られた基本並進ベクトルと原子位置をはりつけてください。

Si_333.harmonic.in

&general
  PREFIX = Si_333
  MODE = suggest
  NAT = 54; NKD = 1
  KD = Si
/

&interaction
  NORDER = 1
/

&cell
  10.2625
  1.500000000000000   1.500000000000000   0.000000000000000 
  1.500000000000000   0.000000000000000   1.500000000000000 
  0.000000000000000   1.500000000000000   1.500000000000000 
/

&cutoff 
  Si-Si None
/

&position
1   0.000000000000000   0.000000000000000   0.000000000000000 
1   0.083333333333333   0.083333333333333   0.083333333333333 
（略）
1   0.750000000000000   0.416666666666667   0.416666666666667 
/

できたらALMを実行します。

$ alm Si_333.harmonic.in > Si_333.harmonic.log

すると Si_333.pattern_HARMONIC というファイルができるので、ALAMODEに付属のPythonスクリプトを使ってQEの入力ファイルを作成します。

$ python -m displace --QE Si_333.pw.in --mag 0.01 --prefix disp -pf Si_333.pattern_HARMONIC

• mag … 原子をどの程度ずらすかで、0.01から0.04程度の値（単位はオングストローム）が推奨されていますので、いろいろ試してください
• prefix … 生成されるQEの入力ファイル名に使われますので、自由に設定してください

今回の場合はQEの入力ファイルが1つだけ(disp1.pw.in)できます。 出てきた入力ファイルをすべてQEで実行します

$ pw.x < disp1.pw.in > disp1.pw.out

• 超格子の計算なので時間がかかります（普通のパソコンだと数時間くらい）。並列計算による実行をお勧めします。

disp1.pw.outに各原子にはたらく力が出力されていますので、これをALAMODEに付属のPythonスクリプトで整理します。

$ python -m extract --QE Si_333.pw.in *.pw.out > DFSET_harmonic

最後にALMを使って原子間力定数を計算します。 Si_333.harminic.inをコピーしてSi.harmonic_opt.inを生成し、次の点を変更します

• &generalフィールドで MODE=optimize に変更
• &optimizeフィールドを作成して、中に DFSET = DFSET_harmonic と記述

Si_333.harmonic_opt.in

&general
  PREFIX = Si_333
  MODE = optimize
  NAT = 54; NKD = 1
  KD = Si
/

&interaction
  NORDER = 1
/

&optimize
  DFSET = DFSET_harmonic
/

&cell
  10.2625
  1.500000000000000   1.500000000000000   0.000000000000000 
  1.500000000000000   0.000000000000000   1.500000000000000 
  0.000000000000000   1.500000000000000   1.500000000000000 
/

&cutoff 
  Si-Si None
/

&position
1   0.000000000000000   0.000000000000000   0.000000000000000 
1   0.083333333333333   0.083333333333333   0.083333333333333 
（略） 
1   0.750000000000000   0.416666666666667   0.416666666666667 
/

できたらALMを実行します

$ alm Si_333.harmonic_opt.in > Si.harmonic_opt.log

Si_333.fcs と Si_333.xml が出力されたら成功です。

Si_333.phband.in

&general
  PREFIX = Si_333
  MODE   = phonons
  FCSXML = Si_333.xml

  NKD = 1; KD = Si
  MASS = 28.0855
/

&cell
  10.2625
  0.500000000000000   0.500000000000000   0.000000000000000 
  0.500000000000000   0.000000000000000   0.500000000000000 
  0.000000000000000   0.500000000000000   0.500000000000000 
/

&kpoint
  1
  G 0.0 0.0 0.0 X 0.5 0.5 0.0 51
  X 0.5 0.5 1.0 G 0.0 0.0 0.0 51
  G 0.0 0.0 0.0 L 0.5 0.5 0.5 51
/

• &cellフィールドは、元の結晶の基本並進ベクトルです。Si.pw.inのものをコピーしてください
• &kpointフィールドの1行目の 1 は、バンド分散の計算を指定しています。2行目以降は分率座標で指定した波数の経路です

$ anphon Si_333.phband.in > Si_333.phband.log

フォノンの分散は Si_333.bands に出力されます

Si_333.phband.in をコピーして Si_333.phdos.in というファイルを作成します

• &generalフィールドでDELTA_E（エネルギーの刻み幅。単位はcm^-1。デフォルトは10でかなり粗い）を指定します。
• &kpointフィールドを状態密度の計算用に書き換えます。

Si_333.phdos.in

&general
  PREFIX = Si_333
  MODE   = phonons
  FCSXML = Si_333.xml

  NKD = 1; KD = Si
  MASS = 28.0855
  
  DELTA_E = 1
/

&cell
  10.2625
  0.500000000000000   0.500000000000000   0.000000000000000 
  0.500000000000000   0.000000000000000   0.500000000000000 
  0.000000000000000   0.500000000000000   0.500000000000000 
/

&kpoint
  2
  40 40 40
/

• &kpointフィールドの1行目の 2 は、状態密度の計算を指定しています。2行目はメッシュの数です。

$ anphon Si_333.phdos.in > Si_333.phdos.log

状態密度は Si_333.dos に出力されます

フォノンを自由ボソンとしたときの熱力学量（比熱、エントロピー、内部エネルギー、自由エネルギー）の温度依存性も Si_333.thermo に出力されます。 温度の刻みはデフォルトで DT = 10 (K)ですので、より細かく取りたい場合は&generalフィールドで DT の値を指定してください。