Thermal Conductivity of a bulk metal through Green-Kubo

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!