1. GW approximation
In DFT most energy eigenvalues don't have a strict physical meaning. This is different for eigenvalues obtained from the GW approximation to many-body perturbation theory. In this tutorial we will use the Spex code to perform single-shot GW calculations on top of converged DFT results to calculate band structures and other quantities. Note that there is also the possibility to perform quasiparticle self-consistent GW calculations by running fleur and Spex repeatedly after each other.
1.1. Obtaining and installing Spex
Download the Spex tarball and copy it to your preferred spex source directory. Unpack the tar with the command 'tar -xf'. In the now created subdirectory invoke the configure script with
MPIFC=mpiifort FCFLAGS=-I/usr/local_rwth/sw/HDF5/1.10.4/intel-intelmpi/include LDFLAGS=-L/usr/local_rwth/sw/HDF5/1.10.4/intel-intelmpi/lib ./configure --with-dft=fleur --enable-mpi=yes
and compile the code with 'make -j4'. You should now have the executables spex.inv (for systems with inversion symmetry) and spex.noinv (without inversion symmetry) in the src subdirectory.
1.2. Si
We use our well-known Si inpgen input
Si bulk &lattice latsys='cF', a0=1.8897269, a=5.43 / 2 14 0.125 0.125 0.125 14 -0.125 -0.125 -0.125
To generate an inp.xml file. Please use the command 'inpgen -gw < your_inpgen_input'. This will include default options in the inp.xml file to easily run fleur together with spex, for example a possibility to read an external kpoint file. Notice that the /calculationSetup/expertModes/@gw parameter is set to 1 now.
Now start with converging the fleur calculation.
Also for Spex we need to set up an input file, which has to have the name 'spex.inp'. A very minimal option is:
BZ 4 4 4 KPT X=1/2*(0,1,1) K=1/8*(3,6,3) #KPTPATH (K,1,X) 40 JOB GW IBZ:(1-10) NBAND 80
where the lines have the following meaning:
| line | meaning |
|---|---|
| BZ 4 4 4 | Size of the kpoint mesh in the Brillouin zone. Note that this kpoint mesh is typically smaller then in DFT, as it is computationally expensive to use a fine kpoint mesh here. |
| KPT X=1/2(0,1,1) K=1/8(3,6,3) | Kpoint labels and position for extra kpoints are given. The Gamma point does not need to be specified as it is the default with the label '1'. Note that the kpoints can be given in squared bracket (cartesian coordinates) or in round brackets for internal coordinates. |
| #KPTPATH (K,1,X) 40 | This line defines the kpoint path for the band structure calculation later on. Therefore it is commented out with a '#'. It defines the path L,Gamma,X together with a step size. A larger number means a finer grid, a more detailed description can be found in the manual. |
| JOB GW IBZ:(1-10) | Defines the job Spex is supposed to do. GW refers to the GW approximation, IBZ refers to all kpoints in the Brillouin zone and (1-10) refers to the band indices for which the quasiparticle correction is evaluated. Note that you could also specify certain kpoints or bands. |
| NBAND 80 | In order to calculate the gw approximation all occupied as well as all unoccupied bands should be taken into account. Since this is impossible one needs to define a cutoff for the number of bands used. This is a parameter that needs to be converged for a precise calculation. |
For a complete description of the spex.inp file please consult the manual.
Now run 'spex.inv -w' so that the 'kpts_gw' file is written (this is done by -w), open the 'inp.xml' once again and add the number of bands to calculate (/calculationSetup/cutoffs/@numbands). This number should be a little higher than the number of bands given in the 'spex.inp' file. Also change the GW input parameter in /calculationSetup/expertModes/@gw to 2. This directs fleur to write out all the wavefunctions for the kpoints given in 'kpts_gw'. Run fleur.
When fleur has finished everything is ready for the actual GW calculation. For this call 'spex.inv | tee spex.out'. This pipes the output to the file 'spex.out'. In that file you will find the output for the GW quasi particle energies per kpoint at the very end. Note that they have a real and an imaginary part which is related to the lifetime of the respective state.
We now want to compute a band structure in order to visualize the difference between the Kohn-Sham eigenenergies (fleur) and the GW energies (spex). One possible way to do so is to perform a GW calculation for each kpoint along a path in the Brillouin zone. For this we need to find coordinates along one path. The Spex feature 'KPTPATH' helps us with this. Uncomment the specific line in the 'spex.inp' file and run spex once again with 'spex.inv -w'. The new file 'qpts' includes all the desired kpoints. Now change your 'spex.inp' file to looks as follows:
BZ 4 4 4 KPT +=(0,0,0) X=1/2*(0,1,1) K=1/8*(3,6,3) #KPTPATH (K,1,X) 40 JOB GW +:(1-10) NBAND 80
The kpoint after the '+' has to be changed to one kpoint from the qpts list, followed by a 'spex.inv -w' run, by a fleur run with 'gw=2', and the Spex calculation. These steps do not need to be performed manually, there is a shell script prepared to do so. It can be found in the folder 'sh' of the spex directory. Please copy the script 'spex.band' to your current direcoty and add the path to the spex executable and to the fleur executable. Now you can run it. You have produced multiple output files from which one needs to extract the relevant values, for that there is a utility compiled called 'spex.extr', it can also be found in the source directory. Please execute it with 'spex.extr g -b -c spex_???.out > bandstr_gw' to get the GW bandstructure and with 'spex.extr k -b spex_???.out > bandstr_ks' to get the DFT band structure.
It is also possible to calculate the dielectric function. To calculate it for the KS system you need to change the 'JOB' line to 'JOB DIELEC 1:{0.0:0.2,0.0005}'. This means calculating the dielectric function from 0.0 Htr to 0.2 Htr with a step size of 0.001 Htr. This is only done for kpoint label '1', meaning for a ´displacement vector of (0,0,0). This is a very good assumption for optical transitions, since in this band picture photons carry a lot of energy with only very little momentum. When you now run spex again you get the files dielec and dielecR. The first one contains the head elements of the (bare) dielectric function and the second one the head element of the renormalized dielectric function . The imaginary part of the 'dielecR' can be identified with an optical absorption spectrum (excluding excitonic effects). We will be using dielecR.
Note: To obtain an EELS spectrum from these values one would have to plot for example in 'gnuplot': 'plot "dielecR" using 1:(-imag(1/(3))))'
Note2: Of course, it is also possible to use SPEX to calculate the dielectric function for the GW system.
Compare the gap extracted from 'dielecR' with the experimental band gap of Si. How large is the difference, why is it there?
2. Exercises
2.1. GaN
GaN crystallizes in Wurtzite (ZnS) structure. Its hexagonal unit cell ('hP' for inpgen) has the lattice parameters a=3.19 Angstrom and c=5.189 Angstrom. It features space group 186 (P63mc). The representative atoms are at (internal coordinates):
| atom species | x | y | z | Wyckoff symbol |
|---|---|---|---|---|
| Ga | 1/3 | 2/3 | 0.0 | 2b |
| N | 1/3 | 2/3 | 0.623 | 2b |
Set up the input for the material and calculate the DFT (PBE) band structure as well as the GW band structure along the path. For the GW calculation use an nband parameter of 200. Compare them. What are the band gaps in both cases? For the GW calculation also calculate the dielectric function and the complex refractive index. According to your calculations: What wavelength does a GaN laser diode have? Look up the experimental band gap of GaN and compare your values to it.
NOTE: M = (0.0,0.5,0.0) ; K = (1/3,1/3,0.0)