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

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
  • bad geometry
  • 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.

Please also see About the estimation of the error in the calculation of thermal conductivity by the Green-Kubo method - #2 by Germain

Thanks for the reply.

I’m not sure I understand what bad geometry means. Is it an error in the data file? This code is a modified version of the one provided in src/PACKAGES/interlayer. with the Green-Kubo algorithm appended.

It can mean one or more of multiple things:

  • inconsistent box dimensions
  • incorrect topology information for the chosen force field or pair style
  • unrelaxed positions with respect to the chosen pair style
  • inconsistent/incompatible units when writing and reading the data file

This should usually result in unexpected behavior (high potential energy, too high/low pressure, buckling) during equilibration.

Please note that LAMMPS examples are primarily meant for demonstrating syntax for specific LAMMPS features, they are not necessarily meant to demonstrate best practices in how to set up simulations in general. When you re-purpose and example input, you need to check it as if you would construct a system yourself from scratch.