# Thermal conductivity calculated by Green-Kubo method

Dear Lammps users,
I am trying to use lammps and Green-Kubo method to study the
thermal transport in thin films. I started with the bulk silicon. I
obtained an average thermal conductivity of around 60 W/m-K at 1000 K
using the SW potentials and a 6*6*6 simulation cell with periodic
boundary conditions, as reported in the literature. The repeatability
of the calculation is very good. I ran 5 simulations and the deviation
was less than 10%. However, when I checked the thermal conductivity
values along the x,y,z directions, I found k11, k22, and k33 were not
equal. Although their average was always around 60 W/m-K, the three
values along different polarizations could hardly reach the same
value, and were not even close. For example, k11, k22, and k33 can be
40, 75 and 68, respectively. I ran several simulations with the same
sampling length and correlation length. It seems that in different
simulations, the magnitudes of the three thermal conductivity
components vary randomly in the range between 30-90 W/m-K while
maintaining almost the same average thermal conductivity. I thought
that the three components should be the same because silicon has a
cubic unit cell. I also ran the example of Green Kubo method (solid
Ar) shown on lammps website and found the similar phenomena. I am
wondering why this phenomena happens. If I want to study the
anisotropies of phonon transport, what can I do to reduce the
uncertainty?

Thanks,

Baoling

This is a good Q - hopefully someone else
will comment. Possibly Reese Jones, who has
a lot of experience with kappa calculations.

Steve

Baoling, I think I’ve seen this also, but I don’t have experience with solids. You can try sampling the autocorrleation more times (for each componet) or adding atoms.

How many atoms do you have?

Tim

Tim,
Thanks. I have used a 6*6*6 simulation cell and each unit cell has
8 atoms. So, there are 1728 atoms in total in the simulation domain.

Baoling

Hi:

I have done a similar calculation for amorphous silica recently with two
different potential parameters. The variation in k11 to K33 was from 2.3
to 1.3 for one set and 1.3 to 1.6 for another potential.

I did notice that the way the integration is carried out in fix
ave/correlate does not work for systems like silica for which the
autocorrelation does not vary monotonically like for example systems
like Argon.

Deena

Baoling, its hard to say if your N is enough, but at least sounds reasonable. I usually compute the HF autocorrelation several thousand times, which I think is standard practice.

Deena- did you mean that “trap” isn’t always summing the variable correctly for the integration?

Tim

Tim:

When you have a highly oscillating autocorrelation "trap" summation is
not working correctly.

In trap I understand that lammps uses trapezoidal rule. In that case I
don't quite understand why division by 2 is conditional during the
summation.

Deena