How to Visualize Kohn-Sham States¶

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.

source tools/excitingjupyter/venv/excitingvenv/bin/activate

0. Before Starting¶

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

1. Example: 3D Plot of Kohn-Sham states in bulk silicon¶

i) Ground-state calculation¶

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 ..

ii) 3D plot: Calculation of KS states¶

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

modify the value of the attribute do to "skip" inside the element groundstate;
add the element 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.

ii) 3D plot: Visualization of KS states¶

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.

How to Visualize 3D-Plots in xsf Format using XCrySDen¶

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:

Press on the main window Tools -> Data Grid. A new window appears.
Press OK: The previous window will disappear and a new window will appear with the control panel of the data grid to be plotted.
In the middle of the left column a box appears, indicating Minimum grid value, Maximum grid value, and Isovalue next to an empty box. In this box the user has to type the isovalue with which the isosurface is to be plotted. As a rule of a thumb, 10% of the Maximum grid value is a reasonable value. Type the isovalue (e.g., 0.00002 in this case) and press the Submit button on the bottom right of the window. The isosurface representing the KS state will now appear on the main window.

An example for the resulting image is the following.

2. Other Types of Plots¶

i) 2D Plot¶

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 ../..

How to Visualize 2D-Plots in xsf Format using XCrySDen¶

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:

Choose BLUE-WHITE-RED in Select color basis;
Choose LOG10 in Select scale function;
Press the Submit button;
Do not worry if you do not see anything appearing: Since the default visualization plane is (x,y), we have to rotate the view in order to properly orient the plot in the XCrySDen window;
To do so, press Display -> Coordinate System in order to visualize the coordinate system;
Orient the view to the (x,z) plane, by pressing the corresponding button on the right column of the XCrySDen main window;
Finally, rotate around the y axis by pressing the button ROT -Y on the right column of the XCrySDen main window, until the x and y arrows of the coordinate system overlap completely.

The corresponding plot should like this:

ii) 1D Plot¶

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.

iii) Plot Other Physical Quantities¶

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`

Exercise¶

Visualize a different KS state. Use 1D, 2D, and 3D visualization tools.
For the example presented in this tutorial, calculate and visualize some of the physical properties reported in the above table. For further details on each properties explore Input Reference.