Exciton Analysis and Visualization

by Dmitrii Nabok, Christian Vorwerk, Olga Turkina, & Caterina Cocchi, for exciting nitrogen

Purpose: In this tutorial, we show how to use advanced visualization tools to analyze excitons, as computed from the solution of the BSE. With the example of lithium fluoride, you will learn how to visualize an exciton wavefunction in real-space and how to identify the main contributions to the excitation from the band structure of the material.

0. Define relevant environment variables

Before starting, be sure that relevant environment variables are already defined as specified in How to Set Environment Variables for Tutorial Scripts. Here, you can find a list of the scripts which are relevant for this tutorial, with a short description:

  • PLOT-dielectric-function.py: Python visualization script for plotting the dielectric function.
  • PLOT-exciton-weights-XX.py: Python visualization scripts for plotting the excitonic weights along a band structure path, where XX stands for LiF or BN.

Furthermore, this tutorial requires the usage of the following visualization packages:

From now on, the symbol $ will indicate the shell prompt.

Requirements: Bash shell.

1. Theoretical background: The electron-hole wavefunction

The two-particle electron-hole wavefunction is calculated from the solution of the eigenvalue problem for the Bethe-Salpeter effective Hamiltonian $H ^{\rm (eff)}$,

(1)
\begin{align} H ^{\rm (eff)}\: | A_{\lambda} \rangle= E_{\lambda}\: |A_{\lambda}\rangle\:, \end{align}

which is introduced and extensively discussed in Excited states from BSE.

Here, we focus specifically on the eigenvectors resulting from the diagonalization, namely $| A_\lambda\rangle = \sum_{v\,c\,\mathbf{k}} A^{\lambda}_{v\,c\,\mathbf{k}} |v\,c\,\mathbf{k}\rangle$, where the coefficients $A^{\lambda}_{v\,c\,\mathbf{k}} = \langle v\,c\,\mathbf{k}|A^{\lambda}\rangle$ enter the expression of the two-particle electron-hole wavefunction

(2)
\begin{align} \Psi_{\lambda}\!(\mathbf{r}_h, \mathbf{r}_e) = \sum_{vc\mathbf{k}}\:A^{\lambda}_{vc\mathbf{k}} \: \psi_{v\mathbf{k}}\!(\mathbf{r}_h) \: \psi^*_{c\mathbf{k}}\!(\mathbf{r}_e), \end{align}

as weighting factor for the transitions between the quasi-particle wavefuctions $\psi_{v\mathbf{k}}\!(\mathbf{r}_h)$ and $\psi_{v\mathbf{k}}\!(\mathbf{r}_e)$. It should be noticed that $\Psi_{\lambda}\!(\mathbf{r}_h, \mathbf{r}_e)$ depends explicitly on the electron and the hole coordinates. As we will show in the following, the two-particle wavefunction can be visualized in real space by fixing either the electron or the hole coordinate.

Another way to analyze the character of the exciton is to adopt a reciprocal-space representation. To do so, it is convenient to introduce the excitonic weights of valence and conduction bands, which are defined respectively as

(3)
\begin{align} w^{\lambda}_{v\mathbf{k}} = \sum_c|A^{\lambda}_{vc\mathbf{k}}|^2 \end{align}

and

(4)
\begin{align} w^{\lambda}_{c\mathbf{k}} = \sum_v|A^{\lambda}_{vc\mathbf{k}}|^2 . \end{align}

In these expressions we recognize again the BSE coefficients $A^{\lambda}_{v\,c\,\mathbf{k}}$ .

2. Absorption spectrum of LiF

We calculate the absorption spectrum of LiF and analyze the character of the electron-hole wavefunction that corresponds to the first and the most intense excitonic peak. A description of the basic setup and parameters to run this calculation can be found in Excited states from BSE. A detailed description of the input parameters can be found in Input Reference. Note that, in order to obtain reliable absorption spectrum, the calculation needs to be converged with respect to several parameters, e.g., the size of the k-mesh (see Excited states from BSE).

Notice that, in any input file in the following, the $EXCITINGROOT variable has to be replaced by its actual value (the path itself) or it can be automatically set by using the command:

$ SETUP-excitingroot.sh
i) Preparation of the input file

Before starting any calculation, create a directory named LiF-exciton and move into it.

$ mkdir LiF-exciton
$ cd LiF-exciton

The file below is an example for the input.
 file which can be used for the LiF ground-state calculation, as well as for the calculations of the optical absorption spectrum of LiF (do not forget to change the path in speciespath). Notice that in the groundstate block the properties element is present, in order to calculate the band structure of the system (for details on this procedure, see Electronic Structure Calculations).

<input>
 
   <title>LiF-BSE-visual</title>
 
   <structure speciespath="$EXCITINGROOT/species/">
      <crystal scale="7.608">
         <basevect>0.5 0.5 0.0</basevect>
         <basevect>0.5 0.0 0.5</basevect>
         <basevect>0.0 0.5 0.5</basevect>
      </crystal>
      <species speciesfile="Li.xml">
         <atom coord="0.0000  0.0000  0.0000" />
      </species>
      <species speciesfile="F.xml">
         <atom coord="0.5000  0.5000  0.5000" />
      </species>
   </structure>
 
   <groundstate
      ngridk="10  10  10"
      xctype="GGA_PBE_SOL"
      gmaxvr="14.0"/>
 
   <properties>
 
      <bandstructure>
         <plot1d>
            <path steps="100">
               <point coord="0.750   0.500   0.250" label="W"    />
               <point coord="0.500   0.500   0.500" label="L"    />
               <point coord="0.000   0.000   0.000" label="GAMMA"/>
               <point coord="0.500   0.500   0.000" label="X"    />
               <point coord="0.750   0.500   0.250" label="W"    />
               <point coord="0.750   0.375   0.375" label="K"    />
            </path>
         </plot1d>
      </bandstructure>
 
   </properties>
 
   <xs 
      xstype="BSE" 
      ngridk="4 4 4" 
      vkloff="0.097 0.273 0.493"
      ngridq="4 4 4"
      nempty="30"
      gqmax="3.0"
      broad="0.007"
      scissor="0.20947"
      tappinfo="true"
      tevout="true">
 
      <energywindow 
         intv="0.0 1.0" 
         points="1200"/>
 
      <screening 
         screentype="full"
         nempty="100"/>
 
      <BSE 
         bsetype="singlet"
         nstlbse="1 5 1 4"/>
 
      <qpointset>
         <qpoint>0.0 0.0 0.0</qpoint>
      </qpointset>
 
      <storeexcitons 
         selectenergy="false"
         MinNumberExcitons="1" 
         MaxNumberExcitons="5"
         MinEnergyExcitons="0" 
         MaxEnergyExcitons="13.8"
         />
 
   </xs>
 
</input>

The first part of the input file above is the same as used in Excited states from BSE. In addition, the new element storeexcitons appears, to enable storing the BSE coefficients. The range for storing the excitons can be defined either as an exciton index or an energy range. When choosing the index range, we need to set the attribute selectenergy="false". In that case, the attributes MinNumberExcitons and MaxNumberExcitons indicate the index of the first and last BSE solution to be stored, respectively. When choosing the energy range we select selectenergy="true" and set the attributes MinEnergyExcitons and MaxEnergyExcitons to the lower and upper energy limit, respectively. The default energy is given in eV. In this example we consider the first 5 solutions of the BSE.

ii) Running calculations

To run the calculation, copy and paste the text of the example above into the file input.xml inside the working directory LiF-exciton and start the calculation by invoking the exciting executable:

$ time excitingser &

The absorption spectra corresponding to the imaginary part of the dielectric function, can be plotted using an xml visualization template, generating agr files readable by xmgrace.

[Absorption spectrum of LiF]

3. Visualization of the exciton wavefunction

i) Preparation of the input file

As a post-processing step, we can now visualize the exciton wavefunction, both in real and in reciprocal space. For this purpose we go back to the directory LiF-exciton and append the following code to the file input.xml inside the xs block:

...
      <plan>
         <doonly task="writeexcitons"/>
         <doonly task="writekpathweights"/>
         <doonly task="excitonWavefunction"/>
      </plan>
 
      <writeexcitons 
         selectenergy="false"
         MinNumberExcitons="1" 
         MaxNumberExcitons="5"
         />
 
      <writekpathweights 
         selectenergy="false"
         MinNumberExcitons="1" 
         MaxNumberExcitons="5"
         />
 
      <excitonPlot epstol="1d-2">
         <exciton lambda="1" fix="hole"/>
         <hole>
            <plot1d>
               <path steps="1">
                  <point coord=" 0.52  0.52  0.52"/>
               </path>
            </plot1d>
         </hole>
         <electron>
            <plot3d>
               <box grid="40 40 40">
                  <origin coord=" -1.0  -1.0  -1.0"/>
                  <point  coord="  2.0  -1.0  -1.0"/>
                  <point  coord=" -1.0   2.0  -1.0"/>
                  <point  coord=" -1.0  -1.0   2.0"/>
               </box>
            </plot3d>
         </electron>
      </excitonPlot>
...

To speed up the calculation, it is recommended to skip the groundstate run, by setting do = "skip" in the element groundstate.

The attribute task = "writeexcitons" allows to print the output files containing the BSE eigenvectors whereas the attribute task = "writekpathweights" projects the excitonic weights on the band structure path given in the file bandstructure.dat that was generated in the first step. These tasks require the presence of corresponding elements writeexcitons and writekpathweights containing the same attributes as the element storeexcitons. The chosen exciton index or energy range must lie within the range given in storeexcitons.

The attribute task = "excitonWavefunction" triggers calculations to generate the real-space representation of the exciton wavefunction without repeating the diagonalization of the BSE Hamiltonian. Since the exciton wavefunction is intrinsically a six-dimensional object, it is convenient to perform the visualization for a fixed position of the hole or of the electron.

More than one exciton can be plotted in the same run, by including the element exciton for each desired solution of the BSE. Each excitation is characterized by an index lambda (reported in the first column of any file EXCITON_BSE-singlet-TDA-BAR_SCR-full_OC*.OUT file in the folder EXCITON) and a flag fix to determine whether the position of the hole or of the electron is fixed in the visualization. Inside the elements hole and electron the real-space visualization volumes are determined, following the conventions explained in How to Visualize Kohn-Sham States.

In our example, the position of the hole is fixed at one point. For this purpose, inside the hole element, the subelement plot1d is included, with steps = "1" inside path. The coordinates of the hole position are given in coord inside point. The chosen values coord = "0.52 0.52 0.52" are slightly shifted from the position of the fluorine atom. When setting the hole (or electron) coordinate, we recommend to avoid the exact atomic position where the electronic wavefunctions typically have nodes, in order to prevent spurious effects in the resulting plot.

Finally, to produce a 3D plot of the electron part of the excitonic wavefunction, in the electron element we include the subelement plot3d with the related attributes following the procedure illustrated in How to Visualize Kohn-Sham States.

Another important aspect to consider when plotting exciton wavefunctions is their effective size. The electron-hole wavefunction determines the probability density of finding, e.g., an electron at a given point for a fixed position of the hole, and, contrary to a Kohn-Sham state, is a non-periodic object. For this reason, it is often necessary to adopt supercells for a meaningful visualization in real space. In this example, we consider for the electron a 3$\times$3$\times$3 supercell surrounding a crystal unit cell. For the origin being at coord = "-1.0 -1.0 -1.0", the supercell is defined by the point elements with coord = "2.0 -1.0 -1.0", coord = "-1.0 2.0 -1.0", and coord = "-1.0 -1.0 2.0". The resolution of the calculated isosurface is controlled by the grid attribute, which should be chosen wisely to reduce the computation time.

Last but not least, the attribute epstol plays an important practical role for reducing the computational efforts. This parameter allows the user to set a tolerance for the BSE eigenvectors $A_{vck}^{\lambda}$, thereby selecting the most relevant transitions responsible for the electron-hole pair formation (see Eq.(4)).

ii) Running calculations

To run the calculation, type once again:

$ time excitingser &

A number of sub-folders and output files are now produced. In the following, we analyze them one by one. X indicates the index $\lambda$ of the BSE solution.

  • The files BEVEC_LAMBDAX.OUT inside the folder BEVEC include the list of $|A^{\lambda}_{v\,c\,\mathbf{k}}|^2$ for each k-point.
  • The files BEVEC_KSUM_LAMBDAX.OUT inside the folder BEVEC_KSUM report $\sum_\mathbf{k}|A^{\lambda}_{v\,c\,\mathbf{k}}|^2$.
  • The files EXEVEC_BSE-singlet-TDA-BAR_SCR-full_QMT001_LAMBDAX.OUT inside the folder EXCITON_EVEC include a list of $|A^{\lambda}_{v\,c\,\mathbf{k}}|$ for each k-point sorted according to their absolute value. Corresponding $Im(A^{\lambda}_{v\,c\,\mathbf{k}})$ and $Re(A^{\lambda}_{v\,c\,\mathbf{k}})$ are included as well.
  • The files WEIGHTS_BSE-singlet-TDA-BAR_SCR-full_QMT001_LAMBDAX.OUT inside the folder KPATHEXC contain the excitonic weights $w^{\lambda}_{v\mathbf{k}}$ and $w^{\lambda}_{c\mathbf{k}}$ for each k-point.
  • The files KPATH_BSE-singlet-TDA-BAR_SCR-full_QMT001_LAMBDAX.OUT inside the folder KPATHEXC contain the excitonic weights $w^{\lambda}_{v\mathbf{k}}$ and $w^{\lambda}_{c\mathbf{k}}$ projected on the band structure path given in the file bandstructure.dat.
  • The files excitonWavefunction-X-Y.xsf and excitonWavefunction-X-Y.cube contain the excitonic wavefunction isosurfaces. Their nomenclature is such that X indicates the sequential order of the exciton elements, and Y is the index of the fixed coordinate (in this case, the hole). Since in our example we plot only the electron distribution of one exciton at a single fixed position of the hole, X = 000001 and Y = 000001.
iii) Real-space visualization and analysis of the exciton wavefunction

We start analyzing the character of the exciton by inspecting the 3D plot of the electron-hole wavefunction. To do so, we can use any visualization program supporting xsf or cube formats, such as VESTA or XCrysDen. For XCrysDen, follow the instructions provided in How to Visualize Kohn-Sham States. To use VESTA, invoke the program by typing in your shell (note that you need VESTA version 3.3.2 or later to be able to reproduce the plot shown below):

$ vesta &

VESTA will open in a separate window. In the menu bar on top press File -> Open… to select the file to be shown. Double-click on excitonWavefunction-000001-000001.xsf and the structure will appear. Tune the isovalue by using the menu bar on the left side. Select Properties -> Isosurfaces and a small window will pop up. Click on the line shown in the white isovalue table in the middle of the window. Then, tune the field Isosurface level by typing 0.01. Visualize the 3$\times$3$\times$3 supercell by selecting Boundary -> New Structure and setting x(min), y(min), z(min) to -1 and x(max), y(max), z(max) to 2. Select the view along the a axis by clicking on a in the top left corner. The resulting plot will resemble the one shown below:

LiF_excitonWavefunction1.png
iv) Reciprocal-space visualization and analysis of the exciton wavefunction

We finally analyze the character of the electron-hole pair in reciprocal space. To do so, we need to visualize the files KPATH_BSE-singlet-TDA-BAR_SCR-full_QMT001_LAMBDAX.OUT containing the information on the valence and conduction bands that mostly contribute to a given exciton in LiF. To visualize the character of the first exciton, type the following lines in your bash shell:

$ cd KPATHEXC/
$ PLOT-exciton-weights-LiF.py KPATH_BSE-singlet-TDA-BAR_SCR-full_QMT001_LAMBDA000001.OUT -10 15 0.5

The first two numbers after the file name define the energy range in eV, the third defines the size of the excitonic weights which are represented by red circles. In the resulting plot you will see the excitonic weights concentrated at the valence-band top and at the conduction-band bottom:

LiF_BSE_bandstructure_exciton_0001.png

The plot is also saved in KPATH_BSE-singlet-TDA-BAR_SCR-full_QMT001_LAMBDA000001.png.

Exercises
  • Give an interpretation of the shape of the calculated electron-hole wavefunction.
  • Modify the input to visualize the exciton wavefunction corresponding to a hole located on a line connecting two neighboring Li and F atoms. Try to understand the sensitivity of the wavefunction to the choice of the origin.
iv) Visualization of core excitations

The procedure shown above can be applied also to the visualization and analysis of excitons in core-level absorption spectra. The elements storeexcitons, writeexcitons, writekpathweights and excitonPlot can be used in this case as well. However, in the reciprocal-space visualization of the excitation it should be noted that only the weights of the electronic part (Eq.(4)), namely of the conduction bands, are plotted. In the following we show how to visualize the electronic part of the core-exciton weights using the example of the spectrum from the B K-edge of cubic BN. For more information on core-level spectra, see X-ray absorption spectra using BSE.

As a first step, repeat the BSE and the band structure calculation for cubic BN in a different folder (you may use the following input file as input.xml).

[Input file to run a calculation of x-ray absorption spectrum of cubic BN]

For visualization purposes we focus on the second exciton of the spectrum which is the most intense one in the lowest-energy region of the spectrum. Following the procedure illustrated in the previous sections, we have to add the following lines to the xs element in order to store and print the weight of the first 5 excitations as well as the corresponding conduction band weights:

...
      <plan>
         <doonly task="bse"/>
         <doonly task="writeexcitons"/>
         <doonly task="writekpathweights"/>
      </plan>
 
      <storeexcitons 
         selectenergy="false"
         MinNumberExcitons="1" 
         MaxNumberExcitons="5"
         />
 
      <writeexcitons 
         selectenergy="false"
         MinNumberExcitons="1" 
         MaxNumberExcitons="5"
         />
 
      <writekpathweights 
         selectenergy="false"
         MinNumberExcitons="1" 
         MaxNumberExcitons="5"
         />
...

Now you can run again exciting for performing the visualization (you may use the following file as input.xml).

[Full input file for exciton visualization in cubic BN]

We run the calculation with the usual command:

$ time excitingser &

We can now visualize the band contributions similar to the example above: 
$ cd KPATHEXC/
$ PLOT-exciton-weights-BN.py KPATH_BSE-singlet-TDA-BAR_SCR-full_QMT001_LAMBDA000002.OUT 0 27 0.7

The resulting image will look like this:

BN_BSE_bandstructure_exciton_0002.png

Literature

  • The discussion of LiF exciton wave function can be found in, e.g., M. Gatti and F. Sottile, Phys. Rev. B 88, 155113 (2013)
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License