Dear All,
I am trying to calculate the thermal conductivity of graphene with LAMMPS. I am using the Muller-Plathe method.
But the value I got is 2 orders larger than the reported value( 4.8410^3 ~ 5.3 10^3 W/mK, Nano Lett., 2008, 8(3), pp902-907). I checked the temperature distribution on my system
and it looks quite linear. I really appreciate if anyone could give me some advice on this. Pleae refer to my code
posted here:
units metal
boundary p p p
atom_style molecular
lattice custom 2.4595 a1 1 0 0 a2 0 1.73203 0 a3 0.0 0.0 10 basis 0 0 0 basis 0.5 0.16666666666666666 0 basis 0.5 0.5 0 basis 0 0.6666666666666666 0
#YL heat flow is in the y direction; z is perpendicular to the single layer graphene; graphene is in the x-y plane
region box block 0 80 0 120 -0.5 0.5
create_box 1 box
create_atoms 1 box
region coldest block INF INF 0 1.5 INF INF
region hottest block INF INF 60 61.5 INF INF
mass 1 12
pair_style tersoff
pair_coeff * * SiC_1994.tersoff C
neighbor 2.0 bin
neigh_modify delay 5
dump snap all cfg 10000 dump.config.*.cfg id type xs ys zs vx vy vz
equilibrium
velocity all create 600 458127641 mom yes rot yes dist gaussian
fix NPT all npt temp 300 300 1 iso 0 0 0.1 drag 1
thermo_style custom step temp etotal vol lx ly lz press pxx pyy pzz
thermo_modify lost warn
thermo 1000
timestep 0.001
run 50000
write_restart restart.NPT50k.equil
reset_timestep 0
unfix NPT
fix temp all temp/berendsen 300 300 0.1
fix nve all nve
variable kb equal 1.3806504e-23
variable ev2J equal 1.60217e-19
compute ke all ke/atom
variable temp atom c_ke*{ev2J}/(1.5*{kb})
fix temp_profile all ave/spatial 1 100000 100000 y lower 1.5 v_temp file temp.profile units lattice #YL for test
#calculate temperature difference***************************
compute hottest_temp all temp/region hottest
compute coldest_temp all temp/region coldest
variable delta_temp equal c_hottest_temp-c_coldest_temp
fix delta_out all ave/time 1 100000 100000 v_delta_temp file delta_temp.dat
fix hottesttemp_out all ave/time 1 100000 100000 c_hottest_temp file hottesttemp.dat
fix coldesttemp_out all ave/time 1 100000 100000 c_coldest_temp file coldesttemp.dat
run 100001
unfix temp
fix heat_swap all thermal/conductivity 10 y 80
fix e_exchange all ave/time 10 1000 100000 f_heat_swap file e_exchange.dat
variable ps2s equal 1.0e-12
variable An2m equal 1.0e-10
variable accumu_thermalK equal f_e_exchange*{ev2J}/(4*0.001*10*{ps2s}196.7624.595*${An2m}*f_delta_out)*511.191
#0.001 is the timestep
#5 is the swap frequency
#196.76 is the width of the graphene layer
#1.42 is the thickness of the graphene layer
#511.191 is the length of the graphene layer (also the heat flux direction)
fix accumu_thermalK_out all ave/time 100000 1 100000 v_accumu_thermalK file thermal_conductivity.dat
the thermal conductivity calculated on step i should be
accumu_thermalK_out (i) minus accumu_thermalK_out (i-1)
run 2000000
Thank you very much.
Best Regards,
Leon