LAMMPS calculated graphene thermal conductivity lower than expected

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
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