[lammps-users] silicon thermal conductivity calculation using Muller-Plathe

Hi all,

I've encountered some issues using fix thermal/conductivity to
calculate the thermal conductivity of silicon. I keep getting a K of
~12 W/K-m when it should be ~150. I searched through the lammps
archives and I found someone else who encountered this problem but he
was using the SW potential and was having units issues (which I don't
think I have). Below is a summary of my input file and my results. I
also attached my complete input file and three images: a plot of the
heat flux v. time, raw temp gradient, and fitted temp gradient data.

General approach: NVT equilibrate @ 300K, NVE to check energy
conservation, 1 ns run to set up temp gradient, another ns run to
capture good data
Box size: 3x3x400 a (w/ a being 5.43 A (I'm using metal units)); this
means my box is 200 nm so I ran a 1 ns equilibration for the time it
takes heat to diffuse across the box (t=L^2/D).
Potential: tersoff using constants from the lammps examples

Thermal gradient: 2.49e9 K/m
Heat flow: 1.58e-7 W
Cross-sectional area: 2.653e-18 m2
Thermal Conductivity: 12.15 W/K-m (heat flow/(2*Area*Thermal Grad))

Thanks in advance for your help,
Zack Cordero

Input file:
#Exploration of different temp fixes

boundary p p p
units metal
echo screen
newton on
log log.si_temp

lattice diamond 5.43
region all prism 0.0 3 0.0 3 0.0 400 0.0 0.0 0.0
create_box 1 all
create_atoms 1 region all
pair_style tersoff
pair_coeff * * Si.tersoff Si
neighbor 2.0 bin
neigh_modify every 1
mass 1 28.0

thermo 100

#nvt equillibration
velocity all create 400 429349 dist gaussian
fix 1 all nvt temp 300 300 .1

timestep .001
run 5000

#nve equillibration
unfix 1
fix 2 all nve
run 5000

#muller-plathe thermal conductivity

#Temp gradient
compute ke all ke/atom
variable temp atom c_ke/(1.5*8.617*10^-5)
fix 3 all ave/spatial 10 100 1000 z lower .02 v_temp file
tmp.profile units reduced

fix 4 all thermal/conductivity 100 z 50

thermo_style custom step temp etotal dt elapsed f_4
log log.warm_up

run 1000000

log log.si_flux
unfix 3
fix 5 all ave/spatial 10 100 1000 z lower .02 v_temp file
tmp.profile units reduced

run 1000000


tempgrad 1.png


Dear Zach,

Did you check what is the average pressure in your calculations? The thermal conductivity may change depending upon if the average value of atomic separation is closer (vice-versa) than their equilibrium value. I am asking as I did not see you running NPT simulation.

Best Regards,

I just looked at the pressure in the NVE part of my run and saw that
it was oscillating around -1600 \pm 100 bars (pretty big
oscillations). So some questions:

This is a huge pressure right?
Why would the system have a negative pressure if I'm heating it up and
I input the correct lattice constant?
Should I use an additional NPT aniso stage w/ 0 external pressure to
let the system relax?
Wouldn't this let the system shrink, decreasing average atomic
separation, thereby increasing the thermal conductivity (like I want)?


You are in right direction. :slight_smile:
Now, Lattice constant is experimental. The potential you are using may not be 100 % agree with the experiment. So, let it relax with with aniso condition.
I have one more suggestion. In order to avoid significant pressure buildup on your thermal source (due to difference of 200 K) between source and sink temperatures, try swapping less number of atoms and less frequently.

Best Regards,

Dear Zachary,
I have the same problem with Si with Tersoff potential. However it was OK with SW and MEAM potentials.
I checked it both MP and Green-Kubo methods. The only reasonable result for SiC with Tersoff potential was
published in work by L Porter, S Yip. But the Tersoff potential was modified compared to what is used in LAMMPS.
German Samolyuk

Hi all,

Since it seems like a lot of people do (or attempt to do) this type of
calculation, I figured I'd keep you all posted on what I've learned so
far. I tried using both an NPT ensemble after my NVE stage and I
tried using an energy min stage at the very start of my run. Neither
change helped very much; I still found a K of around ~12. So I
searched the web and found 2 helpful resources... According to them,
you should see how your thermal conductivity changes as you lengthen
the dimension with the temp gradient along it and then extrapolate to
infinite length.

German, I believe Tuomo recommended this approach to you in a previous
e-mail thread.

They talk about extrapolation here:

and here:


Hi Zack,
Do I understand it right, k=12 is with Trsoff potentials? Did you try MEAM or SW?
It will be interesting to know what your results for this two potentials.

I also think that with dimension of ~ 200 nm, the value of K ~ 12 could be force field issue.