:sequential_nav: next .. _tutorial-basic-Analysis-tools: 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 :ref:`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: .. code-block:: gnuplot 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``: .. code-block:: bash 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`` .. code-block:: bash 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): .. code-block:: bash 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 :ref:`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``: .. code-block:: 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`` .. code-block:: bash 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: .. code-block:: gnuplot 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: .. code-block:: gnuplot 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.