TDDFT Drude term questions

[hongtang2023]

Dear Yambo team,

I am using TDDFT ALDA kernel with Drude term. The Drude term is added by "DrudeWXd= ( 0.0430225165 , 0.012398 ) eV # [Xd] Drude plasmon ".

(a) Without " DrudeWXd= ( 0.0430225165 , 0.012398 ) eV # [Xd] Drude plasmon ", the eel plot shows no peak around 0.043 eV.
(b) However, after adding " DrudeWXd= ( 0.0430225165 , 0.012398 ) eV # [Xd] Drude plasmon ", the calculated eel still shows no peak around 0.043 eV.
Should I expect a plasmon peak around 0.043 eV when I add the " DrudeWXd= ( 0.0430225165 , 0.012398 ) eV " ?

The eel plots from (a) and (b) are almost the same. Is there anything wrong or missing in my input?
Please help me, thank you very much!

Please see my input below:

#
optics # [R] Linear Response optical properties
chi # [R][CHI] Dyson equation for Chi.
tddft # [R][K] Use TDDFT kernel
FFTGvecs= 3 Ry # [FFT] Plane-waves
X_Threads=0 # [OPENMP/X] Number of threads for response functions
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
Chimod= "ALDA" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
FxcGRLc= 2 Ry # [TDDFT] XC-kernel RL size
NGsBlkXd= 2 Ry # [Xd] Response block size
% QpntsRXd
1 | 1 | # [Xd] Transferred momenta
%
% BndsRnXd
1 | 530 | # [Xd] Polarization function bands
%
% EnRngeXd
0.000000 | 0.1500000 | eV # [Xd] Energy range
%
% DmRngeXd
0.100000 | 0.100000 | eV # [Xd] Damping range
%
ETStpsXd= 25 # [Xd] Total Energy steps
DrudeWXd= ( 0.0430225165 , 0.012398 ) eV # [Xd] Drude plasmon
% LongDrXd
1.000000 | 0.000000 | 0.000000 | # [Xd] [cc] Electric Field
%

[Daniele Varsano]

Dear Dr. Hong Tang,
the Drude peak is inserted in the IP response, then the ALDA results comes after the solution of a Dyson equation.
How does compare the two IP spectra? (You should have the IP sepctrum in the same output in column 4).

Best,
Daniele

[hongtang2023]

Dear Daniele,

Thanks so much for helping!
It looks that the column 4 also have no difference.

Please see attached two eel output files (with Drude and without Drude).

Best

[hongtang2023]

Dear Daniele,

Could you please explain a little bit more about
"the Drude peak is inserted in the IP response, then the ALDA results comes after the solution of a Dyson equation.
How does compare the two IP spectra? "

I attached the two eel files (with and without Drude) in the previous response. Please check them for me.
Should I expect the peak at 0.043 eV after I add " DrudeWXd= ( 0.0430225165 , 0.012398 ) eV "? Why did not I see it? Is anything wrong in my input file? see my input file below.

Please help!

Best,

input file:

#
optics # [R] Linear Response optical properties
chi # [R][CHI] Dyson equation for Chi.
tddft # [R][K] Use TDDFT kernel
FFTGvecs= 3 Ry # [FFT] Plane-waves
X_Threads=0 # [OPENMP/X] Number of threads for response functions
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
Chimod= "ALDA" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
FxcGRLc= 2 Ry # [TDDFT] XC-kernel RL size
NGsBlkXd= 2 Ry # [Xd] Response block size
% QpntsRXd
1 | 1 | # [Xd] Transferred momenta
%
% BndsRnXd
1 | 530 | # [Xd] Polarization function bands
%
% EnRngeXd
0.000000 | 0.1500000 | eV # [Xd] Energy range
%
% DmRngeXd
0.100000 | 0.100000 | eV # [Xd] Damping range
%
ETStpsXd= 25 # [Xd] Total Energy steps
DrudeWXd= ( 0.0430225165 , 0.012398 ) eV # [Xd] Drude plasmon
% LongDrXd
1.000000 | 0.000000 | 0.000000 | # [Xd] [cc] Electric Field
%

[Daniele Varsano]

Dear Dr. Hong Tang,
the IP response is not the same, difference are not easy to spot as you have a low resolution, if you plot them in a very small range of energy e.g. [0,0.1] yo ucan appreciate the difference. Anyway, I can see that you have a NaN at zero energy and I believe this should not happen. Developers have worked reecently on the Drude implementation can provide you a more detailed answer.

Best,
Daniele

[DarioALV]

Dear Dr. Hong Tang,

Here I make some clarifications about the use of the Drude model in yambo.

Yambo evaluates intra-band and inter-band contributions at finite q, while in the optical limit q->0 only inter-band contributions are considered.
The inclusion of the Drude model serves then to treat in an approximate way the intra-band contributions in the q->0 limit that are missing only for metals.

The form of the Drude model is: Y_D(w) = (W_D)^2/(w^2 + i*w*g_D)
where W_D and g_D are the real and the imaginary frequencies given by input, ( 0.0430225165 , 0.012398 ) eV in your case.

The Drude model is an additive term to the independent particle (IP) response. As you can see in the equation, it corresponds to a pole
at null frequency w=0, with a residue of (W_D)^2 and a broadening given by g_D. This is exactly what you are getting in columns 4 and 5 of your output.
It can become more evident if you improve the frequency sampling around w=0, as suggested by Daniele.

To answer your question about when to expect a plasmon peak around w = W_D. As explaided by Daniele, the IP response enters in a Dyson equation
when computing the RPA or ALDA responses. This will move intra-band and inter-band poles and change their residues in a non trivial way.
If you consider a given metal for which inter-band contributions are negligible with respect to the intra-band ones, so that the IP response is
given only by a Drude pole, then you may expect to find a plasmon around w = W_D in the RPA/ALDA response. Another possible scenario for the effects
of the Drude tail is to increase the residues of the inter-band poles rather than adding an extra peak.

You can read more about intra-band vs. inter-band and the derivations of the previous result in the following paper and its supplements:
10.1103/PhysRevB.107.155130

Regarding the values of the Drude frequency you are using in your input, I suggest you to make sure that they are appropriate for your system.
To me they seem too small and this could be the reason why you obtain so little effects, basically a small increase in intensity corresponding to the second scenario I described.
It is also possible to evaluate this parameter by integrating e.g. Eq.(17) of 10.1103/PhysRevB.64.195125 or Eq.(13) of 10.1103/PhysRevB.107.155130.
By the way, I made a similar suggestion in another post:
viewtopic.php?t=2347

The NaN in your output could be a numerical problem in your yambo calculation, I do not know well the internal workflow of an ALDA calculation,
for instance I am not sure if there is there could be a problem with the Drude model when forcing the time-ordering to be retarded (GrFnTpXd= "R" in your input).
I guess it should be very similar to RPA. In this regards, could you tell us if you are using smearing in your calculations?
This is my first suspect with numerical problems related to Drude. Anyway, I suggest you to try to increase the broadening to g_D = 0.1 eV to see if the NaN disappears.

I hope I have been helpful.

Best regards,
Dr. Dario A. Leon
Norwegian University of Life Sciences

[hongtang2023]

Dear DarioALV and Daniele,

Thank you very much for inform and suggestions. I am still working on them.

However, I have a simple question right now.

In a TDDFT Alda or LRC calculation, for optical calculation, I input
"
% QpntsRXd
1 | 1 | # [Xd] Transferred momenta
%
"

Is it supposed to be the first Q point or Q= exactly zero?

In my this case, I checked the " r_setup " file, it shows that
"...
IBZ Q-points : 690
BZ Q-points : 2560

K/Q-points units:
rlu = crystal or reduced units; cc = cartesian coordinates; iku = interal k-units

Q [ 1]: 0.000000 0.000000 0.000000 [rlu]
Q [ 2]: 0.000000 0.000000 0.100000 [rlu]
Q [ 3]: 0.000000 0.000000 0.200000 [rlu]
...
"


So, it looks that the first Q point is zero.

Then my further question is that What is the actual value of the first Q point, when the TDDFT code is running?
Is it exactly zero? or some very small numerical value?

For ALDA kernel, there is no need for q in the calculation, since ALDA is fxc(0,0) limit, no q and no frequency dependence.

However, for LRC kernel, fxc is alpha/q**2, it can not be calculated if q is exactly zero, right? So, what is the exact value of the first Q point used in the code?


Then, if the first Q point is assumed to be a very small finite but not zero value, then as you said " Yambo evaluates intra-band and inter-band contributions at finite q, " , so the intra-band contribution is included already? Please help me for this puzzle.

Thank you very much!

Best,
Hong

[Daniele Varsano]

Dear Hong,

the first Q point is zero and when evaluating the response function is meant as the q=0 limit.
All the expressions in this limit are evaluated using a very small, but finite q (q=10^-5).

In this case, the intra-band contribution is not included, as you do not have such a fine grid spacing in your q grid.

Yambo evaluates intra-band and inter-band contributions at finite q

This is meant for any q point of the grid of your sampling, excluding q=0.

Best,

Daniele

[hongtang2023]

Dear DarioALV and Daniele,

Thank you so much! I am making progresses.

1. If I add " DrudeWXd= ( 0.43 , 0.1 ) eV " (still in optical ALDA calculation) , i.e., increase both numbers, there is a sharp peak appearing around 0.6 eV now. However, the "NaN" are still there for the frequency=0 point for both EEL and EPS results. I am try " GrFnTpXd= "T" ", I will update if those "NaN" will disappear or not.

2. Dear Daniele, I remember in another Yambo forum thread "viewtopic.php?t=2347 ", you suggested that " the Drude is a model and needs empirical parameters from input. If they are not available, what you can do is to perform and RPA calculation at the smallest available q vectors,
identify the plasmon frequency and use it in the Drude model, i.e. approximate the q=0 plasmon peak with the one of the smallest available q. " Can we do that for finding the plasmon frequency for my system? or can we do a parabolic fit to get the plasmon frequency with several non-zero q points calculations?

3. Dear DarioALV, in your paper 10.1103/PhysRevB.107.155130, you talked about AL and Na two metals. It seems that for Na, the inter-band contribution is very small for q=0 and the intra-band one (omega_A) is dominant, and omega_A=omega_p, so for optical calculation (q=0) for Na, we have to include the Drude model term with omega_A or omega_p, is this correct? where omega_p is the plasmon frequency.
However, for Al, you said that both inter-band and intra-band contributions are there for q=0. We have omega_A , and omega_E, so, what value will be used for optical calculation (q=0) for AL with the Drude model, omega_A , omega_E or omega_p?

4. Dear DarioALV,, in the middle left on page 7 in 10.1103/PhysRevB.107.155130
" The comparison of ωp and ωA confirms that the experimental
plasmon frequency, ωp, in the case of Na corresponds mainly
to intraband contributions, while for Al there is an important
interband contribution [60], and its use as a Drude intraband
frequency would result in an overestimation of the actual ωA "
What those mean? Do you mean that for Al, we cannot use omega_p in the Drude model for Al? using omega_A?

5. Dear Daniele, the word "RPA " is kind of confusing. in the point 2, do you mean PRA is the "Chimod= "Hartree" in a TDDFT calculation? So if "Chimod= "IP" in a TDDFT calculation, it is a so-called " IP/RPA-NLF" calculation, right?


Please help me!


Best,
Hong

[DarioALV]

Dear Hong,

I'm glad to hear you are making progress.

Here I make more comments on your points:

1. Your results seem reasonable to me, the final peak at 0.6 eV is larger than the intraband frequency in your input 0.43 eV due to interband contributions (the Drude model is only for the intraband part). I'm not sure if the NaN is a bug or not, you could try also to increase the number of k-points and compile the code in double precision to see if it is something numerical.

2. A possible way to estimate the intraband frequency is to follow the recipe in 10.1103/PhysRevB.107.155130. Even without MPA, you can integrate Eq. (13) numerically. Y_E(w) is your response function for q=0, which only account for interband contributions, while Y(w) can be approximated by the smallest q of the grid, which contains both intra and interband contributions. It depends on the system whether or not the plasmon has a parabolic dispersion.

3. As mentioned by Daniele, the only missing term in a standard yambo calculation is the intraband contribution for q=0, both intra and interband terms are included at finite q. The Drude model is a way to introduce the missing term, so it always need to correspond to the intraband frequency omega_A.

> However, for Al, you said that both inter-band and intra-band contributions are there for q=0. We have omega_A , and omega_E, so, what value will be used for optical calculation (q=0) for AL with the Drude model, omega_A , omega_E or omega_p?

What is said in the paper is that in the case of Al the plasmon omega_p has not negligible contributions from both intra and interbands. The correct input for the Drude frequency is still omega_A. It is similar to your example: omega_p > omega_A due to the interband contribution. The only case where one can expect omega_p ~ omega_A is in metals like Na where the interband contributions are negligible.

4. This is related to the previous comments. Your conclusions are correct, the point is that there have been some confusion in the literature, using omega_p as Drude/intraband frequency would result in a further increase of the position of the plasmon.

Best regards,
Dario