Dear all,
I’m trying to calculate the thermal conductivity of a periodic graphene surface using the Green-Kubo strategy with the input card mentioned in the LAMMPS documentation:
https://docs.lammps.org/compute_heat_flux.html
but with some modification to adapt it to my problem:
# LAMMPS input script for thermal conductivity of a periodic graphene sheet
units metal
variable seed equal 24526
variable fileName string 5x5armchair_c.dat
variable T equal 100
variable V equal vol
variable P equal 1.01325 # 1 atm in bar
variable dt equal 0.001 # 1 fs in metal
variable r equal 1000000
variable p equal 200 # correlation length
variable s equal 10 # sample interval
variable d equal $p*$s # dump interval
variable kB equal 1.3806504e-23 # [J/K] Boltzmann
# convert from LAMMPS metal units to SI
variable eV2J equal 96485.31/6.02214e23
variable ps2s equal 1.0e-12
variable A2m equal 1.0e-10
variable convert equal ${eV2J}*${eV2J}/${ps2s}/${A2m}
# setup problem
boundary p p p
atom_style atomic
read_data ${fileName}
mass 1 12.0107 # C
pair_style airebo 3.0
pair_coeff * * ../../../potentials/CH.airebo C # C-C & C-H
neighbor 3 bin
neigh_modify delay 10
# minimization
thermo 50
dump 1 all custom 200 min.lammpstrj id element x y z
dump_modify 1 element C
minimize 1.0e-9 1.0e-10 2000 20000
reset_timestep 0
undump 1
# equilibration and thermalization
timestep ${dt}
thermo $d
velocity all create $T ${seed} mom yes rot yes dist gaussian
fix NPT all npt temp $T $T 100 iso $P $P $d
run 300000
# thermal conductivity calculation, switch to NVE if desired
reset_timestep 0
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 ave running
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}
thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33
dump 1 all custom $d md.lammpstrj id element x y z
dump_modify 1 element C
run $r
variable k equal (v_k11+v_k22+v_k33)/3.0
variable ndens equal count(all)/vol
print "average conductivity: $k[W/mK] @ $T K, ${ndens} /A\^3"
However, I don’t understand why the final value for kappa was 2.47251 [W/mK], while we are waiting for a value higher than this.
We also used this input card to calculate the kappa of a graphene flake into an Ar box and, in this case, we got the expected value, that was 61.9 [W/mK]
for the case of the graphene sheet, I also tested a NVT ensemble, but the result didn’t change. What do you think?
Here are the both .dat, if you want to test it.
graphene sheet:
5x5armchair_c.dat (34.7 KB)
graphene flake:
G586Ar.dat (103.0 KB)
Thnanks a lot!