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.
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
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
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$
*********
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
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?
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!