by Caterina Cocchi, Dmitrii Nabok, & Pasquale Pavone for exciting neon
(Jupyter notebook by Mara Voiculescu & Martin Kuban)
Purpose: In this tutorial you will learn how to visualize Kohn-Sham states (orbitals) with exciting
.
Table of Contents
1. Example: 3D Plot of Kohn-Sham states in bulk silicon
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
In this example, we calculate and visualize the probability density of the highest occupied Kohn-Sham (KS) state at the Γ point of the electronic bandstructure of bulk silicon. To this end, we start by creating a working directory run_KS_silicon and another subdirectory for the ground-state calculation test_GS.
%%bash
mkdir -p run_KS_silicon
cd run_KS_silicon && mkdir -p test_GS
Here, we create an exciting
input file called input.xml corresponding to a ground-state (GS) calculation for bulk silicon, which should appear as the one below.
<input>
<title>Bulk Silicon: Plot example</title>
<structure speciespath="$EXCITINGROOT/species">
<crystal scale="10.26">
<basevect>0.0 0.5 0.5</basevect>
<basevect>0.5 0.0 0.5</basevect>
<basevect>0.5 0.5 0.0</basevect>
</crystal>
<species speciesfile="Si.xml">
<atom coord="0.00 0.00 0.00"></atom>
<atom coord="0.25 0.25 0.25"></atom>
</species>
</structure>
<groundstate
do="fromscratch"
ngridk="8 8 8"
gmaxvr="14"
xctype="LDA_PW">
</groundstate>
</input>
Make sure to set $EXCITINGROOT to the correct exciting
root directory in the speciespath
attribute using the command
%%bash
cd run_KS_silicon/test_GS
python3 -m excitingscripts.setup.excitingroot
cd ../..
Now you can start the calculation by executing the following script
%%bash
cd run_KS_silicon
python3 -m excitingscripts.execute.single -r test_GS
cd ..
In order to calcualte explicitly the KS state in which we are interested, we have to edit the input.xml in the directory test_GS and
do
to "skip"
inside the element groundstate
;properties
after the groundstate
element as shown in the following:...
<properties>
<wfplot>
<kstlist>
<!-- (k-point index) (state index)-->
<pointstatepair>1 4</pointstatepair>
</kstlist>
<!-- WF in the crystal unitcell-->
<plot3d>
<box grid="20 20 20">
<origin coord="-0.5 -0.5 -0.5"/>
<point coord=" 0.5 -0.5 -0.5"/>
<point coord="-0.5 0.5 -0.5"/>
<point coord="-0.5 -0.5 0.5"/>
</box>
</plot3d>
</wfplot>
</properties>
...
The properties
block shown above contains the element wfplot
which allows for the calculation of the KS wavefunctions. The subelement kstlist
defines the quantum numbers of the desidered KS state:
...
<kstlist>
<pointstatepair>1 4</pointstatepair>
</kstlist>
...
In the above case, e.g., the k-point of the KS state is the one in position 1 in the list of k-points reported inside the file KPOINTS.OUT. Furthermore, the band index 4 corresponds to the highest occupied KS state at the k-point specified by the first index in the element pointstatepair
. This can be verified by direct inspection of the file EIGVAL.OUT.
The next important element is plot3d
, which defines the space region where the KS has to be calculated.
...
<plot3d>
<box grid="20 20 20">
<origin coord="-0.5 -0.5 -0.5"/>
<point coord=" 0.5 -0.5 -0.5"/>
<point coord="-0.5 0.5 -0.5"/>
<point coord="-0.5 -0.5 0.5"/>
</box>
</plot3d>
...
The additional element box
is specified by the attribute grid
, which defines the quality of the plot: The larger number of points in this mesh, the more resolved the resulting plot. Be aware that the increase of the number of grid points impacts on the computational costs! By means of the other elements origin
and point
one sets, respectively, the origin and the vectors, expressed in lattice coordinates, defining the space region where the wavefunction is calculated.
Finally, including all changes, the new input.xml will be similar to:
<input>
<title>Bulk Silicon: Plot3D example</title>
<structure speciespath="$EXCITINGROOT/species">
<crystal scale="10.26">
<basevect>0.0 0.5 0.5</basevect>
<basevect>0.5 0.0 0.5</basevect>
<basevect>0.5 0.5 0.0</basevect>
</crystal>
<species speciesfile="Si.xml">
<atom coord="0.00 0.00 0.00"></atom>
<atom coord="0.25 0.25 0.25"></atom>
</species>
</structure>
<groundstate
do="skip"
ngridk="8 8 8"
gmaxvr="14"
xctype="LDA_PW">
</groundstate>
<properties>
<wfplot>
<kstlist>
<pointstatepair>1 4</pointstatepair>
</kstlist>
<plot3d>
<box grid="20 20 20">
<origin coord="-0.5 -0.5 -0.5"/>
<point coord=" 0.5 -0.5 -0.5"/>
<point coord="-0.5 0.5 -0.5"/>
<point coord="-0.5 -0.5 0.5"/>
</box>
</plot3d>
</wfplot>
</properties>
</input>
Now, you can run exciting
again
%%bash
cd run_KS_silicon/test_GS
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ../..
This time the execution of exciting
will produce as a result the wavefunction of the investigated KS state, stored inside the file WF3D.xml.
The file WF3D.xml contains all information necessary for the visualization of the valence-band-top KS state. In order to visualize this state, you have to convert the WF3D.xml file to the xsf format. In order to visualize the underlying structure of the molecule, we transform also the atomic coordinates of the input file into a xsf file:
%%bash
cd run_KS_silicon/test_GS
python3 -m excitingscripts.convert_xml2xsf -f input.xml -d 3D
python3 -m excitingscripts.convert_xml2xsf -f WF3D.xml -d 3D
cd ../..
Finally, we merge the two output xsf files as follows:
%%bash
cd run_KS_silicon/test_GS
cat input.xsf WF3D.xsf > PLOT3D-1-4.xsf
cd ../..
The resulting file PLOT3D-1-4.xsf (the labels 1 and 4 identify the KS state) can be visualized using standard tools like XCrySDen or VESTA.
In order to visualize xsf files using XCrySDen, we have to make sure that XCrySDen is correctly initialized on your workstation. For further instructions and details about initialization, refer to XCrySDen Setup for exciting. Once XCrySDen is initialized, the file PLOT3D-1-4.xsf can be visualized directly by running the following command:
%%bash
cd run_KS_silicon/test_GS
xcrysden --xsf PLOT3D-1-4.xsf >/dev/null 2>&1 &
cd ../..
Now, an XCrySDen window will appear on the screen, showing the atomic structure, which in this example is the bulk silicon. If you want to visualize the isosurface, use the following procedure:
An example for the resulting image is the following.
As a further example of visualization, we will compute and plot a two-dimensional projection of the highest occupied state of bulk Silicon on a crystal plane. For this purpose, we need to modify the properties
block in the input file, by replacing the sub-element plot3d
with plot2d
, as follows:
...
<plot2d>
<parallelogram grid="200 200">
<origin coord="0 0 0"/>
<point coord="0.25 0.25 -0.25"/>
<point coord="0.00 0.00 0.50"/>
</parallelogram>
</plot2d>
...
%%bash
cd run_KS_silicon/test_GS
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ../..
When the plot2d
element is specified, the grid
is defined inside a parallelogram
. To ensure a well resolved visualization, 200 grid points per direction are chosen here. We place the coordinate of the origin
in (0, 0, 0), which coincides with the position of one of the two atoms in the unit cell of bulk silicon. Finally, we project the KS state on a rectangular plane, including on opposite vertices the two atoms in the unit cells, which is spanned by the vectors (1/4, 1/4, -1/4) and (0, 0, 1/2).
After running the calculation, a file WF2D.xml will be printed. Again, we can convert it into xsf format by executing the follow script:
%%bash
cd run_KS_silicon/test_GS
python3 -m excitingscripts.convert_xml2xsf -f WF2D.xml -d 2D
cd ../..
In order to visualize 2D-plots using XCrySDen, we can simply execute the following command:
%%bash
cd run_KS_silicon/test_GS
xcrysden --xsf WF2D.xsf >/dev/null 2>&1 &
cd ../..
An empty XCrySDen window will appear on the screen. As for the case of 3D plots, to visualize the data press Tools -> Data Grid -> OK. A new window appears, containing the setup tools for the visualization of the plot. To obtain the graph shown in the figure below, choose the following parameters in this window:
The corresponding plot should like this:
As a final example in this tutorial, we will calculate the 1D plot of the probability density of the highest occupied KS state along the bond between the two inequivalent atoms in the Si unit cell. Here, we apply the element plot1d
to wfplot
, as done already for plot3d
and plot2d
in the previous sections. In order to calculate this plot, we have to insert the following block in the input file:
...
<plot1d>
<path steps="1000">
<point coord="0.00 0.00 0.00"/>
<point coord="0.25 0.25 0.25"/>
</path>
</plot1d>
...
A similar block is present also in the input file for band-structure calculations (see tutorial Electronic band-structure and density of states). Inside the element plot1d
two or more points are indicated, in order to define the path along which the desired quantity (in the case the probability density of the highest occupied KS state at Γ) is plotted. Since we are interested in visualizing such probability density along the bond between the two atoms in the unit cell, we assign to the element point
the coordinates of the atoms, represented here in lattice units: (0, 0, 0) and (1/4, 1/4, 1/4)). Furthermore, we choose 1000 steps in the point
subelement to ensure a good resolution.
%%bash
cd run_KS_silicon/test_GS
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ../..
After running exciting, a file WF1D.xml is printed. To convert the xml file into a printable (agr) format, we execute the command:
%%bash
cd run_KS_silicon/test_GS
xsltproc $EXCITINGROOT/xml/visualizationtemplates/plot1d2agr.xsl WF1D.xml > PLOT1D-1-4.agr
cd ../..
To visualize the file PLOT1D-1-4.agr, you can execute the script excitingscripts.plot.files
as follows.
%%bash
cd run_KS_silicon/test_GS
python3 -m excitingscripts.plot.files -f PLOT1D-1-4.agr -lx 'Distance [bohr]' -ly 'Probability density [$10^{-3}$]' -x 0 4.44 -ys 1000 -t 'Bulk Silicon' -nl -mtx 5 -mty 5
cd ../..
The resulting plot will look like this:
In this plot, the abscissa indicates the distance along the bond from the atom sitting at the origin. Pay attention that this distance is expressed in Bohr. In the y axis, the probability density of the plotted wave function is indicated. We notice the distribution of the highest occupied KS state along the bond between the atoms. This corresponds exactly to the isosurface and to the contour plot visualized, respectively, as 3D and 2D plots in the previous sections.
The elements plot1d
, plot2d
, and plot3d
can also be used to visualize other quantities than the probability density of a KS state. In the table below some of the main examples are listed:
plot*d |
Elements in properties inside which plot*d can be used |
---|---|
plot1d |
wfplot , chargedensityplot , bandstructure |
plot2d |
wfplot , chargedensityplot , fermisurfaceplot , electricfield , stm |
plot3d |
wfplot , chargedensityplot , fermisurfaceplot , electricfield |