lammps案例：石墨烯热导率模拟计算（EMD方法)

大家好，我是小马老师。
在前面的文章中，介绍了非平衡态下石墨烯的热导率模拟方法，本文介绍第二种热导率模拟方法：使用平衡态分子动力学（EMD）计算热导率。
本文仍然以石墨烯热导率计算为例，以供大家对比参考。
在平衡态下计算热导率，主要计算公式为Green-Kubo。
用到的主要命令为compute heat/flux。
用法为：
compute myFlux all heat/flux myKE myPE myStress
其中，myKE为原子动能，myPE为原子势能，myStress为原子应力。
在使用compute heat/flux命令前，必须提前计算出这三个量的值。
下面给出石墨烯EMD热导模拟代码，代码已经注释。

#模型参数设置
units metal           
dimension 3
boundary p p p        
atom_style atomic
neighbor 0.3 bin
neigh_modify delay 10
timestep 0.001
#读取石墨烯模型文件
read_data gp.data
mass 1 12.0107 
#使用tersoff力场
pair_style tersoff
pair_coeff * * SiC.tersoff C
#热力学信息输出
thermo 10000
thermo_style custom step temp vol press
#温度初始化
velocity all create 300 98989 dist gaussian
fix 1 all npt temp 300 300  0.1 x 0 0 1 y 0 0 1 
run 50000
unfix 1
reset_timestep 0

#Green-Kubo公式中需要的数据
variable A equal lx*ly*1.54      
variable dt equal 0.001          
variable cor equal 10000         
variable sam equal 10            
variable dum equal ${cor}*${sam} 
#单位转换公式
variable kB equal 8.6173324e-5  # eV/K Boltzmann                                                                                               
variable ev2J equal 1.60210e-19 # eV to J/mol                                                                                 
variable A2m equal 1.0e-10      # Angstrom to meter                                                                                           
variable ps2s equal 1.0e-12     # picoseconds to seconds
variable convert equal ${ev2J}/(${A2m}*${ps2s})
#计算热流
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]
variable Jy equal c_flux[2]
variable Jz equal c_flux[3]
fix JJ all ave/correlate ${sam} ${cor} ${dum} c_flux[1] c_flux[2] c_flux[3] type auto file SLG_J.J_t${temp_s}.dat ave running

#计算Green-Kubo公式系数
variable   scale equal (${dt}*${sam})/(${kB}*${temp_s}^2*$A) 
#计算Green-Kubo公式
variable k11 equal (trap(f_JJ[3])*${scale})*${convert}
variable k22 equal (trap(f_JJ[4])*${scale})*${convert}
variable k33 equal (trap(f_JJ[5])*${scale})*${convert}
variable k_total equal (v_k11+v_k22+v_k33)/3.0
fix Ph all phonon 10 10000 0 GAMMA SLG nasr 100 

#启动运算
thermo 100000
thermo_style custom step temp press v_k11 v_k22 v_k33 v_k_total 
fix 1 all nve 
run 1000000
#输出热导率
print "***************************************************** "
print " k11 conductivity : ${k11} [W/mK] @ ${temp_s} K"
print " k22 conductivity : ${k22} [W/mK] @ ${temp_s} K"
print " k33 conductivity : ${k33} [W/mK] @ ${temp_s} K"
print " 3D average conductivity: ${k_total} [W/mK] @ ${temp_s} K"
print "*****************************************************"

更多lammps教程，请关注微信公众号：
lammps加油站