How to Visualize Kohn-Sham States

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

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.

An example for the resulting image is the following (obtained using XCrySDen).

plot3D_xcrysden.png

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
The corresponding plot, obtained with XCrySDen looks like this:
plot2d_log10.png
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:
plot1d-1-4.png

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.
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License