Calculations with SPEX
Note: This tutorial of the CECAM Fleur Workshop 2019 is based on version 05.00 of the Spex code.
Overview
 GW calculations
1.1 Executables and shell scripts
1.2 Silicon
1.3 Gallium nitride
1.4 Antimony telluride: Spinorbit coupling (SOC)
1.5 Sodium  Hubbard U parameters
2.1 Strontium vanadate (SrVO_{3})
2.2 Nickel
1. GW calculations
The GW selfenergy is an approximation to the electronic selfenergy used in manybody perturbation theory. This theory allows excited state properties to be calculated on the basis of a meanfield system, which is usually taken from a previous DFT calculation, e.g., using the localdensity approximation (LDA). One says that LDA is the starting point for the GW calculation.
Thus, in order to carry out a GW calculation, we first have to perform a selfconsistent DFT calculation with the FLEUR code. For the latter, you will use your compiled FLEUR program.
1.1 Installation
As a first step, we are going to install the SPEX code. Copy the file "Tutorial/DAY4/spex05.00.tgz" to your home directory and unpack it:
tar xvzf spex05.00.tgz
Enter the directory "spex05.00":
cd spex05.00
In order to configure the compilation for your system (choice and test of compiler, search for libraries, etc.), you have to run the configure script. The script allows for customization through command line arguments, for example:
./configure withdft=fleur withwan prefix=$HOME
In your Quantum Mobile system, the Wannier90 and HDF5 libraries are, however, not in standard directories. The location of the
Wannier90 library can be given after withwan
. The directories containing the library and include files of HDF5 are specified
using LDFLAGS=L...
and FCFLAGS=I...
:
./configure withdft=fleur withwan=$HOME/codes/wannier903.0.0/ prefix=$HOME LDFLAGS=L/usr/lib/x86_64linuxgnu/hdf5/serial/ FCFLAGS=I/usr/include/hdf5/serial/
A short explanation of the options:
withdft=fleur
: Compile the interface for the new MaX Release 3 of FLEUR (fleur v0.26b is the default),withwan
: Compile with Wannier support (searches for and links in the Wannier90 library),enablempi
: Compile parallel version using MPI 3.1 (not needed for the tutorial),prefix=PREF
: Install executables (and shell scripts) into PREF/bin (defaults to PREF=/usr/local; we choose PREF to be your home directory).
We only have a serial version of the HDF5 library. Therefore, we cannot compile an MPIparallelized version of SPEX.
Compile the source:
make
You have just compiled the executables:
spex.inv
: For systems with inversion symmetry (and without SOC)spex.noinv
: For systems without inversion symmetry (or with SOC)spex.extr
: Extraction utility
These executables can be installed by
make install
Together with the aforementioned executables, we have also installed the following shell scripts
spex
: SPEX launcher scriptspex.band
: For bandstructure calculationsspex.selfc
: For selfconsistent calculations (e.g., QSGW)
The reason for having two executables, one for systems with inversion symmetry (spex.inv
) and the other for systems without (spex.noinv
),
is that many quantities can be declared as realvalued instead of complexvalued arrays in the case of inversion symmetry,
which speeds up the computation considerably.
The executables (spex.inv
or spex.noinv
) can be called directly. Alternatively, the launcher script allows simply running spex
in both cases. You can choose one of the two options. In the following, we will use spex
for simplicity (as in the second option).
The SPEX manual is here. There is also a PDF version of the manual: "spex05.00/docs/spex.pdf".
1.2 Silicon
1.2.1 DFT as the meanfield starting point
Your first GW calculation will be on bulk Silicon. As a first step, we must perform a DFT calculation that creates the necessary input for GW.
Silicon crystallizes in a diamond lattice (fcc lattice with a twoatom basis) with a lattice constant of 5.43 Å (10.26 Bohr).
You can create your own FLEUR input file "inp_Si" as you have already learned it or copy it from the directory "Tutorial/DAY4/inp_Si".
We must tell FLEUR to create some additional files (containing structural information etc.) for the subsequent GW calculation.
This is achieved by running inpgen
with the option gw
:
inpgen gw < inp_Si
(When you investigate the created file inp.xml
, you will see that the option gw
has two effects:
it sets the flag gw
to "1" and declares an optional kpoint file named "kpts_gw".)
Then, perform a DFT calculation for Silicon as usual, which produces a selfconsistent density and effective potential. Make sure that the calculation has converged.
1.2.2 Special kpoint set
As the GW selfenergy is nonlocal, its evaluation involves summations over two kpoints, k and q, and their sum k+q. Thus, we must make sure that the set of kpoints forms an equidistant regular mesh, which contains k+q as an element if k and q are elements of the set. Furthermore, the Γ point (k=0) must be included. (In other words, the kpoint set must form a group under vector addition.) The explicit inclusion of the Γ point is necessary because of the divergence of the Coulomb interaction in the longwavelength limit.
We can now proceed in two ways:
 SPEX creates the new kpoint set (Section 1.2.3). FLEUR runs a oneshot calculation on that new kpoint set (Section 1.2.4). Finally run the GW calculation with SPEX (Section 1.2.5).
 Directly run the GW calculation with SPEX (Section 1.2.5) with an additional line in the input file "spex.inp":
ITERATE SR
.
Obviously, the second procedure is much simpler as it spares us two runs (with SPEX and FLEUR). The first procedure is more flexible, though. It is required, for example, for selfconsistent calculations (QSGW etc.). We leave it up to you, which procedure you choose. You might want to try both.
1.2.3 Create new kpoint set
The new kpoint set can be created in the following way: Create an input file for SPEX, named "spex.inp", with only one entry:
BZ 4 4 4
The entry specifies a 4×4×4 kpoint set. Then, run SPEX with the argument w
spex w
This creates the new kpoint file "kpts_gw".
1.2.4 Run FLEUR on new kpoint set
Now we let FLEUR generate the wavefunctions and energies on this new kpoint set.
Important: This should be a oneshot calculation without selfconsistency. For this purpose, we must modify the FLEUR "inp.xml" file in the following way:

Set the
gw
flag to "2":
<expertModes gw="2" secvar="F"/>

Set
numbands
to "180":
<cutoffs Kmax="3.60000000" Gmax="11.00000000" GmaxXC="9.20000000" numbands="180"/>
Setting gw="2"
lets FLEUR create additional files that are necessary for the GW calculation.
Setting numbands="180"
sets the number of bands (here, 180) returned by the diagonalization routine of FLEUR. This number must comprise occupied states
and a number of unoccupied states, as well, because the GW scheme involves sums over unoccupied states. (The default numbands="0"
would let FLEUR
generate only the occupied states.)
After these changes have been applied to "inp.xml", you can run FLEUR again. This creates the necessary input data for SPEX. Note that you do not have to repeat the
steps described in Section 1.2.3 and 1.2.4 for each SPEX run unless you change the kpoint set (BZ
).
1.2.5 Perform the GW calculation
We have already created the SPEX input file "spex.inp" above. For the GW calculation, we must add a few parameters now. Change your "spex.inp" to look like this:
BZ 4 4 4 KPT X=1/2*(0,1,1) L=1/2*(0,0,1) JOB GW 1:(1,2,5) X:(1,3,5) L:(13,5) NBAND 80
Note: If you have jumped here from Section 1.2.2 (without performing the oneshot FLEUR calculation), you have to add the line
ITERATE SR
to the input file. "SR" stands for "scalarrelativistic".
The second line defines some kpoint labels, here for the X and L point. These labels are used in the job definition (third line):
GW
stands for the GW approximation; then the states (on which the quasiparticle correction is to be applied) are defined
with the syntax "kpoint:(band indices)". The kpoint can be specified by its index (see list of kpoints in the output, e.g, "1" for the Γ point)
or by a kpoint label.
We have chosen to define the band indices in the job definition in such a way that, in the case of degeneracy, only one state out of a group of degenerate states is
calculated. For example, the second, third, and fourth state at Γ (kpoint index 1) are degenerate. You can also include these states with 1:(15)
.
The keyword NBAND
defines the number of bands (occupied and unoccupied) used in the GW calculation. Obviously, this number must not exceed the available number
of bands, which corresponds either to the numbands
parameter of "inp.xml" (180) or, if you use ITERATE
, to the number of LAPW basis functions.
If you omit the NBAND
definition, SPEX takes into account all available states.
Now all necessary files are created, and we can perform the GW calculation with SPEX. When SPEX is run, it writes information to the two standard streams of the shell: standard output (stdout) and standard error (stderr). The first stream is used for the normal output, the second for warning and error messages. Normally, both streams are directed to the screen, but they can also be directed to a file using suitable shell commands. Here are some examples:
spex
: Both output and error streams are written to the screen (not recommended).spex > spex.out
: The output stream is written to the file "spex.out", which you can later open with your favourite editor. The error stream is directed to the screen.spex &> spex.out
: Both streams are written to the file "spex.out".spex  tee spex.out
: The output stream is written to the screen but also to the file "spex.out" (recommended). This allows you to monitor the calculation while it is running. The error stream is written to the screen.spex 2> spex.err  tee spex.out
: The output stream is written to the terminal, but also to the file "spex.out". The error stream is written to the file "spex.err".
You should increase the horizontal size of your terminal a little bit because some lines of the output exceed the standard 80 characters of the terminal.
1.2.6 A longer input file
Before we investigate the output file, let us have a look at the input file again. The one used in the previous section was a minimal one. A typical GW calculation requires a quite large number of parameters. When these parameters are not given in the input file (as above), SPEX uses default parameters or determines reasonable ones from the DFT data. Of course, all parameters can also be defined explicitly in the input file. A longer input file for the above calculation could look like this:
# silicon bulk BZ 4 4 4 KPT X=1/2*(0,1,1) L=1/2*(0,0,1) JOB GW 1:(1,2,5) X:(1,3,5) L:(13,5) NBAND 80 SECTION MBASIS LCUT 4 GCUT 3.6 SELECT 2;3 OPTIMIZE MB 4.0 END SECTION SUSCEP HILBERT 50 30 END SECTION SENERGY MESH 10+2 10.0 CONTINUE END
Note: The empty lines and the twocharacter indentation in the sections are used only for clarity. Everything after #
is interpreted as a comment.
Some explanations:
The section MBASIS
contains the definition of the mixed product basis. The keyword HILBERT
defines the realfrequency mesh for the Hilbert transformation
of the susceptibility (polarization function).
The section SENERGY
defines parameters for the evaluation of the selfenergy. The keyword MESH
defines the imaginaryfrequency mesh.
The keyword CONTINUE
tells SPEX to use analytic continuation (PadÃ© extrapolation) of the selfenergy from the imaginary to the realfrequency axis.
(This is the default.)
1.2.7 Investigating the results in "spex.out"
Have a look at the output file "spex.out". After reading and checking all parameters (routines "getinput" and "checkinput") and setting up the mixed product basis (routine "mixedbasis"), SPEX calculates the Coulomb matrix (routines "structureconstant" and "coulombmatrix"). Then, SPEX evaluates the exchange selfenergy (routines "exchange_core" and "exchange"), which corresponds to a oneshot HartreeFock (HF) calculation, and the HF energies are printed. Then, the correlation selfenergy is constructed in a kpoint loop, which consumes most of the CPU time. For each kpoint the susceptibility (or polarization function ) is calculated (routine "susceptibility"), from which the screened interaction is obtained, and the correlation selfenergy is updated. Finally, the selfenergy is analytically continued to the realfrequency axis and the quasiparticle equation is solved. The results are written in a table:
Bd vxc sigmax sigmac Z KS HF GW lin/dir 1 12.11172 19.10561 6.57060 0.64929 6.20997 13.20387 6.37436 0.72185 1.04269 0.10593 6.37825 0.71237 2 13.21577 14.83300 1.14256 0.76185 5.60804 3.99081 5.24641 0.00147 0.00088 0.00169 5.24744 0.00162 5 11.42949 6.88796 4.16444 0.75863 8.16243 12.70396 8.44847 0.00690 0.00667 0.00489 8.44752 0.00710
The columns are
Bd
: the band index,vxc
: expectation value of , e.g., 12.11172 eV,sigmax
: expectation value of the exchange selfenergy , e.g., 19.10561 eV,sigmac
: expectation value of the correlation selfenergy , e.g., (6.57060+1.04269i) eV (note that the selfenergy is complex),Z
: renormalization factor Z, e.g., (0.649290.10593i) eV,KS
: KohnSham value (from the previous DFT calculation), e.g., 6.20997 eV,HF
: HartreeFock value (oneshot), e.g., 13.20387 eV,GW
: GW quasiparticle values, e.g., (6.37436+0.72185i) eV (solution of a linearized, i.e., approximated, quasiparticle equation) and (6.37852+0.71237i) eV (exact solution of the nonlinear quasiparticle equation). Usually, there is only little difference between the two values.
Look up the KohnSham, HF, and GW direct and indirect gaps (i.e., Γ→Γ and Γ→X) and compare with the experimental values,
3.4 eV and 1.17 eV, respectively. Note that the conduction band minimum is not exactly at the X point, though (see below).
The fourth band is the highest valence band: There are two silicon atoms per unit cell. Each contributes four electrons. So, there are eight electrons in the
unit cell. (You can find this line in the output: Number of valence electrons: 8
.)
Then, dividing by two (two electron spins) gives four.
1.2.8 Contour integration
Instead of analytic continuation, you can use contour integration for evaluating the selfenergy. This is more accurate, but takes more CPU time.
Simply define CONTOUR 0.01
in the file "spex.inp" inside a section SENERGY
:
SECTION SENERGY CONTOUR 0.01 END
(or replace CONTINUE
by CONTOUR 0.01
if you are using the input file of Section 1.2.6). Then rerun the calculation.
(You might want to use a different output file.) The argument "0.01" tells SPEX to calculate the selfenergy at two energies symmetrically around
each KS energy ε (ε0.01 htr and ε+0.01 htr), which enables a linearization of the selfenergy.
Therefore, we only get the "linearized" quasiparticle result.
(We can also obtain both solutions again by specifying a range of values instead of the single argument, e.g., CONTOUR +{0.05:0.05,0.01}
.)
Compare the GW results with the previous calculation.
1.2.9 Band convergence
You can play around with the parameters, e.g., increase the number of bands (keyword NBAND
) to 100 and 120 and check the bandgap convergence.
You can also experiment with other parameters, e.g., those of the mixed product basis.
1.2.10 Denser kpoint sets
If you have time, you can try denser kpoint sets, e.g., 6×6×6.
1.2.11 Band structure calculation
The calculation of a GW band structure is not as straightforward as in the case of DFT, because we cannot arbitrarily add kpoints to the kpoint set.
On the other hand, if we restrict ourselves to the set of kpoints defined by the keyword BZ
, the band structure will have very low k resolution.
Anyway, this will be our first calculation.
1.2.11.1 Keyword PATH
Let us make a band structure calculation for the path from L over Γ to X. In principle, we would have to identify all kpoint indices along this path and define separate jobs in the job definition, but SPEX can do this for you. To define the kpoint path, add the line
KPTPATH (L,1,X)
to "spex.inp". The GW quasiparticle energies are to be calculated on each kpoint along this path.
To this end, use the keyword PATH
in the job definition, replacing the previous definition by
JOB GW PATH:(18)
(To get a somewhat finer band structure, you can also change BZ 4 4 4
to BZ 6 6 6
. If you do this change, you would have to run the oneshot FLEUR
calculation again as explained in the Sections 1.2.3 and 1.2.4, unless you use ITERATE
.)
When you now run SPEX, it will automatically choose the kpoints that lie on the path defined by KPTPATH
and calculate the respective quasiparticle energies.
Afterwards, you can extract the band structure data with the utility "spex.extr"
spex.extr g b spex.out > bandstr
(assuming that the output has been piped into the file "spex.out") and use "gnuplot" or "xmgrace" to plot the data (example Gnuplot: plot "bandstr" with lp
).
As expected, the band structure is not smooth. One solution would be to use denser kpoint sets, but this would entail very expensive calculations.
Note: The absolute quasiparticle energies (y axis) are plotted. The Fermi energy is not subtracted. The valenceband maximum is at around 5.2 eV.
1.2.11.2 Additional kpoints (KPT +...
)
Another possibility is to use a feature that allows us to add arbitrary kpoints, one at a time. Each additional kpoint q must be accompanied with all shifted kpoints q+k. So, for each single kpoint q we must run FLEUR to generate wavefunctions and energies at the shifted kpoint set and then run SPEX. Let us do this for a single kpoint first.
(In principle, running FLEUR can again be avoided by using ITERATE
. However, for simplicity, we will from now on assume that FLEUR is used for the
oneshot calculation.)
We define the conductionband minimum by
KPT +=1/2*(0,0.75,0.75)
The label +
is a special kpoint label for additional kpoints. Then, create the kpoint file (with spex w
) and run FLEUR as explained in Sections 1.2.3 and 1.2.4.
Before running SPEX, the job definition line should be replaced by
JOB GW 1:(15) +:(15)
The bands 15 are calculated for the Γ point and the additional kpoint. After the calculation is done, you can read the fundamental band gap from the output file. It should be smaller by about 0.1 eV than the result for the Γ→X transition calculated above.
Let us now make a full band structure calculation for the same path as above. When you investigate your directory, you should see a file "qpts" that
has been created when we last ran spex w
. This file contains a list of kpoints along the specified path. For each of the kpoints, we would have to run,
in this order, SPEX, FLEUR, and SPEX again. To simplify this task, there is a shell script called "spex.band". This shell script performs all the steps
automatically for you. The script automatically adds the RESTART
option to "spex.inp" (which tells SPEX to calculate the screened interaction only once and to
store it on disc). Specify eight bands with
JOB GW +:(18)
To run the shell script simply type spex.band
. For each kpoint, the script creates a SPEX output file "spex_NUM.out", where NUM is a counter index.
Note: If an output file (e.g., "spex_005.out") already exists in the current directory, the script "spex.band" skips the calculation of the respective kpoint and continues with the next. Therefore, do not forget to remove the files "spex_???.out" when you run "spex.band" again in a later exercise!
From the output files, the band structure data is extracted with
spex.extr g b spex_???.out > bandstr
which writes the data to the file "bandstr". The band structure can then be plotted with "xmgrace" or "gnuplot". Note that the conduction band minimum is, in fact, between the Γ and X point at about k=1/2*(0,0.75,0.75).
If you extract the bandstructure data, instead, with
spex.extr g b c spex_???.out > bandstr
you get an additional column, which contains the imaginary parts of the quasiparticle energies. Remember that these imaginary parts are proportional to the line width and inversely proportional to the excitation lifetime. The line widths can be plotted together with the band structure in gnuplot with
plot "bandstr" using 1:2:(6*abs($3)) with points ps var
You will see that the bands get wider the further they are away from the Fermi energy. (We will come back to this point in Section 1.2.12.)
1.2.11.3 Wannier interpolation
As a third alternative, we can use Wannier interpolation. We have to modify "spex.inp" in the following way
BZ 4 4 4 JOB GW IBZ:(14) KPT X=1/2*(0,1,1) L=1/2*(0,0,1) KPTPATH (L,1,X) SECTION WANNIER ORBITALS 1 4 (sp3) MAXIMIZE INTERPOL END
The job definition lets SPEX calculate quasiparticle energies in the whole irreducible Brillouin zone (IBZ).
The section WANNIER
defines a set of maximally localized Wannier functions (MLWFs) of bonding sp^{3} hybrid orbitals.
These four hybrid orbitals are generated from the four lowest valence bands (first two arguments after ORBITALS
).
Running SPEX will construct the MLWFs and use them to interpolate the band structure.
The band structure data is written to the files "bands0" for the KohnSham energies and "bands1" for the GW quasiparticle energies.
You can use "xmgrace" or "gnuplot" for plotting. Please compare the interpolated band structure stored in "bands1" and the one that we have calculated
explicitly before (file "bandstr"). Note that the energies in the file "bands1" are given in atomic units (1 hartree=27.21138 eV) and must be scaled accordingly.
(In Gnuplot: plot "bands1" u 1:($2*27.21138)
)
When you open the files with an editor, you will see that "bands1" has a column more than "bands0".
This additional column contains the interpolated imaginary parts of the quasiparticle energies.
You can plot the band structure together with the variable point size as described in the previous section.
Since we have defined four Wannier functions, we only get the four valence bands of silicon and no conduction bands.
We will need Wannier functions again when we calculate Hubbard U parameters in Section 2.
1.2.12 Spectral functions
The central quantity in the GW method is the renormalized Green function G, which is related to the selfenergy Σ through the Dyson equation . Every GW program calculates Σ and, then, solves the Dyson equation. The last step is often simplified in such a way that an effective singleparticle equation (the socalled quasiparticle equation) is solved instead, yielding the quasiparticle energies and lifetimes instead of the Green function G, as we have done in this tutorial up to now. One can say that the quasiparticle energies and lifetimes parameterize G. This works well for the quasiparticle peaks but not for manybody states such as plasmon satellites. SPEX can also be used to calculate G directly. The imaginary part of its trace is called the spectral function. To this end, we modify "spex.inp" in the following way:
SECTION SENERGY CONTOUR {2:2,0.01} SPECTRAL END
Here, we recommend to use the contour deformation method (CONTOUR
).
If you now run a GW calculation, you will get a file "spectral" containing spectral functions for each kpoint defined in the job definition. For example,
with
JOB GW +:(18)
you would get the spectral function at the additional kpoint labeled with + (presumably still the kpoint of the conductionband minimum). Of course, you can choose any kpoint (also several kpoints). Have a look at the spectral function (with "gnuplot" or "xmgrace"). The x axis is the energy axis (in Hartree; 1 Hartree = 27.21138 eV). One can imagine the spectral function as a vertical cut through the band structure at the respective kpoint. If you compare with the band structure of silicon, you will see that there is a (quasiparticle) peak in the spectral function for each band at the respective kpoint. The peaks broaden the further they are away from the Fermi energy. One can also observe a broad peak at around 30 eV. This is the plasmon satellite. You should see another broad feature above the Fermi energy, at 2030 eV, which is the plasmon satellite for the inverse photoemission process.
It is also possible to calculate the spectral functions along a path in the Brillouin zone,
to obtain a renormalized band structure, in a similar way as above.
In particular, we would follow the steps of Section 1.2.11.2 or Section 1.2.11.3 with the keyword SPECTRAL
.
(With the procedure of Section 1.2.11.1 one would not get sufficient k resolution.)
In the first case (Section 1.2.11.2), the shell script "spex.band" produces a whole series of files containing spectral functions: "spectral_001", "spectral_002", etc. The corresponding data can be extracted with
spex.extr s spectral_??? > spec
You can use "gnuplot" to plot the resulting renormalized band structure as a colorcoded plot by
set cbrange [0:100] plot "spec" using 1:2:3 palette with lines
Of course, the k resolution is not very high. One can produce a finer kmesh by KPTPATH (L,1,X) 50
(and then running
spex w
to produce the kpoint file "qpts").
The additional parameter sets the number of kpoints along a unit length. The default is 20.
In the second case (Section 1.2.11.3), SPEX immediately produces a spectral file, named "spectralw", which can be plotted in the same way as "spec". Check if you can see the plasmon satellites in the plots!
1.3 Gallium nitride
GaN crystallizes in the Wurtzite (ZnS) structure. Its hexagonal unit cell ("hP" in the input for "inpgen") has the lattice parameters a=3.19 Å and c=5.189 Å with the 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 (or copy it from the directory "Tutorial/DAY4/inp_GaN"), calculate the GW band gap.
It is at the Γ point. You also need to know the band index of the valence band maximum.
Around line 90 of the SPEX output, the number of electrons in the unit cell are given (Number of valence electrons: 36
).
Dividing this number by 2 (for the two spin directions) yields the band index of the valenceband maximum.
So, the band gap is between band 18 and band 19. We are also interested in the 3d levels of Ga. These bands are the first
ten bands (each Ga atom contributes five). Just add the first band to the job definition.
GaN has twice as many atoms as silicon. So, we can reduce the BZ sampling:
BZ 3 3 3
Compare the GW band gap to the DFT one and to experiment (3.39 eV). What wavelength does a GaN laser diode have?
Run a GW calculation for the Ga 3d semicore states (one of the bands 110) and compare the resulting 3d levels with the KohnSham and HF values as well as with experiment (17.1 eV). Note that the semicore level is measured with respect to the valenceband maximum (Band 18 at Γ).
It is also possible to calculate the dielectric function. For this, change the job definition to
JOB DIELEC 1:{0:30eV,0.1eV}
Here, we request a plot of the dielectric function between 0 and 5 eV with a step size of 0.01 eV. (Generally, energy values without a unit are interpreted as values in "Hartree". Alternatively, the values can be given in eV as shown.) This is done for kpoint label "1", meaning for a displacement vector of (0,0,0), corresponding to optical transitions. (The momentum of the photon is much smaller than that of the electrons. The longwavelength limit k=0 is, thus, a good approximation). A SPEX run will produce the files "dielec" and "dielecR", corresponding to the bare and the renormalized dielectric function. The optical absorption spectrum (excluding excitonic effects) can be plotted as a function of frequency ω with ("gnuplot")
plot "dielecR" using 1:3 with lines
For an electronenergy loss spectrum (EELS), we have to plot the imaginary part of the reciprocal value:
plot "dielecR" using 1:(imag(1/($2+{0,1}*$3))) with lines
The literature value of the dielectric constant (ω=0) is about 5.8. The refractive index n is the square root of the dielectric constant. The literature value is n=2.41.
1.4 Antimony telluride: Spinorbit coupling (SOC)
Let us now consider a more advanced system, in which the spinorbit coupling (SOC) plays an important role: Sb_{2}Te_{3}. The input files can be copied from here: "Tutorial/DAY4/Sb2Te3/". These input files are especially prepared for this tutorial (to make a G^{SOC}W^{SOC} calculation possible on a slow system). You do not need to run "inpgen" this time.
Perform a GW calculation without SOC (directory "withoutSOC") as you did for silicon and GaN. Calculate the quasiparticle energies for the states 12 to 15 at the Γ point. The corresponding file "spex.inp" is already prepared. Calculate the value of the direct band gap, Γ→Γ, and compare the KohnSham and GW values. You can determine the band index for the valenceband maximum according to the explanation in the previous section.
Is the GW value larger or smaller than the KS one? The experimental value is under debate (170280 meV).
Now run the GW calculation with SOC (directory "withSOC"). Instead of doing the oneshot calculation with FLEUR (Section 1.2.4), we use
ITERATE FR
in "spex.inp". Here, "FR" stands for "fully relativistic". So, you do not have to rerun FLEUR.
Remember that now the states are doubly degenerate (Kramer's degeneracy). We therefore specify the bands from 23 to 30 in the job definition. Note that the index for the valenceband maximum is doubled as well. (Spin is not a good quantum number anymore.) Calculate again the direct GW band gap and compare with the KohnSham value and with the respective value from the previous calculation (without SOC), as well as with the experiment. Is the GW value now larger or smaller than the KohnSham one? Why?
Look at the states from 23 to 28 (the 6 highest valence bands). Has the spinorbit coupling caused a splitting of states that were otherwise degenerate?
1.5 Sodium
If you still have time, you can make a GW calculation for metallic Na.
Try to calculate a GW band structure from P over Γ to H as explained in Section 1.2.11.2. You can use the kpoint definition
KPT P=1/4*(1,1,1) H=1/2*(1,1,1)
Sodium has four semicore states: 2s and 2p. So, the band indices in the job definition should start with 5. For extracting the band structure data use
spex.extr g b O spex_???.out > bandstr
Compare this band structure to the corresponding HartreeFock band structure (h
instead of g
after spex.extr
).
While the latter is discontinuous at the Fermi energy
(the logarithmic divergence of the effective mass at the Fermi energy cannot be resolved numerically; the calculation yields a discontinuity, instead),
there is no discontinuity in the GW band structure because the correlation selfenergy cancels the discontinuity of the HF term exactly by a term of opposite sign.
Note that the band width is reduced in the GW band structure compared to the KohnSham band structure. (Can be extracted with k
after spex.extr
.)
This reduction is caused by the increased effective mass of the quasiparticles with respect to the bare electrons.
As the electron moves, it has to drag its polarization cloud (Coulomb hole) along, which makes the quasiparticle effectively heavier.
2 Hubbard U parameters
2.1 Strontium vanadate (SrVO_{3})
As a first example, we will calculate the Hubbard U (partially screened Coulomb interaction) parameter for the V t_{2g} orbitals of SrVO_{3} using the cRPA method. You find an input file ("inp_SrVO3") for "inpgen" in the directory "Tutorial/DAY4".
Before we calculate the Hubbard U parameter we want to have a look at the t_{2g} bands, which will form the correlated subspace.
The corresponding band indices are 2123. A basis set for the correlated subspace will be formed by Wannier orbitals.
First, we calculate a band structure with FLEUR: set band="T"
in "inp.xml", then run FLEUR. (You should not forget to set it back to "F"
afterwards.)
The file containing the bandstructure data is called "band.1". Second, we calculate only the bands forming the t_{2g} subspace using
Wannier interpolation, as described in 1.2.11.3.
We want to use the same kpoint path as in FLEUR. Please, look up the special kpoint positions and labels in the file "out.xml" created by FLEUR
and add them to the file "spex.inp" using kpoint labels (KPT
) as described in Section 1.2.5.
You also have to define the kpoint path with KPTPATH
accordingly.
You can use the following input for "spex.inp" but add the KPT
and KPTPATH
lines.
BZ 4 4 4 JOB SECTION WANNIER ORBITALS 21 23 () (t2g) INTERPOL END
The Wannierinterpolated band structure is written to the file "bands0".
You can compare the two band structures in "gnuplot". Make sure you convert Hartree to eV and substract the Fermi energy for the SPEX band structure:
plot "bands.1","bands0" u 1:(($20.354164)*27.21138)
.
A minimal SPEX input file for Hubbard U calculations looks like this:
BZ 3 3 3 JOB SCREENW {0} SECTION SUSCEP HUBBARD END SECTION WANNIER ORBITALS 21 23 () (t2g) END
We have reduced the Brillouinzone sampling to 3×3×3 to speed up the calculation.
The keyword SCREENW
defines the job for calculating the partially screened Coulomb potential U. The expression in the brackets {...} can be
a general frequency range (in the same form as for the dielectric function, Section 1.3). Here, {0} means to calculate the Hubbard U parameter for
the static limit (ω=0). We only need the static U for model Hamiltonians or LDA+U calculations.
In the section SUSCEP
we provide the keyword HUBBARD
to exclude transitions in the correlated subspace, otherwise SPEX calculates
the fully screened Coulomb interaction W.
Wannier functions are constructed internally in the SPEX code.
For this purpose, we have to provide the keyword ORBITALS
and specify the states from which the Wannier functions are to be constructed.
For SrVO_{3}, we construct three V t_{2g} orbitals from the states 17, 18, and 19. The V atom is the second type as specified in the FLEUR input.
Therefore, we use an "empty" definition ()
for the first atom (Sr).
In this calculation, we simply use the firstguess Wannier functions, thus the keyword MAXIMIZE
is excluded.
When SPEX is run, it will calculate the Coulomb matrix U in the mixed product basis and then projects it onto the Wannier basis. You will get some effective parameters at the end of the program run, Kanamori parameters for e_{g} and t_{2g} orbitals (and HubbardHund parameters if full shells are specified). The full interaction matrix is written to the file "spex.cou".
You can run the same calculation without the keyword HUBBARD
. Do the interaction parameters differ from the previous ones? If so, why?
Let us now also include the e_{g} orbitals in the Wannier set:
SECTION WANNIER ORBITALS 21 27 () (t2g,eg) MAXIMIZE END
We have to run the Wannier90 maximal localization procedure in this case. Calculate the Hubbard U parameters and run another calculation with
SECTION WANNIER ORBITALS 21 27 () (t2g) MAXIMIZE END
Compare the Hubbard U parameter for the t_{2g} states! How much do they differ? What do you conclude from the difference?
2.2 Nickel
If you still have time, you can calculate the Hubbard U parameter for Nickel. To be more precise, we want to calculate the U parameter for the 3d states of Ni.
We can mainly proceed as in the previous section. The Brillouinzone sampling (BZ
) should be a bit finer:
BZ 6 6 6
And we need to define the Wannier orbitals (ORBITALS
) differently:
SECTION WANNIER ORBITALS 4 9 (d) END
Why do the band indices start with 4 in this case? (Have a look at the file "inp.xml"!)
As bulk nickel is ferromagnetic and Wannier functions have a formal spin dependence (since they are constructed from wavefunctions), we have to specify the spin indices for the Wannier products explicitly in the job definition. (Usually, the spin dependence is numerically small.)
JOB SCREENW uu{0}
Since we have defined a full shell, the 3d shell, we get HubbardHund parameters (instead of Kanamori parameters).