Thermal Conductivity Questions

Recently, I have been investigating the thermal conductivity methods that LAMMPS offers. From what I have found, there are the Muller-Plathe and Green-Kubo methods which can be used. I am interested in determining the thermal conductivity in Silicon crystal, which I have seen on the lammps mailing list, is not a trivial task.

Thus far, I have only implemented the Muller-Plathe method, as it does not appear to require nearly the same computational time as Green-Kubo. However, I first noticed that the thermal conductivities were fairly small. After reading some more, I found a LAMMPS workshop ppt that talked about the effect of the system size on the results, which sure enough, increases the thermal conductivity. Then they showed how one can extrapolate to infinite length in an attempt to calculate thermal conductivity in bulk materials.

Can anyone offer advice on the validity of this method? I also notice that SW potential gives higher thermal conductivities than tersoff (guessing that is just due to differences in the force field, as no force field is perfect).

Another question: How linear can I really expect the temperature gradient to be? I attached a figure that shows the temperature gradient for a 4x4x96 (lattice units) simulation at 500 K. Unfortunately, despite following the exact method as listed in the PPT, the thermal conductivity I get is about 12 W/mK as opposed to 19.5 W/mK listed at the workshop ppt website. My N value in the fix thermal_conductivity is 50 (looks better with 50 than 100)

(http://lammps.sandia.gov/workshops/Feb10/Dave_Rigby/MD_LammpsWorkshop_Feb2010.pdf)

I am not sure if this could be an artifact of the temperature gradient. Is this best way to measure the gradient to take the point at N/2 and N and find the slope? Clearly, the results are not linear. I only ran the simulation for thermal conductivity for 200 ps (last 100 ps is where I collected data). I know this is not a very long time.

Basically, I am just looking for people who have investigated this sort of thing before and any advice they have. I want to reproduce the thermal conductivity of silicon with tersoff FF, but at 300K, the thermal conductivity that I find is ~3 W/mK with Tersoff and 5x5x200 simulation.

The version of LAMMPS I am using is fairly recent (30 July 2016), so I imagine I am on the right track with everything (not getting any error messages), just finding the thermal conductivity to be much lower than I was expecting.

Thanks,

Jackson

temp_grad.PNG

If I may, I’d like to add to my question.
Using the tersoff Force field to calculate the thermal conductivity of Silicon, I keep getting ~9-12 W/mK which seems way low.

If anyone has experience with this, maybe someone can check my method here and point out if there is a glaring error.

I use the following equation to calculate thermal conductivity:

Thermal Conductivity = Heat Flux / (2 * cross-sectional area * Temperature Gradient)

The heat flux is obtained by the output in thermal conductivity command which outputs a scalar value of the heat transferred. For one case, I get 1286 eV. I then convert this to Joules to get 2e-16 J. Then I divide by the total simulation time which is 200ps (first convert to seconds), then I get a heat flux of 1.64e-6 W.

The “2” in the equation is for the PBCs.

The cross sectional area is just x dimension * y dimension (since thermal conductivity is tested in z direction). This is around 21.7 Angstroms in both dimensions, which gives total area of 474 Angstoms ^2 or 4.74e-18 m^2.

The last thing is just the temperature gradient which I get from my curve which is just the slope. So the value would be in K/Angstrom. I get about 0.95 K/Angstorm which is then 9.5e9 K/m.

Plugging in everything, I get a thermal conductivity of 12.7 W/mK.

I do not expect anyone to check my math here. But is there any number that looks obviously off or is there some part of my method that is incorrect? The temp gradient is obtained from the last 100ps of simulation time, and the first 100ps is just to “warm up.”

I use 200ps for total simulation time because this is the time in which the heat transferred is accumulated.

I should be seeing thermal conductivities closer to 150 W/mK. I even tested the effect on system size and I do not see that much difference. I equilibrated in NPT ensemble for 1ns, then NVT for 1ns, before switching to NVE and thermal conductivity calculation.

Thanks,

Jackson

You could also try one of the other thermal cond methods

in LAMMPS. Section 6.20 highlights them and there

are example scripts in examples/KAPPA. See its README

for details.

Steve

Hi Jackson,

Some answers to your questions below (first and second mails).

best,
Arthur

Materials Design S.A.R.L.

If I may, I’d like to add to my question.
Using the tersoff Force field to calculate the thermal conductivity of Silicon, I keep getting ~9-12 W/mK which seems way low.

If anyone has experience with this, maybe someone can check my method here and point out if there is a glaring error.

I use the following equation to calculate thermal conductivity:

Thermal Conductivity = Heat Flux / (2 * cross-sectional area * Temperature Gradient)

The heat flux is obtained by the output in thermal conductivity command which outputs a scalar value of the heat transferred. For one case, I get 1286 eV. I then convert this to Joules to get 2e-16 J. Then I divide by the total simulation time which is 200ps (first convert to seconds), then I get a heat flux of 1.64e-6 W.

The “2” in the equation is for the PBCs.

The cross sectional area is just x dimension * y dimension (since thermal conductivity is tested in z direction). This is around 21.7 Angstroms in both dimensions, which gives total area of 474 Angstoms ^2 or 4.74e-18 m^2.

The last thing is just the temperature gradient which I get from my curve which is just the slope. So the value would be in K/Angstrom. I get about 0.95 K/Angstorm which is then 9.5e9 K/m.

Plugging in everything, I get a thermal conductivity of 12.7 W/mK.

I do not expect anyone to check my math here. But is there any number that looks obviously off or is there some part of my method that is incorrect? The temp gradient is obtained from the last 100ps of simulation time, and the first 100ps is just to "warm up.”

I use 200ps for total simulation time because this is the time in which the heat transferred is accumulated.

Your procedure seems to be OK, except for the sampling time of the temperature gradient. 100 ps seems really short - did you try increasing that value to say, 2 ns? Are you seeing any changes?

As a very simple and very crude rule of thumbs: the larger the system and its (expected) thermal conductivity, the longer the temperature gradient will take to converge.

I should be seeing thermal conductivities closer to 150 W/mK. I even tested the effect on system size and I do not see that much difference. I equilibrated in NPT ensemble for 1ns, then NVT for 1ns, before switching to NVE and thermal conductivity calculation.

With silicon, you should see a strong system size effect. The value you are getting does not seem that surprising. Try doubling and quadrupling the system size, if possible. I’ll elaborate below concerning size effects.

Thanks,

Jackson

At the end of the day, equilibrium MD (GK) and non-equilibrium MD (MP) require roughly the same computational time.

Size effects can play a significant role depending on the material and the thermodynamic conditions. The finite size of your system directly truncates phonons, according to their wavelength. Increasing the size of your simulation box allows you to represent phonons with larger wavelength. When these long-wavelength modes actively participate in heat transport, this is where you see size effects appearing.

People usually say that size effects are stronger in systems with high thermal conductivity - this is probably not a good correlation. Size effects are stronger in systems where “long-wavelength modes participate in heat transport”, which can be approximated as “long-wavelength propagating modes”. This is why size effects are strong in amorphous Si and weak in amorphous SiO2, both low thermal conductivity materials. See recent work from Larkin and McGaughey: PhysRevB.89.144303.

Inverse extrapolation to infinity is a good solution; however, 1/kappa = f(1/L) is not always linear (see publications from Pawel Keblinski). There is therefore no perfect method allowing to account exactly for size effects in an arbitrary situation.

Here is a graph I’ve adapted from PhysRevB.65.144306, showing 1/kappa = f(1/L) for Si and diamond, for different temperatures, using the SW potential. I strongly recommend reading that paper.

You can see that roughly at your system size (50nm, <=> 2e-2 1/nm), they are obtaining ~ 20 W/m/K at 500 K.

figure3.pdf (35.6 KB)

Should the temperature gradient be virtually symmetric on both sides of the plot? I find that my value at 0 is about 30K higher than at N. If the heat source is at N/2, I would expect them to be about the same.

Jackson

I think it should be roughly symmetric. What is your setup here? Is the heat sink spread between bins 0 and N, and the heat source at the middle of the system?

Arthur

yes, it is. I am a little confused on the thermal conductivity calculation as well.

If I take the temperature gradient between points 1 and 5 and calculate the thermal conductivity. Should not the result be the same as if I take the temperature gradient between point 1 and 2, 2 and 3, 3 and 4, and 4 and 5, calculate the thermal conductivity for each region, and average? These two different methods give different results, though it would seem they should not.

Jackson