Analysis tools¶
Hint
Enter the directories 4ZGNR and/or 10ZGNR.
In this exercise we will cover a few basic/intermediate tools to extract important information from the results of a SIESTA calculation.
The directory contains two subdirectories: one for a thin graphene ribbon with
4 zigzag chains and H-saturated edges (4ZGNR), and another for a thicker 10 chain
ribbon (10ZGNR). If you want to try a faster example, move to the 4ZGNR directory.
The ground state of these ribbons has an antiferromagnetic (AF) spin configuration
where each edge of the ribbon has opposite spin to the other. The DM.ÌnitSpin block
is used to converge the proper AF state, rather than the FM configuration.
The input file ‘zgnr.fdf’ has a few lines that will have to be uncommented during the exercise.
Note
This is a 1D system, which has periodicity only along the x-axis.
However, the %block kgrid_Monkhorst_Pack is set for a very smooth k-sampling
along that direction. This can make the SCF loop quite slow. You can decrease
the number of k-points for a first steps to get a converged DM and then increase
the k-sampling for analysing the results.
Band structure and density of states¶
We first analyse the band-structure for the graphene ribbon and look at the density of states.
To plot the DOS you will use the utility program ‘Eig2DOS’ (see
this how-to for background and installation
notes, if needed). This program reads the file zgnr.EIG (always
produced by SIESTA), which contains all the eigenvalues for each k-point used
to sample the BZ, and the file zgnr.KP, which contains the k-point
sampling information. Useful options to the program (type Eig2DOS
-h for a full list) are the broadening for each state in eV (a value
of the order of the default (0.2 eV) is usually reasonable), the
number of energy points where the DOS will be calculated (200 by
default) and the Emin and Emax of the energy window where the DOS will
be calculated (the default is to compute the DOS in the whole range of
energies available in the EIG file).
For example:
Eig2DOS -k zgnr.KP zgnr.EIG > dos
will compute the DOS in the whole range, using a grid of 200 points, a broadening (“smearing”) of 0.2 eV, and using the k-point info from zgnr.KP. (The k-point file information is important when different k-points have different weights).
Plot the dos using gnuplot, type gnuplot and enter the command:
plot "dos" with lines
The result is a bit confusing, but it is just because the range of energies is very large. Open the zgnr.out file and look for the Fermi energy. EF is also written at the top of the ‘zgnr.EIG’ file.
You can improve the results around the Fermi level (the option ‘-f’ shifts the energies so that the Fermi level is at 0 eV):
Eig2DOS -f -e -3 -E 3 -s 0.05 -k zgnr.KP zgnr.EIG >dos
The input file zgnr.fdf also requests SIESTA to print the information of the
bandstructure along high-symmetry lines in the Brillouin zone (BZ). As this is
a 1D system, the structure is very simple:
WriteKbands .true.
WriteBands .true.
BandLinesScale ReciprocalLatticeVectors
%block Bandlines
1 0.0000000000 0.000000000 0.0000 \Gamma # Gamma is the first point
100 0.5000000000 0.000000000 0.0000 M # 100 points from Gamma and we get to M
%endblock BandLines
Let’s understand these few lines:
The BandLinesScale entry specifies the scale of the k-vectors given in the BandLines block. The possible options are “pi/a” and “ReciprocalLatticeVectors”.
Each line of the BandLines block represents the end-point of a segment to be spanned in the reciprocal space. The comments on the code block above help to familiarize with the format of these k-space segments.
To learn more about the special high-symmetry points and lines of the BZ for this case, you can visit the Bilbao Crystallographic Server (click on ‘Optimized listing of k-vector types…’; then on ‘Brillouin Zone’).
The bands are stored in a the file zgnr.bands.
To plot the band structure you can use the utility program gnubands:
gnubands zgnr.bands > bands.dat
Again, the result covers the whole range of energies, which is too large.
We can improve the results around the Fermi level with a few modifiers
on gnubands
gnubands -G -o bands -F -e -3 -E 3 zgnr.bands
and you can plot with gnuplot, using the bands.gplot
file generated (which has extra information to place and label the
‘ticks’ for each symmetry point in the BZ):
gnuplot -persist bands.gplot
Plotting Projected Density of States¶
There are several utilities that can be used to process the projected DOS,
such as pdosxml, sisl,
or the fmpdos program by Andrei Postnikov. Here we will use mprop,
an utility that is included in the framework for the analysis of the
crystal-overlap populations (COOP/COHP),
which includes the processing of the pDOS as a by-product.
One potential drawback of this approach is that it requires files
storying the wavefunction in the relevant energy range, which for some
systems could be large. The files are generated with the option:
COOP.write T
Hint
Uncomment that line in the input file zgnr.fdf
Once you have run SIESTA with that flag in the fdf,
it will generate a new file zgnr.fullBZ.WFSX. To process the information try this:
ln -sf zgnr.fullBZ.WFSX zgnr.WFSX
mprop [option] pdos
that uses the selection of orbitals’ projection in the file pdos.mpr
Note
In this example, pdos.mpr does the projection on the orbitals of the carbon atoms at the edges of the graphene ribbon. For the 4ZGNR it looks like this:
zgnr # System Label
DOS # ... used for PDOS (or FatBands)
Cedge1 # flag for output file for selected orbitals
2 # index of orbitals selected
Cedge2 # flag for output file for selected orbitals
9 # index of orbitals selected
You could add more lines to plot PDOS over different orbital sets if you want, or
be more specific in your orbital selection. For example: i) projection over the
pz orbitals at edge1 can be obtaining by typing 2_pz in line 4; ii) H_1s will
project over all 1s orbitals in Hydrogen (at both edges).
If you want further information on how to define orbitals, try with mprop -h
You can select the range of energy (eV) in which you plot the DOS, the number of energy points, or the smearing (Gaussian) used:
mprop -n 500 -m -10 -M 3 -s 0.1 pdos
This generates files pdos.*label*.pdos with the projected DOS.
In those files, the first column is the energy (in eV), the following columns contain the
different spin components, and the last column has the spin-sum over the orbitals.
The total DOS (considering all orbitals) is also saved into zgnr.ados.
You can easily plot the projections with gnuplot:
plot "pdos.Cedge1.pdos" u 1:2 w l, "pdos.Cedge1.pdos" u 1:(-column(3)) w l lt 1, \
"pdos.Cedge2.pdos" u 1:2 w l, "pdos.Cedge2.pdos" u 1:(-column(3)) w l lt 3
Fat-bands for Graphene Nanoribbons¶
Now we will use the fat program to compute the projections of
wavefunctions obtained from band-structure calculation with the option
Write.WFS.For.Bands onto specified orbital sets. The resulting
“fatbands” can later be plotted as a weighthed bandstructure in which
contributions of the selected orbitals will be revealed.
Hint
Modify the input file for 4ZGNR, or 10ZGNR, to uncomment the lines requesting the bandstructure calculation
In any case, in the input file zgnr.fdf, you can see a slightly modified band-structure section:
BandLinesScale ReciprocalLatticeVectors
WFS.Write.For.Bands T # For fat-bands analysis
WFS.band.min 10
WFS.band.max 22
%block Bandlines
1 0.000 0.000 0.000 \Gamma # Begin at Gamma
100 0.500 0.000 0.000 X # 60 points to BZ edge
%endblock Bandlines
Note the WFS.* options, requesting the output of wavefunctions for
specific bands (optional). Without them, the program would only generate
a .bands file with band eigenvalues.
Once we run SIESTA with this input file, we will have a
zgnr.bands.WFSX file, which we can process with:
ln -sf zgnr.bands.WFSX zgnr.WFSX
fat [options] fatbands
using the orbital sets in the pdos.mpr file.
Hint
You might need to remove the zgnr.WFSX link if you have used it for the PDOS example.
Note
The input file for fat is similar to the input file for mprop
used to plot the Partial Density of States (PDOS).
Then, a number of .EIGFAT files will be produced. These files contain
both the eigenvalues and the projection weight for each orbital set
in the .mpr file, and can be re-processed by the eigfat2plot program
(in Util/Bands in the Siesta distribution) to produce data files
suitable for plotting by gnuplot
eigfat2plot pdos.Cedge1.EIGFAT > Cedge1.dat
eigfat2plot pdos.Cedge2.EIGFAT > Cedge2.dat
The structure of the output files has 4 columns: k-point index, eigenEnergy, weight, and spin.
The required commands with gnuplot are:
set yrange[0:-8]
plot "Cedge1.dat" using 1:2:(3*column(3)) w p pt 6 ps variable
replot "Cedge1.dat" using 1:2:(3*column(3)) w l palette lw 4
This will plot both spin components simultaneously. To select a specific spin component you can try:
plot "Cedge1.dat" using 1:2:(column(3)*(2-column(4))) w p pt 6 ps variable
replot "Cedge1.dat" using 1:2:(column(3)*(column(4)-1)) w p pt 6 ps variable
Alternatively, users might want to produce their own post-processors.
With these commands you can check how the valence band for spin + is localized at edge1, and the spin - is localized at edge2. For the conduction band, the spin + is localized at edge2, and the spin - is localized at edge1.