Skip to main content

Projected Density of States

Here we continue with our Aluminum example. Often it is needed to know the contribution from each individual atoms and/or each of their orbital contributions. We can achieve that using projwfc.x code. First, we must perform the self consistent field calculation followed by the non-self consistent field calculation with denser k-points.

pw.x < > al_scf.out
pw.x < > al_nscf.out

Then we prepare the input file for projwfc.x:

prefix= 'al',
outdir= '/tmp/',
filpdos= 'al_pdos.dat'

Perform the calculation:

projwfc.x < > al_projwfc.out

Output data format: the DOS values are written in the file {filpdos}.pdos_atm#N(X)_wfc#M(l), where N is atom number, X is atom symbol, M is wfc number, and l=s,p,d,f one file for each atomic wavefunction read from pseudopotential file. The header of file looks like (for spin polarized calculations, we have separate up and down columns):

E  LDOS(E) PDOS_1(E) ... PDOS_{2l+1}(E)

LDOS=m=12l+1PDOSm(E)LDOS = \sum\limits_{m=1}^{2l+1} PDOS_m (E)

PDOSm(E)PDOS_m (E) \rightarrow projected DOS on atomic wfc with component mm.

Orbital order:

  • for l=1l=1:

    • pz (m=0)p_z~(m=0)
    • pxp_x (real combination of m=±1m=\pm 1 with cosine)
    • pyp_y (real combination of m=±1m=\pm 1 with sine)
  • for l=2l=2:

    • dz2 (m=0)d_{z^2}~(m=0)
    • dzxd_{zx} (real combination of m=±1m=\pm 1 with cosine)
    • dzyd_{zy} (real combination of m=±1m=\pm 1 with sine)
    • dx2y2d_{x^2-y^2} (real combination of m=±2m=\pm 2 with cosine)
    • dxyd_{xy} (real combination of m=±2m=\pm 2 with sine)

Let's make our plots:

import matplotlib.pyplot as plt
from matplotlib import rcParamsDefault
import numpy as np
%matplotlib inline

# load data
def data_loader(fname):
import numpy as np

data = np.loadtxt(fname)
energy = data[:, 0]
pdos = data[:, 2]

return energy, pdos

energy, pdos_s = data_loader('../src/al/al_pdos.dat.pdos_atm#1(Al)_wfc#1(s)')
_, pdos_p = data_loader('../src/al/al_pdos.dat.pdos_atm#1(Al)_wfc#2(p)')
_, pdos_tot = data_loader('../src/al/al_pdos.dat.pdos_tot')

# make plots
plt.figure(figsize = (8, 4))
plt.plot(energy, pdos_s, linewidth=0.75, color='#006699', label='s-orbital')
plt.plot(energy, pdos_p, linewidth=0.75, color='r', label='p-orbital')
plt.plot(energy, pdos_tot, linewidth=0.75, color='k', label='total')
plt.xlabel('Energy (eV)')
plt.axvline(x= 7.9421, linewidth=0.5, color='k', linestyle=(0, (8, 10)))
plt.xlim(-5, 27)
plt.ylim(0, )
plt.fill_between(energy, 0, pdos_s, where=(energy < 7.9421), facecolor='#006699', alpha=0.25)
plt.fill_between(energy, 0, pdos_p, where=(energy < 7.9421), facecolor='r', alpha=0.25)
plt.fill_between(energy, 0, pdos_tot, where=(energy < 7.9421), facecolor='k', alpha=0.25)
# plt.text(6.5, 0.52, 'Fermi energy', fontsize= small, rotation=90)

Here is how our projected density of states plot looks like:


We can perform sums of specific atom or orbital contributions using sumpdos.x code if there are multiple ss or pp orbitals:

sumpdos.x *\(Al\)* > atom_Al_tot.dat
sumpdos.x *\(Al\)*\(s\) > atom_Al_s.dat
sumpdos.x *\(Al\)*\(p\) > atom_Al_p.dat