LAMMPSで熱伝導率計算(実践) その3 - こたつむりの備忘録

タイトルが安定しないです

今日もEMDを使った熱伝導率計算、Green-Kubo公式とにらめっこしました

まずはGreen-Kubo公式の表式と、LAMMPSスクリプトとの対応から

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

$\kappa$ : 熱伝導率
$k_B$ : ボルツマン定数
$T$ : 温度
$V$ : 計算モデルの体積
$\vec{J}(t)$ : 時間tにおける熱流束
$\langle \vec{J}(0) \cdot \vec{J}(t) \rangle$ : HCACF(Heat Current Auto-Correlation Function)

## 熱伝導率計算部分

# 熱流束を計算
compute             myKE all ke/atom
compute             myPE all pe/atom
compute             myStress all stress/atom NULL virial
compute             flux all heat/flux myKE myPE myStress
variable            Jx equal c_flux[1]/vol
variable            Jy equal c_flux[2]/vol
variable            Jz equal c_flux[3]/vol

# ミクロカノニカルアンサンブルを設定
fix                 1 all nve
# HCACFを計算
fix                 JJ all ave/correlate $s $p $d &
                    c_flux[1] c_flux[2] c_flux[3] type auto &
                    file profile.heatflux ave running

# xyz各方向の熱伝導率を計算
# trap: HCACF時間積分に相当, scale: 係数部分に相当
variable            scale equal ${convert}*$s*dt/$t/$t/vol/${kB}
variable            k11 equal trap(f_JJ[3])*${scale}
variable            k22 equal trap(f_JJ[4])*${scale}
variable            k33 equal trap(f_JJ[5])*${scale}
variable            kappa equal (v_k11+v_k22+v_k33)/3.0

thermo              $d
thermo_style        custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 v_kappa

run                 1000000  # 1,000,000 steps x 0.8 fs = 800 ps

色々読んでみて、correlation lengthに相当する$pは短すぎるとよくないらしい
dt$s$pが数ps程度になるとよい？

profile.heatfluxにHCACFが出力されているという情報があったのでこう考えたが、間違っている可能性もあります
いまいちHCACFの正体とLAMMPSスクリプトとの対応がわからん。。

fix ave/correlate command — LAMMPS documentation

よく見てみると、trap()関数については以下の記述がありました

The trap function defined for equal-style variables can be used to perform a time integration of this vector of datums, using a trapezoidal rule. This is useful for calculating various quantities which can be derived from time correlation data. If a normalization factor is needed for the time integration, it can be included in the variable formula or via the prefactor keyword.

(自動翻訳)

トラップ用に定義された関数に等しいスタイル変数は、台形法則を用いて、基準点のこのベクトルの時間積分を実行するために使用することができます。これは、時間相関データから導出できるさまざまな量を計算するのに便利です。時間積分に正規化係数が必要な場合は、変数式またはprefactorキーワードに含めることができます。

つまり、fix ave/correlateで計算したものがHCACFに相当し、時間積分がtrap()関数に相当するようですね

ところで、HCACFは時間経過とともに減衰して0に近づき、これが収束時間の目安になるらしいです
LAMMPSの出力からこの減衰を見るために以下の計算をしました

LAMMPSスクリプト

## 熱伝導率計算
compute             myKE all ke/atom
compute             myPE all pe/atom
compute             myStress all stress/atom NULL virial
compute             flux all heat/flux myKE myPE myStress

fix                 1 all nve
fix                 JJ all ave/correlate $s $p $d &
                    c_flux[1] c_flux[2] c_flux[3] type auto &
                    file profile.heatflux ave running

variable            scale equal ${convert}*$s*dt/$t/$t/vol/${kB}
variable            k11 equal trap(f_JJ[3])*${scale}
variable            k22 equal trap(f_JJ[4])*${scale}
variable            k33 equal trap(f_JJ[5])*${scale}
variable            kappa equal (v_k11+v_k22+v_k33)/3.0
variable            h11 equal sum(f_JJ[3])
variable            h22 equal sum(f_JJ[4])
variable            h33 equal sum(f_JJ[5])
variable            HCACF equal (v_h11+v_h22+v_h33)/3.0

thermo              $d
thermo_style        custom step temp v_h11 v_h22 v_h33 v_HCACF v_k11 v_k22 v_k33 v_kappa

run                 10000000

print               "running average conductivity: ${kappa}"

これで標準出力にthermo_styleで指定した熱力学データが出力されます
v_HCACFは累積和になっているので、該当部分を切り出してthermo.datというファイルに保存し、gnuplotで以下のスクリプトを使ってグラフにしました

## normalizeするために最初のデータを参照
# 0          100 -0.014850832 -0.0038024632 -0.00058573814    29083.279    1906.6538    45.242732    10345.059    18.835168    1.2348039  0.029300495    6.6997575
# 10000    98.637033 -0.0018503786  0.017331996 -0.0073324707    12161.638    17045.085    9943.1863     13049.97    11.653242    17.012881     8.619362    12.428495

set logscale x
set term png

set output "hcacf_xyz.png"
oldx = 0
oldy = 0
oldz = 0
plot "thermo.dat" using ($1*0.0008):(dx = ($6-oldx)/(12161.638-29083.279), oldx = $6, dx) w l title "HCACF x", \
  "thermo.dat" using ($1*0.0008):(dy = ($7-oldy)/(17045.085-1906.6538), oldy= $7, dy) w l title "HCACF y", \
  "thermo.dat" using ($1*0.0008):(dz = ($8-oldz)/(9943.1863-45.242732), oldz = $8, dz) w l title "HCACF z"

set output "hcacf.png"
old = 0
plot "thermo.dat" using ($1*0.0008):(dx = ($9-old)/(13049.97-10345.059), old = $9, dx) w l title "HCACF"


set output "kappa.png"
plot "thermo.dat" using ($1*0.0008):13 w l title "kappa"

f:id:justice_vsbr:20180805195145p:plain f:id:justice_vsbr:20180805195512p:plain

HCACFは確かに減衰していますが、文献ではこんなに振動したり負の値になったりせずなめらかに減衰しているので何かがおかしいですね。。
確かめなければいけないパラメータがいくつもあるので骨が折れそうです