Electron-phonon intro: plots of el-ph matrix elements on k-BZ and q-BZ

[sitangshu]

Hi guys,

Many thanks for the yambopy codes. I can now compute all the gkkp matrix elements (expanded in k and q) and the plots as shown here: Electron-phonon intro: plots of el-ph matrix elements on k-BZ and q-BZ.
I can see that in the code, one can change the initial and final electronic band states (scattering) and choose a specific phonon mode. Since the k remains fixed, I understand this as intravalley scattering. :

Code: Select all

# Customly chosen matrix elements
    i_n, i_m = [3,4] #i_n =3 -> valence band, i_n = 4 -> conduction band in 2D-hBN
    i_nu = 3 # LA phonon mode in 2D-hBN at K (ZO mode at Gamma); LO, TO modes are i_nu=4,5
    i_q = 143 # This is the K-point in the hexagonal BZ
    i_k = 143

What I need to compute is the intervalley, where one can change also the k point. say between k to k' in case of single layer mos2. how should i compute that?

Regards
Sitangshu

[palful]

Dear Sitangshu,

The electron phonon matrix element <nk| dV_lq |mk-q> is saved by yambopy in an array with dimension [Nq,Nk,Nl,Nm,Nn].
It is read from the Yambo ndb.elph_*_fragment_* databases, where each fragment corresponds to a value of q.

Here k is the electron momentum, q=k-k' is the momentum transfer (phonon momentum), n and m are electronic band indices and l is the phonon mode index.

Therefore, you have access to all intraband and interband matrix elements. For example, keeping n=m and q=1 gives the intraband value at vanishing momentum transfer. If you want the coupling between K and K' points in MoS2, you have to work with the indices of k and q. For example, you can take the k-index corresponding to K and the q-index such as K-Q=K'.

In order to find out which indices to correspond to which cartesian coordinates, you can plot the annotated k and q grids using yambopy (as seen in the tutorial).

Cheers,
Fulvio

[sirsha]

Dear Sitangshu,

The electron phonon matrix element <nk| dV_lq |mk-q> is saved by yambopy in an array with dimension [Nq,Nk,Nl,Nm,Nn].
It is read from the Yambo ndb.elph_*_fragment_* databases, where each fragment corresponds to a value of q.

Here k is the electron momentum, q=k-k' is the momentum transfer (phonon momentum), n and m are electronic band indices and l is the phonon mode index.

Therefore, you have access to all intraband and interband matrix elements. For example, keeping n=m and q=1 gives the intraband value at vanishing momentum transfer. If you want the coupling between K and K' points in MoS2, you have to work with the indices of k and q. For example, you can take the k-index corresponding to K and the q-index such as K-Q=K'.

In order to find out which indices to correspond to which cartesian coordinates, you can plot the annotated k and q grids using yambopy (as seen in the tutorial).

Cheers,
Fulvio

Hello,

I tried to use elph_plot.py code of yambopy to calculate both intra- and intervalley deformation potentials. As per the code attached (for monolayer MoS2), I use i_n = i_m = 26 as intraband value (CBM occurs at 27th band); i_nu = 1 to choose LA mode of phonon; i_k = 142 (i.e. 'K' point in the BZ plot attached, generated using bz_plot.py) as CBM occurs at the 'K' point. It gives the expanded gkkp matrix elements in q-space (i.e. g_of_q) and it seems that g_of_q is centered at i_k = 142 (i.e. 'K' point in the BZ for electrons). Is it right?

Next using this g_of_q matrix elements I tried to calculate intra- and intervalley deformation potentials for selected phonon modes using the equation that I have attached. Say, for 'Gamma' to 'K' direction (along expanded points 0, 34, 88, 124, 142 in BZ), I calculate M and its slope w.r.t. |q|. This slope near 'Gamma' gives intravalley deformation potentials. Is it right? If yes how to get the intervalley deformation potentials, as in elph_plot.py code I cannot set a particular i_q value, as g_of_q is a complete q-space matrix?

Code: Select all

"""
Tutorial for YamboElectronPhononDB.

Electron-phonon matrix element reading and plotting

EDIT the path below to point to the yambo SAVE_expnd_new folder.
"""
save_path='yambo-with-elph/MoS2'
from yambopy import *
import numpy as np
import matplotlib.pyplot as plt
import math

if __name__ == "__main__":
    """
    Main part of the script
    """
    Kspace_Plot    = False
    Qspace_Plot    = True

    # Customly chosen matrix elements
    i_n, i_m = [26,26]      #i_n =3 -> valence band, i_n = 4 -> conduction band in 2D-hBN
    i_nu = 1                # LA phonon mode in 2D-hBN at K (ZO mode at Gamma); LO, TO modes are i_nu=4,5
    i_q = 142               # This is the K-point in the hexagonal BZ
    i_k = 142

    #                    #
    # Start Yambopy part #
    #                    #

    # Create "lattice" object by reading the ns.db1 database inside the yambo SAVE_expnd_new
    ylat = YamboLatticeDB.from_db_file(filename=save_path+'/SAVE_expnd_new/ns.db1')
    #print(ylat)
    qx = ylat.car_kpoints[:,0]
    qy = ylat.car_kpoints[:,1]

    alat = 5.97229044*(0.528e-10)           # in m
    n = (2*math.pi)/alat
    q = (1e-10)*np.sqrt((n*qx)**2+(n*qy)**2)  # in 1/Ang




    # Create "elphon" object by reading the ndb.elph_gkkp* databases inside the yambo SAVE_expnd_new
    yelph = YamboElectronPhononDB(ylat,folder_gkkp=save_path+'/SAVE_expnd_new',save=save_path+'/SAVE_expnd_new')
    #print(yelph)

    # Print info on how to use this class
    #print(yelph.__doc__)

    # We select a specific K-point to plot |g(Q)| in qspace
    g_of_q = np.abs(yelph.gkkp[:,i_k,i_nu,i_n,i_m]) 

    hbar = (6.62607015e-34)/(2*math.pi)   # in (kg m^2/sec) J-sec
    m_0= (173.29518e-30)*(1.53381e3)   # in kg

    omega = np.zeros(5)                   # in 1/sec
    omega[0] = (1e12)*1.088872    # LA    
    omega[1] = (1e12)*3.399453    # LA   
    omega[2] = (1e12)*5.585689    # LA   
    omega[3] = (1e12)*6.443713    # LA   
    omega[4] = (1e12)*7.042009    # LA   

    M = np.zeros(5)                       # in eV/Ang
    M[0] = (1e-10)*27.2114*g_of_q[0]/math.sqrt(hbar/(2*m_0*omega[0]))
    M[1] = (1e-10)*27.2114*g_of_q[34]/math.sqrt(hbar/(2*m_0*omega[1]))
    M[2] = (1e-10)*27.2114*g_of_q[88]/math.sqrt(hbar/(2*m_0*omega[2]))
    M[3] = (1e-10)*27.2114*g_of_q[124]/math.sqrt(hbar/(2*m_0*omega[3]))
    M[4] = (1e-10)*27.2114*g_of_q[142]/math.sqrt(hbar/(2*m_0*omega[4])) 
    print('D0 (eV/cm)=',M*(1e8))

    q1 = np.zeros(5)
    q1[0] = q[0]
    q1[1] = q[34]
    q1[2] = q[88]
    q1[3] = q[124]
    q1[4] = q[142]

    fit = np.polyfit(q1[0:2+1], M[0:2+1], 1)
    coeff = fit[0]
    intercept = fit[1]
    fit_eq = coeff*q1[0:2+1] + intercept
    plt.plot(q1[0:2+1], fit_eq, color='r')
    plt.plot(q1, M)
    plt.xlabel('q (1/Ang)')
    plt.ylabel('M (eV/Ang)')
    print('D1 (eV)=', coeff)

    
    if Qspace_Plot:
        yelph.plot_elph(g_of_q,s=100,plt_cbar=False,marker='H',cmap='jet')
        yelph.ax.set_title('|g(q)| (Hartree)')
        yelph.ax.set_xlabel('q_x')
        yelph.ax.set_ylabel('q_y')
        plt.show()








    # We select a specific Q-point to plot |g(K)| in kspace
    g_of_k = np.abs(yelph.gkkp[i_q,:,i_nu,i_n,i_m])

    if Kspace_Plot:
        yelph.plot_elph(g_of_k,s=100,plt_cbar=False,marker='H',cmap='jet')
        yelph.ax.set_title('|g(k)| (Hartree)')
        yelph.ax.set_xlabel('k_x')
        yelph.ax.set_ylabel('k_y')
        plt.show()



    #                   #
    # End Yambopy part. #
    #                   # 

Thanks and regards,
Sirsha

[palful]

Dear Sirsha,

I tried to use elph_plot.py code of yambopy to calculate both intra- and intervalley deformation potentials. As per the code attached (for monolayer MoS2), I use i_n = i_m = 26 as intraband value (CBM occurs at 27th band); i_nu = 1 to choose LA mode of phonon; i_k = 142 (i.e. 'K' point in the BZ plot attached, generated using bz_plot.py) as CBM occurs at the 'K' point. It gives the expanded gkkp matrix elements in q-space (i.e. g_of_q) and it seems that g_of_q is centered at i_k = 142 (i.e. 'K' point in the BZ for electrons). Is it right?

It is right. However, in order to make sure that the k and q-point indexing between yambopy and yambo is consistent, I suggest to redo the BZ plot using the option "expand_mode=1" in YamboLatticeDB. This forces yambopy to expand the k-points in the same order as yambo does, e.g.:

Code: Select all

YamboLatticeDB.from_db_file(filename=save_path+'/SAVE/ns.db1',expand_mode=1)

Next using this g_of_q matrix elements I tried to calculate intra- and intervalley deformation potentials for selected phonon modes using the equation that I have attached.

Note that the rescaling by 1/sqrt(2*frequency) in your formula is automatically applied by yambopy (i.e., the "g" are what is plotted). If you want to avoid it and obtain just the "M", you have to read the matrix element with a specific option called "scale_g_with_ph_energies=False", such as:

Code: Select all

YamboElectronPhononDB.read_elph(scale_g_with_ph_energies=False)

If yes how to get the intervalley deformation potentials, as in elph_plot.py code I cannot set a particular i_q value, as g_of_q is a complete q-space matrix?

I am not sure I completely understood the question. In principle, the array yelph.gkkp read by YamboElectronPhononDB has dimensions both for q-space (first axis) and k-space (second axis). You also have the array yelph.car_qpoints which contains the coordinates of the q-points and could be plotted in a similar way as you did for the k-points in bz_plot, so that you can associate a certain q-index to its position. So for instance - assuming that you have all the needed yambo database fragments ndb.elph_gkkp_expanded_fragment_* - you could set any i_q value, and any n and m values.

Cheers,
Fulvio

[rishabhsrsw7]

Hi,

I have been plotting g_of_q for monolayer MoS2. I performed my calculations at 18x18x1 grid and hence have 324 k points. The K point is 322 (bz plot attached), for which I have plotted the intravalley g_of_q for various phonon modes (at the lowest conduction band, i.e. 26).

My i_k is 322, and g_of_q is plotted by keeping the K-point constant. I assumed that changing i_q should give me the plots for intervalley (K to K' if I keep q = 1), but it does not. For i_q 1,2 and even 322, the plots are exactly the same. My question is, how do I plot the g_of_q for intervalley transitions?

Will be very thankful for any insight you could provide.

Regards

[rishabhsrsw7]

Hi,

Awaiting your response.

Thank you and Regards

[rishabhsrsw7]

Hi,

Will be very thankful if you could provide any clarity to this issue.

Regards