Yambopy tutorial: electron-phonon coupling

In this tutorial we will cover the handling of electron-phonon coupling matrix elements by Yambopy. The electron-phonon calculation follows two steps:

1. Quantum Espresso calculation of phonon energies, eigenmodes and variations of the self-consistent potential via ph.x.
2. Electron-phonon matrix element calculation, symmetry expansion and conversion by the LetzElPhC.

For more information on the LetzElPhC code and how you should run step 1., please see the related documentation here.

We can use Yambopy to:

• Run LetzElPhC (both preprocessing and main run) without explicitly writing an input file.
• Convert the resulting databases into Yambo format.
• Analyse the electron-phonon coupling both in LetzElPhC or Yambo format.

In particular we will take a look at:

• Lattice geometry data (Yambo database: ns.db1, Yambopy class: YamboLatticeDB).
• Electron-phonon matrix elements (Yambo databases: ndb.elph_gkkp*, Yambopy class: YamboElectronPhononDB; LetzElPhC databases: ndb.elph, Yambopy class: LetzElphElectronPhononDB).

The scripts of the tutorial, but not the databases, can be found in the yambopy directory:

$cd tutorial/electron_phonon

The full tutorial, including the LetzElPhC and Yambo databases that we will read, can be downloaded and extracted from the yambo website:

 wget https://media.yambo-code.eu/educational/tutorials/files/yambopy_electron_phonon.tar.gz
 tar -xvzf yambopy_electron_phonon.tar.gz
 cd electron_phonon

We will work with monolayer molybdenum disulphide electron-phonon data obtained on a 6x6x1 kpoint grid, including spin-orbit interaction (so we deal with electron-phonon coupling and spinorial wavefunctions). Beware that the parameters are most certainly not converged.

Command line: electron-phonon calculation and databases [OPTIONAL]

Note: you can run this step only if you have compiled LetzElPhC and Quantum ESPRESSO on your machine. If not, please skip to the following sections of the tutorial, as the electron-phonon coupling databases are already provided in the .tar.gz file and you don't need to explicitly recalculate them.

If you run this step, please enter the directory

cd databases

Here you will find:

• A scf input file for MoS2, scf/mos2.scf
• A dvscf input file for MoS2, dvscf/mos2.dvscf
• A Yambo SAVE directory corresponding to a nscf calculation (this has been run already).
• The pseudopotential files to be used for the calculation.

As a first step, you need to run the Quantum ESPRESSO calculations (it is assumed that you are familiar with this). First run the scf calculation inside scf, then copy the mos2.save directory into the dvscf folder and run the dvscf calculation there. When everything is completed correctly, you will have achieved step 1. above.

Now, you will find the phonon results in the dvscf directory. We also have done a non-self-consistent calculation to be used in Yambo, and created the relative SAVE directory. If you check this directory, you will find the ns.wf databases containing the electronic wavefunctions.

Now, we aim to reconstruct the electron-phonon coupling matrix elements from the phonon energies, modes and potential variations (in dvscf) and the Bloch electronic wavefunctions used by Yambo (in SAVE).

In order to do this, we will run the lelphc executable of the LetzElPhC code (this will be step 2. from above). We will run via command line using yambopy, although it will be instructive to have look at the lelphc input files later.

We run in the same directory where the Yambo SAVE is (remember than you can also virtually move it with a symbolic link).

Type:

yambopy l2y

to see the help for the calculation. Have a look at the various flags and their description. For example, if we want to do a serial run of LetzElPhC for bands from [math]\displaystyle{ n_i }[/math] to [math]\displaystyle{ n_f }[/math], we should type:

yambopy l2y -ph path/of/ph_input.in -b n_i n_f

Here [math]\displaystyle{ n_i }[/math] and [math]\displaystyle{ n_f }[/math] are integers representing the initial and final band indices. Please note: in this case, the band index starts from 1 and the the interval of bands read by the code is [math]\displaystyle{ [n_i,n_f] }[/math] including the extrema of the interval.

For our system, we want to do a parallel calculation: 4 MPI tasks on q-point loop (phonon momenta) and 2 MPI tasks on the k-point loop (electron momenta). In addition, we want to explicitly specify the path of the lelphc executable.

We have a final choice to make. We can ask Yambopy to produce the default LetzElPhC electron-phonon databases, or to convert them in Yambo-compatible format. By the default, both will be printed.

• If you don't need the Yambo-style ndb.elph_gkkp* databases you can run yambopy l2y ... with the --no_gkkp option.
• If you don't need the default ndb.elph databases you can run yambopy l2y ... with the --no_lelphc_dbs option.

However, for this tutorial, we will take a look at both databases, so let's run:

yambopy l2y -ph dvscf/mos2.dvscf -b 25 28 -par 4 2 -lelphc path/to/lelphc_exe

We include the last two valence bands of monolayer MoS2 (25 and 26) and the first two conduction bands (27 and 28). If your lelphc executable is in the PATH, you do not need to specify the last flag.

At the end, check your directory: you should find the lelphc.in input file that was run, let's inspect it:

# LetzElPhC input for yambo generated by yambopy
nqpool      = 2
nkpool      = 4
start_bnd   = 25
end_bnd     = 28
save_dir    = ./SAVE
kernel      = dfpt
ph_save_dir = dvscf/ph_save
convention = yambo

Notice the variable convention=yambo: what does it mean? At variance with QE and many other codes, Yambo uses the "backward" momentum transfer convention for electronic scatterings. That is, an electronic transition goes from band [math]\displaystyle{ n }[/math] and momentum [math]\displaystyle{ k-q }[/math] to band [math]\displaystyle{ m }[/math] and momentum [math]\displaystyle{ k }[/math]. In the "forward" momentum transfer convention (the more standard one), the transitions go from [math]\displaystyle{ nk }[/math] to [math]\displaystyle{ mk+q }[/math]. Therefore, this variable ensures that the electron-phonon coupling matrix elements are computed as [math]\displaystyle{ \langle mk|dV|nk-q\rangle }[/math]. However, do not worry: Yambopy can automatically switch between the two conventions when loading the LetzElPhC databases.

You will also find a ndb.elph database that contains the computed electron-phonon matrix elements, in the output format of LetzElPhC.

In addition, if you check the SAVE folder:

ls SAVE/ndb.elph_gkkp*

you will see that yambopy has created the Yambo-compatible electron-phonon databases and placed them inside the SAVE.

If you want to run LetzElPhC directly, without using yambopy, you can refer to its [guide]. Keep in mind that in this case, in order to convert the database to the ndb.elph_gkkp* databases of Yambo, you will then need a couple of lines of python using the Yambopy class ConvertElectronPhononDB in yambopy/letzelph_interface/lelph2y.py.

At the end of this section, you should have produced the same databases that you find in the main electron_phonon folder. You can run the following part either in the directory where you are (by copying the python scripts here), or by going back:

$cd ..

and using the precomputed ones.

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

In this section we will use the script elph_plot.py and read the electron-phonon databases that you generated in the previous section.

In order to read the ndb.elph_gkkp_expanded* databases in python we use the Yambopy class YamboElectronPhononDB, which can be instanced like this:

elph = YamboElectronPhononDB(ylat,folder_gkkp='path/to/elph/folder',save='path/to/SAVE')

(notice that it requires a previous instance of YamboLatticeDB).

Instead, in order to read the ndb.elph databases from LetzElPhC, we use a different class:

elph = LetzElphElectronPhononDB('path/to/ndb.elph')

Now, the elph objects contain phonon frequencies, phonon eigenvectors, q-point information and, of course, the electron-phonon matrix elements [math]\displaystyle{ g_{nm\nu}(k,q) }[/math] where [math]\displaystyle{ n }[/math], [math]\displaystyle{ m }[/math] are electron band states, [math]\displaystyle{ \nu }[/math] is a phonon branch, and [math]\displaystyle{ k }[/math] and [math]\displaystyle{ q }[/math] are the electronic and transfer momenta.

We can print the docstring of the YamboElectronPhononDB class with

print(yelph.__doc__)

and of the LetzElphElectronPhononDB class with

print(lelph.__doc__)

to get an idea of the information stored and of their capabilities.

Now check the elph_plot.py script. Besides choosing which database to load, it also performs two plots:

• Plot of |g(k)| in the k-BZ for selected bands n,m, phonon mode \nu and q [kspace_Plot=True].
• Plot of |g(q)| in the q-BZ for selected bands n,m, phonon mode \nu and k [qspace_Plot=True].

You will see that after loading the database, it first postprocesses the electron-phonon data to obtain a plottable quantity. Indeed, it is not necessarily easy to understand what to plot and what to keep fixed, since the matrix element depends on two momenta and three state indices. In addition, bands or phonon states may be degenerate at certain momenta.

In practice, we set a phonon mode (the eighth mode, this is the branch with symmetry [math]\displaystyle{ A_1^\prime }[/math] at zone center) and either a [math]\displaystyle{ k }[/math] or [math]\displaystyle{ q }[/math] point to keep fixed (in each case, we select the high-symmetry [math]\displaystyle{ K }[/math] point). As for the bands, we consider the top two valence bands. These are mostly degenerate, except around the points [math]\displaystyle{ K }[/math] and [math]\displaystyle{ K^\prime }[/math] where they are split by spin-orbit interaction. We then average over these states before plotting, obtaining the quantity (in k-space below):

[math]\displaystyle{ I^{\mu}_{\lt v\gt q}(\mathbf{k}) = \sqrt{\sum_{v,v^\prime=24,25}|g^{\mu}_{vv^\prime q}(\mathbf{k})|^2} }[/math]

You are invited to play with the electronic, phononic and momentum indices to see what happens! (beware of degenerate states: if not properly accounted for by the sum of the norms, the plots will be phase-dependent and unphysical)