# Graphene lattice conductivity

Dear all,
I am trying to calculate the lattice conductivity using LAMMPS of some BNC 2D models.
To start understanding how the calculation goes, I’m trying to model pristine graphene using the potential reported in this paper: Phys. Rev. B 86, 115410 (2012) - Thermal conductivity of BN-C nanostructures
that says that they manage to reproduce with good accuracy the lattice thermal conductivity of graphene, that should be around 1000W/K.

However, when I try to calculate it, I have values that are way lower (around 150W/Ks). This value is quite close to the one of h-BN which makes me think that this is due to the lack of electron contribution in the LAMMPS calculation. So my question is: since in the paper they claim that they reproduce the experimental value, is my hypothesis correct or something is wrong with the calculation itself?

I report here my input file:

#####################################################

# Variables

#####################################################

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 #2D material
units metal
atom_style full
mass 1 12.0107 # Carbon
pair_style tersoff
pair_coeff * * BNC.tersoff C # chemical
minimize 1e-10 1e-10 10000 100000000

#####################################################

# 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/vol
variable Jy equal c_flux/vol
variable Jz equal c_flux/vol
fix JJ all ave/correlate \$s \$p d & c_flux c_flux c_flux 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){scale} variable k22 equal trap(f_JJ)*{scale}
variable k33 equal trap(f_JJ)
\${scale}
variable k11r equal trap(f_JJ)
variable k22r equal trap(f_JJ)
variable k33r equal trap(f_JJ)
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 100000 #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”

I have not looked at it in too much detail, but here are some stuff I noticed:

1. by `variable k equal (v_k11+v_k22+v_k33)/3` you are seemingly averaging the thermo conductivity of three directions. This is usually good for isotropic materials, but you probably don’t want to do this on 2D or layered materials (where the conductivity heavily depends on the direction).
2. make sure that the total simulation time (100ps in your script) and correlation time (10ps in your script) are long enough. Both should be long enough to get the thermo conductivity converged. Sometimes this can be tricky and I suggest plot the correlation function to check this.
3. the paper you cited used another way to calculate thermal conductivity, so I’m not sure how close the results could be.
1 Like