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”