Some progress in Phys excph

[rreho]

This PR builds upon the work initiated by @muralidhar-nalabothula in #57, and introduces several new features and improvements related to exciton-phonon interactions and phonon-assisted photoluminescence.

Summary of Changes

• New post-processing module: Added a dedicated folder and object for computing phonon-assisted photoluminescence (PL) based on optical property outputs.

• Exciton-phonon coupling: Implemented optimized routines to compute exciton-phonon coupling within yambopy.

• I/O improvements: Fixed data-type inconsistencies and improved reading/writing for several existing classes, notably LetzElPhCode and excitonDB, ensuring precision consistency across modules.

Documentation:

A draft tutorial notebook, exciton_phonon.ipynb, showcasing the new functionality is available.

This notebook can be built as part of the Qurex-Book using make book inside the docs folder (see PR #46).

The book is a work in progress but already includes this new example on phonon-assisted PL for hBN.

• AUTHORS.md: Updated my affiliation.

Next Steps Before Merging

I would like to:

Discuss potential improvements to the structure of folders and classes in yambopy, aiming for cleaner modularity and long-term maintainability.

Perform additional testing with users and various datasets.

Add features related to exciton symmetries.

Improve the exciton_phonon.ipynb tutorial, which is currently a preliminary draft.

Note!!!!

• This PR is built on top of an open PR Tech fixipeps2 #68 . Ideally, first fix-eps2 should be merged, then this one.
  -The interface to Wannier90 is not yet included in this PR. I am holding off on pushing it until further development, but I hope to integrate it in the near future.

[muralidhar-nalabothula]

Hi Fulvio,

Thank you very much for the review. Firstly, I think we should break this PR as this is too large to review. I will reopen my old #57 PR and we can first merge it and later slowly add exph after testing. I also like Riccardo's docs, but better to separate it out.

Dependencies :

• Regarding the dependencies, both are optional, but highly recommended (pykdtree is 2-3 times faster.) Numba is a nice way to speedup.

Degeneracy_finder

• In degeneracy_finder, the loop is there to take care of cases where consecutive values are with in tolerence. for example, consider this array array([1. , 1.001, 1.002, 1.003, 1.004, 1.005, 1.006, 1.007, 1.008, 1.009]), do we treat all of them as generate within 0.001 threshold?. This is debatable, for now I will remove the loop as such cases might be rare near the bandgaps.

Exciton spin

1. Wrong doc string for compute_exc_spin_iqpt. Updated.
2. S_z is hermitian, so should be real, but I would leave it complex.
3. S^2 will always be zeros as the electron and hole "average" cancel each other (This can be seen directly by setting $g_{cc'/vv'}$ in exciton phonon matrix elements to $\delta_{vv'/cc'}$. This implies expectation value of S^2 will not tell you if it is single or triplet.

Exciton wave functions in real space

• Davide changed Q-point to BS_Q in lumen i guess. Maybe we should request him to revert this change.

Regarding the exph implementation, I never had enough time to check/test Riccardo's implementation. I was expecting that they will work as they were directly taken from my PhD scripts repo and modified to use the new yambopy functionality such as wfdb class, new exciton functions etc. Something is not done correctly in that class. I checked with my old PhDscripts, I get similar to yambo figure. Nevertheless, we can also implement the entire exph quickly with the barebones in few lines like the below.

import numpy as np
from yambopy import YamboLatticeDB,YamboWFDB, \
    LetzElphElectronPhononDB,YamboDipolesDB, YamboExcitonDB
from yambopy.bse.exciton_matrix_elements import exciton_X_matelem
from yambopy.bse.rotate_excitonwf import  rotate_exc_wf


path = '.'
bands_range=[24,28]

bsepath =  f'{path}/GW_BSE'
#bsepath =  f'{path}/bse_q0bar_finiteqfull'
savepath = f'{path}/SAVE'
ns_db1 =   f'{savepath}/ns.db1'
ns_wf =    f'{path}/SAVE/ns.wf'
ndb_elph = f'{path}/ndb.elph_s'         # k-q -> k
#ndb_elph = f'{path}/kpq_elph/ndb.elph' # k -> k+q
ndb_dipoles=f'{path}/ndb.dipoles'

# Read lattice
lattice = YamboLatticeDB.from_db_file(filename=ns_db1)
# Read electron-phonon
elph    = LetzElphElectronPhononDB(ndb_elph,read_all=False)
# Read wave functions
wfcs    = YamboWFDB(filename='ns.wf',save=savepath,latdb=lattice,bands_range=bands_range)

## load exe dbs
exdbs = []
for ik in range(wfcs.nkpoints):
    filename = 'ndb.BS_diago_Q%d' % (ik+1)
    excdb = YamboExcitonDB.from_db_file(lattice, filename=filename,
                                             folder=bsepath,
                                             Load_WF=True, neigs=10)
    exdbs.append(excdb)


exph_mat = []
Dmats = wfcs.Dmat()
# compute ex-ph
for iq in range(elph.nq):
    ph_eig, elph_mat = elph.read_iq(iq,convention='standard')
    elph_mat = elph_mat.transpose(1,0,2,4,3)
    #
    idx_BZq = wfcs.kptBZidx(elph.qpoints[iq])
    iq_isymm = lattice.symmetry_indexes[idx_BZq]
    iq_iBZ = lattice.kpoints_indexes[idx_BZq]
    trev  = (iq_isymm >= len(lattice.sym_car) / (1 + int(np.rint(lattice.time_rev))))
    symm_mat_red = lattice.lat@lattice.sym_car[iq_isymm]@np.linalg.inv(lattice.lat)
    #
    exe_iqvec = wfcs.kpts_iBZ[iq_iBZ]#exdbs[iq_iBZ].car_qpoint
    #
    Akq = rotate_exc_wf(exdbs[iq_iBZ].get_Akcv(),symm_mat_red,wfcs.kBZ,exe_iqvec,Dmats[iq_isymm],trev,wfcs.ktree)
    #
    tmp_exph = 0.5 * exciton_X_matelem(np.array([0,0,0]), elph.qpoints[iq], Akq, exdbs[0].get_Akcv(), elph_mat,
                      wfcs.kBZ, contribution='b', diagonal_only=False, ktree=wfcs.ktree)
    ## 0.5 for Ry to Ha.
    exph_mat.append(tmp_exph)


exph_mat = np.array(exph_mat).transpose(0,1,3,2) #(q,nmmodes, initial (Gamma), final (Q))
np.save('Ex-ph.npy',exph_mat)

The above code is consistant with PhD scripts and is gauge invariant. It also does not depend on what convention we use for letz code. Below is what you get from the above script.

So I suggest, we start adding things slowly.

ps. I tried the yambo_ph input, but I got bad values, not sure why (I used lumen ex-ph branch). If possible Could you please see if the above scripts exactly same figure?

Edit: for the above script to work, you need to use PR #57

[palful]

Hi Murali,

Thank you for your comments for reopening the old PR in a more digestible format. I will go ahead with that one for now and the we merge the exc-ph stuff in a subsequent PR (possibly this one).

The branch exc-ph in lumen was recently merged by Davide starting from the one in my yambo fork. The latter is what I used for my tests. I will try to check this and understand what is going on. As for your test, it looks extremely similar to my yambo_ph test (scale and pattern are the same), although there may be some minor differences, to be checked.

Now I go ahead with the other PR.

[rreho]

I believe I fixed the bug relative to the k-space plot pointed out by @palful. Works both in the gamma only case or without.
Some numerics/summation might be involved there, or numerical discrepancies (see below). The colours displayed on each hexagon is not fair, if I were to rescale it would look pr

Passing by, I noted that the way I treat the various energy conversion in LetzElphElectronPhononDB diverged w.r.t to MN's forked yambopy. I tested that both give the same final results, but small numerical discrepancies are there.

I still don't get why you include a 0.5 factor. I do not. Moreover, it appears that the absolute value displayed on the cmap is not exactly equal. Anyone has an idea why the magnitude would be slightly different? I also noted that the raw data values we are dealing with are all quite small and around 1e-5/1e-6. Clearly, in the final result we sum, but still, perhaps small discrepancies are expected.

I would also like to emphasize that I made use of changes done in get_Akcv(), exciton_X_matelem() and exciton_X_matelem that could be another source of small numerical differences w.r.t to MN's or Yambo.

Anyway, whether this fixes the issue raised in the previous comment, I agree with you that we should first merge PR #57. In the meanwhile, I can clean this PR.

[palful]

Thank you Riccardo! Let me know when you finished the cleaning and please re-synchronize this with yambo-code:bug-fixes which now includes PR #57. Then I will go back to testing the small discrepancies.

By the way: what do you mean when you say that there are discrepancies in energy conversion in the letzelphc db class? And which 0.5 factor are you referring to?

[rreho]

The factor 0.5 has been solved, it was simply a energy conversion factor from Rydberg to Hartree that MN needed in his script.

I cleaned the repo, I reproduce MN's plot and introduced a manager for abstract geometry class and IO.
I understand that it might be too involved, therefore you have a working version without such an abstract geometry on this commit 78a266c . I might be the only one that likes it and I know it's quite painful to read (sorry @palful 😂 ). One of the reason I would like to keep it is that I have implemented the symmetries on top of it.

Before merging, we still have four problems:

1. I need to merge bug-fixes here. Gonna do it asap.
2. As it's difficult to disentangle my implementations and yours we are left with few things to check in the merging. In particular, I don't know what to do with dipolesdb.py. You already noted it, and I am not sure what to do with it.
3. It looks to me like yambo_ph disagrees with my and MN python implementation (look orange points). I will test it myself right now and come back here.
4. Discuss whether we should remove from this PR or the repo in general the doc folder that I keep dragging around...

[rreho]

regarding point 1)

• For dependencies, I kept the dependencies as in bug-fixes, but introduced developer and docs options (optional).
• dipolesdb.py must be fixed, I kept the version available in bug-fixes but it's not compatible with my implementation.
• 'exciton_X_matelem' follows https://github.com/muralidhar-nalabothula/yambopy/blob/excitons/yambopy/bse/exciton_matrix_elements.py
• in lelphcdb.py reading of the electron-phonon matrix elements has been updated in bug-fixes. I adapted my routines accordingly, but I need to test it another time. I cannot test until we merge the dipolesdb.py issue.

[rreho]

Fixed shapes of eph_mat_iq with exciton_X_matelem. For testing, I had to reintroduced my own version of dipolesdb.py

[rreho]

@palful I tested the latest version of yambo_ph. I Followed the tutorial https://wiki.yambo-code.eu/wiki/index.php?title=Exciton-phonon_coupling_and_luminescence and got this.

It is still slightly different with respect to what me and MN are getting. I know the reason.
Depending on how you specify exc_in and exc_out you get a different plot (of course). The one providede by @palful both here and in the yambo wiki I believe it has been generated with

exc_in = [0,1,2,3]
exc_out = [1,2,3,4]

[palful]

Hi Riccardo, thanks a lot! Thanks also for your other comment on the lumen PR and checking the tutorial. I'm working on reconfiguring the IO of the dipoles part, I will let you know!

[palful]

Ok, so I updated YamboDipolesDB in the bug-fixes to accept bands_range but it was a real mess. It needed to have compatibility for the DipBandsAll case, with YamboElectronsDB, dipole expansion, k-plots and IP absorption... It took some time to make sure everything worked.

In the process, I changed the class to use the from_db_file method to load the database. Although this makes the update not compatible with users' scripts, on the other hand it made the changes easier and gave the class the same loading method as other classes, avoiding confusion. Now the class can be called as YamboDipolesDB.from_db_file(latdb,filename=f'{dipoles_path}/ndb.dipoles',bands_range=bands_range). At the end of this I will make a new release where this change is noted.

Can you please try to merge this into the PR? Afterwards, I will check the rest!

[rreho]

Impressive work @palful. I was expecting the dipoles to be a nightmare and it required several important decisions on the repo. Thanks!

I merged, adapted my branch to your changes and tested (only) that I reproduce the excph plot above. I also add to change slightly the import orders bsekerneldb.py and dipolesdb.py which I believe was a bug present in bug-fixes.

Let me also state something that has never been mentioned in this PR. We have not included the direct term in the routines related to the computation of luminescence and we have only the phonon-assisted terms.

[palful]

I fixed the import bug you mention, I was using

from yambopy import YamboLatticeDB

instead of the better

from yambopy.dbs.latticedb import YamboLatticeDB

No need to change the order in this way.

For the PL comment, you're right. There are many things that could be added including cumulant expansion, calculation of the exciton renormalization factor, etc. For this PR let us focus on the satellites only, but we can add features later!

[rreho]

I added the direct PL term following formula 17 here https://arxiv.org/pdf/2212.10407 (to be checked).
WARNING I changed the convention in dipolesdb for bands_range to follow Python indexing. This might not be what we want. However, our routines in wfdb and exciton_X_matelem all follow Python indexing.
If @palful wants to support Fortran indexing, we must change accordingly.

[palful]

Hi all, I decided to go with a lighter implementation for various reasons.

First, I was still getting discrepancies in the excph matrix elements and luminescence as coded in this PR.
Second, I think there was a bit too much code: I was not sure that we needed entire new classes in order to compute the various formulas.
In addition, having a class to do everything makes the code more rigid: maybe a user just wants the formula for luminescence, but they have their own set of energies and couplings, for example because they have interpolated them on a denser grid. In the current implementation, it would not be possible for them to use the code.

So in the end I think the simpler the better. And also, with a large addition to the code, if stuff doesn't work, it's very difficult to understand what is going wrong. I have repurposed some of the old scripts by Murali. I will be putting up a different PR soon for you to review!

[rreho]

Hi,
Thank you for your help and feedback. It's ok for me to proceed with a simpler and lighter implementation as I got carried away.

I am surprised about those discrepancies, but I'll wait for the new PR to test there whether everything works.