CP2K Open Source Molecular Dynamics

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== Band structure of WO$_3$ Lattice ======


In this exercise, we will carry out 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|}}
<note important>K-point sampling only support GGA functionals.

Band structure calculations using high-level electronic structure theory, such as hybrid functionals are not available in CP2K.
</note>

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
      UKS
      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 3
         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 plt_band.py>
import numpy as np
import matplotlib.pyplot as plt

file=np.genfromtxt('WO3.bs')

bs = open ('WO3.bs','r')
line = bs.readline().split()
num_points, num_k_points, num_bands = int(line[3]),int(line[6]),int(line[8])
sp=[""]*num_points
i=0
for i in range(num_points):
    line = bs.readline().split()
    sp[i] = line[7]


new=np.resize(file,(num_k_points*2,num_bands,3))
ener=new[:,:,1:2]
ener=np.resize(new[:,:,1:2],(num_k_points*2,num_bands)).transpose()
num_homo = len(new[:,:,2][1][new[:,:,2][1]==1])
num_lumo = len(new[:,:,2][1][new[:,:,2][1]==0])

for i in range(num_homo):
    plt.plot(ener[i,::2],color='k')
for i in range(num_lumo):
    plt.plot(ener[i+num_homo,::2],color='k',linestyle='-')

#plt.ylim([-4,6])
plt.xlim([0,num_k_points-1])
plt.xticks(np.linspace(0,num_k_points-1,num_points),sp)
plt.ylabel("Energy (eV)")
plt.show()

</file>