Bands diagram using VASP and pymatgen

This article presents a python source code in order to plot the bands diagram of graphene calculated using VASP. The plot is done using pymatgen library and RGB coloring adapted from this example. Here is the output obtained with the script :

Graphene bands diagram
Graphene bands diagram

On the band diagram, the contribution of s, (px,py) and pz atomic orbital is map on a RGB scale. Red is associated to a s contribution, green to (px, py) contribution and blue to pz contribution. The same function or procedure can be used to map atomic contributions on the band diagram.

At the end of the article you can find a tarball containing all VASP input/ouput files and the python source code. This calculation is not accurate and is used only for band diagram plotting illustration. For example, you can see a short gap on the DOS which does not exist.

Edit 2016, March 15, update : source code is now python3 and pymatgen version 3.3.5 compatible. You can find the below scripts in an up to date version on github : bandstructureplots

#!/usr/bin/env python
# -*- coding=utf-8 -*-
 
import sys
 
import numpy as np
from numpy import array as npa
 
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
from matplotlib.gridspec import GridSpec
 
import pymatgen as mg
from pymatgen.io.vasp.outputs import Vasprun, Procar
from pymatgen.symmetry.bandstructure import HighSymmKpath
from pymatgen.electronic_structure.core import Spin, Orbital
 
 
def rgbline(ax, k, e, red, green, blue, alpha=1.):
    # creation of segments based on
    # http://nbviewer.ipython.org/urls/raw.github.com/dpsanders/matplotlib-examples/master/colorline.ipynb
    pts = np.array([k, e]).T.reshape(-1, 1, 2)
    seg = np.concatenate([pts[:-1], pts[1:]], axis=1)
 
    nseg = len(k) - 1
    r = [0.5 * (red[i] + red[i + 1]) for i in range(nseg)]
    g = [0.5 * (green[i] + green[i + 1]) for i in range(nseg)]
    b = [0.5 * (blue[i] + blue[i + 1]) for i in range(nseg)]
    a = np.ones(nseg, np.float) * alpha
    lc = LineCollection(seg, colors=list(zip(r, g, b, a)), linewidth=2)
    ax.add_collection(lc)
 
if __name__ == "__main__":
    # read data
    # ---------
 
    # kpoints labels
    path = HighSymmKpath(mg.Structure.from_file("./opt/CONTCAR")).kpath["path"]
    labels = [r"$%s$" % lab for lab in path[0][0:4]]
 
    # bands
    bands = Vasprun("./bands/vasprun.xml").get_band_structure("./bands/KPOINTS", line_mode=True)
 
    # projected bands
    data = Procar("./bands/PROCAR").data
 
    # density of state
    dosrun = Vasprun("./dos/vasprun.xml")
 
    # set up matplotlib plot
    # ----------------------
 
    # general options for plot
    font = {'family': 'serif', 'size': 24}
    plt.rc('font', **font)
 
    # set up 2 graph with aspec ration 2/1
    # plot 1: bands diagram
    # plot 2: Density of State
    gs = GridSpec(1, 2, width_ratios=[2, 1])
    fig = plt.figure(figsize=(11.69, 8.27))
    fig.suptitle("Bands diagram of graphene")
    ax1 = plt.subplot(gs[0])
    ax2 = plt.subplot(gs[1])  # , sharey=ax1)
 
    # set ylim for the plot
    # ---------------------
    emin = 1e100
    emax = -1e100
    for spin in bands.bands.keys():
        for b in range(bands.nb_bands):
            emin = min(emin, min(bands.bands[spin][b]))
            emax = max(emax, max(bands.bands[spin][b]))
 
    emin -= bands.efermi + 1
    emax -= bands.efermi - 1
    ax1.set_ylim(emin, emax)
    ax2.set_ylim(emin, emax)
 
    # Band Diagram
    # ------------
 
    # sum up contribution over carbon atoms
    data = data[Spin.up].sum(axis=2)
 
    # sum up px and py contributions and normalize contributions
    contrib = np.zeros((bands.nb_bands, len(bands.kpoints), 3))
    for b in range(bands.nb_bands):
        for k in range(len(bands.kpoints)):
            sc = data[k][b][Orbital.s.value]**2
            pxpyc = data[k][b][Orbital.px.value]**2 + \
                data[k][b][Orbital.py.value]**2
            pzc = data[k][b][Orbital.pz.value]**2
            tot = sc + pxpyc + pzc
            if tot != 0.0:
                contrib[b, k, 0] = sc / tot
                contrib[b, k, 1] = pxpyc / tot
                contrib[b, k, 2] = pzc / tot
 
    # plot bands using rgb mapping
    for b in range(bands.nb_bands):
        rgbline(ax1,
                range(len(bands.kpoints)),
                [e - bands.efermi for e in bands.bands[Spin.up][b]],
                contrib[b, :, 0],
                contrib[b, :, 1],
                contrib[b, :, 2])
 
    # style
    ax1.set_xlabel("k-points")
    ax1.set_ylabel(r"$E - E_f$   /   eV")
    ax1.grid()
 
    # fermi level at 0
    ax1.hlines(y=0, xmin=0, xmax=len(bands.kpoints), color="k", lw=2)
 
    # labels
    nlabs = len(labels)
    step = len(bands.kpoints) / (nlabs - 1)
    for i, lab in enumerate(labels):
        ax1.vlines(i * step, emin, emax, "k")
    ax1.set_xticks([i * step for i in range(nlabs)])
    ax1.set_xticklabels(labels)
 
    ax1.set_xlim(0, len(bands.kpoints))
 
    # Density of state
    # ----------------
 
    ax2.set_yticklabels([])
    ax2.grid()
    ax2.set_xticks(np.arange(0, 1.5, 0.4))
    ax2.set_xticklabels(np.arange(0, 1.5, 0.4))
    ax2.set_xlim(1e-6, 1.5)
    ax2.hlines(y=0, xmin=0, xmax=1.5, color="k", lw=2)
    ax2.set_xlabel("Density of State")
 
    # s contribution
    ax2.plot(npa(dosrun.pdos[0][Orbital.s][Spin.up]) +
             npa(dosrun.pdos[1][Orbital.s][Spin.up]),
             dosrun.tdos.energies - dosrun.efermi,
             "r-", label="s", linewidth=2)
 
    # px py contribution
    ax2.plot(npa(dosrun.pdos[0][Orbital.px][Spin.up]) +
             npa(dosrun.pdos[1][Orbital.px][Spin.up]) +
             npa(dosrun.pdos[0][Orbital.py][Spin.up]) +
             npa(dosrun.pdos[1][Orbital.py][Spin.up]),
             dosrun.tdos.energies - dosrun.efermi,
             "g-",
             label="(px, py)",
             linewidth=2)
 
    # pz contribution
    ax2.plot(npa(dosrun.pdos[0][Orbital.pz][Spin.up]) +
             npa(dosrun.pdos[1][Orbital.pz][Spin.up]),
             dosrun.tdos.energies - dosrun.efermi,
             "b-", label="pz", linewidth=2)
 
    # total dos
    ax2.fill_between(dosrun.tdos.densities[Spin.up],
                     0,
                     dosrun.tdos.energies - dosrun.efermi,
                     color=(0.7, 0.7, 0.7),
                     facecolor=(0.7, 0.7, 0.7))
 
    ax2.plot(dosrun.tdos.densities[Spin.up],
             dosrun.tdos.energies - dosrun.efermi,
             color=(0.6, 0.6, 0.6),
             label="total DOS")
 
    # plot format style
    # -----------------
 
    ax2.legend(fancybox=True, shadow=True, prop={'size': 18})
    plt.subplots_adjust(wspace=0)
 
    # plt.show()
    plt.savefig(sys.argv[0].strip(".py") + ".pdf", format="pdf")

graphene.zip

Hereafter there are two more examples for Cu and Si. The script may be simpler as all data can be extracted using pymatgen methods.

Si_bands.zip

Cu_bands.zip

Bands diagram of copper

Bands diagram of silicon

39 réflexions au sujet de « Bands diagram using VASP and pymatgen »

  1. Very nice script, but I am getting an error when it tries to load the data.

    Looks something like this:

    : python plot_dos.py 
    Traceback (most recent call last):
      File "plot_dos.py", line 37, in 
        bands = Vasprun("./bands/vasprun.xml").get_band_structure("./bands/KPOINTS", line_mode = True)
      File "/Library/Python/2.7/site-packages/pymatgen/io/vaspio/vasp_output.py", line 443, in get_band_structure
        raise Exception("a band structure along symmetry lines requires a label for each kpoint. "
    Exception: a band structure along symmetry lines requires a label for each kpoint. Check your KPOINTS file
    

    Any ideas how to fix this? I just installed the most current pymatgen software today and am just working through some examples.

    Thanks.

    1. Hi
      I am using the script for MoS2.
      but I have an error that you can see below.
      please help me with it.
      Thx

      Traceback (most recent call last):
      File « ./plot.py », line 85, in
      sc = data[k][b][Orbitals.value]**2
      NameError: name ‘Orbitals’ is not defined

  2. There was a mistake in the KPOINTS file. It must contain a label at the end of each line. Thus it should be :

    calcul des bandes graphene G - M - K - G
    30
    Line-mode
    rec
    0.  0. 0.   ! G
    0.5 0. 0.   ! M
    
    0.5       0.        0. ! M
    0.3333333 0.3333333 0. ! K
      
    0.3333333 0.3333333 0. ! K
    0.        0.        0. ! G
    

    I corrected the graphene.zip tarball and added two more examples.

  3. How to proceed in the case of elements extending to orbital « f »? I see spd_dos localized orbitals « s », « p » and « d » but not how to remain for additional orbital « f » (for example the file Cu).

  4. The script is based on pymatgen methods (http://pymatgen.org/). Indeed, if you look at pymatgen.electronic_structure.dos there is nothing about f orbitals. As I never deal with such heavy elements I do not know how it works for them. You have to try. Load vasprun.xml with pymatgen and look at outputs you get from dos and pdos objects.

    If you succeed in reading the pdos, concerning the colors on the band diagram, You have to choose what you want to see. You cannot have more than 3 quantities which will be mapped on the 3 RGB colors. Thus, again I am not familiar with f elements, you may, for example, map red on s orbitals, green on f orbitals and blue on p and d orbitals… Again you have to try.

  5. Very nice script.

    when I execute this script, I got error message …

    Traceback (most recent call last):
    File « /usr/local/lib/python3.4/site-packages/matplotlib-1.4.3-py3.4-linux-x86_64.egg/matplotlib/colors.py », line 388, in to_rgba_array
    nc = len(c)
    TypeError: object of type ‘zip’ has no len()

    During handling of the above exception, another exception occurred:

    Traceback (most recent call last):
    File « bands_Si.py », line 103, in
    contrib[b,:,2])
    File « bands_Si.py », line 28, in rgbline
    lc = LineCollection(seg, colors=zip(r, g, b, a), linewidth = 2)
    File « /usr/local/lib/python3.4/site-packages/matplotlib-1.4.3-py3.4-linux-x86_64.egg/matplotlib/collections.py », line 1068, in __init__
    colors = mcolors.colorConverter.to_rgba_array(colors)
    File « /usr/local/lib/python3.4/site-packages/matplotlib-1.4.3-py3.4-linux-x86_64.egg/matplotlib/colors.py », line 391, in to_rgba_array
    « Cannot convert argument type %s to rgba array » % type(c))
    ValueError: Cannot convert argument type to rgba array

  6. Hi,
    Thanks for your script, which works quite fine for me. Just one comment, the label positions of band structure should be position-dependent instead of the same distance for all the directions right now.

    Leo

    1. Dear Leo
      Yes it should. I guess there is an attribute or a method of the Kpoints or Band class in pymatgen which should allow to do that but I did not look at this for the moment.
      Germain

  7. This is good code!
    Thank you!

    I have a small question.
    You wrote « e – bands.efermi » in line 114.
    I think dos.efermi is preferable to bands.efermi because bands.efermi depend on the strings of k-points.
    Why did you use bands.efermi?

    1. Kaneko
      Yes it is true that dos.efermi might be more accurate due to a more homogeneous k-points grid. Actually, I used to give dos.efermi as efermi argument of the get_band_structure() method (http://pymatgen.org/pymatgen.io.vasp.html#pymatgen.io.vasp.outputs.Vasprun.get_band_structure). That fix accuracy issue on the fermi level because you manually set it.
      Here, as it is only a graphical example, I do not take care to accuracy. This is probably why I did that in this way.

  8. Hello, Thanks for your script, I am trying to plot s,p,d dos of each atom in a ternary system..
    I tried to modify the script but I couldn’t fine where to chose atom type/number..
    could you tell me please what to add to plot for example « s » dos of Ga and « s » dos of As..

    thank you.

    1. Hello,
      The key point is the contrib tab. You can build a color scale on three quantities of your choice. For example, s for Ga, s for As and p for As. Thus you must compute normalized contributions for each band and each kpoint.
      For another example, you can look at here. Look at input 8 (In [8]) the construction of contrib tab. Here I chose to plot s and p for O and d for Mn.

  9. Hello.
    Thank you for sharing the code.
    I had run the test for the graphene, but I got the following error:
    Python-2.7.8/lib/python2.7/site-packages/matplotlib-1.4.1-py2.7-linux-x86_64.egg/matplotlib/backends/backend_pdf.py », line 418, in __init__
    fh = open(filename, ‘wb’)
    IOError: [Errno 13] Permission denied: ‘/graphene.pdf’
    I guess the code tried to write the output in the root directory. However, I do not familiarize with python and pymatgen so I don’t know how to correct it. Could you please help me with the code?
    I look forward to hearing from you soon. Thank you.

    1. Hello,
      Look at the last line : plt.savefig()
      You have to give as first argument the path to the output file. For example
      plt.savefig("/home/usr/somedir/graphene.pdf")

  10. Is it possible to generate PROCAR into txt file then we can plot it in origin? For example, PROCAR could be shifted to x (symmetry points) and y (energy) columns.

    1. It should be, but there is not a pymatgen method or functions to convert Bands object to txt data. You have to write the file by yourself from Procar object or a band structure object.

  11. The source code is now python3 and pymatgen version 3.3.5 compatible. There were numerous changes in the Procar class in order to support numpy ndarray.

  12. Thanks for updating the script for python3.
    Can you please explain how to select particular ion in the newer version. « dosrun.complete_dos.get_element_spd_dos » does not work for this version.

    1. What do you mean ? It does work but the method return a dict such as :

      Cu = mg.Element("Cu")
      run.complete_dos.get_element_spd_dos(Cu)
      
      {orbitaltype.d: pymatgen.electronic_structure.dos.Dos,
       orbitaltype.p: pymatgen.electronic_structure.dos.Dos,
       orbitaltype.s: pymatgen.electronic_structure.dos.Dos}
      

      The keys of the dict are OrbitalType object and not simply a string of the orbital type such as « S » or « P ».
      Thus, you have to look at pymatgen.electronic_structuer.core.OrbitalType

  13. Thank you for the scripts. It is an excellent piece of code. I am hitting a keyerror while plottinf spd contributions. it says the keys S, P, D does not exist! I also noticed that there are two different functions called in for this job, one is npa and another is spd_dos. In my case the first one isn’t found and the second one is giving me this key error. Please let me know where this is coming from.

    File « Si_bands.py », line 141, in
    ax2.plot(spd_dos[« S »].densities[Spin.up],
    KeyError: ‘S’

    Thank you in advance

    1. It may depends on the pymatgen version you are using. Considering the last version, « S », « P » and « D » keys were deprecated and OrbitalType object are used instead.

      Take a look at pymatgen.electronic_structuer.core.OrbitalType and try the following :

      Cu = mg.Element("Cu")
      spd_dos = run.complete_dos.get_element_spd_dos(Cu)
      print(spd_dos.keys())
      

      I update the script in the archive Si_bands.zip according to the current version of pymatgen.

  14. Hi Germain,
    Thank you for your code!
    I ran the latest version on github and I ran into an error:
    KeyError Traceback (most recent call last)
    /Users/hbliu/plot_bandos_test.py in ()
    83 for b in range(bands.nb_bands):
    84 for k in range(len(bands.kpoints)):
    —> 85 sc = pbands[Spin.up][b][k][name][« s »]**2
    86 pc = pbands[Spin.up][b][k][name][« p »]**2
    87 dc = pbands[Spin.up][b][k][name][« d »]**2

    KeyError:
    Do you know what happens? Thank you so much!

  15. Hi Germain,

    I have run the script successfully. Thank you for this nice script. I have few questions..
    1) In case of spin-polarization, How to plot bands for both spin(spin up and spin down)?

    2) I am interested to look DOS and band structure near Fermi level. How can I specify the energy range in the plot?

    3) Is is possible to plot the band structure with specific orbital(like the DOS of d_xy, d_yz……or t_2g and e_g)?

    Best regards,

    Ashis Kundu

  16. Hello!

    My question is about inconsistency of band and DOS plots of graphene across the Fermi level. Why is there band gap in DOS and how can we overcome this issue or make consistent picture of DOS and band structure? I use vasp.

    1. Hello,
      I agree, it should be consistent.
      I guess it is an accuracy issue because, as I said, I simply did a quick calculation with low accuracy.
      If you perform an accurate calculation I guess it will be ok.
      Germain

  17. Dear all
    Thank you for all you comments.
    I closed the comments on this article in order to encourage you to open issues on the github repository concerning the above scripts.
    To ask questions, report bugs, request features or contribute to these scripts, please open an issue on this page : https://github.com/gVallverdu/bandstructureplots/issues
    The github issues system is better for such discussions.
    see you on github
    Germain

Les commentaires sont fermés.