原文鏈接：http://mdbbs.org/thread-19218-1-1.htmlhtml
看到有很多人在找熱導率計算方面的in文件，我就貢獻三個in文件吧，僅供參考。
同時，附件裏貼出了個人計算結果。EMD的輸出結果（compute heat/flux command+compute tc command的計算結果）中， 「ac.dat」（見附件中的"ac.wmf"）是熱流自相關函數（
我已經修改了compute_tc.cpp，目前輸出的是normalized HCACF，但結果中給出的仍是沒有歸一化的熱流自相關函數，但形狀和歸一化的是同樣的，請注意！
）隨m的變化，"tc.dat"(見附件中的"tc.wmf")是熱導率隨m的變化（m的涵義請參看熱導率計算的Green Kubo離散化公式,見附件"Comparison of atomic-level simulation methods for computing thermal conductivity」中的(9)式），"tc_time.dat"(見附件中的"tc_time.wmf")是熱導率隨時間的變化；NEMD的結果中，"temp.profile"(見附件中的"temp_distribution.wmf")是z方向上的溫度分佈，"thermal_conductivity.dat「(見附件中的"NEMD.wmf")表示熱導率隨時間的變化，二者都包含了fix heat command 和fix thermal/conductivity command的計算結果。
1.用compute heat/flux command+compute tc command獲得熱流自相關函數和熱導率（EMD方法）
# MD simulation of Ar thermal conductivity
# Initialization
unitslj
dimension3
newtonon
boundaryppp
atom_styleatomic
neighbor0.3bin
neigh_modifycheckyes
latticefcc0.844
regionboxblock -44-44-44units lattice
create_box1box
create_atoms1box
mass11.0
velocityallcreate0.71 458127641 mom yesrot yes dist gaussian units box
# LJ potential *********************************************************
pair_stylelj/cut 2.8
pair_coeff111.01.0#LJ parameters for Ar-Ar
fixtemp alltemp/berendsen 0.71 0.71 0.000466
fixnveallnve
thermo_stylecustom step temp etotal vol
thermo_modifylost warn
thermo100
# Run
timestep0.000466
run200000
reset_timestep0
# -------------- Flux calculation in nve ---------------
computemyKE all ke/atom
computemyPE all pe/atom
computemyStress all stress/atom virial
variablefactor_ac equal 1.0
variablefactor_tc equal 1.3806504e-23*sqrt(1.67e-21/6.633e-26)/3.405e-10^2
computejflux all heat/flux myKE myPE myStress
computetc all tc c_thermo_temp c_jflux v_factor_ac v_factor_tc iso first 10000 900000 100000
fixtc_outallave/time111c_tcfiletc_time.dat
thermo_stylecustomsteptemp
restart100000restart.*
run1000000
2. 用fix thermal/conductivity command獲得溫度梯度，進而獲得熱導率(NEMD方法)
# MD simulation of Ar thermal conductivity
# Initialization
unitslj
dimension3
newtonon
boundaryppp
atom_styleatomic
neighbor0.3bin
neigh_modifycheckyes
latticefcc0.844
regionboxblock -44-44-44units lattice
create_box1box
create_atoms1box
regionup1blockINF INFINF INF-0.5-0.25units lattice
regionup2blockINF INFINF INF0.50.75units lattice
regionupunion 2 up1 up2
regiondown1blockINF INFINF INF-3.5-3.25units lattice
regiondown2blockINF INFINF INF3.53.75units lattice
regiondown union 2 down1 down2
mass11.0
velocityallcreate0.71 458127641 mom yesrot yes dist gaussian units box
# Tersoff potential *********************************************************
pair_stylelj/cut 2.8
pair_coeff111.01.0#LJ parameters for Ar-Ar
fixtemp alltemp/berendsen 0.71 0.71 0.0466
fixnveallnve
computekeallke/atom
variabletemp atomc_ke/(1.5*1.0)
fixtemp_profileallave/spatial1100000100000zlower0.25v_tempfiletemp.profileunitslattice
computeup_tempalltemp/region up
computedown_tempalltemp/region down
variabledelta_tempequalc_up_temp-c_down_temp
fixdelta_outallave/time1100000100000v_delta_tempfiledelta_temp.dat
thermo_stylecustom step temp etotal vol
thermo_modifylost warn
thermo100
# Run
timestep0.000466
run100001
unfixtemp
fixheat_swapallthermal/conductivity10z32
fixe_exchangeallave/time1010000100000f_heat_swapfilee_exchange.dat
variablethermal_conductivity equal f_e_exchange/(0.000466*10.0*4.0*f_delta_out)*1.3806504e-23/3.405e-10/3.405e-10*sqrt(1.67e-21/6.633e-26)*6.0/8.0
# 以上variable命令須要特別注意，由於我所模擬的系統，盒子邊長Lx=Ly=Lz，熱導率計算公式通過推導變成爲e_exchange/(4.0*t*L*delta_T)，
# 爲了避免在in文件裏給L賦值，我修改
了fix_thermal_conductivity.cpp文件(見附件)，將e_exchange修改爲了
e_exchange += force->mvv2e * (all[0].value - all[1].value) / (domain->zprd); 同時在end_of_step()
裏添加了一句 「e_exchange = 0.0;「,詳見附件中的fix_thermal
[color=]_conductivity.cpp，這樣所得的e_exchange曲線基本上是一條水平曲線，而不是用原來的fix thermal/conductivity command所獲得的斜向上的曲線，請注意！！！
# 因此纔出現以上variable的表達式。
# 請看明白後再作計算，省得算出錯誤的結果！！！
fixthermal_conductivity_outallave/time1000001100000v_thermal_conductivityfilethermal_conductivity.dat
# Run
run10000000
3. 用fix heat command創建溫度梯度，進而獲得熱導率（NEMD方法）
# MD simulation of Ar thermal conductivity
# Initialization
unitslj
dimension3
newtonon
boundaryppp
atom_styleatomic
neighbor0.3bin
neigh_modifycheckyes
latticefcc0.844
regionboxblock -44-44-44units lattice
create_box1box
create_atoms1box
regionup1blockINF INFINF INF-0.5-0.25units lattice
regionup2blockINF INFINF INF0.50.75units lattice
regionupunion 2 up1 up2
regiondown1blockINF INFINF INF-3.5-3.25units lattice
regiondown2blockINF INFINF INF3.53.75units lattice
regiondown union 2 down1 down2
regionhotblockINF INFINF INF0.00.25units lattice
grouphotregionhot
regioncoldblockINF INFINF INF-4.0-3.75units lattice
groupcoldregioncold
mass11.0
#mass06.633e-26
#epsilon01.67e-21
#sigma03.405e-10
velocityallcreate0.71 458127641 mom yesrot yes dist gaussian units box
# Tersoff potential *********************************************************
pair_stylelj/cut 2.8
pair_coeff111.01.0#LJ parameters for Ar-Ar
fixtemp alltemp/berendsen 0.71 0.71 0.0466
fixnveallnve
computekeallke/atom
variabletemp atomc_ke/(1.5*1.0)
fixtemp_profileallave/spatial1100000100000zlower0.25v_tempfiletemp.profileunitslattice
computeup_tempalltemp/region up
computedown_tempalltemp/region down
variabledelta_tempequalc_up_temp-c_down_temp
fixdelta_outallave/time1100000100000v_delta_tempfiledelta_temp.dat
thermo_stylecustom step temp etotal vol
thermo_modifylost warn
thermo100
# Run
timestep0.000466
run100001
unfixtemp
fixhotallheat150region hot
fixcoldallheat1-50region cold
variablethermal_conductivity equal 50.0*0.5*1.67e-21/3.405e-10/sqrt(6.633e-26/1.67e-21)/((4.0*8.0*8.0*8.0/0.844)^(1.0/3.0)*3.405e-10*2.0*f_delta_out*1.67e-21/1.3806504e-23)*6.0/8.0
fixthermal_conductivity_outallave/time1000001100000v_thermal_conductivityfilethermal_conductivity.dat
# Run
run10000000
dom