Phonon relaxation times at certain q point

[biglinn]

Hi, @DjordjeDangic
I find there are outputs frequency and lifetime in the Phonon_transport_properties_X after calculating thermal conductivity, like this:

Since there are thousands of lines that confuses me, could you explain the structure of all lines?
Because I'd like to calculate the phonon relaxation time at certain q point(such as q=0).

By the way, I see that the total phonon relaxation time contain not only anharmonic scattering but also other scattering in APPENDIX B of
the article 10.1103/PhysRevB.88.045430.
So dose the lifetime in Phonon_transport_properties_X is total phonon relaxation time?

[DjordjeDangic]

Each line is one phonon mode. It has thousands of lines because you have thousands of phonon modes in your system (number of q points times the number of bands).

If you save the ThermalConductivity object (with tc.save_pickle() function where tc is your ThermalConductivity object) you can load it and then find all of your calculated properties. This is shown in the tutorial that is on sscha website. For example, if you need lifetimes at Gamma (change temperature accordingly):

import numpy as np
import cellconstructor as CC
import cellconstructor.ThermalConductivity

tc = CC.ThermalConductivity.load_thermal_conductivity()

temperature = 400.0
key = format(temperature, '.1f')

for iqpt in range(tc.nkpt):
    if(np.linalg.norm(tc.k_points[iqpt]) == 0.0):
        print('Lifetimes in seconds in k point: ', tc.k_points[iqpt])
        print(tc.lifetimes[key][iqpt])

Lifetimes calculated in SSCHA are just due to anharmonicity (just the first term in that reference).

[biglinn]

It really is!
I set kgrid mesh=[8,8,8] in the CC.ThermalConductivity.ThermalConductivity( , and it has 12 phonon bands. So there are 6144(512*12) lines.

But what is the structure of the file?
Is it like that 1-512 lines belong to one band, 513-1024 lines belong to another...
or that 1-12 lines are 12 bands at one point, 13-24 another point... ?

[DjordjeDangic]

The first 12 lines are the first q point (which is usually Gamma). The second 12 lines are from a second q point, etc.

[biglinn]

And I'm still not sure about one thing.
I calculate thermal conductivity by following the tutorial which calculates by potential(such as ML potential).
The two important steps of process is:
1, calculate initial dynamical matrices of corresponding supercell by CC.Phonons.compute_phonons_finite_displacements
2, sample ensembles and calculate energy and forces by compute_ensemble(

Before, I calculated anharmonic dyn following the tutorial https://sscha.eu/Tutorials/tutorial_02_advanced_submission/.
In this workflow, initial dyn come from QE ph.x. Energy and forces are calculated by QE pw.x.

So, I wonder if I can calculate thermal conductivity like this. Because I will use QE to get initial dyn and calculate energy and forces if I research other materials without potential available.
Maybe it's a stupid question, but I think it's best to ask for it.

[DjordjeDangic]

You do not have to use machine learning potential to calculate force constants. You can do everything with DFT. However, keep in mind that to converge results you will need to calculate a large number of configurations and this will not be fast with DFT.

[biglinn]

I get it. As long as I can get SSCHA dyn to calculate 3rd-order IFCs, it's okay.
Thanks!

[biglinn]

I have another question.
The dyn and 3rd-order IFCs are calculated in the same script in the tutorial https://sscha.eu/Tutorials/tutorial_09_TC/.
If I have had calculated dyn, could I directly calculate the 3rd-order IFCs? As follows:

import sscha, sscha.Ensemble, sscha.SchaMinimizer
import sscha.Utilities
import cellconstructor as CC, cellconstructor.Phonons

dyn = CC.Phonons.Phonons("sscha_dyn", 4)
ensemble = sscha.Ensemble.Ensemble(dyn, 0)

dyn_hessian, d3_tensor = ensemble.get_free_energy_hessian(include_v4 = False, get_full_hessian = True, return_d3 = True)
np.save("d3.npy", d3_tensor)

Moreover, before calculating 3rd-order IFCs, there is update weights in the script of toturial:

What role dose it play? And necessary?

[biglinn]

With q-grid X by X by 1 (such as 881, 12x12x1), all q points are in the plane of Gamma point.
And the M-G-K-M path is in the same plane, because of M (0.5, 0.5, 0) and K (0.3333, 0.3333, 0).

Why doesn't the imaginary frequency show at M-G-K-M path of dispersion spectrum?

[biglinn]

I have checked the dynamical matrix files (the first file, corresponding to Gamma point).
The frequencies of first phonon mode are as follows:

Harmonic dyn by QE:
freq ( 1) = -0.162911 [THz] = -5.434134 [cm-1]

SSCHA dyn after one iteration, sscha_dyn_population1_1:
freq ( 1) = -0.00000008 [THz] = -0.00000277 [cm-1]

SSCHA dyn after two iterations, sscha_dyn_population2_1:
freq ( 1) = -0.00000003 [THz] = -0.00000114 [cm-1]

SSCHA dyn after three iterations, sscha_dyn_population3_1:
freq ( 1) = -0.00000011 [THz] = -0.00000358 [cm-1]

SSCHA dyn after four iterations, sscha_dyn_population4_1:
freq ( 1) = 0.00000004 [THz] = 0.00000146 [cm-1]

The imaginary frequency dose exist, just not showed by phonon dispersion.

Although the first mode at Gamma doesn't have imaginary frequency after four iterations, I guess that imaginary frequency may not be eliminated at other area of Brillouin zone. So, perhaps I need more iterations and configurations to get converged SSCHA dyn.

[DjordjeDangic]

As I have told you in one of the previous answers, the best way to check for the imaginary frequencies is to plot the phonon density of states calculated with the same grid you are attempting the thermal conductivity calculation.

Also make sure that the scattering grid is also with only one point in the third direction, ie. 8x8x1.

[biglinn]

Oh, sorry, I missed that answer. Dose SSCHA provide codes to plot phonon density of states?

[biglinn]

I use QE to plot phonon pdos.
I plot pdos of SSCHA dyn with one iteration and with four iterations, as follows:

It doesn't have DOS weight for the negative frequencies. So strange.