This is an old revision of the document!

概要

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

計算手順

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

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

Si.alm0.in をコピーして Si.alm_cubic.in を作成してください。そして、次のようにファイルを編集してください

&generalフィールドのPREFIXをSi333_cubicに変更
&interactionフィールドのNORDERを2に変更
&cutoffフィールドで Si-Si None 7.3 に変更
- これは7.3 bohrの範囲内にあるSi原子間の力のみを計算するという意味です。今の場合は次近接まで考えることになります。
- 原子間の距離は Si.alm0.log 内のINTERACTIONパートでわかります。

Si.alm_cubic.in

&general
  PREFIX = Si333_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
/

&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.alm_cubic.in > Si.alm_cubic.log

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

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

$ python -m displace --QE Si333.pw.in --mag 0.01 --prefix disp -pf Si333_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

めちゃくちゃ時間がかかるので、スクリプトを書いてクラスタ計算機で並列的にやらせるのが良いです。

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

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

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

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

Si.alm1.in

&general
  PREFIX = Si333_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.alm_cubic_optimize.in > Si.alm_cubic_optimize.log

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

熱伝導度

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

Si_phdos.in

&general
  PREFIX = Si333_cubic
  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点から始めましょう
データの解析にはsi333_cubic.resultを使います。このファイルがあれば計算は再スタートになりますので、もしk点を変えて再計算を行う場合は削除するかファイル名を別のものに変更してください。

$ anphon Si_RTA.in > Si_RTA.out

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

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

$ analyze_phonons -m --calc kappa_boundary --size 1.0e+6 si333_cubic.result > si333_boundary_1mm.kl

--sizeでサンプルの大きさ（単位はnm）を指定します

その他、解析できる熱輸送量

フォノンの寿命

$ python -m analyze_phonons --calc cumulative --temp 300 Si333_cubic.result > Si333_tau300K.kl

累積熱伝導度

$ python -m analyze_phonons --calc cumulative --temp 300 --length 10000:5 Si333_cubic.result > Si333_tau300K.kl

熱伝導スペクトル

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

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