EEL from KK transformation computed with BSE

[lilian]

Dear developers,

I am trying to compute electron energy loss spectra using BSE in resonant mode with the Lanczos-Haydock solver. My spectrum is not perfectly converged with respect to the k-grid, but I made sure that there are enough bands to include the energy of the plasmon.
While my absorption spectra (< 15 eV) reproduce well the literature and the experiments with this method, it is not the case for the energy of the plasmon (~20 eV) I obtain.
I noticed this warning message in the output file:

EPS via KK, poles at w>0 reflected also in w<0
[WARNING] [EPS_extend] Possible inaccuracy in the EPS prolungation

I also observed that, plotting the EEL from the o-*eel* file do not give the same result that using the formula: EEL=Im(eps)/(Im(eps)²+Re(eps)²) plotted from the o-*eps*, while it should in my opinion !

My thought is that if the KK transformation is not done properly, the real part of the dielectric function may be badly computed, and thus its second zero (that partly governs the plasmon energy) may be located at the wrong energy.

1) Could you explain me what is the purpose of this warning and the consequences?
2) Do you already observed this kind of discrepancy for the calculation of EELS using o-*eel* and o-*eps* files?
3) Is the KK transformation done in the same way than in the RPA calculation?
4) What could be the cause of the wrong plasmon energy obtained with BSE?

Many thanks in advance!

[Davide Sangalli]

Dear Lilian Khelladi,

Code: Select all

EPS via KK, poles at w>0 reflected also in w<0

This message is there because you are doing a BSE calculation within the "resonant" approximation. In such case eps(w) has poles only for w>0, while it has not for w<0.
Directly extracting the eels with EEL=Im(eps)/(Im(eps)²+Re(eps)²) within "resonant" approximation would lead to a very poor quality eels. Accordingly, the code does the extension and after KK.
If you want to avoid this, simply use in input

Code: Select all

BSEmod="retarded"   # 

The code will do EEL=Im(eps)/(Im(eps)²+Re(eps)²)

Code: Select all

[WARNING] [EPS_extend] Possible inaccuracy in the EPS prolungation

This second message is related on how the KK and to how it is done. There is a check related to the max-value of the imaginary part vs the value at a reference omega (usually w=0 if this is included in your energy axis)
The same procedures are also adopted in G-space with the response function is computed within the resonat approximation (usually not the case)

Best,
D.

[lilian]

Dear Davide,

Thank you very much for your prompt reply.

Following your suggestion, I verified that the EELS obtained from the o-*eel* file in the "resonant" approximation is in excellent agreement with that computed in the "retarded" mode.

My original issue is that, for the system I am studying, the plasmon energy computed within BSE is significantly lower (by a few eV) than the experimental values. From your explanation, I understand that this discrepancy is unlikely to originate from the implementation of the KK transformation, since it should not differ from the RPA case. In addition, as expected from the literature, and confirmed for my system, including the coupling between resonant and anti-resonant terms shifts the plasmon energy to lower values, which further worsens the agreement with experiment.
I also checked that, even if full convergence is not yet achieved, increasing the size of the screening and exchange matrices do not significantly affect the plasmon energy. Similarly, increasing the number of k-points does not seem to raise the plasmon energy either.

Note that I am performing BSE@DFT+U, which in my case is similar to BSE@DFT with a scissor operator.

At this point, I am unsure how to proceed. Have you encountered systems where BSE fails to accurately reproduce plasmon energy? Alternatively, are there other parameters or aspects that you would recommend checking?

Many thanks in advance!

[Davide Sangalli]

Dear Lilian Khelladi,
if I understand correctly, the plasmon from full BSE (transition space including both eh-exchange and eh-direct interaction) is significantly different, lower in energy, compared to the plasmon obtained within RPA.

What do you mean with RPA? A calculation in G-space, or a BSE calculation in transition space with eh-exchange interaction only?
In case you performed RPA in G-space (yambo -o c), I'd try to reproduce the same result in transition space (yambo -o b -k hartree) with coupling on.

The number of bands included is often critical for plasmons

Best,
D.

[lilian]

Dear Davide,

Thank you for your answer.

I computed EELS in four different ways: (a) RPA in G-space, (b) RPA in transition space, (c) full BSE in the resonant approximation, and (d) full BSE including coupling. All calculations are performed on top of DFT+U.

Please note that these results are not fully converged (in particular with respect to the k-point grid). However, I ensured that the BSE calculations include a sufficient number of bands such that the plasmon energy is stable. The system under study is a bulk semiconductor. The corresponding spectra, together with the experimental curve, are shown in the attached figure.

From these results, I observe the following:
1) The RPA results obtained in G-space and in transition space are in excellent agreement, as expected. This agreement could likely be further improved with better convergence.
2) The plasmon energy computed within RPA is higher than the experimental value. Including excitonic effects via BSE shifts the plasmon peak to lower energies, but this correction appears too strong, resulting in a significant underestimation compared to experiment.
3) Including the coupling between resonant and anti-resonant terms in the BSE further lowers the plasmon energy.

Given these observations, I am struggling to understand why the plasmon energy obtained from BSE is in such poor agreement with experiment. I am wondering whether there are aspects of the calculation that I could further improve, or whether this behavior has already been reported for similar systems.

Best regards,

[Davide Sangalli]

Two comments:
1) within RPA, when you work on top of DFT+U the dipoles neglect the term [x,U]. If you have a k-points grid dense enough, you can try the covariant dipoles approach. Jut add

Code: Select all

DipApproach="Covariant"

This will account for the [x,U] term as well. Overall I mostly expect a re-normalization of the intensity. Since you re-scale the experimental data, this is probably not important

2) When moving from RPA to BSE, one should also consider QP corrections as well. Since you neglect QP corrections, it is reasonable that BSE is red-shifting too much the experiment.
However, here you are on top of DFT+U and there are other points to consider
2.a) QP corretions should be defined as <\Sigma_{GW} - Vxc -U > instead the code will neglect the <U> term
2.b) the above mentioned issues with the dipoles, here will also affect the screening which enters the eh-interaction. It is probably a minor effect
2.c) as for the RPA case, there is an overall effect on the intensity, again for the above mentioned issue with the dipoles.

The BSE@DFT+GW, with W=W(RPA@DFT), is a receipt which works fine due to many cancellations. Remember that W enters both in the BSE and in the GW step.

Here you are doing BSE@DFT+U with W(RPA@DFT+U), or eventually you could try, following point 2), BSE@DFT(+U)+GW with W(RPA@DFT+U).

I'd say, it is not guaranteed that everything works fine, and the final outcome is not always predictable.
Similar issues are discussed here, for the case with the PBE0 hybrid functional as a starting point: https://arxiv.org/abs/2212.11500

Best,
D.