首页/文章/ 详情

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

320科技工作室

1年前浏览5118

大家好，我是小马老师。

在前面的文章中，介绍了非平衡态下石墨烯的热导率模拟方法，本文介绍第二种热导率模拟方法：使用平衡态分子动力学（EMD）计算热导率。

本文仍然以石墨烯热导率计算为例，以供大家对比参考。

在平衡态下计算热导率，主要计算公式为Green-Kubo。

用到的主要命令为compute heat/flux。

用法为：

compute myFlux all heat/flux myKE myPE myStress

其中，myKE为原子动能，myPE为原子势能，myStress为原子应力。

在使用compute heat/flux命令前，必须提前计算出这三个量的值。

下面给出石墨烯EMD热导模拟代码，代码已经注释。data文件可自己建模。

案例代码

#模型参数设置units metaldimension 3boundary p p patom_style atomicneighbor 0.3 binneigh_modify delay 10timestep 0.001#读取石墨烯模型文件read_data gp.datamass 1 12.0107 #使用tersoff力场pair_style tersoffpair_coeff * * SiC.tersoff C#热力学信息输出thermo 10000thermo_style custom step temp vol press#温度初始化velocity all create 300 98989 dist gaussianfix 1 all npt temp 300 300  0.1 x 0 0 1 y 0 0 1 run 50000unfix 1reset_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 secondsvariable convert equal ${ev2J}/(${A2m}*${ps2s})#计算热流compute myKE all ke/atomcompute myPE all pe/atomcompute myStress all stress/atom NULL virialcompute 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 100000thermo_style custom step temp press v_k11 v_k22 v_k33 v_k_totalfix 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 "*****************************************************"