Band structure calculation example

As a first step in a band structure calculation, the charge density is calculated (here for fcc-Ag):

#!/usr/bin/env python
from ase.structure import bulk
from espresso import espresso

a=bulk('Ag','fcc',4.0853)
a.calc = espresso(pw=350.,dw=3500.,kpts=(12,12,12),xc='PBE',outdir='agchg')
a.get_potential_energy()
a.calc.save_flev_chg('Ag_efchg.tgz')

save_flev_chg stores both the charge density and also the Fermi level. The charge density defines the DFT-Hamiltonian so we can now calculate the electronic dispersion along lines in the Brillouin zone non-selfconsistently (i.e. at fixed charge density). The stored Fermi level allows us to plot the band structure relative to it.
Here is an example for a non-selfconsistent band structure calculation:

#!/usr/bin/env python
from ase.structure import bulk
from ase.dft.kpoints import ibz_points,get_bandpath
from espresso import espresso

a=bulk('Ag','fcc',4.0853)
a.calc = espresso(pw=350.,dw=3500.,xc='PBE',outdir='agdisp')
a.calc.load_flev_chg('Ag_efchg.tgz')

ip = ibz_points['fcc']
points = ['Gamma','X','W','K','L','Gamma']

bzpath = [ip[p] for p in points]

kpts, x, X = get_bandpath(bzpath, a.cell, npoints=300)

energies = a.calc.calc_bandstructure(kpts)

import pickle
f = open('Agdisp.pickle', 'w')
pickle.dump((points, kpts, x, X, energies), f)
f.close()

Agdisp.pickle contains the band structure, which can be plotted using the following script:

import pickle
import matplotlib.pyplot as plt
s, k, x, X, e = pickle.load(open('Agdisp.pickle'))
symbols = [t.replace('Gamma', '$\Gamma$') for t in s]
emin = -10
emax = 4
plt.figure(figsize=(6, 4))
for n in range(len(e[0])):
    plt.plot(x, e[:, n], 'b-')
plt.plot([X[0],X[-1]], [0,0], 'k--')
for p in X:
    plt.plot([p, p], [emin, emax], 'k-')
plt.xticks(X, symbols)
plt.axis(xmin=0, xmax=X[-1], ymin=emin, ymax=emax)
plt.ylabel('$E-E_{\\rm Fermi}$ [eV]')
plt.savefig('Agdisp.pdf')

The output will look like this:

The bands can be colored according to their atomic orbital character. We therefore calculate the atomic projections of the bands in addition to their energies:

#!/usr/bin/env python
from ase.structure import bulk
from ase.dft.kpoints import ibz_points,get_bandpath
from espresso import espresso

a=bulk('Ag','fcc',4.0853)
a.calc = espresso(pw=350.,dw=3500.,xc='PBE',outdir='agdispproj')
a.calc.load_flev_chg('Ag_efchg.tgz')

ip = ibz_points['fcc']
points = ['Gamma','X','W','K','L','Gamma']

bzpath = [ip[p] for p in points]

kpts, x, X = get_bandpath(bzpath, a.cell, npoints=300)

energies, proj = a.calc.calc_bandstructure(kpts, atomic_projections=True)

import pickle
f = open('Agdispproj.pickle', 'w')
pickle.dump((points, kpts, x, X, energies, proj), f)
f.close()

Agdispproj.pickle contains the band structure and its orbital characters, which can be plotted using the following script:

import pickle
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
import numpy as np
s, k, x, X, e, proj = pickle.load(open('Agdispproj.pickle'))
symbols = [t.replace('Gamma', '$\Gamma$') for t in s]
emin = -10
emax = 4

spd_atom_0 = np.zeros(3, np.float)
#get indices of s-, p- and d-projectors for atom 0
ps = []
pp = []
pdt2g = []
pdeg = []
for i,q in enumerate(proj[0]):
    if q[0]==0: #atom 0;  ==1 for atom 1 etc.
        if q[1]==0:
            ps.append(i)
        elif q[1]==1:
            pp.append(i)
        #one could also distinguish wrt. q[2] (m)
        #to resolve px,py,pz or d3z2-r2,dxz,dyz,dx2-y2,dxy, resp.
        elif q[1]==2:
            if q[2]==1 or q[2]==4:
                pdeg.append(i)
            else:
                pdt2g.append(i)

#sum up absolute projections for s, p, and d on atom 0
nbnd = len(e[0])
nkp = len(k)
projs = np.zeros((nkp,nbnd), np.float)
projp = np.zeros((nkp,nbnd), np.float)
projdt2g = np.zeros((nkp,nbnd), np.float)
projdeg = np.zeros((nkp,nbnd), np.float)
for i in range(nkp):
    for j in ps:
        for b in range(nbnd):
            projs[i][b] += np.abs(proj[1][i][j][b])**2
    for j in pp:
        for b in range(nbnd):
            projp[i][b] += np.abs(proj[1][i][j][b])**2
    for j in pdt2g:
        for b in range(nbnd):
            projdt2g[i][b] += np.abs(proj[1][i][j][b])**2
    for j in pdeg:
        for b in range(nbnd):
            projdeg[i][b] += np.abs(proj[1][i][j][b])**2

#normalize projections for rgb coloring
projsp = projs+projp
maxsp = np.max(projsp)
maxdt2g = np.max(projdt2g)
maxdeg = np.max(projdeg)
invmax = 1./max([maxsp,maxdt2g,maxdeg])
projsp *= invmax
projdt2g *= invmax
projdeg *= invmax

#function to color band structure in rgb color scheme
def rgbline(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=zip(r,g,b,a))
    plt.gca().add_collection(lc)

plt.figure(figsize=(6, 4))
for n in range(len(e[0])):
    rgbline(x, e[:, n], projsp[:,n], projdt2g[:,n], projdeg[:,n])
plt.plot([X[0],X[-1]], [0,0], 'k--')
for p in X:
    plt.plot([p, p], [emin, emax], 'k-')
plt.xticks(X, symbols)
plt.axis(xmin=0, xmax=X[-1], ymin=emin, ymax=emax)
plt.ylabel('$E-E_{\\rm Fermi}$ [eV]')
plt.savefig('Agdispproj.pdf')

The output will look like this: