新知一下
海量新知
5 9 9 1 3 7 0

Visualising the electron density using VASP

学术之友 | 致力于分享和传播理论计算教程 2020/05/18 22:28

注:本文转载于 Ivo Filot 个人博客,更多详情点击文末" 阅读原文 "!

理论计算经验贴收集(长期更新)


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:


新知达人, Visualising the electron density using VASP

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.


新知达人, Visualising the electron density using VASP

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   1Direct  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 = normalENCUT = 400
ISMEAR = 0SIGMA = 0.05
EDIFF = 1E-5EDIFFG = -1E-4NSW = 100
ISIF = 2IBRION = 2ISPIN = 1POTIM = 0.5
IALGO = 38NELMIN = 6NELM = 40
EMIN = -30EMAX = 15NEDOS = 4500
NGX = 150NGY = 150NGZ = 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:

新知达人, Visualising the electron density using VASP

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


新知达人, Visualising the electron density using VASP


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 = mediumENCUT = 400
ISMEAR = 0SIGMA = 0.05
EDIFF = 1E-5EDIFFG = -1E-4NSW = 100
ISIF = 2IBRION = 2ISPIN = 1POTIM = 0.5
IALGO = 38NELMIN = 6NELM = 40
EMIN = -30EMAX = 15NEDOS = 4500
NGX = 150NGY = 150NGZ = 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 (https://github.com/ifilot/edp) . 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. :-)


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.


新知达人, Visualising the electron density using VASP

新知精选

更多新知精选