Table of Contents
概要
- 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