Question about YamboStaticScreeningDB and plot_epsm1

[Kbozier]

Dear all,

Apologies if these are very basic questions (I’m afraid I’m quite new to this area), but I’ve been stuck on the following issue for a little while and would greatly appreciate some guidance.

I have been trying to use Yambo-py to plot the static inverse microscopic dielectric 𝜀^{-1}_{GG} (q, 𝝎=0) for fcc aluminium at the RPA level.
I used the YamboStaticScreeningDB class and attempted to plot using the plot_epsm1 method, which I believe plots the inverse dielectric against the magnitude of q. For ng1=ng2=0 (corresponding to G=[0,0,0]), the plot was very reasonable. However, for all other choices of G vector, I was surprised to see the resulting plot gives a spread of values - and hence, if plotted with a line, a very jagged curve.

yambopy_epsm1_0_0.png

yambopy_epsm1_7_7.png

This is problematic for me, as I'm trying to use the epsm1 read from YamboStaticScreeningDB in Yambo-py to calculate the screened coulomb interaction. However, since the curve is discontinuous for non-zero G, I'm concerned that the final result will then be very sensitive to the choice of q grid and it's not clear to me that it would converge.

Questions:
Question 1: Is there a reason why the database chooses to store epsm1 as a function of the magnitude |q| (rather than the true vector q) that I'm missing?

Question 2: Is this behaviour at non-zero G vectors expected, or is there likely something wrong with my input files (report file attached)?
*I can see why we may expect this, since I’d expect epsm1 to depend on |q+G|, which then depends on the direction of q relative to G.

Question 3: If this behaviour is expected, is it possible in Yambopy to get epsm1 as a function of q (vector) rather than |q| (just the magnitude), as I think this would solve my issues?

Thank you in advance for your help!
Kieran

[palful]

Dear Kieran,

You are indeed right. The reason why you see the jagged plot is because the function _getem1s (which is called by plot_epsm1) automatically returns the eps^-1 values ordered along increasing |q|-values irrespective of the value of G.

Question 1: Is there a reason why the database chooses to store epsm1 as a function of the magnitude |q| (rather than the true vector q) that I'm missing?

Actually this is just the output of the functions that I mentioned above, not how the values are stored in the database. In particular, once you load the database with YamboStaticScreeningDB, say

Code: Select all

yem1s = YamboStaticScreeningDB(save=db_path,em1s=db_path)

You automatically have access to the full thing with the attribute yem1s.X: this is a numpy array with dimensions [Nq, NG, NG] containing all the values of the screening as stored in the yambo database.

The attribute yem1s.X is not exactly the same as eps^-1, though: it's actually the quantity T_{G1,G2}(q)=\sqrt{v_G1(q)} X_{G1,G2}(q) \sqrt{v_G2(q)}. This symmetrized quantity is used by yambo for convenience. The "real" expression for the inverse microscopic screening is eps^{-1}_{G1,G2}(q)=1+v_G1(q)X_{G1,G2}(q) [X being the RPA response function, v being the (cutoffed in 2D) Coulomb potential] which can easily be obtained from T.

Question 2: Is this behaviour at non-zero G vectors expected, or is there likely something wrong with my input files (report file attached)?
*I can see why we may expect this, since I’d expect epsm1 to depend on |q+G|, which then depends on the direction of q relative to G.

Question 3: If this behaviour is expected, is it possible in Yambopy to get epsm1 as a function of q (vector) rather than |q| (just the magnitude), as I think this would solve my issues?

I don't think there's anything wrong. For the plot to appear smooth, you need to plot the eps^-1 values as a function of increasing |q+G| instead of |q| just as you say. You could do it by manipulating yem1s.X directly in your plotting script, or actually you can simply hack the _getem1s function slightly.

In yambopy/dbs/em1sdb.py, you can go to line 328 and replace

Code: Select all

x = [np.linalg.norm(q) for q in self.car_qpoints[indices]]

with

Code: Select all

x = [np.linalg.norm(q+G_vec) for q in self.car_qpoints[indices]]

Where G_vec contains the coordinates of the G-vector in Cartesian coordinates (so not [1,0,0]). If you know the index corresponding to this G-vector, then you can simply use:

Code: Select all

G_vec = self.gvectors[ng1]

I have tried this in a small test case for G-vector [1,0,0] and here's the comparison:

Code: Select all

yem1s.plot_epsm1(ax,ng1=10,ng2=10,ls='-',marker='o',label='Using |q|-ordered getem1s')
plt.plot(x-np.min(x),1.+y.real,ls='-',marker='o',c='red',label='Using |q+G|-ordered getem1s')

test_qplusG.png

It makes sense that this should be the default behaviour of yambopy when plotting finite-G elements of the response function. However, it would only work for the diagonal part where ng1=ng2, otherwise the nonlocality of the function gets in the way of a simple plot.

Cheers,
Fulvio

[Kbozier]

Dear Fulvio,

Thank you so much for your explanation - the yem1s.X numpy array is exactly what I was looking for, and when I went away and made the changes you suggested to _getem1s, the plot I obtained with increasing |q+G| rather than |q| looks much more in line with what I expected.

I’m afraid I did have one final (potentially very basic!) question, relating to transforming between the full Brillouin zone and the irreducible Brillouin zone of q-points.

The yem1s.X database contains the irreducible set, but I was curious to know how one might obtain eps^-1 for a q outside this set (i.e. given I have eps^-1_GG’ (q_IBZ) for all the points in the irreducible Brillouin zone, but what is the process for obtaining eps^-1_GG’ (q) for all points in the full Brillouin zone).

I believe I can use the Expand_kpoints method of the YamboLatticeDB to help me map the general q back to the corresponding q_IBZ, but I’m concerned that I don’t think I can just use eps^-1 _GG’ (q) = eps^-1 _GG’(q_IBZ), because once again of the issue of the direction between q and G.

For example, suppose q_IBZ is along the b1 axis, and q is -q_IBZ (so in the negative b1 direction). In this case, I would expect that if G=[1,0,0] then
eps^-1 _{GG} (q_IBZ) = eps^-1 _ { (-G)(-G) } (q), since then |q_IBZ + G| = |q - G|.

So in summary, the question is: what would be the process for obtaining eps^-1_G1,G2 (q) for all points in the full Brillouin zone from yem1s.X, and would there happen to be any helpful Yambo-py methods I should know about to help me with this?
(alternatively, is it possible to build the screened coulomb matrix using only the irreducible q-points, since in general most scattering vectors q will lie outside the IBZ and require mapping back in?)

Best wishes, and thank you once again for all the help
Kieran

[stevenbos123]

Dear Kieran,

Regarding your question about the YamboStaticScreeningDB, I am happy to provide an answer on behalf of Fulvio if you don't mind.

I’m afraid I did have one final (potentially very basic!) question, relating to transforming between the full Brillouin zone and the irreducible Brillouin zone of q-points.

Within Yambopy, the class YamboEm1sRotate enables users to unfold the em1s database from the irreducible Brillouin zone (IBZ) to the full Brillouin zone (FBZ). You can find it in the file /yambopy/em1s/em1s_rotate.py.

This class takes the em1s database as input and expands the k-points using the Expand_kpoints method. Additionally, it provides the X matrix for the full Brillouin zone through rotations. This approach ensures that the necessary symmetry transformations are applied correctly when working with data outside the IBZ.

Kind regards,
Steven Bos

[Kbozier]

Dear Steven Bos,

Thank you so much - that was exactly the sort of thing I was looking for!

Best wishes and many thanks,
Kieran