Categories
notebook

VASP: Band-structure calculation Tutorial

[toc]

In this tutorial, I will describe how to calculate the Band Structure using VASP, taking Si as an example. In the calculation, I chose a FCC primitive unit cell with 2 atoms in the unit cell. But first lets take some considerations regarding Accurate Band Structure Calculations.

Accurate Band Structure Calculations

A high quality DOS and band-structure calculation might require very fine k-meshes up to 24x24x24 grid points for small unit cells. Particularly in band-structure calculations, where one needs to calculate the eigenvalues along certain high symmetry lines in the BZ (at least 10 k-points are required for reasonably results).

Since, the charge density and the effective potential converge rapidly with increasing number of k-points, it is often helpful to calculate the selfconsistent charge density using a few k-points. In the second step, a non-selfconsistent calculation using the precalculated CHGCAR file from the selfconsistent run (i.e. ICHARG=11) can be performed.

For ICHARG=11 and density functional theory calculations, all k-points become essentially independent, because the charge density and the potential are kept fixed. If necessary it is even possible, to split up the k-points and calculating the eigenvalues individually for each k-point, although with present computing platforms this is rarely required.

Adapted from VASP manual

1st Step

Files needed:

INCAR : Contains all the control parameters.
KPOINTS : Contains the k-point sampling scheme.
POSCAR : Contains information on unit cell and the basis atoms.
POTCAR : Pseudo-Potential File used by VASP.

To calculate the Band structure, we need to first run self-consistently to get the charge density.
INCAR:

Si band-structure
SYSTEM = Si

! Start parameter for this run:
ISTART = 0
       ! Begin from scratch. Does not read existing (if any) wavefunction file (WAVECAR)

NWRITE = 2
      ! This flag determines how much will be written to the file OUTCAR ('verbosity flag').

ICHARG = 2
       ! Initial guess for starting charge-density profile of the whole system.
       ! 2 : Superposition of atomic-charge densities.

ISPIN = 1
      ! Spin polarization keyword. 2 : spin-polarized calculation. 1: default (average).

NWRITE = 2
      ! This flag determines how much will be written to the file OUTCAR ('verbosity flag').

PREC = Accurate
     ! Alternative keyword to define ENCUT and NGX,NGY,NGZ if these are
     ! not explicitly defined by the user in the INCAR file.

! Electronic Relaxation 1
NELM = 90
     ! Maximum number of electronic SC (selfconsistency) steps which may be performed.
NELMIN= 8
     ! Minimum number of electronic SC steps.
NELMDL= 10
     ! Number of non-selfconsistent steps at the beginning.
EDIFF = 0.1E-03
      ! Threshold criterion for termination of SC cycles Units: eV
LREAL = .FALSE.
      ! Real-space projection of SC steps.

! Ionic relaxation
IBRION = 2
       ! Defines ionic motion during structural relaxation, molecular dynamics,
       ! frozen-phonon and vibrational-frequency setups, etc.
       ! Note: 2 : Conjugate-gradient algorithm for structural relaxation.
       ! I consider it a robust algorithm for structural relaxation in
       ! semiconducting systems. May be a bit inefficient close to the minimum,
       ! namely will use a lot of steps to the minimum.

NSW = 0
    ! Sets maximum number of ionic steps. If convergence is reached earlier code
    ! terminates (see below).

EDIFFG = 0.1E-2
       ! Threshold criterion for termination of structural optimization.
       ! Negative sign means that a force tolerance is tested in successive
       ! sc cycles. otherwise, an energy criterion is tested. Units: eV/A.

ISIF = 2
     ! Controls the type of structural optimization needed.
     ! 3 optimizes the cell volume and shape in addition to performing relaxation
     ! for all atoms of the system.
     ! Note: 3 is useful to perform for initial optimization of the primitive bulk cell.
     ! For final supercell defect calculations the keyword 2 is preferred, namely
     ! fix the supercell parameters and relax ionic positions.

POTIM = 0.10
      ! Scaling constant for structural optimization. For relaxation it
      ! scales the ionic displacement for each force step.
      ! Note: From my experience, start with a small value for a completely
      ! unrelaxed configuration, e.g. set POTIM to 0.0002 for instance and then
      ! progressively increase it to speed up convergence. A value of 0.01 is
      ! reasonably safe for medium forces, e.g. ~ 1 eV/Ang..

! DOS related values:
ISMEAR = 0
       ! Defines occupancies (smearing technique) for electronic states.
       ! This keyword is very system specific, e.g. if you are dealing with a
       ! semiconductor/insulator or a metal or a complex system, namely defects
       ! with occupied levels in the gap.
       ! Depends also on the type of calculations one wants to perform, e.g.
       ! DOS, structural relaxation, etc. For instance, for DOS, the tetrahedron
       ! method (ISMEAR = -5) is recommended.

SIGMA = 0.10
      ! Width of smearing (units: eV). Again this is system dependent. Also
      ! will depend upon k-point mesh chosen.

! Electronic relaxation 2 (details)
LWAVE = T
      ! Write WAVECAR file
LCHARG = T
      ! Write CHGCAR file

KPOINTS: (A 4x4x4 mesh is usually good enough!)

Automatic MP mesh
0              ! number of k-points = 0 ->automatic generation scheme
Monkhorst-Pack ! select Monkhorst-Pack (first letter is significant)
4  4  4        ! size of mesh (4x4x4 points along b1, b2, b3)
0. 0. 0.       ! shift of the k-mesh

POSCAR:

Si fcc primite cell
  5.38936		! lattice parameter
 0.5 0.5 0.0
 0.0 0.5 0.5
 0.5 0.0 0.5
  Si
  2
Cartesian
0.0000000000000 0.00000000000 0.0000000000000
0.2500000000000 0.25000000000 0.2500000000000

Once we have all the files, we issue a command to run vasp, such as ‘./vasp>&output&’ (This depends on your system, if a queue system is running, then you need to write a simple script to queue the job).
After this step is completed, we should now have the charge density, which is contained in CHGCAR file. In order to get band structure, we need to do non-self consistent run on each desired k-point, by connecting these information, we can get the E~K dispersion relation, which is the Band structure.

2nd Step

Now we need to get the energy for each k-points based on the charge density we got from the first job.

We need to modify the following parameter in INCAR:

ICHARG=11

ICHARG=11 means reading charge density from CHGCAR and kept constant during the subsequent run.

We also need to modify the KPOINTS file, to specify along which direction (Usually some high symmetry direction) we need VASP to calculate the energy. In VASP 4.6 and higher, there is an easy way of doing this.

k-points along high symmetry lines
 20  ! 20 intersections 
Line-mode
rec
  0.5 0.5 0.5 ! L
  0   0   0   ! gamma
  
  0   0   0   ! gamma
  0.5 0.5 0   ! X
fcc_bz

By specifying the Line-mode, VASP automatically interpolate between the points you want to calculate. The above KPOINTS file instructs VASP to calculate the energy at each k-point between L point and Gamma point, Gamma point and X point (see figure below). Along each line, 20 k-points are calculated.

Note that the coordinates of the k-points can be supplied in cartesian (4th line starts with c or k) or in reciprocal coordinates (4th line starts with r).

Having modified these two files, we rerun the VASP again and we will get all the information we need to plot the band structure (contained in EIGENVAL file).

Output

In EIGENVAL file, there would be lines like this:

  0.2631579E-01  0.2631579E-01  0.0000000E+00  0.2500000E-01
    1       -6.043633
    2        5.995153
    3        6.020739
    4        6.020740
    5        8.618293
    6        8.718029
    7        8.718030
    8        9.779300

The first line is just the k-point, the following 8 lines tell us the energy of 8 bands at that particular energy point. You can visualize these information using any plotting tool. Since you have the Enegy of the bands at each k-point, this is exactly what you want – the Band Structure.

Gnuplot

Now, we need to run a bash script to split the EIGENVAL file so that the eigenvalue of each band 1,2,3… is written in separate files. You could use gnuplot to visualize them.

I write a script like this. Let’s call it ‘band.sh’.

#!/bin/bash
# a shell script band to split EIGENVAL file
mkdir tmp
for c in {1..8}
do
        grep '  '$c'      ' EIGENVAL>|tmp/band$c.dat
done

cd tmp
for d in {1..8}
do
        awk 'BEGIN{FS=OFS="      "}{$1="";sub("      ","")}1' band$d.dat > out$d.dat
done

sed -i -e '$a\ ' *
sed -i -e '$a\ ' *

rm band*
rename out out0 out?.dat
rename out out0 out??.dat
rename out out0 out???.dat

cat *.dat > ../bandst.dat
cd ..
rm -r tmp

In my case I have 8 bands. We must change the script according to the number of calculated bands.
You need to be careful about the space in grep command. Double check if my space is the same with your EIGENVAL file. If not, just copy from your file and paste it in ‘ ‘

Now, how to run a shell script in Linux. Just type:

$ chmod 755 band.sh
$ ./band.sh

You will see the file ‘bandst.dat’ in your folder.

Now, how to plot these data using gnuplot?
Hit:

$ gnuplot

You enter the environment of Gnuplot and you could type some commands like these:

gnuplot > reset
gnuplot > plot 'bandst.dat' index 1:8 with lines

Change index parameter according to the bands you want to plot.

bandst_si

The resulting plot is something like this:

You can now play with the style and appearance.