lammps计算热导率大体可以分为三种方法：EMD、NEMD、rNEMD。
本文介绍非平衡模拟（NEMD）热导率的计算方法。
NEMD下热导率计算公式为:
k=dQ*dZ/dt
k：热导率
dQ:能量通量。
dT/dZ:温度梯度。
非平衡态下计算热导率的原理是在材料的两端或者中间设置热源和冷源。
在热源输入热量，在冷源抽取等量的热量，在这个过程中计算材料内部会有热量的传递，产生温度梯度。
在lammps模拟中得到能量通量和温度梯度，代入公式即可求得热导率。
热量的输入方式有多种，如fix langevin、fix heat、fix ehex。
本文以官方案例in.heat为例，介绍fix heat命令注入热源并计算热导率的方法。
#模拟体系在xyz三个方向的尺寸
variable x equal 10
variable y equal 10
variable z equal 20
#晶格
variable rho equal 0.6
#温度
variable t equal 1.35
#截断半径
variable rc equal 2.5
#模型参数，lj单位
units lj
atom_style atomic
#生成box、填充原子
lattice fcc ${rho}
region box block 0 $x 0 $y 0 $z
create_box 1 box
create_atoms 1 box
mass 1 1.0
#温度初始化
velocity all create $t 87287
#力场参数设置，lj力场
pair_style lj/cut ${rc}
pair_coeff 1 1 1.0 1.0
#近邻列表定义
neighbor 0.3 bin
neigh_modify delay 0 every 1
#设置热源、冷源区域
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
#nvt下进行平衡模拟
fix 1 all nvt temp $t $t 0.5
thermo 100
run 1000
#调整温度达到设定温度
velocity all scale $t
unfix 1
#设定nve系综
fix 1 all nve
#热源处持续输入热量
fix hot all heat 1 100.0 region hot
#冷源处持续抽走热量
fix cold all heat 1 -100.0 region cold
#设定热力学参数输出
thermo_style custom step temp c_Thot c_Tcold
thermo 1000
#运行10000步
run 10000
#计算原子动能
compute ke all ke/atom
#根据原子动能计算原子温度
variable temp atom c_ke/1.5
#沿z方向进行切块，把原子划入到不同的块中
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.heat
#计算温差
variable tdiff equal f_2[11][3]-f_2[1][3]
#计算平均温差
fix ave all ave/time 1 1 1000 v_tdiff ave running start 13000
#输出平均温差
thermo_style custom step temp c_Thot c_Tcold v_tdiff f_ave
#运行20000步
run 20000
dQ可通过fix heat命令中注入的热量计算。
温度梯度中的dT可从log文件中得到温差，dZ为z方向尺度。
得到这三个数据代入公式即可求得热导率。
官方案例单位为lj，如果换成其他材料，请注意单位转换！
END
不定期分享lammps小技巧，敬请扫码关注微信公众号：lammps加油站。
