View source for Modena 2025 : Yambopy part 2

In this tutorial we will first plot and interpolate GW band structures calculated with YAMBO, then deal with the analysis of excitons from a BSE calculation.

Make sure you have completed [[Part 1|Modena 2025: Yambopy part 1]] and have the relevant databases in your '''$SCRATCH''' directory. 
 cd $SCRATCH
 cd YAMBOPY_TUTORIALS

= GW calculations =


==Tutorial 1: GW bands==

Yambopy can be used either to run Yambo and QE calculations, or to analyse the results of QE and Yambo by dealing with their generated databases. This is done with a variety of classes included in the qepy (for QE) or yambopy (for Yambo) modules. 
In the case of Yambo GW quasi-particle calculations, we can use the yambopy class '''YamboQPDB''' to read the database produced by the simulation.
Enter in the folder
 cd databases_qepy/gw-bands

We can use this to find the scissor operator, plot the GW bands and to interpolate the GW bands on a smoother k-path. The example runs by typing:

 python plot-qp.py

See below for an explanation of the tutorial. As usual, we can import the qepy and yambopy libraries:
  from yambopy import YamboLatticeDB,YamboQPDB
  from qepy import Path

We define the k-points path.
  npoints = 10
  path = Path([ [[  0.0,  0.0,  0.0],'$\Gamma$'],
                [[  0.5,  0.0,  0.0],'M'],
                [[1./3.,1./3.,  0.0],'K'],
                [[  0.0,  0.0,  0.0],'$\Gamma$']], [int(npoints*2),int(npoints),int(sqrt(5)*npoints)] )

Importantly, the number of points is a free choice used only in case of interpolation. We can increase the variable npoints as we wish, it just means that the interpolation step will take more time. In order to analyse GW results we need to have the file related to the basic data of our Yambo calculation ('''SAVE/ns.db1''') and the netcdf file with the quasi-particle results ('''ndb.QP'''). We load the data calling the respective classes:

 # Read Lattice information from SAVE
 lat  = YamboSaveDB.from_db_file(folder='SAVE',filename='ns.db1')
 # Read QP database
 ydb  = YamboQPDB.from_db(filename='ndb.QP',folder='qp-gw')
(as you know by now, in the yambopy module each class is specialised to read a specific Yambo database)

The first option is to plot the energies and calculate the ideal Kohn-Sham to GW scissor operator. We need to select the index of the top valence band:

 n_top_vb = 4
 ydb.plot_scissor_ax(ax,n_top_vb)

Yambopy displays the fitting and also the data of the slope of each fitting. Notice that this is also a test if the GW calculations are running well. '''If the dependence is not approaching linear you should double-check your results!'''

[[File:Figure 6-slope-scissor.png|400px|center]]

In this case the slope is:

 valence bands:
 slope:     1.0515886598785766
 conduction bands:
 slope:     1.026524081134514
 scissor list (shift,c,v) [eV,adim,adim]: [1.8985204833551723, 1.026524081134514, 1.0515886598785766]

In addition to the scissor operator, we can plot the GW (and DFT) band structure along the path. The first choice would be to plot the actual GW calculations, without interpolation, to check that the results are meaningful (or not). This plot is independent of the number of k-points ('''npoints''') that we want to put in the interpolation. The class '''YamboQPDB''' finds the calculated points that belong to the path and plots them. Be aware that if we use coarse grids the class would not find any point and the function will not work.

 ks_bs_0, qp_bs_0 = ydb.get_bs_path(lat,path)
 ks_bs_0.plot_ax(ax,legend=True,c_bands='r',label='KS')
 qp_bs_0.plot_ax(ax,legend=True,c_bands='b',label='QP-GW')

[[File:Figure 7-GW-band-structure-non-interpolated.png|400px|center]]

The interpolation of the DFT and GW band structures looks similar:

 ks_bs, qp_bs = ydb.interpolate(lat,path,what='QP+KS',lpratio=20)
 ks_bs.plot_ax(ax,legend=True,c_bands='r',label='KS')
 qp_bs.plot_ax(ax,legend=True,c_bands='b',label='QP-GW')

The '''lpratio''' can be increased if the interpolation does not work as well as intended. The SKW interpolation scheme is the same one implemented in abipy (the python interface for the Abinit DFT code).

[[File:Figure 8-GW-band-structure-interpolated.png|400px|center]]

Finally, we can compare the calculated GW eigenvalues with the interpolation.

[[File:Figure 8-GW-band-structure-comparison.png|400px|center]]

= Excitons =
Let us now move to the directory where the exciton data for the next tutorials are stored:
 cd ../../databases_yambopy

== Tutorial 2: read and sort exciton data ==
In this section we will use the script <code>exc_read.py</code> and read the exciton data resulting from a BSE calculation in order to familiarise with their structure. This time we are going to read the <code>ndb.BS_diago_Q*</code> databases which correspond to BSE calculations at various momenta Q (Q=1 being the zero-momentum optical absorption limit - notice that python indexing starts from 0, while the databases are counted starting from 1).

In order to read them in python we use the Yambopy class <code>YamboExcitonDB</code>, which can be instanced like this:
 yexc = YamboExcitonDB.from_db_file(ylat,filename=f'path/to/BSE/databases/ndb.BS_diago_Q{i}')
(where i is the number of BS_diago_Qi file. Notice that it requires a previous instance of <code>YamboLatticeDB</code>).

Now the <code>yexc</code> object contains information about the exciton energies (<code>yexc.eigenvalues</code>), wave functions (<code>yexc.eigenvectors</code>) and intensities (<code>yexc.l_residuals</code>, <code>yexc.r_residuals</code>) for Qpoint i_Q. You can start by keeping <code>Q1</code> which is the optical limit and then check other momenta.

The exciton wave function components A_{vck}^{\lambda Q} (\lambda being the exciton index) are stored as A_{t}^{\lambda} in each Q-dependent database, with t being the transition index grouping the single-particle indices vck. Therefore, <code>yexc.eigenvectors</code> is an array with (N_excitons,N_transitions) dimension. In order to make the connection t->kcv, the table <code>yexc.table</code> provides k,v,c and spin indices corresponding to each transition.

In order to get a clearer picture, take a look at the script <code>exc_read.py</code>. Run it and try to understand the outputs.
  python exc_read.py

=== Tutorial 3: plot exciton weights on k-BZ and (interpolated) band structure ===
In this section we will use the script <code>exc_kspace_plot.py</code> and work with the quantities defined in the previous section.

We will focus on the i_Q=0 case (the optical limit).
We want to check the relative importance of the various (vk->ck) electronic transitions in forming a certain excitonic state. 

Since many excitonic states are degenerate, we first need to tell Yambopy which states we want to check. For example, in monolayer hexagonal boron nitride, the lowest-bound exciton - forming the main peak in the BSE absorption spectrum - is doubly degenerate. In order to analyse it, in the <code>exc_kspace_plot.py</code> script we set
 states = [1,2]
(we can also set [3,4] to check the second excitonic state, and then [5] for the non-degenerate third state, and so on; we can check degeneracies by printing the energies using the script <code>exc_read.py</code>).

In order to make a band structure plot, it is necessary to read the band energies with the <code>YamboElectronsDB</code> class
 yel = YamboElectronsDB.from_db_file(folder='path/to/SAVE')
and specify the path in the BZ using crystal coordinates via the <code>Path</code> function:   
 npoints = 20
 path = Path([ [[  0.0,  0.0,  0.0],'$\Gamma$'],
               [[  0.5,  0.0,  0.0],'M'],
               [[1./3.,1./3.,  0.0],'K'],
               [[  0.0,  0.0,  0.0],'$\Gamma$']],
               [int(npoints*2),int(npoints),int(sqrt(5)*npoints)] )

Now let's take a closer look at the <code>exc_kspace_plot.py</code> script. You will see that it performs three plots:
* Plot of |A(k)| in the k-BZ for the selected exciton state. You can see how most contributions come from the area around K and the MK direction [<code>Kspace_Plot=True</code>]. 
* Plot of |A(v,c,k)| on the electronic band structure for the selected exciton state. Here you can better see how the top valence and bottom conduction around MK are mostly involved [<code>Bands_Plot=True</code>].
* Plot of interpolated |A(v,c,k)| on the interpolated electronic band structure for the selected exciton state [<code>Bands_Plot_Interpolated=True</code>].

As always, try to play with the parameters, for example by changing the excitonic states to analyse.

[[File:Exc bands1 hbn.jpg|YamboExcitonDB plot from yambopy tutorial|300px| Uninterpolated exciton weights on top of band structure]][[File:Exc bands2 hbn.jpg|YamboExcitonDB plot from yambopy tutorial|300px]]

=== Tutorial 4: plot BSE optical absorption with custom options ===
In this section we will use the script <code>exc_abs_plot.py</code> and work with the quantities defined in the previous sections.

The class <code>YamboExcitonDB</code> also enables us, of course, to read the dielectric function eps(w) and print absorption and related spectra in python.

The advantages over using the standard human-readable yambo outputs are mainly two:
# We can freely select energy ranges, energy steps, peak broadenings and number of excitonic states included without having to rerun the <code>yambo</code> executable every time. This reduces the chance of mistakenly editing the yambo input and allows for greater freedom since these plots can now be done in your local machine without having to transfer the output data from a remote cluster.
# The dielectric function constructed from netCDF database variables has a much higher precision (number of significant digits) than the human readable output, which may be important if you want to perform additional operations in python (e.g., computing refractive index, reflectivity and other optical spectra from it, performing finite-difference derivatives, etc).

We can read eps(w) as a <code>numpy</code> array with the specified options as (energy values are in eV):
 w, eps = yexc.get_chi(nexcitons='all',emin=0,emax=7,estep=0.005,broad=0.04)
Or, you can directly prepare a plot of its imaginary part with
 w, chi = yexc.plot_chi_ax(ax,nexcitons='all',emin=0,emax=7,estep=0.005,broad=0.04)
where <code>ax</code> is a previously defined <code>Axes</code> object from <code>matplotlib</code>.

In order to see how this works, let's take a look at the <code>exc_abs_plot.py</code> script.
You will see that it performs the two operations described above:
* Plotting the BSE absorption spectrum [<code>Absorption_Plot=True</code>].
* Reading the complex dielectric function.

Try to customise the plots by changing the <code>nexcitons</code>, <code>emin</code>, <code>emax</code>, <code>estep</code>, <code>broad</code> parameters.

[[File:Exc abs hbn.jpg|YamboExcitonDB plot from yambopy tutorial|500px]]

=Links=
* Back to [[Modena 2025 : Yambopy part 1]]
* Back to [[Modena 2025]]