Density of States
In this example we will cover how to use ase-espresso to calculate density of states.
The following script generates and saves the density of states.
#######################################################
## Tutorial by: Hongliang Xin
## Topic: Density of States in ase-espresso
##
## dependencies: ase, espresso, pickle
##
## Generates the data files necessary to analyze
## the density of states for a Pd(111) slab
#######################################################
#!/usr/bin/env python
from ase import io, Atom, Atoms
from ase.constraints import FixAtoms
from ase.optimize import QuasiNewton
from ase.lattice.surface import *
from ase.visualize import view
from ase.lattice.general_surface import surface
from espresso import espresso
import os
# Surface
a0 = 3.99663379586
bulk = Atoms('Pd4',
scaled_positions=[(0, 0, 0),
(0.5, 0.5, 0),
(0.5, 0, 0.5),
(0, 0.5, 0.5)],
cell=[a0, a0, a0],
pbc=True)
atoms = surface(bulk,(1, 1, 1),4,0.0)
atoms.center(vacuum=10.0, axis=2)
io.write('InitialGeom.traj',atoms)
print atoms
print 'Unit cell array \n', atoms.get_cell()
atoms.set_initial_magnetic_moments([1.0 for a in atoms])
print 'Initial magnetic moment array is', atoms.get_initial_magnetic_moments()
# calculator
calc=espresso(pw=500,
dw=5000,
kpts=(6,6,1),
dipole={'status':True},
outdir='./',
xc='BEEF',
nbands=-30,
spinpol=True,
occupations = 'smearing', # 'smearing', 'fixed', 'tetrahedra'
smearing = 'fd',
sigma = 0.1,
output = {'avoidio':False,
'removewf':True,
'wf_collect':False},
convergence = {'energy':1e-6*13.6,
'mixing':0.1,
'maxsteps':500,
'diag':'david'},
nosym=True,
psppath='/nfs/slac/g/suncatfs/sw/external/esp-psp/v2')
if os.path.exists('qn.traj') and os.path.getsize('qn.traj')!=0:
atoms=io.read('qn.traj',index=-1)
atoms.set_calculator(calc)
constraint = FixAtoms(mask=[(atom.position[2]-min(atoms.positions[:,2]))<1.5*a0/sqrt(3.0) for atom in atoms])atoms.set_constraint(constraint)
print 'Constraint array is ', constraint
#view(atoms)
from ase.io.trajectory import PickleTrajectory
Traj = PickleTrajectory('qn.traj','a',atoms)
dyn = QuasiNewton(atoms,trajectory=Traj,logfile='qn.log',restart='qn.pckl')
dyn.run(fmax=0.05)
if True:
dos = calc.calc_pdos(nscf=True,
kpts=(12,12,1),
Emin=-15.0,
Emax=15.0,
ngauss=0,
sigma=0.2,
DeltaE=0.01,
tetrahedra=False,
slab=True)
import pickle
f = open('dos.pickle', 'w')
pickle.dump(dos, f)
f.close()
if False:
dos = calc.calc_pdos(nscf=True,
kpts=(24,24,1),
Emin=-15.0,
Emax=15.0,
DeltaE=0.01,
tetrahedra=True,
slab=True)
import pickle
f = open('dos_tetra.pickle', 'w')
pickle.dump(dos, f)
f.close()The following script will plot the density of d-states saved in the pickle file.
#######################################################
## Tutorial by: Hongliang Xin
## Topic: Plotting the density of states data from ase-espresso
##
## dependencies: pickle, matplotlib, numpy
##
## Plots the d density of states of a surface Pd atom in Pd(111) slab
#######################################################
#!/usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt
from pylab import *
import os, sys, pickle
def read_dos(dir,tetrahedra=False):
import pickle
try:
if tetrahedra==True:
f = open(dir + '/dos_tetra.pickle')
energies, dos, pdos = pickle.load(f)
f.close()
else:
f = open(dir + '/dos.pickle')
energies, dos, pdos = pickle.load(f)
f.close()
except:
print "No Density of States DATA Found."
sys.exit(1)
return energies, dos, pdos
rcParams['figure.figsize'] = 6*1.67323,4*1.67323
rcParams['ps.useafm'] = True
plt.rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
rcParams['pdf.fonttype'] = 42
matplotlib.rc('xtick.major', size=6)
matplotlib.rc('xtick.minor', size=3)
matplotlib.rc('ytick.major', size=6)
matplotlib.rc('ytick.minor', size=3)
matplotlib.rc('lines', markeredgewidth=0.5*2)
matplotlib.rc('font', size=12*2)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlabel('Energy (eV)')
ax.set_ylabel('PDOS (Arb. Unit)')
num = 12
energies,dos,pdos = read_dos('./')
resov_dos = pdos[num]['d'][0] + pdos[num]['d'][1]
ax.plot(energies,resov_dos,
color='k',
linestyle='-',
label='Pd d DOS')
ax.axis([-10.,10.,0,4.])
plt.xticks([-10.0,-5.0,0.0,5.,10.0])
plt.yticks([0,1,2,3,4])
ax.minorticks_on()
leg=plt.legend(loc=1,prop={'size':24},numpoints=1)
leg.draw_frame(False)
left = 0.125 # the left side of the subplots of the figure
right = 0.95 # the right side of the subplots of the figure
bottom = 0.2 # the bottom of the subplots of the figure
top = 0.95 # the top of the subplots of the figure
wspace = 0.35 # the amount of width reserved for blank space between subplots
hspace = 0.35 # the amount of height reserved for white space between subplots
subplots_adjust(left=left, bottom=bottom, right=right, top=top, wspace=wspace, hspace=hspace)
plt.savefig('dos.png', format='png')
plt.show()Here is the figure generated from the script, showing the d-states.
The pdos variable in the pickle files holds all projected density of states and each channel can be accessed in the following way:
pdos[atom_nr][band][channel]
where atom_nr (int) selects a particular atom, band selects a band 's', 'p', 'd' and channel the component of that band.
Note: The tetrahedra method can be used by turn on the flag in the pdos calculator as shown above. It usually requires more k-points. It is also very important to know the numbering scheme of different channels of density of states.
1. For spin-unpolarized calculations:
for l=1 (p-orbitals):
0 sum
1 pz
2 px
3 py
for l=2 (d-orbitals):
0 sum
1 dz2
2 dzx
3 dzy
4 dx2-y2
5 dxy
2. For spin-polarized calculations:
for l=1 (p-orbitals):
0 sum up
1 sum down
2 pz up
3 pz down
4 px up
5 px down
6 py up
7 py down
for l=2 (d-orbitals):
0 sum up
1 sum down
2 dz2 up
3 dz2 down
4 dzx up
5 dzx down
6 dzy up
7 dzy down
8 dx2-y2 up
9 dx2-y2 down
10 dxy up
11 dxy down