# 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
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!

Looking at your system I’am not sure the command:

``````fix          NPT all npt temp \$T \$T 100 iso \$P \$P \$d
``````

Is really adapted to equilibrate your system. Did you check the equilibrated average temperature and pressure tensor? I suggest you (re)read the `iso` keyword and `Pdamp` and `Tdamp` parameters entries of the manual.

Also when computing k, do you really need to take the 3 dimensions and the full volume of the simulation box into account for an isolated graphene sheet?

Why do we multiply sample interval (s) in this line?