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.

Three different charged states of zirconium phosphate
Three different charged states of zirconium phosphate

Happy computing,

Levi

28 comments Add yours
  1. Hi. I’m new in python so the are questions.
    In what directory place and launch Bands.py? Working directory that contents output files of QE ?
    How/ Where to determine names of datafile, symmetryfile?

  2. Hi, Thanks for this nice tutorial,

    I am trying to compile you python script but actually not plot is showin up? Any idea what could it be?

    1. Hi Zack,

      I am not entirely sure. Do you get an output of *.png or *.pdf in the current working directory? If you are running on a machine that does not have an X-Windows system, it might not show up directly.

      The other, and more obvious, problem would be with matplotlib or numpy. I recommend installing anaconda3 to run all my scripts. See: https://www.continuum.io/downloads

      Cheers,

      Levi

      1. Hi Levi,

        Thanks for taking time to answer my question.

        Actually I am working with ubuntu and when I run the code it just doesn’t create anything folder in the working directory. At the begining I though it was a problem of my file. But when I downlead your files I got the same problem.

        Regards,

        Zakaria

        1. HI Levi or Zack

          did you ever get this solved for Ubuntu. I even tried it on windows and it didn’t work. I am about to try it on anaconda on ubuntu. Please any help would be great

  3. Dear Sir,
    I want to Label M,K,G brillouin zone points in x axis in the dotted line regions, How can I do that.Please tell me the necessary commands.

  4. Hi Levi
    I am having problems running bands.x.

    I have finished running pw.x bands.out. This was done on a cluster where you provide the input file and the .pbs file is constructed.

    Now i want to continue with bands.x.
    This is the way i did it but it did not work.
    I created this file called bands.in

    &bands
    prefix = ‘ZnO_Co_C’
    filband = ‘bands.dat’
    no_overlap=.true.
    /

    and placed in the same folder where I have run scf and bands.
    I then tried to run bands.x by using bands.x in the .pbs file and bands.in (above) but it did not .

    Can you please help?

    Thanks
    George

  5. Dear colleagues
    Thank for your detail tutorial on band structure calculation . I am asking how did you add the vdw_corr=’dft-d in quantum espresso caculation

  6. Hello dear Levi Lentz
    Thank you for you helpfull tutorial. i run successfully 1. pw.x relax.out , 2. pw.x nscf.out , 3. pw.x bands.out , 4. bands.x bandx.out. but i cant plot data!
    can you help me how to call that two function?
    i have problem in subplot section. i call as fallow:

    symmetryfile=Symmetries(‘bandx.out’)
    fig, axs = plt.subplots(1,1,figsize=(9, 4.5), tight_layout=True)
    bndplot(‘Bandx.dat.gnu’,2.3439,symmetryfile’,axs,fig)
    plt.show()
    thank for your helping.
    sumi.

    1. Hello dear David Barcene

      The term “high-symmetry” exists in the file bandx.out (The result of input (for example bandx.in)

      & BANDS
      prefix = “1”
      outdir = “.”
      filband = “Bandx.dat”
      no_overlap = .true.
      /

      generate two OUTPUT files: bandx.dat and bandx.out.
      It should be noted that I performed as follows:
      Path-toPath/bin/bands.x bandx.out

      So the block in Python script:

      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

      it is to collect the results “high-symmetry” in the file: bandx.out

          1. Sorry I could not write in the comment the lower and upper mathematical signs, respectively.

            Path-toPath/bin /bands.x , bandx.in , bandx.out

            – replace the two commas (,) by the lower and upper mathematical signs, respectively.

          2. Don’t worry, thanks for the information, I actually got successful in doing my code and generating the figures I needed.

  7. Can you post a simple code just to plot one band structure calculation and then give an example for plotting different band structure plots as you have shown here?

  8. I see your code converts the data in the gnu file from Ry to eV. I don’t know if this is a change from the older versions of quantum espresso, but the data in the gnu file generated by bands.x appears to already be in eV. I had to take out the conversion to get plots with your program to match the output of plotbands.x, which asks for the range of the plot in eV.

Leave a Reply

Your email address will not be published. Required fields are marked *