# Calculated graphene thermal conductivity lower than expected

Hi all,

I am trying to simulate the thermal conductivity of a square graphene sheet (~3000 atoms), but my results are an order of magnitude lower than the expected values (~150W/mK vs 4000W/mK). Does anyone have any insights as to why this is happening? 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 p
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

There are lots of possible causes:

• typos in the input or mistakes in converting values
• too short a simulation to get converged data
• inconsistent use of force field settings
• unsuitable reference data (you need to compare to simulation data using the exact same setup, settings, and force field).

Neither of those are LAMMPS specific.