前回

justice-vsbr.hatenablog.com

GaNで熱伝導率計算をしてみよう

SWポテンシャルの下で構造最適化したw-GaNに対して、EMDで熱伝導率を解析してみます。

とりあえず実行

## 変数設定
variable            p       equal 200     # correlation length
variable            s       equal 10      # sample interval
variable            d       equal $p*$s   # dump interval
variable            t       equal 100     # temperature
variable            tstep   equal 0.001   # timestep
variable            kB      equal 1.3806504e-23
variable            eV2J    equal 1.6021766e-19
variable            ps2s    equal 1.0000000e-12
variable            ang2m   equal 1.0000000e-10
variable            convert equal ${eV2J}*${eV2J}/${ps2s}/${ang2m}

## 基本設定
units               metal
dimension           3
boundary            p p p
atom_style          atomic
newton              on
timestep            ${tstep}

## 構造データ読み込み、初期速度設定
read_data           GaN.sw-optimized
velocity            all     create  $t 87287 mom yes rot yes dist gaussian

## ペアポテンシャル
pair_style          sw
pair_coeff          * *     GaN.sw      Ga  N

## 熱平衡状態計算
fix                 1 all nvt temp $t $t 0.01 drag 0.2
thermo              1000
run                 10000
velocity            all scale $t
unfix               1

## 熱伝導率計算：略

#=> kappa = -15.4177052995385 (W/mK)

🤔🤔🤔

(c.f.100KにおけるGaNの熱伝導率の実験値は100(W/mK)のオーダー)

とりあえずMDのtimestepを1fs~10fsで変えながら計算してみたところ、以下のようになりました(9,10fsの場合はtimestep大きすぎて系が爆発したので計算できず)

f:id:justice_vsbr:20180803003617p:plain

収束も何もあったもんじゃないですね笑
一応4fsではそれらしい値を示していますが偶然だと思います

もちろん他にもいじるべきパラメータがたくさんあるので、どこをどうしたらちゃんと計算できるのか、これから調べます
特に、Green-Kubo公式を解くために熱流束の自己相関関数のアンサンブル平均(もう既にナンノコッチャですが)を計算しますが、この計算におけるパラメータ(p,s,d)がどういう役割を果たすのかわかってないです

読み物

https://arxiv.org/pdf/1005.4478.pdf

熱伝導率を求めるGreen-Kubo公式の表式は以下の通り

$\kappa = \dfrac{1}{2k_{B}T^{2}V} \int_0^{\infty} dt \langle \vec{J}(0) \cdot \vec{J}(t) \rangle$

係数はいいとして、重要なのは0~∞の時間積分と、被積分関数である熱流束の自己相関関数のアンサンブル平均です

シミュレーションでは有限時間しか計算できないので上の式を厳密には計算できず、なるべく長い時間(タイムステップ数)で計算することで真の値に収束させる方針になります

この論文では金属であるSiやGeを対象としているため、デバイ温度(この温度以上では電子による熱伝導への影響を無視できるという指標)以上の1000Kを設定温度にしていますが、GaNは半導体でありもともと格子振動だけが熱キャリアと考えられます
HCACFが素早く減衰するのは低波長フォノンの寄与が大きいことを示すらしいですが、存在するフォノンの波長は計算モデルの大きさに制限されます。論文中では計算モデルを大きくしても同じように減衰したことから、計算モデルの大きさとタイムステップ数はここで十分、としています
各アンサンブル平均におけるHCACFの標準偏差と平均からnumerical errorを見てみると、ある時点(tc)から急激に大きくなります。よってあまりに長時間計算してもこのnumerical errorが累積してしまいよくない。さらにこのtcは初期速度や計算モデルサイズにより変化し、決めうちでタイムステップ数を決めるのも難しい…
有限サイズ効果もあるため、計算モデルサイズを少しずつ大きくし、サチるところを探しています

読み取れた計算条件は

タイムステップ: 0.8fs
平衡状態計算: 10⁵ steps
熱伝導率計算: 2 x 10⁶ steps
- サンプリング: every 1 step
- correlation length: 10³
- 試行回数: 8

ひとまずここまで
進捗があり次第更新します

こたつむりの備忘録

研究・技術的なことをメモしていきます

LAMMPSで熱伝導率計算(実践) - EMD

前回

GaNで熱伝導率計算をしてみよう

とりあえず実行

読み物