Help with thermal conductivity

Hello Timothy,

Your answer was very helpful, but unfortunately I have not been able to solve my problem completely yet. Hence this email.

So, I started with the Argon example and calculated the ACF by adding the 3 outputs in the ‘J0Jt.dat’ file written after every Nrepeat steps. The ACF decays to zero very quickly and the thermal conductivity value is consistent (0.28 W/mK). I played with the Nrepeat and Nevery parameters and plotted them (attached).

With this background I started to work on Silicon bulk case (6X6X6). Since it is equilibrium MD, I equilibriated the structure and made sure that kinetic energy, potential energy and pressure are minimized at 500 K and used the equilibriated restart file for kappa calculations. Since I am using Stillinger-Waber potential I changed the units according to metal units (for Argon It was done corresponding to the real units).

However, I am still not able to get the ACF decay to zero. The thermal conductivity values are changing with changing Nrepeat and Nevery parameters (figures attached) in the fix ave/correlate command. The ACF profiles do not match with the ones available in literature. The thermal conductivity values are very different than published work. I am not able to figure out the mistake in my approach. For last few days, I have searched the mailing list thoroughly but could not figure out the the problem. I am attaching my input script, log file and the J0Jt.dat file.

Please help me out.



in.txt (1.74 KB)

log.lammps (45.6 KB)




Hi Sankha,

Hi Sankha

Silicon is the second most difficult system for calculating the thermal conductivity (diamond is the most difficult). The issue issue is the very high thermal conductivity arising from the very regular, rigid structure. This translates into very long phonon scattering lengths. In practice this translates into a strong dependency on the system size and the need to run quite long calculations to get long enough correlation times. It also manifests itself in a very large isotope effect, with the thermal conductivity of an isotopically pure system at low temperatures being several times larger than silicon with the natural isotopic mix.

With patience you can get reasonable results, but it is quite some work. I’d suggest that work your way up to silicon and pay attention to the isotopic masses since most experiments are on natural isotopic mixtures. The lower the thermal conductivity the easier the calculations will be, so argon is a very easy case. You might start with germanium (I think Stillinger-Weber covers that) since it has a lower thermal conductivity, or a random alloy of Si-Ge would be even easier place to start since the thermal conductivity is in the range of 1-10 W/mK rather than a few hundred up to about 4000 W/mK for pure silicon, depending on the temperature. The Ioffe Institute has a nice website ( with experimental data for the alloys, etc.

I previously answered a very similar question on the mailing list ( and quoted some values that I had calculated using the version of the ACF that we originally implemented in LAMMPS for use in MedeA. The version of our code that ended up in LAMMPS as fix ave/correlate lost some of the averaging that we had implemented in the ACF, and we also work very hard outside LAMMPS to combine the ACF’s over a large number of time scales into a single, combined ACF for the analysis, which may be needed for the systems with high thermal conductivity.

Silicon Germanium ~Si43Ge (2.3%) Si43Ge5(10.4%) 25%
N Atoms Supercell Lambda (W/m.k) Lambda (W/m.K) Lambda (W/m.K) Lambda (W/m.K) Lambda (W/m.K)

64 2x2x2 40.9 12.7 9.5 (3.1%)
512 4x4x4 102.3 29.3 8.5 (2.4%)
4096 8x8x8 153.2 33.8 10.2 (2.3%) 2.6 (10.4%) 1.9
32768 16x16x16 170 52.4

experiment 130 58

These calculations are for the natural isotopic mix, sampling the ACF for 5 ns using a 4 fs step with our forcefield (pcff+) rather than Stillinger-Weber and result in ACF’s like this for the 8x8x8 cell:

Where the grey envelope is an estimate of the error. The sharp discontinuities are apparently an artifact of the way that we combine multiple ACF’s to produce the overall ACF. If we fit the initial rise to the plateau and ignore the subsequent wandering around:


As to why the calculated results do not match experiment for silicon, I am not sure. It could of course be the forcefield, but I do wonder about the experiment. As I pointed out, even the small differences of the isotopic mass make a large difference to the thermal conductivity because they help scatter the phonons. The mean free path for the phonons before scattering is quite large up to microns or perhaps mm in diamond, so it may be that “perfect” silicon crystals still have enough defects to significantly lower the thermal conductivity from that of a truly perfect crystal. I should try running one of these again with a single vacancy and see the effect is. :slight_smile:

Hopefully this helps you see at least what the target is. The good news is that you are probably not interested in the pure silicon, and anything else, including your nanowire, will be easier!