by Caterina Cocchi, Dmitrii Nabok, & Pasquale Pavone for exciting oxygen
Purpose: In this tutorial you will learn how to visualize Kohn-Sham states (orbitals) with exciting.
Table of Contents
|
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.
- CONVERT-xml2xsf.sh: (Bash) shell script for converting xml to xsf files.
- PLOT-files.py: Python visualization tool (see The python script "PLOT-files.py").
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
- 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
$ time exciting_smp
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.
$ CONVERT-xml2xsf.sh 3D WF3D.xml
In order to visualize the underlying structure of the molecule, we transform the atomic coordinates of the input file into a xsf file:
$ CONVERT-xml2xsf.sh 3D input.xml
Finally, we merge the two output xsf 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).
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:
$ CONVERT-xml2xsf.sh 2D WF2D.xml PLOT2D-1-4.xsf
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 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 path subelement to ensure a good resolution.
After running exciting, a file WF1D.xml is printed. To convert the xml file into a printable (agr) format, we execute the script:
$ xsltproc $EXCITINGROOT/xml/visualizationtemplates/plot1d2agr.xsl WF1D.xml > PLOT1D-1-4.agr
To visualize the file PLOT1D-1-4.agr you can use the script PLOT-files.py (for a detailed description of the script arguments see The python script "PLOT-files.py") as follows
$ PLOT-files.py -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
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:
|
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 band-structure and density of states.
- 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.