Some questions about the dielectric functions calculated in YAMBO

[HaozheLi]

Dear developers,

I have some questions about how to get the meaningful dielectric functions when using YAMBO to do the GW-BSE calculation for some 2D materials. I have already known how to use RIM and CUT settings to treat low-dimensional systems, and have successfully got the 2D polarizability alpha. Below are my questions:

1. Is the output file o-GW.eps_q1_diago_bse of GW-BSE calculation containing the 2D dielectric functions \epsilon? Does this file is produced by using the formula \epsilon(q) = 1 + 2 \pi * \alpha_2D * |q|. If not this case, what is this file?

2. How to get the optical conductivity or absorbance using 2D polarizability alpha? I note that the tutorial https://www.yambo-code.eu/wiki/index.ph ... al_systems written,

For this reason yambo produces another file which contains the macroscopic polarizability (called here 'o-2D_WR_WC.alpha_q1_diago_bse') which is a well defined quantity strictly related to optical conductivity or absorbance.

Does this mean we can directly use the imagine part of alpha to describe optical conductivity or absorbance? Is there any reference formula for me to get it?

3. If I use yambopy to get the BSE optical spectra using yambopy/tutorial/databases_yambopy/exc_abs_plot.py, what the obtained \chi is? In the tail of the definition of python function get_chi, the author comment that

chi = 1. + chi*cofactor #We are actually computing the epsilon, not the chi.

So I get the epsilon instead of polarizability?

Please give me some advice. Thanks in advance.

[Daniele Varsano]

Dear HaozheLi ,

the o.eps* file contains the macroscopic epsilon calculated as a 3D system, but when the coulomb potential is 2D it behaves as 1/q in the long wavelength limit, and you have that its imaginary part is zero and its real part is equal to 1, as it should for a 2D system (note that it is not strictly (1,0), because a finite q0=1e^-5 is used for the long wavelength limit). Moreover, it tends to that limit for infinite volume.
It is produced as output for historical reason, and you can avoid to look at it.

The polarizability alpha is essentially proportional to the absorption of the system and does not depend on the volume of the supercell.

About the question on the yambo-py output, you will get an answer from the developers soon.

Best,
Daniele

[HaozheLi]

Dear Daniele,

Thanks for your reply!

Do you mean that when we deal with a 2D material with cutoff and RIM, we can not use the o.eps* file to describe the dielectric function anymore? And we can not derive the absorption coefficient from it using the formula Absorb = \sqrt(2) * \omega /c * \sqrt (\sqrt (Im^2 + Re^2) - Re)) like we dealing with 3D dielectric function, right?

Now, we can directly use the imaginary part of alpha to describe the absorption of the system, right?

Best,
Haozhe

[HaozheLi]

Dear Daniele,

I also have some questions about calculating the lifetime of excitons. To get this, the most important thing is the exciton transition dipole. I have read a lot posts in this forum, especially these of Quxiao, like viewtopic.php?t=1947. However, I'm still a little confused about the detail way to get the exciton transition dipole.

1. Firstly, according to my understanding, the oscillator strengths are contained in o*exc*sorted file (second column, but normalized to 1) . So to get the non-normalized one, should I multiply residuals^2 or residuals?

2. Secondly, I want to check that whether I am right to use mu_S^2 = (non-normalized oscillator strength ) * V * N_k / (4 * pi * HA2EV * n_spin) to get the modulus square of exciton dipoles. Here mu_S is the exciton transition dipole, V is the volume of my 2D unit cell (containing the vacuum layer), N_k is the total number of k-mesh (like 31x31x1), n_spin is 2 (without considering SOC and spin polarization). By the way, according to Daniele'r reply, the unit of oscillator strength (before normalization) is eV, why? Should not it be Bohr^2 if the oscillator strength = the modulus square of the exciton dipole and the unit of dipole is Bohr in yambo? That is the most confusing thing for me.

Please give me some advice. Thanks in advance.

Best,
Haozhe

[Daniele Varsano]

Dear Haozhe,

1) Yes, you should look at alpha and discard epsilon.
2) The strength reported in the output file is
Residuals(:) = real(BS_R(:)*conjg(BS_R(:)),SP)*real(spin_occ,SP)/(2._SP*pi)**3*d3k_factor*4._SP*pi/q_norm*HA2EV
see /ypp/excitons/excitons_sort_and_report.F

and BS_R is defined in
src/bse/K_diago_hermitian_residuals.F,
and essentially are \sum_cvk Acvk*dipole_cvk sqrt(Eck-E_vk)

They are normalized to its maximum that is also reported in the output file, so you need to multiply the strength by that renormalization factor.
It seems to me that the expression you report is correct, besides a qo^2 factor (q_norm). This is the long wavelength limit momentum, which is set to the square of q0_def_norm=1.E-5_SP.

In this way, mu_S^2 should be Bohr^2.

Anyway, please double-check along the code in the indicated subroutine that all this is correct.

Alternatively, using yambo-py you can read separately the BSE eigenvector Acvk and the dipoles_cvk and build the mu_S.

Best,
Daniele

[HaozheLi]

Dear Daniele,

Thanks for your help!
According to your reply here and before (viewtopic.php?t=1742&hilit=oscillator&s ... d&start=10), the residuals in o*exc*sorted file are Res = |BS_R|**2 * n_spin *4*pi *HA2EV / (Nk * V * |q_0|**2), right?
So the exciton transition dipole I want to get is \mu_s = |BS_R|/|q_0|, which has a unit of length, right?
Then I can easily get that using formula: mu_S^2 = (|BS_R|/|q_0|)**2 = (non-normalized oscillator strength ) * V * N_k / (4 * pi * HA2EV * n_spin).

Best,
Haozhe