======= Projected density of states for WO$_3$ =======

In this exercise, we will carry out Density Of States(DOS) and band structure calculation using K-point sampling for Cubic lattice WO$_3$. The reference DOS and band structure you can find in [[http://pubs.acs.org/doi/abs/10.1021/cm3032225|this paper]]

{{:exercises:2017_uzh_cmest:wo3.jpeg?1200|}}

====== Getting the PDOS ======

In the following exercise we are going to look at the density of states of <chem>WO3</chem>:

Similar to the previous exercise we write the coordinates in term of the unit cell:

<code - WO3-PDOS.inp>
&GLOBAL
   PROJECT WO3-pdos
   RUN_TYPE ENERGY
   PRINT_LEVEL MEDIUM
&END GLOBAL

&FORCE_EVAL
   METHOD Quickstep
   &DFT
      BASIS_SET_FILE_NAME  BASIS_MOLOPT
      POTENTIAL_FILE_NAME  POTENTIAL

      &POISSON
         PERIODIC XYZ
      &END POISSON

      &SCF
         SCF_GUESS ATOMIC
         EPS_SCF 1.0E-6
         MAX_SCF 300
         ADDED_MOS 100
         &DIAGONALIZATION
            ALGORITHM STANDARD
            EPS_ADAPT 0.01
         &END DIAGONALIZATION
         &SMEAR  ON
            METHOD FERMI_DIRAC
            ELECTRONIC_TEMPERATURE [K] 300
         &END SMEAR

         &MIXING
            METHOD BROYDEN_MIXING
            ALPHA 0.2
            BETA 1.5
            NBROYDEN 8
         &END MIXING

      &END SCF
      &XC
         &XC_FUNCTIONAL PBE
         &END XC_FUNCTIONAL
      &END XC
      &PRINT
         &PDOS
            # print all projected DOS available:
            NLUMO -1
            # split the density by quantum number:
            COMPONENTS
         &END
      &END PRINT
   &END DFT

   &SUBSYS
      &CELL
         ABC [angstrom] 3.810000 3.810000 3.810000
         PERIODIC XYZ
         MULTIPLE_UNIT_CELL 2 2 2
      &END CELL
      &TOPOLOGY
         MULTIPLE_UNIT_CELL 2 2 2 
      &END TOPOLOGY
      &COORD
         SCALED
         W 0.0 0.0 0.0
         O 0.5 0.0 0.0
         O 0.0 0.5 0.0
         O 0.0 0.0 0.5
      &END
      &KIND W
         ELEMENT W
         BASIS_SET DZVP-MOLOPT-SR-GTH
         POTENTIAL GTH-PBE
      &END KIND
      &KIND O
         ELEMENT O
         BASIS_SET DZVP-MOLOPT-SR-GTH
         POTENTIAL GTH-PBE
      &END KIND

   &END SUBSYS

&END FORCE_EVAL

</code>

The replication of the unit cell is necessary since the program samples only at the $\Gamma$ point unless instructed otherwise and we will otherwise do get a meaningful sampling of the density of states (e.g. the grid over the Brillouin Zone will be too coarse). Another option (which we will look into in the next exercise) is to sample over k-points instead.

What you will get in addition to the output file is a file named ''WO3_pdos-k1-1.pdos'' (to be precise, you will get one such file per atom kind but here we only have one, carbon) with a content similar to:

<code>
# Projected DOS for atomic kind W at iteration step i = 0, E(Fermi) =     0.144475 a.u.
#     MO Eigenvalue [a.u.]      Occupation                 s                py                pz                px               d-2               d-1                d0               d+1               d+2               f-3               f-2               f-1                f0               f+1               f+2               f+3
       1         -2.621088        2.000000        0.87115225        0.00000000        0.00000000        0.00000000        0.00000000        0.00000000        0.00000256        0.00000000        0.00000256        0.00000000        0.00000000        0.00000000        0.00000000        0.00000000        0.00000000        0.00000000
       2         -2.621080        2.000000        0.85006340        0.00131878        0.00166813        0.00134861        0.00000000        0.00000000        0.00515023        0.00000000        0.00441289        0.00006412        0.00000000        0.00003847        0.00012977        0.00003934        0.00000000        0.00006557

[...]
</code>

The columns correspond to the orbitals present in the basis set (hence //projected// DOS). You would now do a convolution plot using a gaussian to get a smooth DOS

To the convoluted DOS, you may want to check [[http://wiki.wpi.edu/deskinsgroup/Density_of_States| this website]]. Here, it is provided two Python script to do the convolution. 
Download two files [[http://wiki.wpi.edu/images//images/b/b7/Pdos.txt|pdos.py]] and [[http://wiki.wpi.edu/images//images/4/41/Get-smearing-pdos.txt| get-smearing-pdos.py]] to your folder. And execute the Python script using 
 <code>python get-smearing-pdos.py file.pdos</code>

Alternatively, you could also use the [[https://raw.githubusercontent.com/dev-zero/cp2k-tools/master/scripts/cp2k_pdos.py|Python script]] developed by Tiziano: <code>python cp2k_pdos.py -h</code>

<note important>Different $\sigma$ values give you different convolution, which mean the lineshape is different. A reasonable $\sigma$ value is required to get a good PDOS plot. When visualize the PDOS, only energy region close to the Fermi level is interesting. One need to adapt the xrange properly.</note>

Please also note the unit of the energy, it is in $E_h$. When looking at DOS plots you may want to convert it to Electronvolt instead. In the convolution program, this has been done in the code.

While some of the new options to help with convergence are of numerical nature, [[howto:static_calculation#adding_smearing|the smearing is not]].


<note>
  * Repeat the above calculation for the different multiple cells 3x3x3, 4x4x4
  * Get the Total DOS and PDOS of O$_2p$ and W$_5d$ orbitals and compare to the literature value.
  * Do you see why it is necessary to do the unit cell replication?  
  * What is the value of WO$_3$ band gap? Compare the plots for 3x3x3 and 4x4x4.
  * .. which state ($s$, $p_x$, ..) is mainly responsible for that?
  * Change the $\sigma$ (or $width$) value in convolution program, what effect does it have on the PDOS plot ? 
</note>