Exciton Analysis and Visualization

by Dmitrii Nabok, Christian Vorwerk, & Caterina Cocchi, for exciting carbon

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. In particular, the relevant script for this tutorial is

• bandstructure_exciton.py: Python script for the visualization of bandstructure contribution to the exciton.

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}\rangleenter 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 coefficientsA^{\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 spectra, the calculation needs to be converged with respect to several parameters, such as 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.xml 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"
scissor="0.20947"
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 MinNumberExcitons="1" MaxNumberExcitons="5"/>

<writeexcitons MinNumberExcitons="1" MaxNumberExcitons="5"/>

</xs>

</input>


The first part of the input file above is the same used in Excited states from BSE. In addition, the new elements storeexcitons and writeexcitons appear, to enable storing and printing the BSE coefficients. The attributes MinNumberExcitons and MaxNumberExcitons indicate in both case the index of the first and last BSE solution to be stored, respectively. In this example we consider only 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 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. ### 2. 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 append the following code inside the xs block: ... <plan> <doonly task="writebevec"/> <doonly task="excitonWavefunction"/> </plan> <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". The attribute task = "writebevec" allows to print the output files containing the BSE eigenvectors. 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_BSEsinglet_*.OUT file) 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 submelement 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 to 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 output files are now produced. In the following, we analyze them one by one.

• The files BEVEC_LAMBDAXXXX.OUT include the list of $|A^{\lambda}_{v\,c\,\mathbf{k}}|^2$ for each k-point. XXXX is the index $\lambda$ of the BSE solution.
• The files BEVEC_KSUM_LAMBDAXXXX.OUT report $\sum_\mathbf{k}|A^{\lambda}_{v\,c\,\mathbf{k}}|^2$. XXXX is the index $\lambda$ of the BSE solution.
• The file bandstructure.dat contain the information on the band structure, that are needed for the analysis of the exciton in reciprocal space.
• The files exciton_evec_XXXX.dat contain the information on the exciton weights, $A^{\lambda}_{v\,c\,\mathbf{k}}$ that are needed for the analysis of the exciton in reciprocal space. XXXX is the index $\lambda$ of the BSE solution.
• The files excitonWavefunction-XX-YYYYYY.xsf and excitonWavefunction-XX-YYYYYY.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 = 01 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

$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-01-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 white isovalue table in the midddle of the window. Then, tune the field Isosurface level by typing 0.01. The resulting plot will resemble the one shown below: ##### 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 post-process the files bandstructure.dat and exciton_evec_XXXX.dat with the python script bandstructure_exciton.py. To visualize the character of the first exciton, type the following line in your bash shell: $ bandstructure_exciton.py bandstructure.dat exciton_evec_0001.dat > plot_exciton_0001.dat


The resulting file plot_exciton_0001.dat contains the information on the valence and conduction bands that mostly contribute to the lowest-energy exciton in LiF. This file can be opened with xmgrace, in order to do this, follow the instructions given below.

The final result is plot like the following. ##### 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, 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 calculation for cubic BN, (you may use the following input file as input.xml).

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, first we calculate the band structure. To adapt the calculation to cubic BN, we change the element properties in the input file to:

...
<properties>
<bandstructure>
<plot1d>
<path steps="100">
<point coord="0.500   0.000   0.500" label="X"    />
<point coord="0.500   0.250   0.750" 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.000   0.500" label="X"    />
</path>
</plot1d>
</bandstructure>
</properties>
...


In addition, we have to add the following lines to the xs element in order to store and print the weight of the first 5 excitations:

...
<plan>
</plan>

<storeexcitons MinNumberExcitons="1" MaxNumberExcitons="5"/>
<writeexcitons MinNumberExcitons="1" MaxNumberExcitons="5"/>
...


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

We run the calculation with the usual command:

$time excitingser &  Once the calculation is finished, we use the files bandstructure.dat and exciton_evec_0002.dat to visualize the band contributions. To do so, we run the corresponding script by typing the following bash command line: $ bandstructure_exciton.py bandstructure.dat exciton_evec_0002.dat xas > plot_exciton_0002.dat


As you can see the keyword xas has been added here in order to enable the visualization of the contributions to the conduction states only. We can now visualize the band contributions similar to the example above. If you want to use xmgrace, you can download here the parameter file bands_BN_exciton.par.

The resulting image will look like this: ### Literature

• The discussion of LiF exciton wave function can found in, e.g., M. Gatti and F. Sottile, Phys. Rev. B 88, 155113 (2013)