LAMMPSで熱伝導率計算
LAMMPS
LAMMPS Molecular Dynamics Simulator
C++で書かれた高速で動作する分子動力学法(molecular dynamics: MD)計算パッケージ
詳しい説明は省きます(自分も理解してない部分が多い)
熱伝導率解析
MDで熱伝導率を解析する手法は主に２つに分類されます
非平衡MD(non-equilibrium MD: NEMD):
系に高温領域と低温領域を設け、その間に発生する熱流束と温度勾配から、フーリエ則を用いて熱伝導率を解析平衡MD(equilibrium MD: EMD):
グリーン–久保公式 - Wikipedia
Green-Kubo公式は平衡状態における物理量の"揺らぎ"から輸送係数(今回は熱伝導率)を求めるもの
正直に言うとこの"平衡"と"非平衡"の違い、完全には理解してません。。
ただ、熱力学的平衡 - Wikipediaによると
熱力学的に非平衡 (non-equilibrium) であるとは、上記の熱的、力学的、化学的平衡のいずれかが満たされていない状態であり、系に物質またはエネルギーの正味の流れ、あるいは相転移などが生じる。またこのような非平衡状態は不安定であるため別の状態へ転移するが、転移速度が極めて遅いために不安定な状態が維持される場合、この状態を準安定状態という。
とある。
NEMDはシミュレーション中高温領域で加熱&低温領域で放熱してエネルギーの流れがあるという点で非平衡、EMDは断熱状態であるという点で平衡、という理解でいいんでしょうかね？
サンプルコードとその解説
LAMMPSのexample/KAPPA/ディレクトリには熱伝導率解析を行うためのサンプルコードとその実行ログが入ってます
6. How-to discussions — LAMMPS documentation
READMEをざっと訳すと、
このディレクトリにはレナードジョーンズ流体に対して5つの方法(NEMD4つ、EMD1つ)で熱伝導率を解析するためのスクリプトが入っている
ただし、あくまでサンプルであって、十分な平衡状態・精度が保証されているものではない
固体に対しても同様に適用可能
in.langevin: ランジュバン動力学を用いて高温/低温領域の温度を制御(NEMD)
in.heat: HEXアルゴリズムを用いて高温/低温領域の熱を制御(NEMD)
in.ehex: eHEXアルゴリズムを用いて(ry
in.mp: Muller-Plathe法を用いて運動エネルギーが最大/最小の粒子同士のエネルギーを交換(NEMD)
in.heatflux: Green-Kubo公式を用いる(EMD)
熱伝導率を計算するためには、NEMDの場合計算後にログファイルから熱流束と温度勾配を計算する必要があり、EMDの場合はスクリプトの中でそのまま計算、表示することができる
LAMMPSはコマンドが大量にあり、スクリプトをパッと見てもさっぱりわからないので、上にある5つの手法それぞれのサンプルスクリプトを可能な範囲で解説していきます。間違っていたら指摘してください。。
in.langevin
## 変数設定 # 単位格子を各軸方向に何倍するか variable x equal 10 variable y equal 10 variable z equal 20 # レナードジョーンズ流体のパラメータ(単位系はlj) variable rho equal 0.6 # 密度 variable t equal 1.35 # 初期温度 variable rc equal 2.5 # カットオフ半径 variable tlo equal 1.0 # 低温領域温度 variable thi equal 1.70 # 高温領域温度 ## 初期設定 units lj # 単位系 atom_style atomic # 原子が持つ属性情報を指定(atomicはデフォルト) lattice fcc ${rho} # 単位格子を指定 region box block 0 $x 0 $y 0 $z # 単位格子をx,y,z倍した領域を作成 & "box"とタグ付け create_box 1 box # "box"を計算モデルに設定 create_atoms 1 box # 計算モデルにlatticeで指定した構造に従い原子を配置 mass 1 1.0 # type1の原子の質量を指定 velocity all create $t 87287 # 温度t&シード値87287の下でランダムに速度を設定 pair_style lj/cut ${rc} # 原子間ポテンシャルとしてレナードジョーンズポテンシャルを指定 pair_coeff 1 1 1.0 1.0 # ポテンシャルに用いるパラメータ neighbor 0.3 bin # 近接原子リストを作成するアルゴリズムを指定 neigh_modify delay 0 every 1 # 近接原子リストをいつ更新するかを指定 ## 高温/低温領域作成(この段階ではまだ温度設定はしない) # z方向に熱流束を発生させるため、xy平面方向に広がるような領域"hot", "cold"を作成 # それぞれの領域の温度を計算して"Thot", "Tcold"とタグ付け region hot block INF INF INF INF 0 1 region cold block INF INF INF INF 10 11 compute Thot all temp/region hot compute Tcold all temp/region cold ## 第一段階：全体の温度を$tで均一にして平衡状態になるまでMD # NVTアンサンブル, 能勢-フーバーサーモスタット, 温度はtで一定 fix 1 all nvt temp $t $t 0.5 thermo 100 # 熱動力データ(温度や圧力など)を100ステップごとにロギング run 1000 # 1000ステップMD計算を実施 velocity all scale $t # 現在の温度を計算し、温度tになるように速度をスケーリング unfix 1 # 上のfixコマンドの設定を破棄 ## 第二段階：用意した高温/低温領域に温度を設定して再度平衡状態になるまでMD fix 1 all nve # NVEアンサンブルはミクロカノニカルアンサンブルであり、全体が断熱状態となる # 高温/低温領域にランジュバンサーモスタットを設定, tally yesで熱浴と授受したエネルギーを監視 fix hot all langevin ${thi} ${thi} 1.0 59804 tally yes fix cold all langevin ${tlo} ${tlo} 1.0 287859 tally yes # fixとcomputeをタグを使って紐付け(hot-Thot, cold-Tcold) fix_modify hot temp Thot fix_modify cold temp Tcold variable tdiff equal c_Thot-c_Tcold # c_XXは、"compute XX"で計算した値に置換される # ロギングする数値を設定 f_XX, v_XXは同様にfixコマンド、variableコマンドで計算した値に置換される thermo_style custom step temp c_Thot c_Tcold f_hot f_cold v_tdiff thermo 1000 run 10000 ## 第三段階：熱伝導率計算 compute ke all ke/atom # 運動エネルギー(kinetic energy)を計算 variable temp atom c_ke/1.5 # atom styleを指定することで、ログファイルに出力する度に原子ごとの平均を計算してくれる # 高温/低温領域をリセット fix hot all langevin ${thi} ${thi} 1.0 59804 tally yes fix cold all langevin ${tlo} ${tlo} 1.0 287859 tally yes fix_modify hot temp Thot fix_modify cold temp Tcold # 高温/低温領域の実際の温度差(v_tdiff)の累積時間平均を計算 # この設定では、10, 20, ..., 1000ステップの平均値を1000ステップ目における値とする fix ave all ave/time 10 100 1000 v_tdiff ave running thermo_style custom step temp c_Thot c_Tcold f_hot f_cold v_tdiff f_ave # z軸に沿って原子を20個(=1/0.05)のグループに分割 compute layers all chunk/atom bin/1d z lower 0.05 units reduced # グループごとにv_tempを計算、平均をとってprofile.langevinに出力 fix 2 all ave/chunk 10 100 1000 layers v_temp file profile.langevin run 20000
in.heat
NEMDなので流れはlangevinと同じ
主な違いはfixコマンドにある
tdfifの時間平均の仕方がlangevinと違う理由は不明
## 第二段階 # 毎ステップ, それぞれのregionに対して100の熱量(単位はunitsコマンドに従う)を加算/減算 fix hot all heat 1 100.0 region hot fix cold all heat 1 -100.0 region cold ## 第三段階 fix ave all ave/time 1 1 1000 v_tdiff ave running start 13000
in.ehex
in.heatとほぼ同じで、fix heatコマンドがfix ehexコマンドになっているだけ
in.mp
第二段階からの流れが少し違う 少し複雑なアルゴリズムになっているので、詳細は以下のリンク参照
fix thermal/conductivity command — LAMMPS documentation
## 第二段階 compute ke all ke/atom variable temp atom c_ke/1.5 fix 1 all nve compute layers all chunk/atom bin/1d z lower 0.05 units reduced fix 2 all ave/chunk 10 100 1000 layers v_temp file profile.mp # モデルをz方向に20分割し、Muller-Platheアルゴリズムに従って10ステップごとに1つの原子ペアのz方向の運動エネルギーを交換 fix 3 all thermal/conductivity 10 z 20 variable tdiff equal f_2[11][3]-f_2[1][3] thermo_style custom step temp epair etotal f_3 v_tdiff thermo 1000 run 20000 # 他のNEMDより多めのステップを取る ## 第三段階 # リセット fix 3 all thermal/conductivity 10 z 20 fix ave all ave/time 1 1 1000 v_tdiff ave running thermo_style custom step temp epair etotal f_3 v_tdiff f_ave run 20000
in.heatflux
以下のリンクも参照
compute heat/flux command — LAMMPS documentation
## 変数設定 # x,y,z,rho,t,rc,tlo,thiはNEMDと同じなので略 # Green-Kubo公式では熱流束 variable p equal 200 # 熱流束の自己相関を計算するときの長さ variable s equal 10 # 熱流束をサンプリングする間隔 variable d equal $p*$s # 出力間隔 ## 初期設定 # NEMDと同じなので略 ## 第一段階：熱平衡状態になるまで緩和 # NEMDと同じなので略 ## 第二段階：熱伝導率計算 reset_timestep 0 # タイムステップカウントをリセット # KE:運動エネルギー、PE:ポテンシャルエネルギー、Stress:圧力を計算 # 上記3つを使って熱流束fluxを計算 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 # 計算した熱流束を元にその時間相関を計算し、その平均を取る # type autoを指定することで自己相関を計算 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 $s*dt/$t/$t/vol variable k11 equal trap(f_JJ[3])*${scale} variable k22 equal trap(f_JJ[4])*${scale} variable k33 equal trap(f_JJ[5])*${scale} thermo $d thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 run 100000 # 精度よく求めるためには長時間のシミュレーションをする必要がある variable kappa equal (v_k11+v_k22+v_k33)/3.0 print "running average conductivity: ${kappa}"
注意事項
格子定数がそのポテンシャルに対して最適かどうか
上記のNEMD/EMD計算中は計算モデルの大きさは固定されていますが、
原子間ポテンシャルに対して大きすぎたり小さすぎたりすると、圧力値がおかしなことになり負の熱伝導率が出力されたりします(EMDで確認)
LAMMPSで構造最適化も可能なので、事前に行っておく必要があります
熱伝導率計算に入る前に平衡状態になっているかどうか確認すること
サンプルでは流体を取り扱っている(流体と固体どう区別するのかは謎)ので少ないタイムステップで平衡状態になるが、そのまま流用すると平衡状態にならないまま次に進んでしまうことがあります
サンプルスクリプトに準備段階におけるロギングも仕込まれているので、これを元に温度やエネルギーの変化をよく見ておく必要があります
単位系を確認
compute heat/flux command — LAMMPS documentation
リンク先ページ下部サンプルスクリプトにあるように、unitsコマンドで指定したLAMMPSにおける単位系からSI単位系に変換するような工夫が必要
サンプルではunits realとしているため、それに応じて変換している
# 略 # convert from LAMMPS real units to SI variable kB equal 1.3806504e-23 # [J/K] Boltzmann variable kCal2J equal 4186.0/6.02214e23 variable A2m equal 1.0e-10 variable fs2s equal 1.0e-15 variable convert equal ${kCal2J}*${kCal2J}/${fs2s}/${A2m} ... variable scale equal ${convert}/${kB}/$T/$T/$V*$s*${dt} variable k11 equal trap(f_JJ[3])*${scale} variable k22 equal trap(f_JJ[4])*${scale} variable k33 equal trap(f_JJ[5])*${scale}
