上一篇我們介紹了，本篇我們介紹如何由LAMMPS實現EMD 方法計算熱導率。
分子動力學模擬計算熱導率優缺點
相比其他常用計算熱導率方法（常用計算方法的簡介請參見），分子動力學計算方法能夠計算更大的體系，計算原子數通常在數百到數十萬個范圍，尺度在納米到微米區間。分子動力學模擬包含了聲子與聲子之間的各階散射，能夠更加真實地模擬聲子間的相互作用。但是該方法需要用到經典作用勢，計算熱導率的準確程度與所用的勢函數類型及參數密切相關。該方法也不適合用于低溫，通常要求體系的溫度要高于模擬物質的德拜溫度，如果溫度過低，需要考慮對溫度進行量子修正。
平衡態分子動力學
分子動力學計算熱導率的常用方法有平衡態分子動力學（EquilibriumMolecularDynamics,EMD）和非平衡態分子動力學(NonequilibriumMolecularDynamics,NEMD)。兩種體系在計算過程中都處于穩態（體系的溫度分布不隨時間發生變化），其中，前者不存在溫度梯度，通過原子間的微熱流來實現熱量交換，繼而由統計力學原理來得到熱導率。后者存在一個穩定的溫度梯度，類似于實驗中直觀的熱導率測量方法，所以又被稱為直接法（Direct method）,我們在后面的文章會詳細介紹采用NEMD計算熱導率的方法，在這里略去不表。
EMD的優勢&劣勢：尺寸效應相對較小，理論上采用周期性邊界條件可以實現計算無限大的體系。但是該方法要得到收斂的熱流需要較長的弛豫時間，需要比較大的計算量。
具體地，EMD方法計算熱導率由如下的Green-Kubo公式得到：
其中，
是熱導率張量，
別代表
三個方向，V是模擬體系的體積，
是玻爾茲曼常數，T表示模擬體系溫度，
表示熱流自相關函數，尖括號表示對熱流自相關函數進行系綜平均，t表示熱流自相關函數的積分上限，τ表示相關時間。對于三維各項同性材料，熱導率張量的非對角張量一定為零，我們通常會計算三個方向的熱導率數值
，然后取平均得到最終的熱導率數值。對于一維結構（單鏈、納米線等）只需要對某一個方向的熱流自相關函數進行積分，對二維結構，可以分別計算兩個面內兩個方向的熱導率然后取平均。
計算時需要對時間步長、熱流取樣間隔、相關時間步、總的記錄熱流時間步等參數做收斂性測試。通常最終的熱導率計算結果需要使用不同的初始速度計算多次，然后取平均值。
EMD方法計算實例
以下給出LAMMPS采用EMD方法計算固體氬晶格熱導率實例
[注]：滑動屏幕可以查看完整代碼
1# sample LAMMPS input script for thermal conductivity of liquid LJ`
2# Green-Kubo method via compute heat/flux and fix ave/correlate`
3# settings`
4variable ? ?x equal 10 ?`#模擬體系X方向長度`
5variable ? ?y equal 10 ?`#模擬體系Y方向長度`
6variable ? ?z equal 10 ?`#模擬體系Z方向長度`
7variable ? ?rho equal 0.6 `#固體氬晶格常數`
8variable ?t equal 1.35 `#模擬體系溫度`
9variable ? ?rc equal 2.5 ? ?`#截斷半徑`
10variable ? ?p equal 200 ? ? `# 熱流相關長度`
11variable ? ?s equal 10 ? ? ?`# 熱流取樣間隔`
12variable ? ?d equal $p*$s ? `# 輸出熱流自相關函數間隔`
13`# setup problem`
14units ? ? ? ?lj ?`#模擬體系采用的` **單位**`（具體各個物理量的單位查LAMMPS手冊）`
15timestep ? ? 0.005 `#時間步長0.005tau`
16atom_style ? ?atomic `#定義模擬體系原子所具有的屬性`
17lattice ? ? ? ?fcc ${rho} ? ? `#定義晶格類型和晶格常數`
18region ? ? ? ?box block 0 $x 0 $y 0 $z `#模擬體系大小`
19create_box ? ?1 box ? `#創建模擬體系box`
20create_atoms ? ?1 box ? `#創建模擬體系原子`
21mass ? ? ? ?1 1.0 ? `#給定原子質量`
22velocity ? ?all create $t 87287 `#給定特定溫度下的初始速度，其中87287是隨機數，可以任意取值`
23pair_style ? ?lj/cut ${rc} ? `#定義原子間作用勢和截斷半徑`
24pair_coeff ? ?1 1 1.0 1.0 `#勢函數參數值`
25neighbor ? ?0.3 bin `#定義近鄰原子`
26neigh_modify ? ?delay 0 every 1 `#創建近鄰原子列表`
27# 1st equilibration run
28fix ? ? ? ?1 all nvt temp $t $t 0.5 ?`#NVT系綜弛豫`
29thermo ? ? ? ?100 ? ? ? ? ? ? `#每100步輸出一次默認參數`
30run ? ? ? ?1000 ? ? ? ? ? ? ? ?`#跑1000步`
31velocity ? ?all scale $t ? ? ? ? ? `#調整體系溫度以達到給定值`
32unfix ? ? ? ?1
33# thermal conductivity calculation
34reset_timestep ?0 ? ? ? ? ? ? ? ?`#重新設置時間步為0`
35compute ? ? ? ? myKE all ke/atom ? ?`#計算每個原子的動能`
36compute ? ? ? ? myPE all pe/atom ? ?`#計算每個原子的勢能`
37compute ? ? ? ? myStress all stress/atom NULL virial ? ?`#計算每個原子的應力張量`
38compute ? ? ? ? flux all heat/flux myKE myPE myStress `#計算熱流`
39fix ? ? ? ? ? ? ? ?1 all nve ? `#NVE系綜下記錄熱流`
40fix ? ? ? ? ? ? JJ all ave/correlate $s $p $d & `#計算熱流自相關函數并輸出到指定文件`
41 ? ? ? ? ? ? ? ?c_flux[1] c_flux[2] c_flux[3] type auto &
42 ? ? ? ? ? ?file profile.heatflux ave running
43variable ? ? ? ?scale equal $s*dt/$t/$t/vol
44variable ? ? ? ?k11 equal trap(f_JJ[3])*${scale} `#對熱流自相關函數做積分`
45variable ? ? ? ?k22 equal trap(f_JJ[4])*${scale}
46variable ? ? ? ?k33 equal trap(f_JJ[5])*${scale}
47thermo ? ? ? ? ? ?$d `#輸出變量間隔`
48thermo_style ? ?custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 `#自定義輸出物理量類型`
49run ? ? ? ? ? ? 100000 ? ?`#跑100000步，可以得到10組熱流數據以及熱導率數據，采用最后一組數據`
50variable ? ? ? ?kappa equal (v_k11+v_k22+v_k33)/3.0 ? ?`#對三個方向熱導率做平均`
51print ? ? ? ? ? "running average conductivity: ${kappa}" `#在屏幕上輸出熱導率數值`
下圖為筆者采用上述LAMMPS腳本計算得到的10組熱流自相關函數和熱導率結果：
本文由邵成和鮑華編輯
原創文章，作者：計算搬磚工程師，如若轉載，請注明來源華算科技，注明出處：http://www.zzhhcy.com/index.php/2024/03/30/ad31e39a3e/