by Dmitrii Nabok, Christian Vorwerk, Keith Gilmore, Olga Turkina, & Caterina Cocchi for exciting neon
(Jupyter notebook by Mara Voiculescu and Martin Kuban)
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.
Table of Contents
1. Theoretical background: The electron-hole wavefunction
Read the following paragraphs before starting with the rest of this tutorial!
Before running any Jupyter tutorials, please refer to the 00_before_starting.md
document on how to correctly set up the environment. This only needs to be done once. After which, the venv can be (re)activated from exciting
's root directory:
source tools/excitingjupyter/venv/excitingvenv/bin/activate
Furthermore, this tutorial requires the usage of the following visualization packages:
which is introduced and extensively discussed in the tutorial 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
$$\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
$$\begin{align} w^{\lambda}_{v\mathbf{k}} = \sum_c|A^{\lambda}_{vc\mathbf{k}}|^2 \end{align}$$and
$$\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}}$.
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.
As a first step, create a directory named run_LiF_exciton.
%%bash
mkdir -p run_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 the tutorial 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="F.xml">
<atom coord="0.5000 0.5000 0.5000" />
</species>
<species speciesfile="Li.xml">
<atom coord="0.0000 0.0000 0.0000" />
</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.
To run the calculation, save the text of the example above into the file input.xml inside the working directory run_LiF_exciton.
Before starting the calculation, make sure to set the speciespath
:
%%bash
cd run_LiF_exciton
python3 -m excitingscripts.setup.excitingroot
cd ..
Then run the calculation:
%%bash
cd run_LiF_exciton
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ..
The calculation may take a few minutes. When it finishes, we can plot the absorption spectrum just as in the previous tutorial.
%%bash
cd run_LiF_exciton
cp EPSILON/EPSILON_BSE-singlet-TDA-BAR_SCR-full_OC11.OUT singlet-TDA
python3 -m excitingscripts.plot.files -f singlet-TDA -lx 'Energy [eV]' -ly 'Im $\varepsilon_M$' -t 'Macroscopic dielectric function' -g -rc -cy 3 -x 0 27 -nl
cd ..
The result should look like the following.
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. 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.
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 run_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 ground-state run by inserting the attribute do="skip"
inside the groundstate
element.
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 the tutorial 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×3×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.
%%bash
cd run_LiF_exciton
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ..
A number of sub-folders and output files are now produced. In the following, we analyze them one by one. X indicates the index λ of the BSE solution.
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.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 the tutorial 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):
%%bash
cd run_LiF_exciton
vesta
cd ..
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×3×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:
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, execute the following commands:
%%bash
cd run_LiF_exciton
cp KPATHEXC/KPATH_BSE-singlet-TDA-BAR_SCR-full_QMT001_LAMBDA000001.OUT LiF-exciton-1
python3 -m excitingscripts.plot.exciton_weights LiF LiF-exciton-1 -10 15 0.5
cd ..
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.
In this plot you see the excitonic weights concentrated at the valence-band top and at the conduction-band bottom.
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, 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 the tutorial 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:
%%bash
mkdir -p run_BN_exciton
You may use the following input file as input.xml for the initialization run.
<input>
<title>Cubic boron nitride</title>
<structure autormt="true" speciespath="$EXCITINGROOT/species/">
<crystal scale="6.83136">
<basevect>0.5 0.0 0.5</basevect>
<basevect>0.0 0.5 0.5</basevect>
<basevect>0.5 0.5 0.0</basevect>
</crystal>
<species speciesfile="B.xml" >
<atom coord="0.00 0.00 0.00" />
</species>
<species speciesfile="N.xml">
<atom coord="0.25 0.25 0.25" />
</species>
</structure>
<groundstate
ngridk="4 4 4"
xctype="GGA_PBE_SOL"
gmaxvr="14.0" />
<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>
<xs
xstype="BSE"
ngridk="2 2 2"
vkloff="0.05 0.15 0.25"
ngridq="2 2 2"
nempty="30"
gqmax="3.0"
broad="0.018"
tevout="true" >
<energywindow
intv="13.4 16.5"
points="1500" />
<screening
screentype="full"
nempty="100" />
<BSE
xas="true"
xasspecies="2"
xasatom="1"
xasedge="K"
bsetype="singlet"
nstlxas="1 15" />
<qpointset>
<qpoint>0.0 0.0 0.0</qpoint>
</qpointset>
</xs>
</input>
After updating the speciespath
with
%%bash
cd run_BN_exciton
python3 -m excitingscripts.setup.excitingroot
cd ..
run exciting
with the usual command:
%%bash
cd run_BN_exciton
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ..
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"
/>
...
Now, run exciting
again to generate files for the visualization with the usual command:
%%bash
cd run_BN_exciton
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ..
We can now visualize the band contributions similar to the example above:
%%bash
cd run_BN_exciton
cp KPATHEXC/KPATH_BSE-singlet-TDA-BAR_SCR-full_QMT001_LAMBDA000002.OUT BN-exciton-2
python3 -m excitingscripts.plot.exciton_weights BN BN-exciton-2 0 27 0.7
cd ..
The resulting image (BN-exciton-2.png) will look like this: