原文连接：http://mdbbs.org/thread-19218-1-1.html
看到有不少人在找热导率计算方面的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