2015-04-28

# 1. Introduction

Let us start by stating what we would like to obtain. We wish to visualise the electron density of the molecular orbitals of CO using VASP. We wish to create a surface cut rather than generating a 3D image containing an isosurface, as the former more clearly conveys the electron density distribution. In this short tutorial, I will show you how to create the image below: Figure 1: The electron density plots of the molecular orbitals of CO. The yellow line in the lower-left image is added to indicate the nodal plane.

The molecular orbitals of CO are formed from the atomic orbitals of C and O. The molecular orbitals of CO can be identified from the molecular orbital diagram of CO as shown below. Each of these orbitals have specific symmetry characteristics and these symmetry characteristics give rise to the name of the molecular orbitals. The molecular orbital lowest in energy, which is a non-bonding molecular orbital that only contains core electrons, is termed the 1σ orbital. The σ denotes the symmetry characteristics of the orbital. σ-symmetry means that the orbital can be rotated arbitrarily around the bonding axis of C and O. Hence, there are an infinite number of mirror planes parallel to this rotation axis. The four molecular orbitals lowest in energy all have this particular symmetry. The fifth molecular orbital however, has π- symmetry and is hence denoted as the 1π molecular orbital. π-symmetric orbitals have a nodal plane on the bonding axis and hence only have a C2 rotation operation. Consequently, they also have only two mirror planes parallel to this bonding axis. The electron density plots shown above are in agreement with these symmetry characteristics. Figure 2: The occupied orbitals of CO, that are constructed from the valence atomic orbitals of C and O, are in order of low to high energy: 3σ, 4σ, 1π and 5σ.

# 2. VASP calculations

In order to construct the electron density, we first have to perform a series of VASP calculations.

• Perform a structure optimisation
• Perform a single point calculation and calculate the DOS for the final state (printed in DOSCAR)
• Identify the molecular orbitals using the DOS plot (select the "energy range")
• Perform for each identified molecular orbital another calculation where the partial charge corresponding to the energy range is saved (printed in PARCHG)

## 2.1 POSCAR and INCAR files

Our POSCAR files looks as follows:


CO
1.0000000000000000
10.0000000000000000    0.0000000000000000    0.0000000000000000
0.0000000000000000   10.0000000000000000    0.0000000000000000
0.0000000000000000    0.0000000000000000   10.0000000000000000
1   1
Direct
0.5000000000000000  0.5000000000000000  0.500000
0.5000000000000000  0.5000000000000000  0.617000



and the INCAR file as:


SYSTEM = CO

LREAL=.FALSE.
LWAVE=  .TRUE.
LCHARG=  .TRUE.

PREC = normal
ENCUT = 400

ISMEAR = 0
SIGMA = 0.05

EDIFF = 1E-5
EDIFFG = -1E-4
NSW = 100

ISIF = 2
IBRION = 2
ISPIN = 1
POTIM = 0.5

IALGO = 38
NELMIN = 6
NELM = 40

EMIN = -30
EMAX = 15
NEDOS = 4500

NGX = 150
NGY = 150
NGZ = 150



The KPOINTS file specifies the gamma-point (1x1x1 k-points) and the POTCAR file corresponds to one with the potentials of C and O (in that order).

I am not going to discuss all the directives in detail (the VASP manual is really good at that), but I will explain a couple of directives. The INCAR file tells the program to optimise the geometry (IBRION = 2). The WAVECAR and CHGCAR files will be generated and a DOSCAR file will be made with a resolution of 4500 data points on the interval [-30 eV;15 eV]. The mesh of the CHGCAR is 150x150x150 grid points.

After the optimisation run, change NSW to 0 and set ISTART = 1. Copy the CONTCAR to POSCAR and start the calculation again. (this way, we have the DOS of the final state and not an averaged DOS). The generated DOSCAR file has 6 header lines and then 4500 lines containing the DOS and the integrated DOS as a function of the energy. To extract that part (and ignore the header lines) and save it DOS-total, we run:


sed -n '7,4506p' DOSCAR > DOS_total



This file can be readily visualised using your favourite program (Excel, Origin, GnuPlot, your pick...). I chose Origin and made the following graph: Figure 3: The density of states (DOS) and integrated DOS of CO. Each of the peaks corresponds to a molecular orbital.

In this graph, I have depicted the molecular orbitals 3σ, 4σ, 1π and 5σ. The corresponding energy intervals (in eV) are (roughly):

 3σ 4σ -30 -27 -15 -13 -13 -10 -10 -5

For each of the energy intervals, we are going to run VASP again with the following example INCAR file. This INCAR file is for the energy interval [-30; -27.5], but you can chose any interval you like. The LPARD directive instructs the program to generate a PARCHG file of the partial charge corresponding to the wavefunctions between -30 and -27.5 eV.


SYSTEM = CO

LREAL=.FALSE.
LWAVE=  .TRUE.
LCHARG=  .TRUE.

PREC = medium
ENCUT = 400

ISMEAR = 0
SIGMA = 0.05

EDIFF = 1E-5
EDIFFG = -1E-4
NSW = 100

ISIF = 2
IBRION = 2
ISPIN = 1
POTIM = 0.5

IALGO = 38
NELMIN = 6
NELM = 40

EMIN = -30
EMAX = 15
NEDOS = 4500

NGX = 150
NGY = 150
NGZ = 150

LPARD = .TRUE.
EINT = -30 -27.5



The resulting PARCHG files are stored separately and will be used for the visualisation in the next section.

# 3. Visualising electron density

To visualise the electron density, we are going to use the EDP program, which is freely available via its Github repository. You can compile the program yourself or obtain it from the repository (if you are running Debian).

## 3.1 Visualising electron density with EDP

To use the mathematical terms, the electron density is a scalar field. In other words, it is a scalar value that depends on the position in the unit cell. Although the electron density is a smooth function, it is stored at particular grid points in the CHGCAR and PARCHG file. The resolution of that grid can be set by specifying NGX, NGY and NGZ in the INCAR file. The purpose of EDP is to plot the electron density on a so-called cutting plane. The end result is therefore a heat map (as shown above). In order to specify the cutting plane, the user is required to provide two vectors and one point. In our example, CO is placed along the z-axis at the centre of a cubic unit cell of 10x10x10 angstrom. The centre of the unit cell is at position (5.0, 5.0, 5.0). If we want to project the electron density on an (xz)- plane cutting through the point (5.0, 5.0, 5.0), then the command is:


./bin/edp -i path/to/CHGCAR -v 1,0,0 -w 0,0,1 -p 5,5,5 -s 100 -o test.png



The tags -v and -w specify the vectors and -p specifies the point. The -s tag is used to set the resolution of the final image. The value of 100 means that the final image has 100px per angstrom. Finally, the -i and -o tags designate the input (a CHGCAR of PARCHG file) and output files (an .png image file).

Maybe you are wondering how the values are generated that do not lie exactly at the grid points specified in the CHGCAR file. To plot such values, a trilinear interpolation scheme is used. A detailed description of that algorithm is given here.

## 3.2 Plotting the orbitals

Plotting the four binding orbitals has now become a simple matter of running the program for each of the PARCHG files corresponding to the molecular orbitals. To exemplify:


./bin/edp -i path/to/PARCHG_3s -v 1,0,0 -w 0,0,1 -p 5,5,5 -s 100 -o 3s.png
./bin/edp -i path/to/PARCHG_4s -v 1,0,0 -w 0,0,1 -p 5,5,5 -s 100 -o 4s.png
./bin/edp -i path/to/PARCHG_1p -v 1,0,0 -w 0,0,1 -p 5,5,5 -s 100 -o 1p.png
./bin/edp -i path/to/PARCHG_5s -v 1,0,0 -w 0,0,1 -p 5,5,5 -s 100 -o 5s.png



Finally, you can collect all images and create a nice compilation in your favourite photo-editing program (such as what I did in Figure 1).

## 3.3 Concluding remarks

Obviously, you can use the above procedure to visualise any kind of molecule. Especially simple molecules such as N2 or H2 are easy to visualise. Less obvious is that you can also use the above procedure to visualise the molecular orbitals of a substrate adsorbed on a catalytic surface or the electron density difference between CO in its molecular state and C and O in their atomic states. It all depends on the way the CHGCAR and PARCHG files are created. I hope this tutorial was informative. If you have used this tutorial and found it useful for your work, please drop a line or send me an e-mail. Especially if you have made some nice pictures. :-)

Bart Zijlstra has made a very nice analysis of the results of the electron density visualisation routine as a function of the number of k-points and the smearing value. The result can be seen below. Increasing the number of K-points shows a more detailed result, yet there is no significant change in the characteristic features. #### Drop a line

Question:
What is the answer to One + Seven? zdenko_dot_kolaric_at_siol_dot_net
2016-04-27 19:51:00
Please for something similar for SO2 and SO3.
Best Regards
Zdenko

Question:
What is the answer to One + Three? jkshenton_at_gmail_dot_com
2016-07-14 11:01:39
Hi! Nice work with edp! I've generated lots of nice pictures using it, but I have one question:

How are the saturation levels chosen? I mean, I'd like to plot a colour bar next to it in some units (say eV/bohr^3). How would I go about doing this?

Thanks
Kane

Question:
What is the answer to Ten + Seven? ivo_at_ivofilot_dot_nl
2017-12-31 21:29:24
This feels a bit like necroposting. Nevertheless, I have implemented a legend into EDP. You can enable the legend using the "-l" directive and it will show which colors represent which values. It also mentions the units (which are in eV/A^3).

Question:
What is the answer to Four + Two?

52722
9
28-07-2014
11689
0
30-06-2016
10391
1
22-03-2015
5896
0
29-12-2017
5670
0
27-12-2017