# Calculation of thermal conductivity

Hello.
I try to calculate thermal conductivity using two methods: Green-Kubo and Muller-Plathe. For the same systems I have two different values for thermal conductivity. (Muller-Plathe: 25 W/mK, Green-Kubo: 0.25 W/mK). Can you help me to find mistakes in my scripts?
My script for Muller-Plathe is:

units metal
boundary p p p
atom_style atomic
lattice diamond 5.43
region main block -3 3 -3 3 0 6
create_box 1 main
region myfirst block -3 3 -3 3 0 6
create_atoms 1 region myfirst
group atom1 region myfirst
set group atom1 type 1
pair_style meam
pair_coeff * * library.meam Si92 NULL Si92
timestep 0.0001
#equilibration
velocity all create 500.0 4928459 dist gaussian
fix 1 all nvt temp 500 500 0.002
thermo 1
run 10000
#Muller-Plathe conductivity calculation
compute 1 all ke/atom
variable temp atom c_1/(1.5*8.617332478e-5)
fix 2 all ave/spatial 100 100 10000 z lower 6 v_temp file tempprofile.txt units box
fix tc all thermal/conductivity 100 z 6
fix 3 all ave/time 100 100 10000 f_tc file thermal_conductivity.txt
thermo_style custom step temp etotal ke vol
run 1000000

Thermal conductivity is 25 W/m*K .

My script for Green-Kubo (manual.pdf p.287) is:
# Sample LAMMPS input script for thermal conductivity of solid Ar

units real
variable T equal 500
variable V equal vol
variable dt equal 1.0
variable p equal 200 # correlation length
variable s equal 10 # sample interval
variable d equal \$p*\$s # dump interval
# convert from LAMMPS real units to SI

variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable kCal2J equal 4186.0/6.02214e23
variable A2m equal 1.0e-10
variable fs2s equal 1.0e-15
variable convert equal \{kCal2J\}\*{kCal2J}/\{fs2s\}/{A2m}

# setup problem

dimension 3
boundary p p p
lattice diamond 5.43
region box block -3 3 -3 3 0 6
create_box 1 box
create_atoms 1 box
pair_style meam
pair_coeff * * library.meam Si92 NULL Si92
timestep \${dt}
thermo \$d
# equilibration and thermalization
velocity all create \$T 102486 mom yes rot yes dist gaussian
fix NVT all nvt temp \$T \$T 10 drag 0.2
run 100000

# thermal conductivity calculation, switch to NVE if desired
#unfix NVT
#fix NVE all nve

reset_timestep 0
compute myKE all ke/atom
compute myPE all pe/atom
compute myStress all stress/atom 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 &amp; &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;c\_flux\[1\] c\_flux\[2\] c\_flux\[3\] type auto file J0Jt\.dat ave running &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;variable scale equal {convert}/\${kB}/\$T/\$T/\$V*s\*{dt}
variable k11 equal trap(f_JJ[3])*\{scale\} &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;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
run 100000
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"

Thermal conductivity is 0.25 W/m*K .

Sorry - no one is likely to run two complex
input scripts for you and figure out why they give
different answers. One obvious issue to
be careful of is that you are running one in
real units, the other in metal units, so the
output of the two would have to be
transformed carefully to compare them.

Steve