LAMMPS calculated graphene thermal conductivity lower than expected

Question

I have been trying to calculate the thermal conductivity of a square graphene sheet of 3200 atoms, but I get a k value of 150W/mK, which is many times smaller than the published 4000W/mK. I am fairly new to LAMMPS and I'm not sure where I'm going wrong since it isn't a coding error. Does anyone have any insights? Thanks in advance

#SLG GK
#####################################################
# Variable Definition ###############################
#####################################################
variable        damp equal 0.5
variable        seed equal 102486 #just a random number
variable        T equal 300 # Temperature
variable        V equal vol #volume
variable        P equal 1.01325 # 1 atm in bar 
variable        dt equal 0.001 # 1 fs in metal (timestep)
variable        r equal 100000 #??
variable        p equal 10000    # correlation length should be 10000
variable        s equal 1        # sample interval should be 1
variable        d equal $p*$s    # dump interval
variable        kB equal 1.3806504e-23    # [J/K] Boltzmann constant

# Convert from LAMMPS metal units to SI
variable        eV2J equal  1.60218e-19 # eV to Joule
variable        ps2s equal 1.0e-12 # picosec to sec
variable        A2m equal 1.0e-10 #angstrom to m
variable        convert equal ${eV2J}*${eV2J}/${ps2s}/${A2m} #

#####################################################
# Initializing simulation box #######################
#####################################################
dimension       3
boundary        p p f
units           metal
atom_style      full
read_data       SLG2.data
mass            1 12.0107   # Carbon
pair_style      airebo 2 
pair_coeff      * *     CH.airebo C # chemical
neighbor        2.0 bin
neigh_modify    every 1
neigh_modify    delay 0
neigh_modify    check yes
minimize 1e-10 1e-10  10000 100000000
#delete_atoms overlap 0.1 all all

#####################################################
# 1st equilibration run #############################
######################################################
#At the beginning of each simulation, a system is optimized and then equilibrated
#with NPT(isothermal-isobaric) ensemble under 300K temperature and zero pressure.
#All of the simulations are performed by time step of 1fs with periodic boundary condition(PBC)
#NPT 2.5ns 2500st, NVT 1ns 1000st, NVE 1ns 1000st
timestep        ${dt}
thermo          100
velocity        all create 300 12345 mom yes rot yes dist gaussian
fix             1 all npt temp 300.0 300.0 0.2 iso 0.0 0.0 0.1
thermo_style    custom step temp press
run             2500
unfix           1
fix             1 all nvt temp 300.0 300.0 0.2
thermo_style    custom step temp press
run             1000
unfix           1
fix             1 all nve
thermo_style    custom step temp press
run             1000
#unfix          1
reset_timestep  0
#####################################################
# Green Kubo Method #################################
#####################################################

thermo          $d
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]/vol
variable        Jy equal c_flux[2]/vol
variable        Jz equal c_flux[3]/vol
fix             JJ all ave/correlate $s $p $d & 
                c_flux[1] c_flux[2] c_flux[3] type auto file J0Jt.dat ave running #Nevery Nrepeat Nfreq
variable        scale equal ${convert}/${kB}/$T/$T/$V*$s*${dt}
variable        k11 equal trap(f_JJ[3])*${scale}
variable        k22 equal trap(f_JJ[4])*${scale}
variable        k33 equal trap(f_JJ[5])*${scale}
variable        k11r equal trap(f_JJ[3])
variable        k22r equal trap(f_JJ[4])
variable        k33r equal trap(f_JJ[5])
thermo_style    custom step temp v_V v_convert v_k11r v_k22r v_k33r v_scale v_Jx v_Jy v_Jz v_k11 v_k22 v_k33 
dump myDump1 all atom 100 SLG-GK.lammpstrj
run             50000 #should be 100000

variable    k equal (v_k11+v_k22+v_k33)/3
variable    ndens equal count(all)/vol
print           "average conductivity: $k[W/mK] @ $T K, ${ndens} /A\^3"

print       "SIMULATION DONE"



#heatflux = W/m2=J/sm2=x eV/psA2

#V/kBT2= m3K/Js = A3K/eVK2

#conversion factor = A2m / ps2s*K

P.S. I am also trying to simulate the thermal conductivity of stacked graphene on hBN

Answer 1

Your result is probably correct for this model.

I'm actually not surprised that the thermal conductivity is much too low for graphene. The Airebo potential, being a classical atomistic potential, doesn't represent electronic degrees of freedom, so any contribution of conduction electrons to the thermal conductivity is not represented. I'd expect then that the thermal conductivity of this model of graphene would be the same order of magnitude as the conductivity of its insulator sibling, hexagonal boron nitride. Indeed, hexagonal boron nitride has a thermal conductivity of a few hundreds of W/(m K).

Figure 2 of this paper gives values similar to yours for the Airebo potential: