Visualising the electron density of the binding orbitals of the CO molecule using VASP

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)

If you know how to do this, skip to the next section, else continue reading.

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

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

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

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. Detailed compilation instructions are provided there, but will be repeated in short here.

3.1 Compiling EDP

To compile EDP, you need to have the development libraries of pcrecpp, tclap and Cairo installed. In general, these libraries are available via the package manager of your distribution. To compile the program, you simply need to type make. The executable is stored in the bin folder. If the program does not compile, you might need to tweak the Makefile.

3.2 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.3 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.4 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. :-)

4.Comments

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.

If you have questions or comments, feel free to drop a line! Like what you read? Share this page with your friends and colleagues.

Comments

Question:
What is the answer to Six + Four?
Please answer with a whole number, i.e. 2, 3, 5, 8,...
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 Two + Eight?
Please answer with a whole number, i.e. 2, 3, 5, 8,...
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 Nine + Three?
Please answer with a whole number, i.e. 2, 3, 5, 8,...