Questions regarding Phys. Rev. B 90, 134312 (2014) - "Role of anharmonic phonon scattering in the spectrally decomposed thermal conductance at planar interfaces"

[thejashc]

Q1: I am interested to know how the interatomic force constants in Equation 10 and 11 were calculated. Were they calculated analytically or numerically? Also, given the LJ interactions are a spherically symmetric potential I would assume that the kij in Equation 10 is a diagonal matrix?
What about gamma in Equation 11?

In sum, calculating the harmonic forces between atoms i and j comes down to calculating their position per timestep, and calculating only once the interatomic force constants (kij and gammaij).

Q2: Is the F_i\alpha term in Equation 10 different from the interatomic force Fij? Is F_i\alpha the sum of F_ij over all j ?

Q3: The framework mentioned in this article should work as long as the forces between dissimilar materials are two-body interactions, correct?
In the article APPLIED PHYSICS LETTERS 103, 051602 (2013) the Si-Ge are modelled using SW potentials which are intrinsically three-body potentials. The framework like you mentioned is perhaps not applicable to potentials which are not two-body.

Sääskilahti et al. - 2014 - Role of anharmonic phonon scattering in the spectrally decomposed thermal conductance at planar inte.pdf

[ksaaskil]

Hi Thejas, thank you for your interest in the paper! It's been more than ten years since I wrote it with my colleagues, so I'll try my best to answer your questions.

Q1: I am interested to know how the interatomic force constants in Equation 10 and 11 were calculated. Were they calculated analytically or numerically? Also, given the LJ interactions are a spherically symmetric potential I would assume that the kij in Equation 10 is a diagonal matrix? What about gamma in Equation 11?

They were calculated numerically using the code from this repository similar to fc.py. I cannot immediately say if the LJ interaction matrix is diagonal in the coordinate space, could be.

Q2: Is the F_i\alpha term in Equation 10 different from the interatomic force Fij? Is F_i\alpha the sum of F_ij over all j ?

Yes and yes.

Q3: The framework mentioned in this article should work as long as the forces between dissimilar materials are two-body interactions, correct? In the article APPLIED PHYSICS LETTERS 103, 051602 (2013) the Si-Ge are modelled using SW potentials which are intrinsically three-body potentials. The framework like you mentioned is perhaps not applicable to potentials which are not two-body.

The framework was expanded to many-body potentials in the paper https://journals.aps.org/prb/abstract/10.1103/PhysRevB.91.115426

Also see https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.094301

[thejashc]

Many thanks, Kimmo. I am trying the calculations and will come back to this thread.
In the meanwhile, I do have a follow-up question.

They were calculated numerically using the code from this repository similar to fc.py. I cannot immediately say if the LJ interaction matrix is diagonal in the coordinate space, could be.

Would you happen to remember/know if computing the force constants using the code provided would give you the same result as the second partial derivative of LJ potential ? I would expect that the two routes are identical.

Yes and yes.
The framework was expanded to many-body potentials in the paper https://journals.aps.org/prb/abstract/10.1103/PhysRevB.91.115426

Also see https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.094301

Thanks for clarifying.

[ksaaskil]

Many thanks, Kimmo. I am trying the calculations and will come back to this thread. In the meanwhile, I do have a follow-up question.

They were calculated numerically using the code from this repository similar to fc.py. I cannot immediately say if the LJ interaction matrix is diagonal in the coordinate space, could be.

Would you happen to remember/know if computing the force constants using the code provided would give you the same result as the second partial derivative of LJ potential ? I would expect that the two routes are identical.

Yes, assuming you evaluate the second derivative at the average position of the atoms.

You may also want to take a look at my paper from Phys. Rev. E where we calculated the spectral heat current decomposition for a solid-liquid interface. Correct approach depends on the problem you're trying to solve.

Yes and yes.
The framework was expanded to many-body potentials in the paper https://journals.aps.org/prb/abstract/10.1103/PhysRevB.91.115426
Also see https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.094301

Thanks for clarifying.

[thejashc]

Yes, assuming you evaluate the second derivative at the average position of the atoms.

I can imagine that the analytical route can get messy very quickly if the interfacial region is large and interaction cutoff is also large, since it requires summing over all j's. This would ultimately reduce to a computer calculation. But nonetheless, thanks for the input.

You may also want to take a look at my paper from Phys. Rev. E where we calculated the spectral heat current decomposition for a solid-liquid interface. Correct approach depends on the problem you're trying to solve.

Indeed, I am familiar with your PRE article since I am looking at solid-liquid interfaces. But, your PRE article focusses on a full sdhc calculation. I am particularly interested in applying your (in)elastic decomposition of sdhc for solid-liquid interfaces.

[thejashc]

Thanks for sharing your scripts, Kimmo.

I’ve now tried to reproduce one case from your PRB paper - the 30 K Ar/heavy-Ar interface. I adapted your workflow into C++ to integrate it with my existing setup.

I’m quite close to your results (features are quite similar), but I still see some systematic differences in the elastic spectra. However, I’ve verified that the total conductance at T = 30 K agrees with your reported value to within 1%.

One possible source of the discrepancy could be the interatomic force constants. I extracted the Kij values from a single restart file using the finite-displacement approach, with a displacement of 0.001 Å. I suspect the Kij values may depend on the configuration used to compute them.

From your scripts, I had the impression that the force constants were computed from a single configuration. Just to confirm, did you average the Kij over multiple configurations, or were they obtained from a single snapshot?

[ksaaskil]

Hi, the harmonic approximation is made around the average position of the atoms. So we calculated the average positions from the simulation, confirmed the displacements were small and then calculated force constants at the average positions using the finite displacement method.

You can also try and compare to the Green's function method at low temperature and see you get the same result.

[thejashc]

I've now calculated the KIj's from the averaged atom positions. I've also calculated the spectra from the analytical Kij's. The analytical Kij's are identical to the Kij's from the FDM technique. I've also tried to verify that using just plain compute forces from lammps is identical to force/tally compute which calculates forces between specific groups of atoms. They too are identical.

Now, the agreement with Fig. 4a of your manuscript is much better, though some minor discrepancies exist. I am not sure why these discrepancies still exist.

Like you mentioned, its perhaps best to do some calculations at 1 K where the spectra can be compared with the Green's Function approach.

[thejashc]

Just to be doubly sure, your manuscript mentions that the fcc lattice constant changes between 1.5496 sigma_LJ at T=0 K to 1.58 sigma_LJ at T=40 K. When constructing the geometry for T=30 K, I assumed it can be linearly interpolated between the sigma's at 0 and 40 K cases. If the lattice constant is slightly different, I would expect the final spectra would also be slightly different.

[thejashc]

The 1K case compares very well with your manuscript. I am not sure what happens at the 30 K case.
But this is good enough, I think.

Thanks a lot for your help, Kimmo.