Density of States Calculation

Quantum EspressoLevi Lentz7

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
module load espresso/openmpi-intel
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.

7 comments Add yours

  1. 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

    Reply

    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

      Reply
  2. Pingback: Another pDOS Plotter – Levi Lentz's Blog

  3. 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 7: module: command not found
    ./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$
    *********

    Reply

  4. 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

    Reply

  5. 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?

    Reply

  6. Dear Levi
    Thank you so much! Can you kindly tell me how to carry out DOS using hybrid functionals. NSCF calculations do not account for hybrid functionals!

    Reply

Leave a Reply

Your email address will not be published. Required fields are marked *