====== 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|}}
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.
To get the band structure for WO3, only a few changes are required compared to the previous example for [[PDOS|calculating the PDOS]]:
&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
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.
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//.
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:
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
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:
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
[...]
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.
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.
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.
gnuplot>set yrange [-8:14]
gnuplot>plot for [i=4:23] "WO3.bs.set-1.csv" u 0:i w l t ""
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()