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.84*10^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.76*24.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