

In [ ]:

    
from __future__ import print_function
import sisl
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

Analyzing output from Siesta makes analysis much easier since many things are intrinsically enabled through sisl, i.e. orbital numbering etc. may easily be handled with sisl.

In this example we will use Siesta to calculate the heavy part, i.e. the projected density of states for a fine $k$-point sampling. The formula used to calculate the projected density of states is: \begin{align} \mathrm{PDOS}_\nu(E) &= \int \mathrm d k \sum_i \psi^*_{i,k,\nu} [\mathbf S(k) | \psi_{i,k}\rangle]_\nu D(E-\epsilon_i), \\ \mathrm{DOS}(E) &= \int \mathrm d k \sum_i D(E-\epsilon_i) = \sum_\nu \mathrm{PDOS}_\nu(E), \end{align} with $D(\Delta E)$ being the distribution function (in Siesta this is the Gaussian distribution function). Finally we will also compare with the same calculation using sisl to ensure that correctness of both methods. Note that sisl is (for now) not MPI parallelized and thus for large systems it may be way more efficient to use Siesta to calculate PDOS, but use sisl to post-process the PDOS data.

The system will again be graphene.



In [ ]:

    
graphene = sisl.geom.graphene(1.44)
graphene.write('STRUCT.fdf')
graphene.write('STRUCT.xyz')

To tell Siesta to calculate the projected DOS one has to add these flags to the input (if you want you can play with the numbers):



In [ ]:

    
open('PDOS.fdf', 'w').write("""
# K-point sampling
%block PDOS.kgrid.MonkhorstPack
  63  0 0
   0 63 0
   0  0 1
%endblock
# Energy grid to calculate the PDOS on
# E_min = -20 eV
# E_max = 15 eV
# Gaussian width: \sigma 2 ^{1/2} = 0.2 eV
# Number of points: 3500 (spacing of 0.01 eV)
%block Projected.DensityOfStates
-20.00 15.00 0.200 3500 eV
%endblock
""");

Please open PDOS.fdf and RUN.fdf.

Why do you see two different $k$-point samplings?

Now run Siesta to calculate the electronic structure and the PDOS:

 siesta RUN.fdf | tee RUN.out

Subsequently we will read in the PDOS information and post-process it.
First we will data found in the siesta.PDOS.xml file, 3 quantities will be returned:

The geometry (as found in the XML file), this also contains all atomic species etc.
The energy grid used
Orbital and energy resolved DOS

Before analyzing, please examine the data returned by viewing the shapes and printing the geometry.



In [ ]:

    
geometry, E, PDOS = sisl.get_sile('siesta.PDOS.xml').read_data()



In [ ]:

    
# Examine output by printing shapes and geometry

Exercises

Plot the total DOS (invoke Eq. 2 in the top paragraph)
Plot the total PDOS on the first atom
HINT: Sum all orbitals on the first atom
Figure out which orbital contributes to the Dirac cone by plotting the PDOS for all orbitals
Read in the Hamiltonian (as done in S 1) and create a Monkhorst-Pack grid as the one in Siesta input. Use this snippet to calculate the PDOS using sisl
```
mp = sisl.MonkhorstPack(H, [nx, ny, 1])
pdos = mp.asaverage().PDOS(E, eta=True)
```
Search the API documentation for the MonkhorstPack.asaverage method and figure out what it does. Note, there are other MonkhorstPack.as* methods, these are all extremely usefull when calculating a large quantity of data in a Brillouin-zone object.
Compare the Siesta PDOS with sisl PDOS, why are they different?
HINT: $\sigma$. Check the API for the PDOS method.



In [ ]:

    
plt.plot(E, PDOS.sum(0), label='DOS')
plt.xlabel(r'$E - E_F$ [eV]')
plt.ylabel(r'DOS [1/eV]')

# Add more plots for individual PDOS


# Add legend(s)
plt.legend();