Difference between revisions of "How to analyse excitons - ICTP 2022 school"

From The Yambo Project
Jump to navigation Jump to search
 
(31 intermediate revisions by 4 users not shown)
Line 1: Line 1:
In this tutorial you will learn (for 3D-hBN)  how to:
In this tutorial you will learn how to:
* analyze a BSE optical spectrum in terms of excitonic eigenvectors and eigenvalues
* analyze an optical spectrum obtained from BSE in terms of excitonic energies and composition
* look at the spatial distribution of the exciton
* look at the spatial distribution of the exciton
* For a similar tutorial focusing on 2D-hBN and including finite-momentum BSE, go to [[How to analyse excitons|this page]].
 
This is demonstrated for 3D hBN, for which you should have obtained the needed data from previous calculations.
 
'''Warning''': for this demonstration, we did not use converged parameters (in particular they are not converged for the k-grid). Note that one needs to converge carefully all parameters to obtain scientific relevant results. This will be the topic of [[How to choose the input parameters | one of the next tutorials]].  
 
   
   
== Prerequisites ==
== Prerequisites ==
'''Previous modules'''
'''Previous modules'''
* You must have completed the [[Calculating optical spectra including excitonic effects: a step-by-step guide]] tutorial
* You must have completed the [[Calculating optical spectra including excitonic effects: a step-by-step guide]] tutorial.
'''You will need''':
'''You will need''':
* The <code>SAVE</code> databases for 3D hBN  
* The <code>SAVE</code> databases for 3D hBN  
* The <code>ndb.BS_diago*</code> databases for 3D hBN
* The <code>3D_BSE</code> directory containing the <code>ndb.BS_diago*</code> databases for 3D hBN
* <code>ypp </code> executable
* <code>ypp </code> executable
* <code>xcrysden</code> executable
* <code>xcrysden</code> or <code>VESTA</code> executables
* <code>gnuplot</code> executable
* <code>gnuplot</code> executable
'''
'''


All the databases required for this tutorial should be in the <code>YAMBO_TUTORIALS/hBN/</code> directory.
Recall that in order to plot excitonic wave functions, you need to have done a BSE calculation with either the <code>diago</code> or <code>slepc</code> solvers and explicitly set the flag <code>WRbsWF</code> to print the wave function databases to disk.
<!--
== YAMBO calculations ==
== YAMBO calculations ==
If you have completed the tutorials on 3D hBN you should have all the databases required to do this tutorial in your <code>YAMBO_TUTORIALS/hBN/SAVE</code> and <code>YAMBO_TUTORIALS/hBN/3D_BSE</code> directories.
Enter the <code>YAMBO_TUTORIALS/hBN/</code> directory. In there you have the <code>SAVE</code> directory


  $ ls ./SAVE
  $ ls ./SAVE
  ndb.gops ndb.kindx ns.db1 ns.kb_pp_pwscf_fragment_1 ....
  ndb.gops ndb.kindx ns.db1 ns.kb_pp_pwscf_fragment_1 ....
and the <code>3D_BSE</code> directory


  $ ls ./3D_BSE
  $ ls ./3D_BSE
  ndb.BS_Q1_CPU_0 ndb.BS_diago_Q1 ndb.dip_iR_and_P_fragment_1 ndb.em1s_fragment_1 ...
  ndb.BS_Q1_CPU_0 ndb.BS_diago_Q1 ndb.dip_iR_and_P_fragment_1 ndb.em1s_fragment_1 ...
!-->
==List the excitonic energies==
First, we will obtain the list of all the exciton energies (that is the eigenvalues of the two-particle Hamiltonian) sorted both by energy and strength. The exciton strengths or intensities are defined from the expression of the optical absorption spectrum as
[[File:strengh.png|none|x150px|]]
That is, they are given by the linear combination of the square of the dipole transition matrix elements between electron-hole pairs, where A<sup>&lambda;</sup><sub>''eh''</sub>  is the exciton composition in terms of electron-hole pairs.


==Sort the excitonic eigenvalues==
To read this information from the Yambo databases produced in [[Bethe-Salpeter solver: diagonalization| this previous step]], we use the Yambo pre- and post-processing utility <code>ypp</code>.
We are now going to use <code>ypp</code> to post-process our BSE results.
 
First, type
Type:


  $ ypp -J 3D_BSE -e s 1
  $ ypp -J 3D_BSE -e s 1


which means that we are sorting the excitons for the q-index = 1 (optical limit q=0).
The new generated file ''o-3D_BSE.exc_qpt1_E_sorted '' (along with ''o-3D_BSE.exc_qpt1_I_sorted '') reports
This instructs the code to list the excitons (<code>-e s</code>) in the database ''3D_BSE'' (<code>-J 3D_BSE</code>) for the q-index = 1 (optical limit q=0).  
the energies of the excitons and their Dipole Oscillator Strengths sorted by Energy (Intensity).
 
[[File:strengh.png|none|x120px|]]
The results are output in the files ''o-3D_BSE.exc_qpt1_E_sorted'' and ''o-3D_BSE.exc_qpt1_I_sorted''.
Open the first file and look inside. The first exciton is doubly degenerate at the energy of 3.54 eV, and the second one - also doubly degenerate - has the highest strength (normalized to 1). The latter exciton is the one appearing as the major peak in the optical absorption plots of the [[Bethe-Salpeter solver: diagonalization|previous tutorials]].
They report the energies of the excitons and their strengths. The strengths are normalised to the largest strength. Thus in the list, the brightest exciton has strength 1.
In ''o-3D_BSE.exc_qpt1_E_sorted'' the exciton energies are sorted by energy  (lower to higher), in ''o-3D_BSE.exc_qpt1_I_sorted'' by the oscillator strength (larger to smaller).
 
Inspect ''o-3D_BSE.exc_qpt1_E_sorted'':  
 
The exciton with the lowest energy (4.98 eV) is doubly degenerate. The strength is ~0, thus this exciton is dark, does not contribute to the absorption spectrum.  


You can also make a plot:
The second-lowest exciton is also doubly degenerate and it is the brightest (one of the components has strength 1, the other ~0.03). The latter exciton is the one responsible for the largest peak in the optical absorption plots of the [[Bethe-Salpeter solver: diagonalization|previous tutorials]].
 
To better visualize these results, we plot the exciton strength (normalized to 1, arbitrary units) versus the energy (in eV):
  $ gnuplot
  $ gnuplot
  gnuplot> set style line 2 lc rgb 'black' pt 7  # circle
  gnuplot> set style line 2 lc rgb 'black' pt 7  # circle
  gnuplot> plot 'o-3D_BSE.exc_qpt1_E_sorted' with points ls 2 title 'Strengths'
  gnuplot> plot 'o-3D_BSE.exc_qpt1_E_sorted' with points ls 2 title 'Strengths'


[[File:strenght.png|none|600px|]]
[[File:Exc weights3DhBN2.png|none|600px]]


Attention: the convergence of these results with different k-points grids is mandatory!
== Calculate the exciton oscillator strength and amplitude ==


== Calculate the exciton oscillator strenght and amplitude ==
Once inspected the list of excitons sorted by energies and strengths, we will look at the composition of some of these excitons in terms of single-particle states (electron-hole pairs).


We can now analyze the excitons in terms of single-particle states: which electronic transitions are the most relevant? In order to do that, create the appropriate input with
Type:  
  $ ypp -F ypp_AMPL.in -J 3D_BSE -e a 1  
  $ ypp -F ypp_AMPL.in -J 3D_BSE -e a 1  


Suppose you wish to analyze the first 4 excitons (i.e., the first 2 doubly-degenerate excitons). Then, change this line as:
This instructs the code to create the the input file ''ypp_AMPL.in'' (<code>-F ypp_AMPL.in</code>) to analyze the composition of the excitons (<code>-e a</code>) in the database ''3D_BSE'' (<code>-J 3D_BSE</code>) for the q-index = 1 (optical limit q=0).
 
You can now edit the input file ''ypp_AMPL.in''. From the list of excitons, we should have got an idea of which excitons we wish to analyze. The ordering from the list of exciton sorted by energies in ''o-3D_BSE.exc_qpt1_E_sorted'' is used to identify the excitons to analyse.
 
For example, to analyse the 4 lowest-energy excitons (i.e., the first 2 doubly-degenerate excitons of which the first is dark and the second is bright) change this line as:
  States= "1 - 4"              # Index of the BS state(s)
  States= "1 - 4"              # Index of the BS state(s)


Close the input and run ypp
Close the input file and run <code>ypp</code> again by typing


  $ ypp -F ypp_AMPL.in -J 3D_BSE
  $ ypp -F ypp_AMPL.in -J 3D_BSE
This carries out the analysis for the specified excitons. The results are output in ''o-3D_BSE.exc_qpt1_amplitude_at_*'' and  ''o-3D_BSE.exc_qpt1_weights_at_*'', where ''*'' identifies the exciton according to the order in ''o-3D_BSE.exc_qpt1_E_sorted''.


  $ ls  o*exc*at*
  $ ls  o*exc*at*
  o-3D_BSE.exc_qpt1_amplitude_at_1 o-3D_BSE.exc_qpt1_weights_at_1 ...
  o-3D_BSE.exc_qpt1_amplitude_at_1 o-3D_BSE.exc_qpt1_weights_at_1 ...
   
   
For an exciton  <math>|\lambda></math> , ''o-3D_BSE.exc_qpt1_weights_at_*'' report the Weights
For an exciton  <math>|\lambda></math> , ''o-3D_BSE.exc_qpt1_weights_at_*'' reports the weights defined as
[[File:Weights.png|none|x60px|]]
[[File:Weights.png|none|x60px|]]
and ''o-3D_BSE.exc_qpt1_amplitude_**'' report the amplitudes  
and ''o-3D_BSE.exc_qpt1_amplitude_* '' reports the amplitudes defined as
[[File:Ampl.png|none|x70px|]]
[[File:Ampl.png|none|x70px|]]


Open the file ''o-3D_BSE.exc_weights_at_3''
 
Inspect, for example, the file ''o-3D_BSE.exc_weights_at_3'' (which we learned from ''o-3D_BSE.exc_qpt1_E_sorted'' has the largest strength)


  # Band_V          Band_C          Kv-q ibz        Symm_kv        Kc ibz          Symm_kc        Weight          Energy
  # Band_V          Band_C          Kv-q ibz        Symm_kv        Kc ibz          Symm_kc        Weight          Energy
Line 78: Line 112:
   8.00000000      9.00000000      13.0000000      1.00000000      13.0000000      1.00000000    0.738964081E-1  4.81087065   
   8.00000000      9.00000000      13.0000000      1.00000000      13.0000000      1.00000000    0.738964081E-1  4.81087065   


The third exciton is mostly composed only single-particle transitions from VBM to CBM at points H (last k-point of the grid, number 14) and K (number 13) of the hexagonal Brillouin zone. All the contributions weighing less than 5% are not shown by default. Recall from the previous analysis that exciton states 3 and 4 (degenerate) are the optically active ones.
The third exciton is mostly composed of single-particle transitions from VBM to CBM at point H (last k-point of the grid, number 14) of the 3D hexagonal Brillouin zone, with contributions also coming from point K (number 13). All the contributions weighing less than 5% are not shown by default.  


[[File:Amplitude_plot.png|none|600px|]]
Recall from the previous analysis that of the first 4 exciton states, exciton states 3 and 4 (degenerate) are the optically active ones. We plot then the amplitude of the lowest energy bright exciton:
$ paste o-3D_BSE.exc_qpt1_amplitude_at_3 o-3D_BSE.exc_qpt1_amplitude_at_4 > o-3D_BSE.exc_qpt1_amplitude_at_3_4
$ gnuplot
gnuplot> set xlabel 'Energy (eV)'
gnuplot> set ylabel 'Amplitude'
gnuplot> p 'o-3D_BSE.exc_qpt1_amplitude_at_3_4' u 1:($2+$4)/2 w l t 'Bright exciton'
 
[[File:ExcAmp3DhBN.png|none|600px]]
 
Note that while this plot is related to the peaks in the optical absorption spectrum, it does not take into account the dipole matrix elements and the peaks appear at the energies of the single-particle transitions rather than the exciton energies.


== Plot the exciton spatial distribution ==
== Plot the exciton spatial distribution ==


To see the spatial character of the exciton YPP writes the exciton spatial distribution, in other words the probability to find the electron somewhere in the space when the hole is fixed in a give position.
A further analysis we carry out is on the spatial distribution of the exciton, specifically of the probability of finding the electron at a certain position r when the hole is fixed in a position r'.
Different output formats can be selected and 1D,2D,3D plots done.
 
Create the input and change the size of the cell where to see the exciton.
Create the input file ``ypp_WF.in`` using the command below and change the various parameters as shown below:
Note that  If the k-grid of the BSE simulation is a NxNx1  the exciton has an induced fictitious periodicity every Nx Nx1 Cell  of the simulation.
 
For hBN-2D this is not a problem because the exciton is strongly localized but in other systems with  more delocalized excitons to look at the real exciton size it is necessary to use
  $ ypp -F ypp_WF.in -J 3D_BSE -e w
very large k-grids in the BSE
  $ ypp -F ypp_WF.in -J 2D_WR_WC  -e w 1


  excitons                    # [R] Excitons
  excitons                    # [R] Excitons
  wavefunction                # [R] Wavefunction
  wavefunction                # [R] Wavefunction
  Format= "x"                  # Output format [(c)ube/(g)nuplot/(x)crysden]
  Format= "x"                  # Output format [(c)ube/(g)nuplot/(x)crysden]
  Direction= "12"              # [rlu] [1/2/3] for 1d or [12/13/23] for 2d [123] for 3D
  Direction= "123"              # [rlu] [1/2/3] for 1d or [12/13/23] for 2d [123] for 3D
  FFTGvecs=  3951       RL   # [FFT] Plane-waves
  FFTGvecs=  30       Ry   # [FFT] Plane-waves
  States= "1 - 1"              # Index of the BS state(s)
  States= "3 - 3"              # Index of the BS state(s)
  Degen_Step=  0.0100  eV    # Maximum energy separation of two degenerate states
  Degen_Step=  0.0100  eV    # Maximum energy separation of two degenerate states
  % Cells
  % Cells
Line 103: Line 144:
  %
  %
  % Hole
  % Hole
  2.4    | 1.400    | 0.00    |        # [cc] Hole position in unit cell
  2.35800028 | 1.36139178 | 7.08835602 # [cc] Hole position in unit cell (positive)
 
  %
 
Close the input and run ypp
  $ ypp -F ypp_WF.in -J 2D_WR_WC
 
$ xcrysden --xsf o-2D_WR_WC.exc_2d_1.xsf


[[File:exc_BN2D.png|none|400px|]]
<code>Cells</code>, is the size of the cell where the exciton will be visualised.
If the k-grid used for the BSE calculations is XxYxZ, then the exciton has an induced fictitious periodicity for every XxYxZ cells of the simulation.
For hBN, this is not a problem because the first bright exciton is strongly localized, but in other systems, where excitons are more delocalized, one must use very large k-grids for the BSE calculation.


<!--
<code>Degen_Step</code>, tells the maximum energy separation of two degenerate states. Since this value is small enough, the code will automatically recognize excitons 3 and 4 are degenerate and merge them: you can check this in the log file.
Or alternatively
$ xcrysden --sushi
[[File:SUSHI.JPG|none|400px|]]
-->


== Plot electron/hole average density (only in Yambo 5.x)==
Plots can be in 1, 2 or 3D (<code>Direction</code>) and the output format is chosen with <code>Format</code>.


Another way to analyze excitons, it is the possibility to plot the average electron/hole densities defined as:
A non-trivial input parameter is the hole position. The best position for the hole is where the valence electrons contributing to the exciton are localised. In the case of hBN, they are on the nitrogen atoms. We can obtain the Cartesian coordinates of these atoms from our previously generated ''o-3D_BSE.exc_qpt1_E_sorted'' file:
# Atom 1 with Z 5 [cc]:  2.35800028      1.36139178      0.00000000 
# Atom 2 with Z 5 [cc]: -2.35800028    -1.36139178      6.08835602 
# Atom 1 with Z 7 [cc]: -2.35800028    -1.36139178      0.00000000 
# Atom 2 with Z 7 [cc]:  2.35800028      1.36139178      6.08835602
We take the positive values of the nitrogen (Z=7, last line) and shift the vertical position of the hole from 6.088 to 7.088 so that the hole does not end up in the center of the atom (and therefore, possibly in a node of the wave function). 


[[File:Electron hole density.png|400px]]
Close the input file and run <code>ypp</code> by typing
$ ypp -F ypp_WF.in -J 3D_BSE


to generate the corresponding input just type
After the calculation is completed, we can visualise the output with xcrysden or with VESTA (if you are running this tutorial on a remote cluster, you should download the *.xsf file using <code>scp</code> or <code>rsync</code>):


  ypp -F ypp_WF.in -e w -avehole
  $ xcrysden --xsf o-3D_BSE.exc_qpt1_3d_3.xsf


and choose the exciton you want to plot. The electron/hole average densities correspond to generalized valence/conduction orbitals for a given exciton. They are interesting in particular for molecular crystals because they allow distinguishing charge-transfer versus Frenkel excitons, from the relative position of the electron/hole densities.
$ VESTA o-3D_BSE.exc_qpt1_3d_3.xsf


If you want see an example of hole/electron density of excitons please have a look to ref. <ref>Strongly Bound Excitons in Metal-Organic Framework MOF-5: A Many-Body Perturbation Theory Study, A. R. Kshirsagar et al., [https://doi.org/10.26434/chemrxiv.14034917.v1 preprint ChemRxiv]</ref>.
[[File:ExcWF3DhBN.jpg|none|600px]]


== Interpolate exciton dispersion (only in Yambo 5.x)==
We notice that the electron is completely confined on the boron atoms and on the same layer of the hole. For comparison, see for example ref. <ref>Huge Excitonic Effects in Layered Hexagonal Boron Nitride, B. Arnaud et al., [https://arxiv.org/abs/cond-mat/0503390 preprint ArXiv]</ref>.


Starting from Yambo 5.x it is possible to calculate excitons at finite momentum '''q'''. The momentum can be specified with the variable
== Plot electron/hole average density (only in Yambo 5.x)==
%BSEQptR
  iq_start | iq_end |                            # [BSK] Transferred momenta range
%
where iq_start and iq_end is the first and last index of momentum in the irreducible Brillouin zone (IBZ). You can find the list of q-vectors in the r_setup report.
If you calculate the BSE for all q-points of the IBZ, then you can interpolate  exciton dispersion along any direction of the full Brillouin zone(BZ) using the
interpolation scheme of ref. <ref>Warren E. Pickett, Henry Krakauer, and Philip B. Allen  [https://journals.aps.org/prb/abstract/10.1103/PhysRevB.38.2721 PRB 38 p2721 (1988)]</ref>, by doing:


ypp -e i
Another way to analyze the spatial character of excitons is to plot the average electron/hole densities defined as:


and then setting:
[[File:Electron hole density.png|400px]]


excitons                        # [R] Excitonic properties
To generate the input file ypp_WF.in type
interpolate                      # [R] Interpolate
States= "1 - 4"                  # Index of the BS state(s)
INTERP_mode= "BOLTZ"                # Interpolation mode (NN=nearest point, BOLTZ=boltztrap aproach)
% INTERP_Grid
  -1 |-1 |-1 |                            # Interpolation BZ Grid
%
#PrtDOS                          # Print Exciton Density of States
BANDS_steps= 100                # Number of divisions
cooIn= "rlu"                    # Points coordinates (in) cc/rlu/iku/alat
cooOut= "rlu"                    # Points coordinates (out) cc/rlu/iku/alat
%BANDS_kpts                      # K points of the bands circuit
0 | 0 | 0 |
0.333333333333333 | 0.33333333333333 | 0 |
%


Running ypp, it will interpolate the first 4 excitons energies on 100 points along with the line Gamma->K. If you plot the output file "o.excitons_interpolated",
ypp -F ypp_WF.in -e w -avehole
you can compare the interpolated results with the original points calculated in the BSE, that are available in the report "r_excitons_interpolate", and you get something like:


[[File:Excitons dispersion.png|center | Yambo exciton dispersion]]
and choose the exciton you want to plot according to the list of energy-sorted excitons. The electron/hole average densities correspond to generalized valence/conduction orbitals for a given exciton. They are interesting in particular for molecular crystals because they allow distinguishing charge-transfer versus Frenkel excitons, from the relative position of the electron/hole densities.


If you want to see an example of hole/electron density of excitons please have a look at Ref. <ref>Strongly Bound Excitons in Metal-Organic Framework MOF-5: A Many-Body Perturbation Theory Study, A. R. Kshirsagar et al., [https://doi.org/10.26434/chemrxiv.14034917.v1 preprint ChemRxiv]</ref>.


It is also possible to interpolate excitons on a given regular grid by setting <code>INTERP_Grid</code>, for example in the figure below we interpolate the dispersion of the lowest exciton in MoS<sub>2</sub> calculated on a 39x39x1 and interpolated in a 117x117x1 grid in all the Brillouin zone:
== Summary ==
From this tutorial you've learned:
* How sort excitonic states by energy and intensity
* Analyse their composition in reciprocal space and in terms of single-particle transitions
* Visualize the exciton wave function in real space


[[File:Screenshot from 2021-03-24 11-38-46.png|center| 700px|Yambo tutorial image]]
=References =
<references />


Finally if you turn on the flag <code>PrtDOS</code> the code calculate the excitonic density of states on the excitons interpolated on the regular grid given by  <code>INTERP_Grid</code>.
==Navigate==
* [[Calculating optical spectra including excitonic effects: a step-by-step guide]] tutorial
* Back to [[Rome 2023#Tutorials]]
* Back to [[ICTP 2022]]
* [[BSE hBN Yambo Virtual 2021 version|Back to CECAM 2021 tutorial page for the BSE]]
* [[Tutorials|Back to tutorials menu]]
* [[Modules|Back to technical modules menu]]
* Similar tutorial for 2D-hBN and including finite-momentum BSE: [[How to analyse excitons|this page]].

Latest revision as of 16:42, 19 May 2023

In this tutorial you will learn how to:

  • analyze an optical spectrum obtained from BSE in terms of excitonic energies and composition
  • look at the spatial distribution of the exciton

This is demonstrated for 3D hBN, for which you should have obtained the needed data from previous calculations.

Warning: for this demonstration, we did not use converged parameters (in particular they are not converged for the k-grid). Note that one needs to converge carefully all parameters to obtain scientific relevant results. This will be the topic of one of the next tutorials.


Prerequisites

Previous modules

You will need:

  • The SAVE databases for 3D hBN
  • The 3D_BSE directory containing the ndb.BS_diago* databases for 3D hBN
  • ypp executable
  • xcrysden or VESTA executables
  • gnuplot executable

All the databases required for this tutorial should be in the YAMBO_TUTORIALS/hBN/ directory.

Recall that in order to plot excitonic wave functions, you need to have done a BSE calculation with either the diago or slepc solvers and explicitly set the flag WRbsWF to print the wave function databases to disk.


List the excitonic energies

First, we will obtain the list of all the exciton energies (that is the eigenvalues of the two-particle Hamiltonian) sorted both by energy and strength. The exciton strengths or intensities are defined from the expression of the optical absorption spectrum as

Strengh.png

That is, they are given by the linear combination of the square of the dipole transition matrix elements between electron-hole pairs, where Aλeh is the exciton composition in terms of electron-hole pairs.

To read this information from the Yambo databases produced in this previous step, we use the Yambo pre- and post-processing utility ypp.

Type:

$ ypp -J 3D_BSE -e s 1


This instructs the code to list the excitons (-e s) in the database 3D_BSE (-J 3D_BSE) for the q-index = 1 (optical limit q=0).

The results are output in the files o-3D_BSE.exc_qpt1_E_sorted and o-3D_BSE.exc_qpt1_I_sorted. They report the energies of the excitons and their strengths. The strengths are normalised to the largest strength. Thus in the list, the brightest exciton has strength 1. In o-3D_BSE.exc_qpt1_E_sorted the exciton energies are sorted by energy (lower to higher), in o-3D_BSE.exc_qpt1_I_sorted by the oscillator strength (larger to smaller).

Inspect o-3D_BSE.exc_qpt1_E_sorted:

The exciton with the lowest energy (4.98 eV) is doubly degenerate. The strength is ~0, thus this exciton is dark, does not contribute to the absorption spectrum.

The second-lowest exciton is also doubly degenerate and it is the brightest (one of the components has strength 1, the other ~0.03). The latter exciton is the one responsible for the largest peak in the optical absorption plots of the previous tutorials.

To better visualize these results, we plot the exciton strength (normalized to 1, arbitrary units) versus the energy (in eV):

$ gnuplot
gnuplot> set style line 2 lc rgb 'black' pt 7   # circle
gnuplot> plot 'o-3D_BSE.exc_qpt1_E_sorted' with points ls 2 title 'Strengths'
Exc weights3DhBN2.png

Calculate the exciton oscillator strength and amplitude

Once inspected the list of excitons sorted by energies and strengths, we will look at the composition of some of these excitons in terms of single-particle states (electron-hole pairs).

Type:

$ ypp -F ypp_AMPL.in -J 3D_BSE -e a 1 

This instructs the code to create the the input file ypp_AMPL.in (-F ypp_AMPL.in) to analyze the composition of the excitons (-e a) in the database 3D_BSE (-J 3D_BSE) for the q-index = 1 (optical limit q=0).

You can now edit the input file ypp_AMPL.in. From the list of excitons, we should have got an idea of which excitons we wish to analyze. The ordering from the list of exciton sorted by energies in o-3D_BSE.exc_qpt1_E_sorted is used to identify the excitons to analyse.

For example, to analyse the 4 lowest-energy excitons (i.e., the first 2 doubly-degenerate excitons of which the first is dark and the second is bright) change this line as:

States= "1 - 4"              # Index of the BS state(s)

Close the input file and run ypp again by typing

$ ypp -F ypp_AMPL.in -J 3D_BSE

This carries out the analysis for the specified excitons. The results are output in o-3D_BSE.exc_qpt1_amplitude_at_* and o-3D_BSE.exc_qpt1_weights_at_*, where * identifies the exciton according to the order in o-3D_BSE.exc_qpt1_E_sorted.

$ ls  o*exc*at*
o-3D_BSE.exc_qpt1_amplitude_at_1 o-3D_BSE.exc_qpt1_weights_at_1 ...

For an exciton [math]\displaystyle{ |\lambda\gt }[/math] , o-3D_BSE.exc_qpt1_weights_at_* reports the weights defined as

Weights.png

and o-3D_BSE.exc_qpt1_amplitude_* reports the amplitudes defined as

Ampl.png


Inspect, for example, the file o-3D_BSE.exc_weights_at_3 (which we learned from o-3D_BSE.exc_qpt1_E_sorted has the largest strength)

# Band_V          Band_C          Kv-q ibz        Symm_kv         Kc ibz          Symm_kc         Weight          Energy
#    
  7.00000000      10.0000000      14.0000000      2.00000000      14.0000000      2.00000000     0.395135850      4.35248947   
  7.00000000      10.0000000      14.0000000      1.00000000      14.0000000      1.00000000     0.394993663      4.35248947   
  8.00000000      9.00000000      14.0000000      2.00000000      14.0000000      2.00000000     0.391943455      4.35241365   
  8.00000000      9.00000000      14.0000000      1.00000000      14.0000000      1.00000000     0.391800284      4.35241365   
  7.00000000      10.0000000      13.0000000      2.00000000      13.0000000      2.00000000     0.745555162E-1   4.81094742   
  7.00000000      10.0000000      13.0000000      1.00000000      13.0000000      1.00000000     0.745274872E-1   4.81094742   
  8.00000000      9.00000000      13.0000000      2.00000000      13.0000000      2.00000000     0.739243180E-1   4.81087065   
  8.00000000      9.00000000      13.0000000      1.00000000      13.0000000      1.00000000     0.738964081E-1   4.81087065   

The third exciton is mostly composed of single-particle transitions from VBM to CBM at point H (last k-point of the grid, number 14) of the 3D hexagonal Brillouin zone, with contributions also coming from point K (number 13). All the contributions weighing less than 5% are not shown by default.

Recall from the previous analysis that of the first 4 exciton states, exciton states 3 and 4 (degenerate) are the optically active ones. We plot then the amplitude of the lowest energy bright exciton:

$ paste o-3D_BSE.exc_qpt1_amplitude_at_3 o-3D_BSE.exc_qpt1_amplitude_at_4 > o-3D_BSE.exc_qpt1_amplitude_at_3_4
$ gnuplot
gnuplot> set xlabel 'Energy (eV)' 
gnuplot> set ylabel 'Amplitude' 
gnuplot> p 'o-3D_BSE.exc_qpt1_amplitude_at_3_4' u 1:($2+$4)/2 w l t 'Bright exciton'
ExcAmp3DhBN.png

Note that while this plot is related to the peaks in the optical absorption spectrum, it does not take into account the dipole matrix elements and the peaks appear at the energies of the single-particle transitions rather than the exciton energies.

Plot the exciton spatial distribution

A further analysis we carry out is on the spatial distribution of the exciton, specifically of the probability of finding the electron at a certain position r when the hole is fixed in a position r'.

Create the input file ``ypp_WF.in`` using the command below and change the various parameters as shown below:

$ ypp -F ypp_WF.in -J 3D_BSE -e w
excitons                     # [R] Excitons
wavefunction                 # [R] Wavefunction
Format= "x"                  # Output format [(c)ube/(g)nuplot/(x)crysden]
Direction= "123"               # [rlu] [1/2/3] for 1d or [12/13/23] for 2d [123] for 3D
FFTGvecs=  30        Ry    # [FFT] Plane-waves
States= "3 - 3"              # Index of the BS state(s)
Degen_Step=   0.0100   eV    # Maximum energy separation of two degenerate states
% Cells
 5 | 5 | 1 |                             # Number of cell repetitions in each direction (odd or 1)
%
% Hole
2.35800028 | 1.36139178 | 7.08835602 # [cc] Hole position in unit cell (positive)
%

Cells, is the size of the cell where the exciton will be visualised. If the k-grid used for the BSE calculations is XxYxZ, then the exciton has an induced fictitious periodicity for every XxYxZ cells of the simulation. For hBN, this is not a problem because the first bright exciton is strongly localized, but in other systems, where excitons are more delocalized, one must use very large k-grids for the BSE calculation.

Degen_Step, tells the maximum energy separation of two degenerate states. Since this value is small enough, the code will automatically recognize excitons 3 and 4 are degenerate and merge them: you can check this in the log file.

Plots can be in 1, 2 or 3D (Direction) and the output format is chosen with Format.

A non-trivial input parameter is the hole position. The best position for the hole is where the valence electrons contributing to the exciton are localised. In the case of hBN, they are on the nitrogen atoms. We can obtain the Cartesian coordinates of these atoms from our previously generated o-3D_BSE.exc_qpt1_E_sorted file:

# Atom 1 with Z 5 [cc]:  2.35800028      1.36139178      0.00000000   
# Atom 2 with Z 5 [cc]: -2.35800028     -1.36139178      6.08835602   
# Atom 1 with Z 7 [cc]: -2.35800028     -1.36139178      0.00000000   
# Atom 2 with Z 7 [cc]:  2.35800028      1.36139178      6.08835602 

We take the positive values of the nitrogen (Z=7, last line) and shift the vertical position of the hole from 6.088 to 7.088 so that the hole does not end up in the center of the atom (and therefore, possibly in a node of the wave function).

Close the input file and run ypp by typing

$ ypp -F ypp_WF.in -J 3D_BSE

After the calculation is completed, we can visualise the output with xcrysden or with VESTA (if you are running this tutorial on a remote cluster, you should download the *.xsf file using scp or rsync):

$ xcrysden --xsf o-3D_BSE.exc_qpt1_3d_3.xsf
$ VESTA o-3D_BSE.exc_qpt1_3d_3.xsf
ExcWF3DhBN.jpg

We notice that the electron is completely confined on the boron atoms and on the same layer of the hole. For comparison, see for example ref. [1].

Plot electron/hole average density (only in Yambo 5.x)

Another way to analyze the spatial character of excitons is to plot the average electron/hole densities defined as:

Electron hole density.png

To generate the input file ypp_WF.in type

ypp -F ypp_WF.in -e w -avehole

and choose the exciton you want to plot according to the list of energy-sorted excitons. The electron/hole average densities correspond to generalized valence/conduction orbitals for a given exciton. They are interesting in particular for molecular crystals because they allow distinguishing charge-transfer versus Frenkel excitons, from the relative position of the electron/hole densities.

If you want to see an example of hole/electron density of excitons please have a look at Ref. [2].

Summary

From this tutorial you've learned:

  • How sort excitonic states by energy and intensity
  • Analyse their composition in reciprocal space and in terms of single-particle transitions
  • Visualize the exciton wave function in real space

References

  1. Huge Excitonic Effects in Layered Hexagonal Boron Nitride, B. Arnaud et al., preprint ArXiv
  2. Strongly Bound Excitons in Metal-Organic Framework MOF-5: A Many-Body Perturbation Theory Study, A. R. Kshirsagar et al., preprint ChemRxiv

Navigate