Yambo Community Forum

Posted: **Sat Sep 07, 2024 8:05 am**

Dear guys,

Can you please check the probable unit ambiguity in the line 214-215 of the module elphondb.py for plotting the elph gkkp elements. The line is:

Code: Select all

ph_E = self.ph_energies[iq,inu]/ha2ev # Put back the energies in Hartree units
        g[iq,:,inu,:,:] = dvscf[iq,:,inu,:,:]/np.sqrt(2.*ph_E)

The issue is that the plot of intravalley gkkp for MoS2 monolayer we obtained for TA branch, seems fine when compared to PRB 87, 115418 (2013), but the numeric values are wrong.

Upon inspection, we found that the term g[iq,:,inu,:,:] should be in Hartree.
Therefore, the elements dvscf[iq,:,inu,:,:] should be in Ha/angstrom. But, I cant find the units of dvscf[iq,:,inu,:,:].
I can see only the denominator scaling factor which is in sqrt(Hartree).

Should it be this:

Code: Select all

g[iq,:,inu,:,:] = self.alat*(dvscf[iq,:,inu,:,:]/np.sqrt(2.*ph_E))*(ha2eV)**(3/2)/(2*np.pi)

Addidtionally the q_x and q_y values seems not to be in 2pi/a.
Can you please assist us?

Regards,
Sitangshu

Posted: **Mon Sep 09, 2024 12:10 pm**

Dear Sitangshu,

Thank you for pointing this out. We have tested many times the units of the electron-phonon matrix elements across different systems and we are pretty sure they are correct. The full |g| matrix element has energy units, and since |g| = |dvscf|/sqrt(E_ph), that means that |dvscf| has units of energy^{3/2}. In Yambo, the standard energy units are hartree so the |dvscf| is given in hartree^{3/2}, and in order to obtain the correct |g| in hartree we need to temporarily convert back the phonon energies, since they by default they are given in eV like any other energy eigenvalue in yambopy.

Now, looking at Figure 6(a) of the paper you cite, it looks like they plot the value of |g| in eV. Indeed, I made a test and tried to plot that quantity by converting the gkkp array in eV (parameters: 12x12x1 k/q-grid plus spin-orbit so I have i_k=143, then for the TA mode I took i_nu=1 and for bottom conduction band i_c1=26):

Code: Select all

g_of_q = np.abs(yelph.gkkp[:,i_k,i_nu,i_c1,i_c1])*ha2ev

By doing this I get values pretty similar to the 2013 paper (see figure below).

As for the kpoint units: you are right, they are different from the quantum espresso "cart. coord. in units 2pi/alat". I also used to get confused all the time, so there is a function in YamboLatticeDB called get_units_info() as a reminder about reciprocal-space unit conventions between QE and Yambo. It always prints this output:

Yambo cartesian units [cc in yambo]:
:: self.car_kpoints*2.*pi

QE cartesian unists [cart. coord. in units 2pi/alat in QE:
:: self.car_kpoints*self.alat[0]

Internal yambo units [iku]:
:: self.iku_kpoints

Reduced coordinates [rlu in yambo, cryst. coord. in QE]:
:: self.red_kpoints

From which you can see that in Yambo, the Cartesian coordinates are different from QE because they are not in units of 2pi/alat. In yambopy, the Cartesian coordinates are in units of 2pi but not in units of 1/alat. So, in order to get the same as they use in the 2013 paper, we need to rescale like this: ylat.car_kpoints*ylat.alat[0].

Below an example on how one could do the rescaling while still using the function plot_elph from YamboElectronPhononDB (but of course it may be more convenient to do your own easily editable plotting script, this is taken from tutorial/databases_yambopy/elph_plot.py):

Code: Select all

    if Qspace_Plot:
        from yambopy.plot.plotting import BZ_hexagon   # We need to rescale the BZ border to QE cc in 2pi/a
        QE_cc_kpoints = ylat.car_kpoints*ylat.alat[0]  # Rescale kpoints to QE cc in 2pi/a
        yelph.rlat = yelph.rlat*ylat.alat[0]           # Rescale rec. lattice vectors to QE cc in 2pi/a
        yelph.plot_elph(g_of_q,kcoords=QE_cc_kpoints,s=100,plt_cbar=True,marker='H',cmap='hot') # Give explicitly rescaled kpts
        yelph.ax.patches[-1].remove()                  # Remove default BZ borders in previous units
        scaled_BZ = BZ_hexagon(yelph.rlat,center=(0.,0.),orientation=np.radians(30),color='white',linewidth=2) # New BZ borders
        yelph.ax.add_patch(scaled_BZ)                  # Add new BZ borders
        yelph.ax.set_title('|g(q)| (eV)')
        yelph.ax.set_xlabel('q_x')
        yelph.ax.set_ylabel('q_y')
        plt.show()
        plt.savefig('test_g_q_sitangshu_b%d.png'%i_c1,dpi=150)

Below is what I get in my test with updated units for |g| and qpoints. This is just a test with a coarse grid and I'm using SOC which probably wasn't used in the 2013 paper (plus different pseudopotential, lattice parameter optimization, etc). All in all, it looks pretty similar to their plot to me. But maybe I missed something, in that case please don't hesitate to let me know!

test_g_q_sitangshu_b2.png

Best,
Fulvio

Posted: **Tue Sep 10, 2024 8:30 am**

Many thanks Fulvio for your nice posts. Yes, we are now getting the plots.
In this paper PRB 85, 115317, (2012), Fig. 4 (top panel) I can see |M_{q\lambda}| plots in eV/ang. Can you suggest how can I plot them?

Regards,
Sitangshu

Posted: **Mon Oct 07, 2024 5:38 pm**

Dear Sitangshu,

The M-matrix elements in that paper should correspond to the matrix elements without the frequency prefactor, just as they are saved in the database. In order to obtain them, you should call the read_elph function of the ElectronPhononDB object like this:

Code: Select all

yelph.read_elph(scale_g_with_ph_energies=False)

In this way, the array yelph.gkkp should now contain the M-elements. As for the units, I think they depend on the dimensional prefactor, you can check if the results are the same, if not it should just be a scaling factor.

Cheers,
Fulvio

Yambo Community Forum

Unit ambiguity in elphondb.py

Unit ambiguity in elphondb.py

Re: Unit ambiguity in elphondb.py

Re: Unit ambiguity in elphondb.py

Re: Unit ambiguity in elphondb.py