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σ.
In order to construct the electron density, we first have to perform a series of VASP calculations.
If you know how to do this, skip to the next section, else continue reading.
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
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):
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.
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).
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
-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
-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.
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).
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.