How to Visualize Kohn-Sham States

by Caterina Cocchi, Dmitrii Nabok, & Pasquale Pavone for exciting nitrogen

Purpose: In this tutorial you will learn how to visualize Kohn-Sham states (orbitals) with exciting.

---

0. Define relevant environment variables

Read the following paragraphs before starting with the rest of this tutorial!

Before starting, be sure that relevant environment variables are already defined as specified in How to set environment variables for tutorials scripts. Here, you can find a list of the scripts which are relevant for this tutorial, with a short description.

• EXECUTE-single.sh: (Bash) shell script for running a single exciting calculation.

From now on the symbol $ will indicate the shell prompt.

Requirements: Bash shell and lxml library.

Important note: All input parameters are in atomic units!

---

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 KS-silicon and we move inside it.

$ mkdir KS-silicon
$ cd KS-silicon

Here, we create an exciting (xml) 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

$ SETUP-excitingroot.sh

To start the calculation, you may run the script EXECUTE-single.sh.

$ EXECUTE-single.sh test-GS

ii) 3D plot: Calculation of KS states

After you have completed the ground-state run, move to the directory test-GS, where the results are stored.

$ cd test-GS

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

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

$ time excitingser

This time the execution of exciting will produce as a result the wavefunction of the investigated KS state, stored inside the file WF3D.xml.

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

$ xsltproc $EXCITINGVISUAL/plot3d2xsf.xsl WF3D.xml > WF3D.xsf

In order to visualize the underlying structure of the molecule, we transform the atomic coordinates of the input file into a xsf file:

$ xsltproc $EXCITINGCONVERT/xmlinput2xsf.xsl input.xml > input.xsf

Finally, we merge the two files as follows:

$ cat input.xsf WF3D.xsf > PLOT3D-1-4.xsf

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 typing on the shell:

$ xcrysden --xsf PLOT3D-1-4.xsf >/dev/null 2>&1 &

Now, a 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:

1. Press on the main window Tools -> Data Grid. A new window appears.
2. Press OK: The previous window will disappear and a new window will appear with the control panel of the data grid to be plotted.
3. 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 (obtained using XCrySDen).

---

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

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:

$ xsltproc $EXCITINGVISUAL/plot2d2xsf.xsl WF2D.xml > PLOT2D-1-4.xsf

[How to visualize 2D-plots in xsf format using XCrySDen]

In order to visualize 2D-plots using XCrySDen, we can simply type on the shell:

$ xcrysden --xsf PLOT2D-1-4.xsf >/dev/null 2>&1 &

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:

1. choose BLUE-WHITE-RED in Select color basis;
2. choose LOG10 in Select scale function;
3. press Submit button;
4. 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;
5. to do so, press Display -> Coordinate System in order to visualize the coordinate system;
6. orient the view to the (x,z) plane, by pressing the corresponding button on the right column of the XCrySDen main window;
7. 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, obtained with XCrySDen looks like this:

This plot represents a projection of the highest occupied KS state, already visualized as a 3D isosurface in the previous section, on a plane containing the two inequivalent atoms in the unit cell. The red region in the middle of the plot indicates the localization of the probability density of the highest occupied KS state along the bond between the two Si atoms. Notice that an intense region of probability density for this state can be observed in the vicinity of the atoms themselves, which are placed at the bottom right and top left vertices of the visualized plane.

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 Electronic-structure calculations). 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 path subelement to ensure a good resolution.

After running exciting, a file WF1D.xml is printed. To convert the xml file into the agr format, which can be read by xmgrace, we execute the script:

$ xsltproc $EXCITINGVISUAL/plot1d2agr.xsl WF1D.xml > PLOT1D-1-4.agr

To visualize the file PLOT1D-1-4.agr with xmgrace, just type:

$ xmgrace PLOT1D-1-4.agr >/dev/null 2>&1 &

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 elemets 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:

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.

---

4. Links to other visualization tutorials

• The visualization of the desity of states (DOS) and of the electronic bandstructure of a crystal is discussed in Electronic-structure calculations.
• In the last section of How to run calculations for simple molecules you can learn how to calculate and visualize molecular orbitals, with explicit application to the HOMO and LUMO states.
• In the application tutorial Fermi surface plot you can learn how to calculate and visualize the Fermi surface for metallic systems.
• Many other specific visualization scripts are used in the other exciting Tutorials.