CP2K Open Source Molecular Dynamics

User Tools

Site Tools

Trace:
exercises:2017_uzh_cmest:pdos

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revisionPrevious revision
Next revision		Previous revision
exercises:2017_uzh_cmest:pdos [2017/10/17 15:46] – jglanexercises:2017_uzh_cmest:pdos [2020/08/21 10:15] (current) – external edit 127.0.0.1
Line 1: Line 1:
-======= Projected density of states for graphene and h-BN =======+======= Projected density of states and Band structure for WO$_3$ =======
  
-In the following exercise we are going to look at the density of states of WO$_3$:+In this exercise, you 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: Similar to the previous exercise we write the coordinates in term of the unit cell:
Line 110: Line 116:
  <code>python get-smearing-pdos.py file.pdos</code>  <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.+Alternatively, you could also use the [[https://raw.githubusercontent.com/dev-zero/cp2k-tools/master/scripts/cp2k_pdos.py|Python script]] developed by Tiziano, the parser bug was fixed already.
  
-<note warning>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>+<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. 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.
Line 118: Line 124:
 While some of the new options to help with convergence are of numerical nature, [[howto:static_calculation#adding_smearing|the smearing is not]]. 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   * Repeat the above calculation for the different multiple cells 3x3x3, 4x4x4
-  * Do you see why it is necessary to do the unit cell replication? Hints: does WO$_3$ have a band gap? Compare the plots for 3x3x3 and 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?   * .. which state ($s$, $p_x$, ..) is mainly responsible for that?
 +  * Change the $\sigma$ value in convolution program, and determine a reasonable value for the PDOS plot
 +</note>
 +
 +====== Getting the band structure of WO$_3$ Lattice ======
 +
 +To get the band structure for <chem>WO3</chem>, only a few changes are required compared to the previous example for [[PDOS|calculating the PDOS]]:
 + 
 +<code - WO3-bs.inp>
 +&GLOBAL
 +   PROJECT WO3-kp-bs
 +   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
 +      &QS
 +         EXTRAPOLATION USE_GUESS ! required for K-Point sampling
 +      &END QS
 +      &SCF
 +         SCF_GUESS ATOMIC
 +         EPS_SCF 1.0E-6
 +         MAX_SCF 300
 +
 +         ADDED_MOS 2
 +         &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
 +      &KPOINTS
 +         SCHEME MONKHORST-PACK 3 3 1
 +         WAVEFUNCTIONS REAL
 +         SYMMETRY .FALSE.
 +         FULL_GRID .FALSE.
 +         PARALLEL_GROUP_SIZE -1
 +      &END KPOINTS
 +      &PRINT
 +         &BAND_STRUCTURE
 +            ADDED_MOS 2
 +            FILE_NAME WO3.bs
 +            &KPOINT_SET
 +               UNITS B_VECTOR
 +               SPECIAL_POINT ???   #GAMA
 +               SPECIAL_POINT ???   #X
 +               SPECIAL_POINT ???   #M
 +               SPECIAL_POINT ???   #GAMA
 +               SPECIAL_POINT ???   #R
 +               SPECIAL_POINT ???   #M
 +               NPOINTS ???
 +            &END
 +         &END BAND_STRUCTURE
 +      &END PRINT
 +   &END DFT
 +
 +   &SUBSYS
 +      &CELL
 +         ABC [angstrom] 3.810000 3.810000 3.810000
 +         PERIODIC XYZ
 +         MULTIPLE_UNIT_CELL 1 1 1
 +      &END CELL
 +      &TOPOLOGY
 +         MULTIPLE_UNIT_CELL 1 1 1
 +      &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>
 +
 +<note important>At present, it is not possible to get the projected density of states when doing a K-Point calculation. The special points should be given in terms of the b-vectors.</note>
 +
 +Some notes on the input file:
 +  * By specifying the ''KPOINT'' section you are enabling the K-Point calculation.
 +  * While you could specify the K-Points directly, we are using the Monkhorst-Pack scheme [(http://journals.aps.org/prb/abstract/10.1103/PhysRevB.13.5188)] to generate them. The numbers following the keyword ''MONKHORST-PACK'' specify the tiling of the brillouin zone.
 +  * After the basic calculation, CP2K calculates the energies along certain lines, denoted as ''KPOINT_SET'' (when you check [[https://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/PRINT/BAND_STRUCTURE/KPOINT_SET.html|the documentation]] you will note that this section can be repeated).
 +  * The keyword ''NPOINTS'' specifies how many points (in the addition to the starting point) should be sampled between two special points.
 +  * The ''SPECIAL_POINT'' keyword is used to specify the start-, mid- and endpoints of the line. Those points usually denote special points in the reciprocal lattice/unit cell, like the $\Gamma$, $M$ or $K$ point. You can find the definition for these in the appendix section of [[http://www.sciencedirect.com/science/article/pii/S0927025610002697|this paper]]. This keyword can also be specified multiple times, making it possible to directly get the band structure for a complete //path//.
 +
 +<note tip>You are encouraged to use [[ http://tools.materialscloud.org/seekpath|SeeK-path Tool]] when doing the sampling via K-Points to get a skeleton input file for CP2K with the important paths in the reciprocal space. Give SeeK-path the following as your XYZ file and specify a simple cubic cell with the lattice constant $a$ as specified below as well:
 +<code - WO3-cubic.xyz>
 +4
 +WO3; a=3.810000
 + W    0.000000    0.000000    0.000000
 + O    1.905000    0.000000    0.000000
 + O    0.000000    1.905000    0.000000
 + O    0.000000    0.000000    1.905000
 +</code>
 +</note>
 +
 +
 +Now, when you run this input file you will get in addition the the output file, a file named ''WO3.bs'' which will look similar to the following:
 +
 +<code>
 + SET:       1                 TOTAL POINTS:      26
 +   POINT   1                     ********    ********    ********
 +   POINT   2                     ********    ********    ********
 +   POINT   3                     ********    ********    ********
 +   POINT   4                     ********    ********    ********
 +   POINT   5                     ********    ********    ********
 +   POINT   6                     ********    ********    ********
 +       Nr.    1    Spin 1        K-Point  0.00000000  0.00000000  0.00000000
 +               20
 +           -73.66652408    -38.53370023    -37.80464132    -37.79327769
 +           -16.71308703    -16.11075946    -16.02553853     -1.43495530
 +            -1.34739188     -1.33357408      0.37912017      0.38948689
 +             0.39582882      0.40030859      0.46965212      0.47418816
 +             2.60728842      2.62105342      3.16044140      6.99806305
 +       Nr.    2    Spin 1        K-Point  0.00000000  0.10000000  0.00000000
 +               20
 +           -73.66647294    -38.53337818    -37.80859042    -37.79536623
 +           -16.67479677    -16.09554462    -15.96731960     -1.68492873
 +            -1.44087258     -1.34318045      0.09257368      0.13769271
 +             0.21643888      0.38447849      0.44179002      0.45757924
 +             2.61768501      3.02368022      3.51828287      7.06644645
 +
 +[...]
 +</code>
 +
 +For each set there is a block named ''SET'' with the special points listed as ''POINT'', followed by sub-blocks for each K-Point containing the energies for each MO.
 +
 +<note>
 +Your tasks:
 +
 +  * Lookup the special points for the $\Gamma$, $X$,$M$,$R$ points in the [[http://pubs.acs.org/doi/abs/10.1021/cm3032225|mentioned paper]] (make sure you choose the right lattice). Calculate and plot the band structure for WO3 lattice along $\Gamma$, $X$,$M$,$\Gamma$,$R$,$M$ (you are free to decide whether to use multiple K-Point sets are multiple special points in a single set). Mark the special points. Choose an appropriate number of interpolation points to get a smooth plot.
 +  * Compare your plot with plots from literature. What is different?
 +  * How many orbital energies do you get and why? Try to change the input to get more unoccupied orbitals.
 +</note>
 +
 +
 +To convert the band structure file to a file which can be plotted directly, you can use the script ''cp2k_bs2csv.py'' from below, which when passed a band structure file ''WO3.bs'' as an argument will write files ''WO3.bs-set-1.csv'' for each set containing the K-Point coordinates and the energies in one line.
 +
 +To plot the ''WO3.bs-set-1.csv'' file, you can either load it into $MATLAB$ or use $GNUPLOT$ command line.
 +<code>gnuplot>set yrange [-8:14]
 +gnuplot>plot for [i=4:23] "WO3.bs.set-1.csv" u 0:i w l t ""
 + </code>
 +<file python cp2k_bs2csv.py>
 +#!/usr/bin/env python
 +"""
 +Convert the CP2K band structure output to CSV files
 +"""
 +
 +import re
 +import argparse
 +
 +SET_MATCH = re.compile(r'''
 +[ ]*
 +  SET: [ ]* (?P<setnr>\d+) [ ]*
 +  TOTAL [ ] POINTS: [ ]* (?P<totalpoints>\d+) [ ]*
 +  \n
 +(?P<content>
 +  [\s\S]*?(?=\n.*?[ ] SET|$)  # match everything until next 'SET' or EOL
 +)
 +''', re.VERBOSE)
 +
 +SPOINTS_MATCH = re.compile(r'''
 +[ ]*
 +  POINT [ ]+ (?P<pointnr>\d+) [ ]+ (?P<a>\S+) [ ]+ (?P<b>\S+) [ ]+ (?P<c>\S+)
 +''', re.VERBOSE)
 +
 +POINTS_MATCH = re.compile(r'''
 +[ ]*
 +  Nr\. [ ]+ (?P<nr>\d+) [ ]+
 +  Spin [ ]+ (?P<spin>\d+) [ ]+
 +  K-Point [ ]+ (?P<a>\S+) [ ]+ (?P<b>\S+) [ ]+ (?P<c>\S+) [ ]*
 +  \n
 +[ ]* (?P<npoints>\d+) [ ]* \n
 +(?P<values>
 +  [\s\S]*?(?=\n.*?[ ] Nr|$)  # match everything until next 'Nr.' or EOL
 +)
 +''', re.VERBOSE)
 +
 +if __name__ == '__main__':
 +    parser = argparse.ArgumentParser(description=__doc__)
 +    parser.add_argument('bsfilename', metavar='bandstructure-file', type=str,
 +                        help="the band structure file generated by CP2K")
 +
 +    args = parser.parse_args()
 +
 +    with open(args.bsfilename, 'r') as fhandle:
 +        for kpoint_set in SET_MATCH.finditer(fhandle.read()):
 +            filename = "{}.set-{}.csv".format(args.bsfilename,
 +                                              kpoint_set.group('setnr'))
 +            set_content = kpoint_set.group('content')
 +
 +            with open(filename, 'w') as csvout:
 +                print(("writing point set {}"
 +                       " (total number of k-points: {totalpoints})"
 +                       .format(filename, **kpoint_set.groupdict())))
 +
 +                print("  with the following special points:")
 +                for point in SPOINTS_MATCH.finditer(set_content):
 +                    print("  {pointnr}: {a}/{b}/{c}".format(
 +                        **point.groupdict()))
 +
 +                for point in POINTS_MATCH.finditer(set_content):
 +                    results = point.groupdict()
 +                    results['values'] = " ".join(results['values'].split())
 +                    csvout.write("{a} {b} {c} {values}\n".format(**results))
  
 +</file>
exercises/2017_uzh_cmest/pdos.1508255201.txt.gz · Last modified: 2020/08/21 10:15 (external edit)