概要

Siを例に熱伝導度などの熱輸送関係の量を計算します。
熱伝導の計算にはフォノンの3次の非調和項が必要で、計算が非常に重いです。クラスタ計算機の使用をおすすめします。
事前にフォノンの分散と状態密度のうち、原子間力定数の計算まで済ませておいてください。

計算手順

3次の原子間力定数の計算

有限温度の計算で重要な寄与がある3次の原子間力定数の計算を行います

フォノンの分散と状態密度で作成した Si_333.harmonic.in をコピーして Si_333.cubic.in を作成してください。そして、次のようにファイルを編集してください

&generalフィールドのPREFIXをSi_333_cubicに変更
&interactionフィールドのNORDERを2に変更（NORDER+1次の原子間力が考慮されます）
&cutoffフィールドで Si-Si None 7.3 に変更
- これは2次の項はカットオフを設けません（None）が、3次の項は7.3 bohrの範囲内にあるSi原子間の力のみを計算するという意味です。今の場合は次近接まで考えることになります。
- 原子間の距離は Si_333.harmonic.log 内のINTERACTIONパートに出力されています。

Si_333.cubic.in

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

&interaction
  NORDER = 2
/

&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 7.3
/

&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.cubic.in > Si_333.cubic.log

Si_333_cubic.pattern_HARMONIC, Si333_cubic.pattern_ANHARM3 の2つのファイルができます。
Si_333_cubic.pattern_HARMONIC は調和項の原子間力を求めるときに使った Si_333.pattern_HARMONIC と同一です。

ALAMODEに付属のPythonスクリプトを使ってQEの入力ファイルを作成します。

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

今回の場合はQEの入力ファイル(disp01.pw.in, disp02.pw.in, …)が20個できます。出てきた入力ファイルをすべてQEで実行します

$ pw.x < disp01.pw.in > disp01.pw.out
$ pw.x < disp02.pw.in > disp02.pw.out
（以下略）
$ pw.x < disp20.pw.in > disp20.pw.out

普通のパソコンですと1つ計算するのに半日かかるので、スクリプトを書いてクラスタ計算機で並列的にやらせるのが良いです。

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

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

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

&generalフィールドでMODE=optimizeに変更
&optimizeフィールドを作成して、中に DFSET = DFSET_cubic と FC2XML = Si333.xml と記述（Si333.xmlは調和項の原子間力で得たファイル）

Si_333.cubic_opt.in

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

&interaction
  NORDER = 2
/

&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 7.3
/

&optimize
  DFSET = DFSET_cubic
  FC2XML = Si333.xml
/

&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.cubic_opt.in > Si_333.cubic_opt.log

Si_333_cubic.fcs, Si_333_cubic.xmlが出力されたら成功です。

熱伝導度

以下の Si_333.rta.in というファイルを作成します

Si_333.rta.in

&general
  PREFIX = Si_333
  MODE   = RTA
  FCSXML = Si333_cubic.xml

  NKD = 1; KD = Si
  MASS = 28.0855

  DT = 1.0
/

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

&kpoint
  2
  10 10 10
/

計算は調和項のときに比べて重いので、最初は少ないk点から始めましょう
データの解析には Si_333.result を使います。このファイルがあれば計算は再スタートになりますので、もしk点を変えて再計算を行う場合は削除するかファイル名を別のものに変更してください。

$ anphon Si.rta.in > Si.rta.out

結果は Si_333.kl に出力されます

ちなみに、この結果だと絶対零度で熱伝導度が発散するという物理的におかしな結果になります。低温では端の効果が大きいので、これを取り込んだ解析を行います

$ python -m analyze_phonons --calc kappa_boundary --size 1.0e+6 Si_333.result > Si_333_boundary_1mm.kl

--sizeでサンプルの大きさ（単位はnm）を指定します
他にもいろいろ解析できるので、詳細はマニュアルを読んでください

熱伝導スペクトル

Si_333.rta.in をコピーして Si.rta2.in というファイルを作成します。以下の点を編集します。

&generalフィールドに EMIN = 0; EMAX = 550; DELTA_E = 1.0 を追加
&analysisフィールドを作成して、中に KAPPA_SPEC = 1 と記述

Si_333.rta2.in

&general
  PREFIX = Si_333
  MODE   = RTA
  FCSXML = Si333_cubic.xml

  NKD = 1; KD = Si
  MASS = 28.0855

  EMIN = 0; EMAX = 550; DELTA_E = 1.0
/

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

&kpoint
  2
  30 30 30
/

&analysis
 KAPPA_SPEC = 1
/

$ anphon Si_333.rta2.in > Si_333.rta2.out

スペクトルの計算なのでk点数がたくさん必要です。収束は必ず確認しましょう。

結果は Si_333.kl_spec に出力されます。さまざまな温度のスペクトルがまとめて出力されるので、特定の温度（例えば300K）の結果を抜き出したいときは次のようにします（チュートリアルのページにあった方法）

$ awk '{if ($1 == 300.0) print $0}' Si_333.kl_spec > Si_333_300K.kl_spec

Table of Contents

概要

計算手順

3次の原子間力定数の計算

熱伝導度

熱伝導スペクトル