# Plotting density of states of Fe(BCC) using Quantum ESPRESSO

I have run the DFT calculation and plotted DOS of Fe(BCC) in Quantum ESPRESSO. I have attached the plot: (click on image to see clearly)

But It's a wrong plot as it is not matched with the referenced DOS plots of Fe(BCC). Fermi energy is around 18 eV from the calculation. Now, I am stuck about how to get the right plot. I have already converged lattice constant and magnetic moment with plane wave numbers. I have also attached the PWSCF input files to perform the calculation.Now I have following questions:

1)How can I plot DOS of Fe(BCC) rightly(shape of the curve is completely wrong)?
2)Is there any changes which I have to do in my input files?
3)Since it is a ferromagnetic system, DOS at Fermi level should be relatively low for spin down electrons, which is not appearing in the graph. Why is so?

Input file for SCF calculation:

&CONTROL
title = 'Magnetization of Fe' ,
calculation = 'scf' ,
outdir = '.' ,
pseudo_dir = '.' ,
prefix = 'Fe' ,
etot_conv_thr = 1.0D-6 ,
forc_conv_thr = 1.0D-6 ,
tprnfor = .true. ,
tstress = .true. ,
/
&SYSTEM
ibrav = 3,
A = 2.84197 ,
nat = 1,
ntyp = 1,
nbnd = 30,
ecutwfc = 100 ,
ecutrho = 1000 ,
occupations = 'smearing' ,
degauss = 0.05 ,
smearing = 'gaussian' ,
nspin = 2 ,
starting_magnetization(1) = 0.1,
/
&ELECTRONS
conv_thr = 1.0D-6 ,
/
&IONS
/
&CELL
/
ATOMIC_SPECIES
Fe   55.84500  Fe.pbe-spn-kjpaw_psl.0.2.1.UPF
ATOMIC_POSITIONS alat
Fe      0.000000000    0.000000000    0.000000000
K_POINTS automatic
10 10 10   1 1 1


Input file for nscf calculation:

&CONTROL
title = 'Magnetization of Fe' ,
calculation = 'nscf' ,
outdir = '.' ,
pseudo_dir = '.' ,
prefix = 'Fe' ,
etot_conv_thr = 1.0D-6 ,
forc_conv_thr = 1.0D-6 ,
tprnfor = .true. ,
tstress = .true. ,
/
&SYSTEM
ibrav = 3,
A = 2.84197 ,
nat = 1,
ntyp = 1,
nbnd = 30,
ecutwfc = 100 ,
ecutrho = 1000 ,
occupations = 'smearing' ,
degauss = 0.05 ,
smearing = 'gaussian' ,
nspin = 2 ,
starting_magnetization(1) = 0.1,
/
&ELECTRONS
conv_thr = 1.0D-6 ,
/
&IONS
/
&CELL
/
ATOMIC_SPECIES
Fe   55.84500  Fe.pbe-spn-kjpaw_psl.0.2.1.UPF
ATOMIC_POSITIONS alat
Fe      0.000000000    0.000000000    0.000000000
K_POINTS automatic
20 20 20   1 1 1

• Can you clarify what is wrong? I dont look at DOS plots often so maybe it's obvious. – Tristan Maxson Oct 23 '20 at 21:44
• I agree with @TristanMaxson, you should clarify what is wrong with the calculation. – ProfM Oct 23 '20 at 21:51
• @UjjawalM. if this is a ferromagnetic system, have you checked that you have a non-zero magnetic moment at the end of the DFT calculation? If you don't then the problem isn't with the DOS part of the calculation, but with the DFT calculation itself. – ProfM Oct 24 '20 at 7:28
• I am not familiar with QE input files, But the appearance of your DOS may be dominated by a broadening of the DOS at every calculated point. Try to find out which parameter is responsible for this and reduce the broadening. Also please consider a more reasonable energy range, e.g., +-10 eV around the Fermi level. I think your DOS does not look too bad. The resolution is just very low. – Gregor Michalicek Oct 24 '20 at 18:52
• @Ujjawal M.: From the documentation of the QE input file it is not clear to me whether degauss fulfills this purpose or if it only affects the occupations. But you may try out "tetrahedra" instead of "smearing" for you occupations parameter. (I don't know whether you have to perform further adaptions then.) – Gregor Michalicek Oct 25 '20 at 8:43

After some testing, I've found that the discrepancy isn't from pseudopotentials, number of atoms in the unit cell, or any other of your calculation parameters. I tested a few things like k-point mesh, energy cutoff, and swapping pseudopotentials. I also found you could sometimes converge to a non-magnetic ground state with certain combinations unless you set the starting magnetization to a higher value than the 0.1 you use in your calculation (0.4 worked). I also reduced the number of bands in your calculation input file. Is there any reason you need so many in a system with 16 electrons? 1.3 * (1/2 # of electrons) usually works fine. At most I've ever set is # bands = # electrons.

First, pseudopotentials. Your reference calculation uses a LOW accuracy ultrasoft pseudopotential, from PSLibrary 1.0.0. The target energy cutoff is around 45 Ry. Your calculation uses PAW, also from PSLibrary (but the older, more recommended version from the SSSP library), and is the higher-accuracy version with larger cutoff (typically over 75-80 Ry). This shouldn't really explain the large discrepancy in DOS in itself, unless something was seriously wrong with the pseudopotential.

Next, the unit cell. You used the primitive BCC cell, while the reference uses the conventional cell. QE uses symmetry to reduce the computational workload, and these two approaches should be completely equivalent for a ferromagnetic system where you only need one type of atom.

I did some quick test calculations on my laptop. The issue comes from the nscf calculation in the reference. There is likely something unconverged in the reference calculation you used. From my replication test, the nscf Fermi energy difference vs. the scf calculation, as well as the cell pressure of the final vc-relax verification scf step being larger than the minimum converged value in the final bfgs step, both indicate that the calculations are possibly not converged with respect to k-points and/or energy cutoff.

You can see in this plot, the weird DOS comes up after the reference's nscf calculation. Your scf+nscf calculation is in agreement with the reference's scf calculation. I'm not sure why the plot on that website seems to look more like the scf result rather than the nscf, but I ran the input files directly from how they were provided on that blog post.

The lesson: don't trust a calculation just because it's been posted online in a tutorial. I think it's just for instruction, not rigorously tested for convergence and possible issues stemming from that.

• Your answer has clarified all my doubts.This time, I have converged carefully, parameters like kinetic energy cutoff, k-points mesh, optimised lattice constant, starting magnetization(0.6) and then plotted DOS using occupations as tetrahedra and yes it worked for me now. I think before I was having problem with gaussian broadening. Thanks for the lesson, too. – Ujjawal M. Oct 26 '20 at 17:59
• I have one more question, if I am performing spin polarization calculation in a crystal with different atomic species (e.g. Co2MnSi), do I have to set starting magnetization for all the different atomic species seperately? – Ujjawal M. Oct 26 '20 at 18:07
• @UjjawalM. Did you use "tetrahedra" for the occupations? It is designed to work with semiconductors, not with metal. Try to use MV cold smearing, it is more "stable" than the gaussian method. Using the gaussian occupation requires checking for the convergency against the smearing value. – Anibal Bezerra Oct 26 '20 at 19:23
• @UjjawalM. Basically yes. I set starting_magnetization per atomic species based on the magnetic ordering I'm trying to converge to. Keep in mind it's an initial condition to nudge things in the right direction, it's not a constraint (like tot_magnetization). – Kevin J. M. Oct 26 '20 at 21:08
• @Anibal Bezerra, If I used 'cold' smearing, do I have to converge gaussian broadening again by changing 'degauss' unlike the 'tetrahedra' method that does it by itself? – Ujjawal M. Oct 27 '20 at 7:58

I have noticed the error, I believe. The example you link in the comments has two Fe atoms whereas your model only has one. These expected results seem fairly consistent with an antiferromagnetic/ferrimagnetic system I believe, which will not be possible to model in a single atom representation.

Please try to make a 2 atom unit cell and see if this fixes your problems. This is a very surface level guess, but we can refine this answer if this turns out not to be the problem. Also with such a small cell, maybe you can avoid the NSCF run altogether and just run a SCF calculation at 20x20x20 with your optimized geometry.

Here is an exercise from the GPAW documentation as some additional reading.

Suggestion 2:

Looking at the example you linked, I see they are using an ultrasoft potential and you are using a PAW potential. Reading on the mailing list I see this can maybe have trouble with AFM structures, but as you note this should converge to an FM structure. Can you check in the final output what the magnetization actually is? Also try an ultrasoft potential.

• @Tristan Maxson, Thanks for your response! I have again plotted DOS curve with 2 atoms in a unit cell but still I am getting the same plot as above. I have edited my question accordingly. What should I do now? – Ujjawal M. Oct 24 '20 at 5:47
• @UjjawalM. See new edit – Tristan Maxson Oct 24 '20 at 17:40
• @Tristan Maxson, Thanks for your response, total magnetization is 2.11 Bohr mag/cell for the primitive unit cell(1 atom) and 4.22 Bohr mag/cell for the conventional unit cell (2 atoms). Can you tell me how ultrasoft potential is differed by PAW potentials? – Ujjawal M. Oct 25 '20 at 8:44
• I don't actually know much about the difference personally, I just know they differ between your example and your input. – Tristan Maxson Oct 25 '20 at 11:45
• @UjjawalM. Actually, your current plot looks mostly correct? You might want to flip it so both are overlaid to be able to directly compare but I definitely see a shift in peaks as expected across the fermi level that looks similar. If that was from the PAW potentials it seems to have solved it (any other differences might be from broadening). – Tristan Maxson Oct 25 '20 at 12:08