Large difference in Kappa_P (Lattice thermal conductivity including ph-el scattering) when using symmetry and not using symmetry.

[AnishDas-2025]

Dear experts,

I am using Phoebe v1.0 (installed from the instruction given in the manual)

I am calculating Kappa_P (Lattice thermal conductivity) including ph-el scattering. The q-grid is well converged for 300K from the only ph-ph calculation. I am using "iterative BTE" (also checked with "variational" method) with 1) useSymmetries = false, and 2) useSymmetries = true. However, I found a large difference in the values of Kappa_P when calculated using symmetry and not using symmetry. The large difference also exists in the RTA level. I think the values should be same irrespective of symmetry. I also checked increasing the K point (keeping q grid fixed) to the commensurate grid of q and also both k and q points. But still there are large difference in the values for using symmetry and not using symmetry. I inspected the Ph-Ph and Ph-el linewidths for both the cases. Ph-Ph linewidths are same for both the cases; however, they are not same for the Ph-el scattering. Ph-el scattering using symmetry looks pretty weird.

In that case, which values should I accept, using symmetry or not using symmetry? Actually, for this case we don't have any experimental data reported in the literature.

Below are attached the input file, Kappa_P values for both symmetry on and off, Ph-Ph and Ph-el linewidths.

I appreciate for your valuable suggestion
Thanks Anish

Input File

appName = "phononTransport"
sumRuleFC2 = "crystal"
phFC2FileName = "espresso.fc"
phFC3FileName = "FORCE_CONSTANTS_3RD"
windowType = "population"
qMesh = [22,22,22]
kMesh = [22,22,22] # Also checked for 44,44,44
withIsotopeScattering = true
electronH0Name = "wannier_tb.dat"
elphFileName = "phoebe.elph.hdf5"
temperatures = [300.]
chemicalPotentials = [16.23]
smearingMethod = "adaptiveGaussian"

useSymmetries = false # Or useSymmetries = true
scatteringMatrixInMemory=true
enforcePositiveSemiDefinite = false

solverBTE = ["iterative"]
maxIterationsBTE = 200
convergenceThresholdBTE = 1.0e-6
wsVecFileName = "wsvec.dat"

Kappa_P values

Not using symmetry: 58.8879 (xx), 59.8071 (yy), 117.731 (zz)
Using symmetry: 75.1392 (xx), 75.1453 (yy), 149.093 (zz)

[jcoulter12]

Hi @AnishDas-2025,

Thanks for reporting this -- I'm doing some tests of my own now, so stay tuned for an update.
I would say this isn't that surprising to me, but I'll come up with a good explanation of this for you and get back to you. In the meantime if you've tried a higher kmesh sampling, can you show me that or report what you find for the lattice thermal conductivity?

One other thing that might be possible just for testing the phel contribution is to lower your window width a bit -- phel scattering typically involves states very close to mu (basically, scattering directly across the Fermi surface), as it's got a term f(E') - f(E) ~ df/dE in the scattering rate.

My suspicion is that it's somewhat sensitive to interpolation error in the elph matrix elements, and also in my experience is very sensitive to convergence of the k-mesh. ph-ph calculations are much easier than el-ph ones, on both issues.

More soon -- and thanks again for the detailed report,
Jenny

[AnishDas-2025]

Dear Dr. Coulter,
Thank you for your valuable suggestion and sorry for the late response as I was performing the TEST you mentioned. Your suggestion always led me to better understand the working method of the code.

I used adaptive smearing and q = 22x22x22. Here are my observations (Please let me know if I am wrong anywhere),

1. For 300 K, I increased the K point sampling from K = q to 5q in case of symmetry imposed, and q to 3q for symmetry off. In both cases, the Kappa_P_Trace/3 values lie around 81 to 76 (within 6.5%) respectively, for K=5q (with symmetry) and k = 2q (no symmetry), figure attached below.

Does it mean that, if we want to impose symmetry then we need very dense K point sampling (but faster computation), while if we don’t consider symmetry then less dense K point sampling is enough (but slower computation)

2. From the perspective of linewidths for the symmetry ON case, obviously I got better scattering rate and less spurious effect for dense sampling case, K = 5q. Also, the Kappa_P_Trace/3 values converging towards dense sampling. Figure attached below.

For the non-symmetry case, similar like to symmetry imposed one, dense sampling of K mesh shows better linewidths. But interestingly, changing from k = q (Kappa_P_Trace/3 = 78.81) to k = 2q (75.57) does not lead much change in the value of Kappa_P_Trace/3, but the line widths for both are quite different, specifically below 1E-4. Figure attached.

3. I also compared the linewidths between Symmetry_On_K=5q and Symmetry_Off_K=2q. Here I also find that the linewidths quality of Symmetry Off case is better than the Symmetry On case. Figure attached.

4. Now, for the window selection part. I compared two cases with Symmetry ON ,1) K = 5q, and 2) K = q. Both for dn/dT and df/dT < 1.0e-10, and dn/dT and df/dT < 1.0e-5. I did not compare the Symmetry OFF case.

For the K= 5q case, no change in KAPPA_P for both the windows (< 1.0e-10 and < 1.0e-5). Also, there is no change in the linewidths. Figure attached.

For the K = q case, the KAPP_P values are unconverged. No change in KAPPA_P for both the windows (< 1.0e-10 and < 1.0e-5). But the linewidth for the small window is slightly better than the large window. Figure attached. So, I think as long as I have dense K sampling, I can use less window.

So, in a nutshell, is the choice of “Symmetry = OFF” is a better choice for Ph-el scattering? I also faced this problem in my last work, then I shifted to the “Symmetry = OFF” then check the convergence. We got excellent agreement with the experiment. But in my previous work, I was lucky enough that the system was a bit complex and I don’t need ultra dense grid, so I could perform with “less dense + symmetry OFF” grid. But in this case, I may need to go with a very dense grid even at 300 K. I may reduce the window, but a dense grid without symmetry is always challenging.

Thanks Anish

[jcoulter12]

Hi Anish,

Thank you for sharing this detailed calculation.
You wrote such a nice, detailed message -- let me try to summarize to make sure I understood all your points:

a) The calculations become closer when they are properly converged. This can be seen from the no-sym k = 2xq vs sym k = 5xq calculation.

b) Small linewidths appear in some cases when the window is extended
I can explain why the linewidths are different when the window is bigger, this is very simple and physical -- you see that the large linewidths (the important ones) are the same. The very small linewidths, which will contribute almost nothing to transport because they are on the edge of the population window cutoff (phel scattering has a f(E') - f(E) ~ df/dE term in it), show up when you make the window bigger. This is all ok -- this is why it doesn't change your value of kappa. In the case that's poorly converged in k-points, you get very few accepted ph-el scattering processes and therefore very tiny linewidths (as in your final plot).

c) Now, most importantly, as you point out -- the no sym case looks nicer than the sym case does at the same lower value of k.
For this one, I want to ask a couple follow up questions to make sure I get the right answer -- first, are you doing just an RTA calculation alone here? Can I see the input file?
I've gone through my implementation for the RTA case and it seems ok to me, but I want to be sure for you.

Best,
Jenny

[AnishDas-2025]

Dear Dr. Coulter,
Thank you for your valuable suggestion.

a) Yes, that is true. sym k = 5xq and no-sym k = 2xq are definitely closer, but with little more small linewidths appearing in the sym-ON case. I think this is the reason, I am getting little larger Kappa_P for sym k = 5xq (81.17197 W/mK) than the no-sym k = 2xq (75.56947 W/mK). But symmetry-ON dense grid calculation is very faster than the symmetry-OFF less dense grid.

b) Yes, the window selection is clear to me now. Thats why for the low T (~ 50 K), which requires very dense grid (because df/dE is very small at low T and to maintain sufficient grid points in that small Fermi surface), I reduced the window to include denser grid.

c) I am doing the iterative BTE calculation because in this system the in-scattering might be important, but I guess the lifetime/linewidths are printed in RTA only; rta_phel_relaxation_times.json.

Sure, the input file is attached below.

Input File

appName = "phononTransport"
sumRuleFC2 = "crystal"
phFC2FileName = "espresso.fc"
phFC3FileName = "FORCE_CONSTANTS_3RD"
windowType = "population"
qMesh = [22,22,22]
kMesh = [22,22,22] # Also checked for 44,44,44, 66,66,66, 88,88,88 etc.
withIsotopeScattering = true
electronH0Name = "wannier_tb.dat"
elphFileName = "phoebe.elph.hdf5"
temperatures = [300.]
chemicalPotentials = [16.23]
smearingMethod = "adaptiveGaussian"

useSymmetries = true # Or useSymmetries = false
scatteringMatrixInMemory=true
enforcePositiveSemiDefinite = false

solverBTE = ["iterative"]
maxIterationsBTE = 200
convergenceThresholdBTE = 1.0e-6
wsVecFileName = "wsvec.dat"

[jcoulter12]

OK, this is very helpful to understand -- there's something specific I want to check about the iterative case with symmetries, I'll get back to you later today or sometime tomorrow about this.

Thanks again for the detailed reports --
Jenny

[AnishDas-2025]

Thank you for your support.

[jcoulter12]

Just wanted to write to give you a small update because it's been longer than I initially said -- I've done a bunch of work to see if I can improve this, it needs just a little more time before I'm sure, but I would speculate that even with what I've attempted to upgrade so far, it doesn't make a huge change in the convergence -- this may just be the nature of this kind of scattering calculation.

One thing I did want to note -- you're showing the trace of the thermal conductivity, but did you look at the off diagonal elements as well? I suspect those are much smaller in the case with symmetry, as this often happens and it's a positive thing. The transport tensor becoming more "diagonal" is not uncommon to see with/without symmetries in any transport calculation. This may make up some of the difference between your sym/no sym calculations... curious if this is true.

I'll update again soon,
Jenny

[jcoulter12]

Hi @AnishDas-2025,

Again providing an update to keep you posted -- I believe I've figured out something that should improve/fix this and I have implemented it. It seems to work, but I will do final tests and push the fix for you over this weekend.

Thanks very much for your patience and for reporting this behavior. The related issue was extremely subtle, and took a bit to debug properly.

Best,
Jenny

[jcoulter12]

Ok, I have merged a PR to develop that I hope is going to help you.
Thank you again both for reporting behavior that didn't seem right, and for being very patient with me while I looked into it!

Best,
Jenny

[AnishDas-2025]

Dear Dr. Coulter,
Firstly, I would like to apologise for being very late to reply because last week I was traveling too much and could not manage a stable network to reply.

Secondly, thank you very much for your effort to fix the issue. How to install it? From scratch? Or "git pull && cd build && make phoebe"

1. The present system, I am working on possesses large acoustic-optical gap in phonon, so I guess there will not be any off-diagonal contribution in phonon-only transport. I checked it using "windows=nothing".

2. In other hand, the system is also metalic with some entangled bands near the Ef. So there might be some contribution of off-diagonal terms in the phonon-electron and only-electron transport. I didn't thoroughly checked it, but certainly I will do at high T after I return to office. I also guess, the off-diagonal terms should not have large effect for this system. In that case, yes it is true that transport tensor will be diagonal. Though, I will check once the contribution of off-diagonal with and without symmetry.

3. Actually, I noticed that, this problem with symmetry and without symmetry is more pronounced at low T (50K in my case). At high T both are OK, apart from some low line-widths, appearing with sym-ON. This should not be problem now, as you fixed it.

4. Some other queries:
   i) I read in the manual that "SOC" is not supported yet, but I saw a tag "hasSpinOrbit" in the output. So, does SOC work now?

   ii) I saw a discussion over separating Umklapp and Normal process with "outputUNTimes = true". But, in my case, when I use this keyword then no seperate keys regarding the U and N process are generated in the "rta_phph_relaxation_times.json" file. I checked this feature two months ago, not recently though.

   iii) The present system I am working on, has some particular feature in its optical phonon bands. In practical, I am expecting that AAO and AOO process might have some significant contribution apart from the AAA process. In that case, may you please suggest me that "which portion of the code" should I focus on to seperately print the AAA, AAO, AOO, and OOO processes. Though I am quite new in the transport field, but may be I can try to understand that how the code works.

Thank you again to fix the issue.
Anish

[jcoulter12]

Hi Anish,

No worries about reply time, just wanted to keep you up to date earlier because things were taking longer than I anticipated. Everyone is busy. :)

You can install the update by going into your Phoebe directory just by pulling changes down from github. You can do this by typing (as you suggested):

git pull 
cd build
make -j 2 

For your questions:
1+2) It's hard to know without trying when beyond-BTE contributions will be important -- I would suggest doing as you've done and just making a couple checks to be sure. For phonons, it seems there are some number systems where people have observed beyond BTE effects to be important, but in my experience with electron transport, I've only seen a few systems where Wigner transport has appeared important -- these have been materials with linear topological band crossings near Ef.

You can see two examples here (+ of course, plenty of correlated electron materials where ee interactions result in spectral bands, which are likely out of scope of the BTE even with Wigner contributions):
https://doi.org/10.1103/PhysRevX.15.021054
https://www.sciencedirect.com/science/article/pii/S2542529321000730?via%3Dihub

3. Actually, I noticed that, this problem with symmetry and without symmetry is more pronounced at low T (50K in my case). At high T both are OK, apart from some low line-widths, appearing with sym-ON. This should not be problem now, as you fixed it.

I've made improvements to the phel scattering code -- but you may still find some difference between symmetry on and off. It's just as I mentioned earlier -- you will likely see the transport tensor become more diagonal (off diagonal elements smaller). This is a good thing about using symmetry, and I would encourage you to use it where you can. Basically, in the interpolation procedures used to calculate matrix elements and scattering rates from them, there can be a little noise. Symmetry is helping to reduce this.

4i) I read in the manual that "SOC" is not supported yet, but I saw a tag "hasSpinOrbit" in the output. So, does SOC work now?

Unfortunately, not for QE. Basically, Phoebe is internally completely ready to handle calculations with SOC. However, I need to be very careful to interface properly with the DFPT el-ph calculation of QE, handling information about the gauge of the wavefunction, to finish this interface. It's just tricky and time consuming. Meanwhile, there are other codes to supply el-ph in development in the community -- to my knowledge, a finite difference el-ph code working with QE and also one to interface with VASP. I'm thinking it's better to wait to interface to one of those, as it will hopefully be simpler.

In the meantime, I have another DFT package interface which is new and not yet publically documented in Phoebe, but should work with SOC. If you wanted to use this other package (one called JDFTx) I could explain how. I would call this interface a beta feature right now, but I use it personally in my research

ii) I saw a discussion over separating Umklapp and Normal process with "outputUNTimes = true". But, in my case, when I use this keyword then no seperate keys regarding the U and N process are generated in the "rta_phph_relaxation_times.json" file. I checked this feature two months ago, not recently though.

I think a warning is printed that UN output is only implemented with the RTA right now, so it might not be in output if you run with the full matrix. Currently also, I am providing the UN rates for all the processes summed up, and this will appear in the json file containing the total rates rta_ph_relaxation_times.json.
A way to get separate rates currently -- you can run calculations which only perform the calculation of each rate, so that the total rate output in rta_ph_relaxation_times.json is just phel / phph, etc. If you suggest improvements to this, I'll certainly consider them.

iii) The present system I am working on, has some particular feature in its optical phonon bands. In practical, I am expecting that AAO and AOO process might have some significant contribution apart from the AAA process. In that case, may you please suggest me that "which portion of the code" should I focus on to seperately print the AAA, AAO, AOO, and OOO processes. Though I am quite new in the transport field, but may be I can try to understand that how the code works.

This is a bit tricky to separate out, even for me personally. In order to look at that information, one would probably need to print the entire scattering matrix to file and look at specific elements, but that file is very very large and because the data itself is stored in parallel during the calculation, it's non-trivial to write it to HDF5. One thing you can do to get a look at this plot:
https://phoebe.readthedocs.io/en/develop/tutorials/lifetimes.html

And it will show you for each phonon mode on the bands what scattering rate they have in the RTA.
Right now, I think this plot won't work with phel scattering, though. If this becomes necessary, let me know and I'll see if I can find time to unlock it, though it may take me some time, it still it helps to know what users need.

Another strategy might be to set a windowType = "energy" to select out just the branches you are interested in. This would get you the scattering contribution associated with that energy range of modes from all intermediate phonons -- basically something similar to the above plot, but for the full BZ rather than a path, allowing you to make a lifetimes vs energy plot as you have above (also working for phel processes, if they are important).

Overall, it's hard to pick out specific phonon modes because there's not a simple way to trace one phonon across the BZ, due to branch crossings changing the ordering of phonon modes. If I think of another way to do this, I'll let you know.

Let me know if I missed any questions,
Jenny