Density of States Calculation

In my line of research, an important tool that we use is the so-called density of states of a system. This is a very important tool for detecting and predicting such properties of a system such as the band gap and other properties. One thing I hate doing is having to remember how to do something. So here is a quick and dirty guide on how to calculate a density of states and then a projected density of states for a given system.

The first step is to have your structure fully relaxed. This has/will be addressed in a different post. Assuming that you have your final structure, there are four calculations that you must run: SCF, NSCF, DOS, and PDOS. In that order. If you are running in the same directory that you did your relax calculation, you can skip the SCF calculation to save time.

Below is the overall script to do all four calculations at once:

```#!/bin/bash
#PBS -N big_solar
#PBS -q batch
#PBS -l nodes=3:ppn=16
source /opt/intel/bin/ifortvars.sh intel64
cd \$PBS_O_WORKDIR
for i in 32; do
cat >solar\$i.in <<!
&control
calculation='scf',
restart_mode='from_scratch',
pseudo_dir = '.',
outdir='.',
prefix='solar',
/
&system
ibrav = 0,
nat = 44,
ntyp = 7,
ecutwfc = 55,
occupations='smearing',
degauss=0.005
input_dft='pbe'
nbnd = 130,
/
&electrons
conv_thr = 1.0e-9
/
&ions
ion_dynamics = 'bfgs'
/
&cell
cell_factor = 10
/
ATOMIC_SPECIES
Ti 47.867 Ti.pbe-sp-van_ak.UPF
O 32.00 O.pbe-rrkjus.UPF
C 12.00 C.pbe-rrkjus.UPF
N 14.00 N.pbe-rrkjus.UPF
H 1.00 H.pbe-rrkjus.UPF
Zr 91.224 Zr.pbe-nsp-van.UPF
S 32.065 S.blyp-van_ak.UPF
ATOMIC_POSITIONS {angstrom}
O 0.201052626 0.117197688 -2.412592252
O 3.284433490 0.128313559 -1.885528864
Ti 3.477028511 -1.371724909 -0.465777813
O 3.367055504 0.117679575 0.638257284
O 1.291476196 -1.379145081 -0.345371479
O 5.420202800 -1.355334499 -0.425299042
Ti 0.115451674 0.113809527 -0.565510518
N 0.073173746 0.032772649 2.001178917
C 1.091955189 0.257194887 2.847475514
C 0.957760874 0.203553204 4.230604641
C -0.286425158 -0.079190634 4.812606085
C -1.335150931 -0.310556683 3.912120124
C -1.114203749 -0.261271424 2.551310465
H 2.053731208 0.468278680 2.375957869
H 1.832812082 0.384723971 4.856552022
H -2.344594941 -0.520738780 4.258310603
H -1.942177672 -0.418961814 1.863974609
S 0.616222323 0.148669389 7.489049423
C -0.549585563 -0.150398897 6.238441034
C -1.768550490 -0.477279017 6.815767535
C -1.757322166 -0.477728022 8.203263296
C -0.527117054 -0.150594118 8.758553901
H -2.657754138 -0.701823036 6.235615304
H -2.633027166 -0.700520683 8.804588172
N 0.142062313 0.028913595 12.949777886
C -1.055124926 -0.216218341 12.443936671
C -1.290774078 -0.275019025 11.083029833
C -0.236429059 -0.076894305 10.180286855
C 1.022857586 0.193996105 10.739973760
C 1.171500558 0.242401856 12.124726073
H -1.905943225 -0.331262720 13.154607926
H -2.307150591 -0.472162299 10.752096645
H 1.890303626 0.371536119 10.100426101
H 2.135699105 0.429119712 12.613559132
O 1.325503235 -1.477854440 15.331959399
Zr 3.397991810 -1.487862763 15.595233668
O 5.476566344 -1.472387174 15.596623898
Zr 0.134495781 0.018318748 16.019557986
O 3.469139781 -0.003731377 14.151748393
O 3.171735574 -0.013688428 16.937523690
O 0.155681673 0.287055234 17.956260162
H -0.045891335 -0.735517777 -2.839244983
H 2.579292516 0.112310278 -2.560928846
H 0.125466118 1.082981469 18.523758919
K_POINTS (automatic)
2 4 1 1 1 1
CELL_PARAMETERS { angstrom }
6.68600 0.00000 0.00000
0.00000 2.95800 0.00000
0.00000 0.00000 \$i
!
echo "a= \$i" ; mpirun pw.x -in solar\$i.in | tee solar\$i.out
done

for i in 32; do
cat >solar\$i.in <<!
&control
calculation='nscf',
restart_mode='from_scratch',
pseudo_dir = '.',
outdir='.',
prefix='solar',
/
&system
ibrav = 0,
nat = 44,
ntyp = 7,
ecutwfc = 55,
occupations='smearing',
degauss=0.005
input_dft='pbe'
nbnd = 130,
/
&electrons
conv_thr = 1.0e-9
/
&ions
ion_dynamics = 'bfgs'
/
&cell
cell_factor = 10
/
ATOMIC_SPECIES
Ti 47.867 Ti.pbe-sp-van_ak.UPF
O 32.00 O.pbe-rrkjus.UPF
C 12.00 C.pbe-rrkjus.UPF
N 14.00 N.pbe-rrkjus.UPF
H 1.00 H.pbe-rrkjus.UPF
Zr 91.224 Zr.pbe-nsp-van.UPF
S 32.065 S.blyp-van_ak.UPF
ATOMIC_POSITIONS {angstrom}
O 0.201052626 0.117197688 -2.412592252
O 3.284433490 0.128313559 -1.885528864
Ti 3.477028511 -1.371724909 -0.465777813
O 3.367055504 0.117679575 0.638257284
O 1.291476196 -1.379145081 -0.345371479
O 5.420202800 -1.355334499 -0.425299042
Ti 0.115451674 0.113809527 -0.565510518
N 0.073173746 0.032772649 2.001178917
C 1.091955189 0.257194887 2.847475514
C 0.957760874 0.203553204 4.230604641
C -0.286425158 -0.079190634 4.812606085
C -1.335150931 -0.310556683 3.912120124
C -1.114203749 -0.261271424 2.551310465
H 2.053731208 0.468278680 2.375957869
H 1.832812082 0.384723971 4.856552022
H -2.344594941 -0.520738780 4.258310603
H -1.942177672 -0.418961814 1.863974609
S 0.616222323 0.148669389 7.489049423
C -0.549585563 -0.150398897 6.238441034
C -1.768550490 -0.477279017 6.815767535
C -1.757322166 -0.477728022 8.203263296
C -0.527117054 -0.150594118 8.758553901
H -2.657754138 -0.701823036 6.235615304
H -2.633027166 -0.700520683 8.804588172
N 0.142062313 0.028913595 12.949777886
C -1.055124926 -0.216218341 12.443936671
C -1.290774078 -0.275019025 11.083029833
C -0.236429059 -0.076894305 10.180286855
C 1.022857586 0.193996105 10.739973760
C 1.171500558 0.242401856 12.124726073
H -1.905943225 -0.331262720 13.154607926
H -2.307150591 -0.472162299 10.752096645
H 1.890303626 0.371536119 10.100426101
H 2.135699105 0.429119712 12.613559132
O 1.325503235 -1.477854440 15.331959399
Zr 3.397991810 -1.487862763 15.595233668
O 5.476566344 -1.472387174 15.596623898
Zr 0.134495781 0.018318748 16.019557986
O 3.469139781 -0.003731377 14.151748393
O 3.171735574 -0.013688428 16.937523690
O 0.155681673 0.287055234 17.956260162
H -0.045891335 -0.735517777 -2.839244983
H 2.579292516 0.112310278 -2.560928846
H 0.125466118 1.082981469 18.523758919
K_POINTS (automatic)
2 4 1 1 1 1
CELL_PARAMETERS { angstrom }
6.68600 0.00000 0.00000
0.00000 2.95800 0.00000
0.00000 0.00000 \$i
!
echo "a= \$i" ; mpirun pw.x -in solar\$i.in | tee solar\$i.out
done

cat > bigsolar.dos.in << EOF
&dos
outdir='.'
prefix='solar'
fildos='solar.dos'
Emin=-20, Emax=25.0, DeltaE=0.1
/

EOF

mpirun dos.x < bigsolar.dos.in | tee bigsolar.dos.out

cat > bigsolar.pdos.in << EOF
&projwfc
outdir='.'
prefix='solar'
io_choice='both'
Emin=-20, Emax=25.0, DeltaE=0.1
ngauss=1, degauss=0.02
/
EOF

mpirun projwfc.x < bigsolar.pdos.in > bigsolar.pdos.out

```

What is nice about this method is you do not have to worry about the number of processors that you used in your SCF calculation. Unless you specify wf_collect (http://www.quantum-espresso.org/wp-content/uploads/Doc/INPUT_PW.html#id1486037) in the control block, every calculation that reuses data from the SCF/NSCF calculations must be run on the same number of processors.

Of note in this is the actual DOS and PDOS calculations. The Emin and Emax should be more or less centered around the fermi level. As a rule of thumb, I add around 20eV to either side of the fermi level to get a good idea of what the DOS looks like, especially around the fermi level. The fermi level is outputted in your scf/relax calculation that should have been completed prior to running a dos calculation:

```
llentz@head:~/solar_grid/dos\$ cat solar.out | grep Fermi

```

In my case, the Fermi level was ~1.3eV for the above system, hence the given emin and emax values. It is of note that it does not add that much computation time to increase the values, you just  will end up with a much longer graph in the end as the states far away from the fermi level are normally of no consequence.

Once this script is complete, you will end up with a bunch of files like this:

```solar.pdos_atm#12(C)_wfc#2(p)
solar.pdos_atm#13(C)_wfc#1(s)
solar.pdos_atm#13(C)_wfc#2(p)
solar.pdos_atm#14(H)_wfc#1(s)
solar.pdos_atm#15(H)_wfc#1(s)
solar.pdos_atm#16(H)_wfc#1(s)
solar.pdos_atm#17(H)_wfc#1(s)
solar.pdos_atm#18(S)_wfc#1(s)
solar.pdos_atm#18(S)_wfc#2(p)
solar.pdos_atm#19(C)_wfc#1(s)
solar.pdos_atm#19(C)_wfc#2(p)
solar.pdos_atm#1(O)_wfc#1(s)
solar.pdos_atm#1(O)_wfc#2(p)
solar.pdos_atm#20(C)_wfc#1(s)
solar.pdos_atm#20(C)_wfc#2(p)
solar.pdos_atm#21(C)_wfc#1(s)
solar.pdos_atm#21(C)_wfc#2(p)
solar.pdos_atm#22(C)_wfc#1(s)
solar.pdos_atm#22(C)_wfc#2(p)
solar.pdos_atm#23(H)_wfc#1(s)
```

The files you will be most interested in will be the *.pdos_* files. We will discuss them in a later blog post about posting them.

1. Samira says:

Dear Levi,

Thank you in advance for your useful notes. I just started learning quantum espresso and your blog is so helpful for me. However, I am wondering about two points. occupations=’smearing’ in NSCF calculation. (I used to think it should be ‘tetrahedra’ !) and the next thing is “io_choice=both” in PDOS . (I couldn’t find it in Projwfc input data. what does it mean?) Would you please it clear to me?

Thank you and best wishes
Samira

1. Hello Samira,

As far as occupations=’smearing’ this is because the systems I am working with have a lot of electrons, and makes the convergence a bit better. If you are simply after a very accurate DOS, tetrahedra also works. I would avoid fixed occupations whenever possible.

I also just looked at the projwfc.x input file and it appears that it now ignores the io_choice input parameter. I cant find something readily available for what it means, but I imagine that it has to do with input/output of the files from the projwfc.x program.

Let me know if you have more questions,

Levi

2. N.khossossi says:

Dear LEVI LENTZ,

Thanks sir, Your tutorials are so beneficial for me. I’m new in python/QE_package so the are questions.
I am trying to compile python script but actually not result is showing up? Any idea what could it be?
i have in my terminal as result :
*********
khossossi@nabil:~/Desktop\$ ./QE.py
./QE.py: line 212: warning: here-document at line 9 delimited by end-of-file (wanted `!’)
./QE.py: line 213: syntax error: unexpected end of file
khossossi@nabil:~/Desktop\$
*********

3. Viet says:

Dear Levi,
Thank you for your post, This is truly useful for a beginner.
I took a confusing that the Fermi level in “smearing” and “tetrahedra” is really different. In Tetrahedra looks more exactly.
And, compared with the band structure calculation that used “smearing” or “fixed”. The Fermi level in DOS (used by “tetrahedra”) is different in Band structure. So, which one is better for Fermi level consideration? and Why do not use Tetrahedra in Band structure calculation?

Thank you.
Viet

4. Samira says:

Dear Levi,
Do you have any script for calculating the d-band center? I have search a lot, but I don’t think there is such a thing in QE? Am I right?