Band Diagram Tutorial for Quantum Espresso

Hello,

In this tutorial we will go over how to perform a band structure calculation in Quantum Espresso and then how to plot it in python. This tutorial assumes you use the following things:

1. Quantum Espresso 5.3+
2. xcrysden
3. python3.4+
4. Matplotlib
5. Numpy

The python requirements will only be if you want to use the plotting tools I have generated to plot my own personal data. This tutorial should also work for a range of Quantum Espresso (QE) versions as well, however they changed the way post-processing takes input files around 5.x so you may get some errors if you use them verbatim. If you do experience errors, feel free to email me if you cannot find proper documentation online for your QE version.

All the files from this tutorial can be downloaded here

Band Structure Calculation

The bandstructure calculation requires four major steps: a relax calculation, a dense nscf (or scf if you are adventurous) k-point calculation, a bands calculation, finally followed post-processing with bands.x and plotband.x if you desire. I will go over these, however I will be plotting the bandstructure with a custom python script that you can customize on your own.

In this tutorial I will be using an example from my ongoing research: zirconium phosphate. The unit cell is about the same size of silicon so you should be able to run it on your laptop, but it is a bit more interesting as it has three different elements and a whole lot of electrons. It also is not well studied, so the band structure can’t look too bad.

Relax and dense nscf calculations

The first step is to fully relax your system. As an example, here is the relax.in file that I have used:

 &control
    calculation='relax',
    restart_mode='from_scratch',
    pseudo_dir = '.',
    outdir='.',
    prefix='1',
    wf_collect=.true.
    verbosity='high'
    forc_conv_thr = 0.01D-3
    etot_conv_thr = 1D-6
    wf_collect = .true.
 /
 &system
    ibrav = 0,
    nat = 11,
    ntyp = 4,
    ecutwfc = 50,
    occupations='smearing',
    degauss=0.005
    vdw_corr='dft-d'
    tot_charge=0.5
 /
 &electrons
    conv_thr = 1.0e-9
    electron_maxstep = 150,
    mixing_beta=0.3
    diagonalization='cg'
 /
 &ions
    ion_dynamics = 'bfgs'
 /
 &cell
    cell_factor = 10
    cell_dynamics='damp-w'
/
ATOMIC_SPECIES
P 30.97 P.pbe-n-van.UPF
Zr 91.224 Zr.pbe-nsp-van.UPF
O 16 O.pbe-rrkjus.UPF
H 1 H.pbe-rrkjus.UPF
ATOMIC_POSITIONS {crystal}
Zr       0.000381444   0.000263039  -0.000011516
P        0.333734416   0.666926994   0.300249262
P        0.667031911   0.333607637   0.699728761
O        0.139956220   0.782172447   0.215796553
O        0.218465107   0.357889207   0.215824072
O        0.642777110   0.860716177   0.215852669
O        0.860780811   0.218340246   0.784215740
O        0.782240518   0.642635199   0.784154708
O        0.357982949   0.139830740   0.784122021
H        0.333724901   0.666936549   0.551898251
H        0.667099278   0.333619644   0.448080371
K_POINTS {automatic}
3 3 3 0 0 0
CELL_PARAMETERS {angstrom}
5.417202237 7.58679E-06 -1.5423E-05
-2.708594549  4.691425372 4.98917E-06
-1.58786E-05  -3.23754E-06  5.575997826

Once this calculation finished, I used my existing ScfCreator.py to generate the nscf file. The only note here is that I actually modified the k-point grid to be 12 12 12 0 0 0 in the nscf calculation to be super dense. This will make your bands calculation converge much smoother and should prevent issues down the road.

Bands Calculation

The bread and butter of a band structure calculation is actually in the bands calculation. This is also performed with pw.x, with the input of calculation=’bands’ in lieu of relax or nscf. The trick here is you want to pick your k-path in reciprocal space. Normally the path you choose can tell you about the unique physics of your system. For a general calculation such as this, I simply go over all the high-symmetry points in k-space. This gives a qualitative view of the bandstructure and can guide future work on different paths in k-space. Generally I use a cool little program called xcrysden to select this path. Below are some pictures to guide you on selecting this path:

You can see in the last slide that I went over a range of k-space that selected all the unique high-symmetry points. Remember, right now all we care about is the qualitative look of the band structure. If you were to be looking into this for rigorous scientific applications, you would need to be more careful about how you choose these points. At any rate, once you have the output of xcrysden (here I made 300 points along the kpath that we have selected), you should end up with a bands.in file such as this:

 &control
    calculation='bands',
    restart_mode='from_scratch',
    pseudo_dir = '.',
    outdir='.',
    prefix='1',
    wf_collect=.true.
    verbosity='high'
    forc_conv_thr = 0.01D-3
    etot_conv_thr = 1D-6
    wf_collect = .true.
 /
 &system
    ibrav = 0,
    nat = 11,
    ntyp = 4,
    ecutwfc = 50,
    occupations='smearing',
    degauss=0.005
    vdw_corr='dft-d'
    tot_charge=0.5
 /
 &electrons
    conv_thr = 1.0e-9
    electron_maxstep = 150,
    mixing_beta=0.3
    diagonalization='cg'
 /
 &ions
    ion_dynamics = 'bfgs'
 /
 &cell
    cell_factor = 10
    cell_dynamics='damp-w'
/
ATOMIC_SPECIES
P 30.97 P.pbe-n-van.UPF
Zr 91.224 Zr.pbe-nsp-van.UPF
O 16 O.pbe-rrkjus.UPF
H 1 H.pbe-rrkjus.UPF
ATOMIC_POSITIONS {crystal}
Zr       0.000381444   0.000263039  -0.000011516
P        0.333734416   0.666926994   0.300249262
P        0.667031911   0.333607637   0.699728761
O        0.139956220   0.782172447   0.215796553
O        0.218465107   0.357889207   0.215824072
O        0.642777110   0.860716177   0.215852669
O        0.860780811   0.218340246   0.784215740
O        0.782240518   0.642635199   0.784154708
O        0.357982949   0.139830740   0.784122021
H        0.333724901   0.666936549   0.551898251
H        0.667099278   0.333619644   0.448080371
CELL_PARAMETERS {angstrom}
5.417202237 7.58679E-06 -1.5423E-05
-2.708594549  4.691425372 4.98917E-06
-1.58786E-05  -3.23754E-06  5.575997826
 K_POINTS crystal
         300
        0.0000000000    0.0000000000    0.0000000000    1.0
        0.0142857143    0.0000000000    0.0000000000    1.0
        0.0285714286    0.0000000000    0.0000000000    1.0
        0.0428571429    0.0000000000    0.0000000000    1.0

This has been truncated for length as under the K_POINTS card I would have 300 lines. If you would like to see the whole file, please download the tarball that is provided with this post. Once you have this, running pw.x < bands.in > bands.out will perform a bandstructure calculation along the given k-path.

Running band.x

The final step to get the plottable bands from the output of QE is to run the post processing command bands.x. This simply extracts the high symmetry points and the band energies into plottable files for analysis. The input file is below:

&BANDS
  prefix="1"
  outdir="."
  filband="Bandx.dat"
/

Running this as bands.x < bandx.in > bandx.out will create all the files you need to plot your band structure. I have created a python script to help with the plotting. However, you can just use plotbands.x to readily generate a view-able plot.

plotbands.x

I am not a huge fan of it, but QE does provide another post-processing tool to plot out the information generated by bands.x. You can run plotbands.x without an input file as it will guide you through all the required information for a plot. However, if you would like to use an input file, you can use the following as an guide:

Bandx.dat !output of bands.x
-1.9546, 4.0454 !Energy range of plot in Emin Emax
plotbands.xmgr !output for xmgrace plotting
plotbands.ps ! output for a ps plot
1.0454  ! The Fermi level (marks it on the plot)
1.0, 1.0454 ! The Estep and Reference Energy

Plotting the bands the Python Way

As I mentioned, I really dont like how QE plots the bands. Not that this is plotted wrong, but the plots are a bit ugly. So in this section I will show you how to plot the bands with the matplotlib and numpy.

The brunt of the processing is converting the bandx.dat.gnu data into x-y data that plots properly. To do this, you need to go band-by-band and extract the data into energy-positional data. This is relatively straight forward as the data in bandx.dat.gnu is formatted in a predicable manner. It stores the lowest energy band first then the next lowest then the next and on and on until you have gone over the entire band structure. I have written up a python function that can easily extract this into band-wise data that can then be plotted. You can download the fully annotated file here.

#!/usr/bin/env python
# This was written by Levi Lentz for the Kolpak Group at MIT
# This is distributed under the MIT license
import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import matplotlib.gridspec as gs
import sys

#This function extracts the high symmetry points from the output of bandx.out
def Symmetries(fstring): 
  f = open(fstring,'r')
  x = np.zeros(0)
  for i in f:
    if "high-symmetry" in i:
      x = np.append(x,float(i.split()[-1]))
  f.close()
  return x
# This function takes in the datafile, the fermi energy, the symmetry file, a subplot, and the label
# It then extracts the band data, and plots the bands, the fermi energy in red, and the high symmetry points
def bndplot(datafile,fermi,symmetryfile,subplot,label):
  z = np.loadtxt(datafile) #This loads the bandx.dat.gnu file
  x = np.unique(z[:,0]) #This is all the unique x-points
  bands = []
  bndl = len(z[z[:,0]==x[1]]) #This gives the number of bands in the calculation
  Fermi = float(fermi)
  axis = [min(x),max(x),Fermi - 4, Fermi + 4]
  for i in range(0,bndl):
    bands.append(np.zeros([len(x),2])) #This is where we storre the bands
  for i in range(0,len(x)):
    sel = z[z[:,0] == x[i]]  #Here is the energies for a given x
    test = []
    for j in range(0,bndl): #This separates it out into a single band
      bands[j][i][0] = x[i]
      bands[j][i][1] = np.multiply(sel[j][1],13.605698066)
  for i in bands: #Here we plots the bands
    subplot.plot(i[:,0],i[:,1],color="black")
  temp = Symmetries(symmetryfile)
  for j in temp: #This is the high symmetry lines
    x1 = [j,j]
    x2 = [axis[2],axis[3]]
    subplot.plot(x1,x2,'--',lw=0.55,color='black',alpha=0.75)
  subplot.plot([min(x),max(x)],[Fermi,Fermi],color='red',)
  subplot.set_xticklabels([])
  subplot.set_ylim([axis[2],axis[3]])
  subplot.set_xlim([axis[0],axis[1]])
  subplot.text((axis[1]-axis[0])/2.0,axis[3]+0.2,label,va='center',ha='center',fontsize=20)

You will need to extend this function for your own plotting needs. By extending this function I was able to make the following plot for a study we are doing on zirconium phosphate.

Happy computing,

Levi