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

read_data graphene.data

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