热导率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 "*****************************************************"