Exciton Analysis and Visualization

by Dmitrii Nabok, Christian Vorwerk, Keith Gilmore, Olga Turkina, & Caterina Cocchi, for exciting oxygen

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-files.py: Python visualization tool for multiple files with the same name in different directories (for more details, see The python script "PLOT-files.py").
  • 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:

The symbol $ will indicate the shell prompt.

Requirements: Bash shell, python, XCrystDen, VESTA.

Prerequisites: Excited states from BSE.

1. Theoretical background: The electron-hole wavefunction

☛ The electron-hole wavefunction

2. Absorption spectrum of LiF

You should have already calculated the absorption spectrum of LiF in the tutorial Excited states from BSE. Here, we repeat this calculation and analyse the character of the electron-hole wavefunction that corresponds to the first and the most intense excitonic peak.

i) Preparation of the input file

Create a directory named LiF-exciton and move into it.

$ mkdir LiF-exciton
$ cd LiF-exciton

Below is an example input file for the LiF ground-state calculation, as well as for the calculations of the optical absorption spectrum of LiF. Notice that within the groundstate block the properties element is present in order to calculate the band structure of the system (for details, see Electronic band-structure and density of states).

<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="8 8 8"
      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="3 3 3" 
      vkloff="0.097 0.273 0.493"
      ngridq="3 3 3"
      nempty="30"
      gqmax="2.5"
      broad="0.007"
      scissor="0.20947"
      tappinfo="true"
      tevout="true">
 
      <energywindow 
         intv="0.0 1.0" 
         points="1200"/>
 
      <screening 
         screentype="full"
         nempty="60"/>
 
      <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 the calculation

To run the calculation, copy and paste the text of the example above into the file input.xml inside the working directory LiF-exciton.

Before starting the calculation, make sure to set the species path:

$ SETUP-excitingroot.sh

Then run the calculation:

$ time exciting_smp

The calculation may take a few minutes. When it finishes, we can plot the absorption spectrum just as in the previous tutorial.

$ cp EPSILON/EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT singlet-TDA
$ PLOT-files.py -f singlet-TDA  -lx 'Energy [eV]'  -ly 'Im $\varepsilon_M$'  -t 'Macroscopic dielectric function'  -g  -rc  -cy 3  -x 0 27  -nl
epsilon_Im.png

The result is again stored in PLOT.png and PLOT.eps.

3. Visualization of the exciton wavefunction

To analyze the excitonic states we inspect the output files contained in the folder EXCITON. This directory contains the solutions of the diagonalization of the BSE. Specifically, the files EXCITON_BSE-singlet-TDA-BAR_SCR-full_OCXX.OUT, where XX=11,22,33 include the relevant information about the excitation energies and intensities. Each file has six columns. After the progressive index referring to the excitation number, we find the excitation energies, exciton binding energies, oscillator strengths, and finally the real and the imaginary part of the transition coefficients from which the oscillator strengths are computed (see also SAD-2009). An important quantity is the exciton binding energy, which is defined as the difference between the lowest excitation energy and the band gap. This quantity is printed in the file header and labeled as E_shift. Negative values of the exciton binding energies are associated with bound excitons.

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 stay in directory LiF-exciton.

Now, append the following to the file input.xml inside the xs block:

...
      <plan>
         <doonly task="writebevec"/>
         <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 printing of 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 exciting_smp

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.0037, then press OK. Visualize the 3$\times$3$\times$3 supercell by selecting Boundary and changing in the window which will appear XCrySDen XSF file to XCrySDen XSF structure data. In the same window, set x(min), y(min), z(min) to -1 and x(max), y(max), z(max) to 2, then press again OK. Select the view along the a axis by clicking on a in the top left corner. Finally unclick Show session. The resulting plot will resemble the one shown below:

LiF_excitonWavefunction.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
$ cp KPATH_BSE-singlet-TDA-BAR_SCR-full_QMT001_LAMBDA000001.OUT LiF-exciton-1
$ PLOT-exciton-weights-LiF.py LiF-exciton-1 -10 15 0.5

In the last command line above, 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. The plot is saved in LiF-exciton-1.png, which is shown below.

LiF-exciton-1.png

In this plot you see the excitonic weights concentrated at the valence-band top and at the conduction-band bottom.

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:

$ cd ../..
$ mkdir BN-exciton
$ cd BN-exciton

You may use the following input file as input.xml for the initialization run.

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

After updating the speciespath with

$ SETUP-excitingroot.sh

run exciting with the usual command:

$ time exciting_smp

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 should 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="writebevec"/>
         <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"
         />
...

The full input file for producing visualization files is given in the following.

☛ Input file for producing visualization files

Now, run exciting again to generate files for the visualization with the usual command:

$ time exciting_smp

We can now visualize the band contributions similar to the example above:

$ cd KPATHEXC
$ cp KPATH_BSE-singlet-TDA-BAR_SCR-full_QMT001_LAMBDA000002.OUT BN-exciton-2
$ PLOT-exciton-weights-BN.py BN-exciton-2 0 27 0.7

The resulting image (BN-exciton-2.png) will look like this:

BN-exciton-2.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