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:

- Quantum Espresso 5.3+
- xcrysden
- python3.4+
- Matplotlib
- 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

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?

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?

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

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

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.

Hi ANINDYA BOSE

did you run successfully bands.py. I am also facing same problem as ZACK facing ?

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

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

Hi Levi

could you do some tutorial on band unfolding, or give me some suggestion on how I do it.

Thank you!

Nada

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.

Has the issue resolved now?

Dear Levi Lenz:

What file do you refer to as ”Symmetry file”?