How to obtain the quasiparticle band structure of a bulk material: hBN
In this tutorial you will learn how to:
 Calculate quasiparticle corrections in the HartreeFock approximation
 Calculate quasiparticle corrections in the GW approximation
 Choose the input parameters for a meaningful converged calculation
 Plot a band structure including quasiparticle corrections
We will use bulk hBN as an example system. Before starting, you need to obtain the appropriate tarball. See instructions on the main tutorials page.
We strongly recommend that you first complete the First steps: a walk through from DFT to optical properties tutorial.
The aim of the present tutorial is to obtain quasiparticle correction to energy levels using manybody perturbation theory (MBPT).
The general nonlinear quasiparticle equation reads:
As a first step we want to evaluate the self energy Σ entering in the quasiparticle equation. In the GW approach the selfenergy can be separated into two components: a static term called the exchange selfenergy (Σx), and a dynamical term (energy dependent) called the correlation selfenergy (Σc):
We will treat these two terms separately and demonstrate how to set the most important variables for calculating each term. In practice we will compute the quasiparticle corrections to the one particle KohnSham eigenvalues obtained through a DFT calculation.
The steps are the following:
Contents
 1 Step 1: The Exchange Self Energy or HF quasiparticle correction
 2 Step 2: The Correlation SelfEnergy and Quasiparticle Energies
 3 Step 3: Interpolating Band Structures
 4 Step 4: Summary of the convergence parameters
 5 Step 5: Eigenvalue only selfconsistent GW0
 6 Step 6: Exercise: Convergence with respect K points
 7 References
Step 1: The Exchange Self Energy or HF quasiparticle correction
We start by evaluating the exchange SelfEnergy and the corresponding Quasiparticle energies (HartreeFock energies). Follow the module on Hartree Fock and then return to this tutorial How to obtain the quasiparticle band structure of a bulk material: hBN.
Step 2: The Correlation SelfEnergy and Quasiparticle Energies
Once we have calculated the exchange part, we next turn our attention to the more demanding dynamical part. The correlation part of the selfenergy in a plane wave representation reads:
In the expression for the correlation self energy, we have (1) a summation over bands, (2) an integral over the Brillouin Zone, and (3) a sum over the G vectors. In contrast with the case of Σx, the summation over bands extends over all bands (including the unoccupied ones), and so convergence tests are needed. Another important difference is that the Coulomb interaction is now screened so a fundamental ingredient is the evaluation of the dynamical dielectric matrix. The expression for the dielectric matrix, calculated at the RPA level and including local field effects, has been already treated in the Local fields tutorial.
In the following, we will see two ways to take into account the dynamical effects. First, we will see how to set the proper parameters to obtain a model dielectric function based on a widely used approximation, which models the energy dependence of each component of the dielectric matrix with a single pole function. Secondly, we will see how to perform calculations by evaluating the dielectric matrix on a regular grid of frequencies.
Once the correlation part of the selfenergy is calculated, we will check the convergence of the different parameters with respect to some final quantity, such as the gap.
After computing the frequency dependent selfenergy, we will discover that in order to solve the quasiparticle equation we will need to know its value at the value of the quasiparticle itself. In the following, unless explicitly stated, we will solve the nonlinear quasiparticle equation at first order, by expanding the selfenergy around the KohnSham eigenvalue. In this way the quasiparticle equation reads:
where the normalization factor Z is defined as:
The Plasmon Pole approximation
As stated above, the basic idea of the plasmonpole approximation is to approximate the frequency dependence of the dielectric matrix with a single pole function of the form:
The two parameters R_{GG'} and Ω_{GG'} are obtained by a fit (for each component), after having calculated the RPA dielectric matrix at two given frequencies. Yambo calculates the dielectric matrix in the static limit ( ω=0) and at a user defined frequency called the plasmonpole frequency (ω=iωp). Such an approximation has the big computational advantage of calculating the dielectric matrix for only two frequencies and leads to an analytical expression for the frequency integral of the correlation selfenergy.
Input file generation
Let's start by building up the input file for a GW/PPA calculation, including the calculation of the exchange selfenergy. From yambo H
you should understand that the correct option is yambo x p p g n
:
$ cd YAMBO_TUTORIALS/hBN/YAMBO $ yambo x p p g n F gw_ppa.in
Let's modify the input file in the following way:
HF_and_locXC # [R XX] HartreeFock Selfenergy and Vxc ppa # [R Xp] Plasmon Pole Approximation gw0 # [R GW] GoWo Quasiparticle energy levels em1d # [R Xd] Dynamical Inverse Dielectric Matrix EXXRLvcs = 40 Ry # [XX] Exchange RL components Chimod= "Hartree" # [X] IP/Hartree/ALDA/LRC/BSfxc % BndsRnXp 1  10  # [Xp] Polarization function bands % NGsBlkXp= 1 RY # [Xp] Response block size % LongDrXp 1.000000  1.000000  1.000000  # [Xp] [cc] Electric Field % PPAPntXp = 27.21138 eV # [Xp] PPA imaginary energy %GbndRnge 1  40  # [GW] G[W] bands range % GDamping= 0.10000 eV # [GW] G[W] damping dScStep= 0.10000 eV # [GW] Energy step to evaluate Z factors DysSolver= "n" # [GW] Dyson Equation solver ("n","s","g") %QPkrange # [GW] QP generalized Kpoint/Band indices 7 7 8 9 %
Brief explanation of some settings:
 Similar to the Hartree Fock study, we will concentrate on the convergence of the direct gap of the system. Hence we select the last occupied (8) and first unoccupied (9) bands for kpoint number 7 in the
QPkrange
variable.  We also keep
EXXRLvcs
at its converged value of 40 Ry as obtained in the Hartree Fock tutorial.  For the moment we keep fixed the plasmon pole energy
PPAPntXp
at its default value (=1 Hartree).  We keep fixed the direction of the electric field for the evaluation of the dielectric matrix to a nonspecific value:
LongDrXp
=(1,1,1).  Later we will study convergence with respect to
GbndRnge
,NGsBlkXp
, andBndsRnXp
; for now just set them to the values indicated.
Understanding the output
Let's look at the typical Yambo output. Run Yambo with an appropriate J
flag:
$ yambo F gw_ppa.in J 10b_1Ry
In the standard output you can recognise the different steps of the calculations: calculation of the screening matrix (evaluation of the non interacting and interacting response), calculation of the exchange selfenergy, and finally the calculation of the correlation selfenergy and quasiparticle energies. Moreover information on memory usage and execution time are reported:
... <> [05] Dynamic Dielectric Matrix (PPA) ... <08s> Xo@q[3] ######################################## [100%] 03s(E) 03s(X) <08s> X@q[3] ######################################## [100%] (E) (X) ... <43s> [06] Bare local and nonlocal ExchangeCorrelation <43s> EXS ######################################## [100%] (E) (X) ... <43s> [07] Dyson equation: Newton solver <43s> [07.01] G0W0 (W PPA) ... <45s> G0W0 (W PPA) ######################################## [100%] (E) (X) ... <45s> [07.02] QP properties and I/O <45s> [08] Game Over & Game summary
Let's have a look at the report and output. The output file o10b_1Ry.qp contains (for each band and kpoint that we indicated in the input file) the values of the bare KS eigenvalue, its GW correction and the correlation part of the self energy:
# # Kpoint Band Eo EEo ScEo # 7.000000 8.000000 0.411876 0.563191 2.322443 7.000000 9.000000 3.877976 2.413822 2.232241
In the header you can see the details of the calculations, for instance it reports that NGsBlkXp
=1 Ry corresponds to 5 Gvectors:
# X G`s [used]: 5
Other information can be found in the report file r10b_1Ry_em1d_ppa_HF_and_locXC_gw0, such as the renormalization factor defined above, the value of the exchange selfenergy (nonlocal XC) and of the DFT exchangecorrelation potential (local XC):
[07.02] QP properties and I/O ============================= Legend (energies in eV):  B : Band  Eo : bare energy  E : QP energy  Z : Renormalization factor  So : Sc(Eo)  S : Sc(E)  dSp: Sc derivative precision  lXC: Starting Local XC (DFT) nlXC: Starting nonLocal XC (HF) QP [eV] @ K [7] (iku): 0.000000 0.500000 0.000000 B=8 Eo= 0.41 E= 0.98 EEo= 0.56 Re(Z)=0.81 Im(Z)=.2368E2 nlXC=19.13 lXC=16.11 So= 2.322 B=9 Eo= 3.88 E= 6.29 EEo= 2.41 Re(Z)=0.83 Im(Z)=.2016E2 nlXC=5.536 lXC=10.67 So=2.232
Extended information can be also found in the output activating in the input the ExtendOut
flag. This is an optional flag that can be activated by adding the V qp
verbosity option when building the input file. The plasmon pole screening, exchange selfenergy and the quasiparticle energies are also saved in databases so that they can be reused in further runs:
$ ls ./10b_1Ry ndb.pp ndb.pp_fragment_1 ... ndb.HF_and_locXC ndb.QP
Convergence tests for a quasi particle calculation
Now we can check the convergence of the different variables entering in the expression of the correlation part of the self energy.
First, we focus on the parameter governing the screening matrix you have already seen in the RPA/IP section. In contrast to the calculation of the RPA/IP dielectric function, where you considered either the optical limit or a finite q response (EELS), here the dielectric matrix will be calculated for all qpoints determined by the choice of kpoint sampling.
The parameters that need to be converged can be understood by looking at expression of the dielectric matrix:
where χ_{GG'} is given by
and χ^{0}_{GG'} is given by

NGsBlkXp
: The dimension of the microscopic inverse matrix, related to Local fields 
BndsRnXp
: The sum on bands (c,v) in the independent particle χ^{0}_{GG'}
Converging Screening Parameters
Here we will check the convergence of the gap starting from the variables controlling the screening reported above: the bands employed to build the RPA response function BndsRnXp
and the number of blocks (G,G') of the dielectric matrix ε^{1}_{G,G'} NGsBlkXp
.
In the next section we will study convergence with respect to the sum over states summation (sum over m in the Σ_{c} expression); here let's fix GbndRnge
to a reasonable value (40 Ry).
Let's build a series of input files differing by the values of bands and block sizes in χ^{0}_{GG'} considering BndsRnXp
in the range 1050 (upper limit) and NGsBlkXp
in the range 1 to 5 Ry. To do this by hand, file by file, open the gw_ppa.in file in an editor and change to:
NGsBlkXp = 2 Ry
while leaving the rest untouched. Repeat for 3 Ry, 4 Ry etc. Next, for each NGsBlkXp
change to:
% BndsRnXp 1  20  # [Xp] Polarization function bands %
and repeat for 30, 40 and so on. Give a different name to each file: gw_ppa_Xb_YRy.in with X=10,20,30,40 and Y=1,2,3,4,5 Ry.
This is obviously quite tedious. However, you can automate both the input construction and code execution using bash or python scripts (indeed later you will learn how to use the the Yambopythontool for this task). For now, you can use the simple generate_inputs_1.sh bash script to generate the required input files (after copying the script you need to type $ chmod +x name_of_the_script.sh
).
Finally launch the calculations:
$ yambo F gw_ppa_10b_1Ry.in J 10b_1Ry $ yambo F gw_ppa_10b_2Ry.in J 10b_2Ry ... $ yambo F gw_ppa_40b_5Ry.in J 40b_5Ry
Once the jobs are terminated we can collect the quasiparticle energies for fixed BndsRnXp
in different files named e.g. 10b.dat, 20b.dat etc. for plotting, by putting in separate columns: the energy cutoff; the size of the G blocks; the quasiparticle energy of the valence band; and that of the conduction band.
To do this e.g. for the 10b.dat file you can type:
$ cat o10b*  grep "8.000"  awk '{print $3+$4}'
to parse the valence band quasiparticle energies and
$ cat o10b*  grep "9.000"  awk '{print $3+$4}'
for the conduction band; and put all the data in the 10b.dat files. As there are many files to process you can use the parse_qps.sh script to create the 10b.dat file and edit the script changing the number of BndsRnXp
for the other output files.
Once we have collected all the quasiparticle values we can plot the gap, or the valence and conduction band energies separately, as a function of the block size or energy cutoff:
$ gnuplot gnuplot> p "10b.dat" u 1:3 w lp t " Valence BndsRnXp=10", "20b.dat" u 1:3 w lp t "Valence BndsRnXp=20",.. gnuplot> p "10b.dat" u 1:4 w lp t " Conduction BndsRnXp=10", "20b.dat" u 1:4 w lp t "Conduction BndsRnXp=20",.. gnuplot> p "10b.dat" u 1:($4$3) w lp t " Gap BndsRnXp=10", "20b.dat" u 1:($4$3) w lp t "gap BndsRnXp=20",..
or both using e.g. the ppa_gap.gnu gnuplot script:
Looking at the plot we can see that:
 The two parameters are not totally independent (see e.g. the valence quasiparticle convergence) and this is the reason why we converged them simultaneously
 The gap (energy difference) converge faster than single quasiparticle state
 The convergence criteria depends on the degree of accuracy we look for, but considering the approximations behind the calculations (plasmonpole etc.), it is not always a good idea to enforce too strict a criteria.
 Even if not totally converged we can consider that the upper limit of
BndsRnXp
=30, andNGsBlkXp
=3Ry are reasonable parameters.
Converging the sum over states in the correlation selfenergy
From now on we will keep fixed these parameters and will perform a convergence study on the sum over state summation in the correlation selfenergy (Σc) GbndRnge
.
In order to use the screening previously calculated we can copy the plasmon pole parameters saved in the 30b_3Ry
directory in the SAVE
directory. In this way the screening will be read by Yambo and not calculated again:
$ cp ./30b_3Ry/ndb_pp* ./SAVE/.
(Note: you may have to delete these files before running the BSE tutorials)
In order to use the database we have to be sure to have the same plasmonpole parameters in our input files.
Edit gw_ppa_30b_3Ry.in and modify GbndRnge
to order to have a number of bands in the range from 10 to 80 inside different files named gw_ppa_Gbnd10.in, gw_ppa_Gbnd20.in etc. You can also run the the generate_inputs_2.sh bash script to generate the required input files.
Next, launch yambo for each input:
$ yambo F gw_ppa_Gbnd10.in J Gbnd10 $ yambo F gw_ppa_Gbnd20.in J Gbnd20 ...
and as done before we can inspect the obtained quasiparticle energies:
$ grep 8.00 oGbnd*  awk '{print $4+$5}'
for the valence bands, and
$ grep 9.00 oGbnd*  awk '{print $4+$5}'
for the conduction band.
Collect the results in a text file Gbnd_conv.dat containing: the number of Gbnd, the valence energy, and the conduction energy. Now, as done before we can plot the valence and conduction quasiparticle levels separately as well as the gap, as a function of the number of bands used in the summation:
$ gnuplot gnuplot> p "Gbnd_conv.dat" u 1:2 w lp lt 7 t "Valence" gnuplot> p "Gbnd_conv.dat" u 1:3 w lp lt 7 t "Conduction" gnuplot> p "Gbnd_conv.dat" u 1:($3$2) w lp lt 7 t "Gap"
Inspecting the plot we can see that:
 The convergence with respect to
GbndRange
is rather slow and many bands are needed to get converged results.  As observed above the gap (energy difference) converges faster than the single quasiparticle state energies.
Accelerating the sum over states convergence in the correlation selfenergy
In general the convergence with respect the sum over states can be very cumbersome. Here we show how it can be mitigated by using a technique developed by F. Bruneval and X. Gonze ^{[1]}. The basic idea relies in replacing the eigenenergies of the states that are not treated explicitly by a common energy, and take into account all the states, which are not explicitly included in the calculation through the closure relation. To apply this technique in Yambo we need to activate the optional terminator variable GTermKind. We can activate it by adding the quasiparticle verbosity to the command line when generating the input file:
$ yambo p p g n F gw_ppa_Gbnd10.in V qp
and you can edit the input file by setting:
GTermKind= "BG"
or simply you can add by hand this line in all the other gw_ppa_GbndX.in input files. Now we can repeat the same calculations
$ yambo F gw_ppa_Gbnd10.in J Gbnd10_term $ yambo F gw_ppa_Gbnd20.in J Gbnd20_term ...
and collect the new results:
$ grep 8.000 *term.qp  awk '{print $4+$5}' $ grep 9.000 *term.qp  awk '{print $4+$5}'
in a new file called Gbnd_conv_terminator.dat. Now we can plot the same quantities as before by looking at the effect of having introduced the terminator.
gnuplot gnuplot> p "Gbnd_conv.dat" u 1:2 w lp t "Valence", "Gbnd_conv_terminator.dat" u 1:2 w lp t "Valence with terminator" gnuplot> p "Gbnd_conv.dat" u 1:3 w lp t "Conduction", "Gbnd_conv_terminator.dat" u 1:3 w lp t "Conduction with terminator" gnuplot> p "Gbnd_conv.dat" u 1:($3$2) w lp t "Gap", "Gbnd_conv_terminator.dat" u 1:($3$2) w lp t "Gap with terminator"
Valence band energy wrt
GbndRange
with and without the terminator technique ^{[1]}Conduction band energy wrt
GbndRange
with and without the terminator technique^{[1]}
We can see that:
 The convergence with respect to
GbndRange
is accelerated in particular for the single quasiparticle states.  From the plot above we can see that
GbndRange
=40 is sufficient to have a converged gap
Beyond the plasmon pole approximation: a full frequency approach  real axis integration
All the calculations performed up to now were based on the plasmon pole approximation (PPA). Now we remove this approximation by evaluating numerically the frequency integral in the expression of the correlation selfenergy (Σc). To this aim we need to evaluate the screening for a number of frequencies in the interval depending on the electronhole energy difference (energy difference between empty and occupied states) entering in the sum over states. Let's build the input file for a full frequency calculation by simply typing:
$ yambo F gw_ff.in g n
and we can set the variable studied up to now at their converged value:
%BndsRnXd 1  30  % NGsBlkXd=3 Ry %GbndRange 1  40  % EXXRLvcs=40 Ry % LongDrXd 1.000000  1.000000  1.000000  # [Xd] [cc] Electric Field %
and next vary in different files (name them gw_ff10.in etc.) the number of frequencies we evaluate for the screened coulomb potential. This is done by setting the values of ETStpsXd=10 , 20, 50 , 100, 150, 200, 250. You can also run the the generate_inputs_3.sh bash script to generate the required input files.
Next launch yambo:
$ yambo F gw_ff10.in J ff10 $ yambo F gw_ff50.in J ff50
...
Clearly as we are evaluating the screening for a large number of frequencies these calculations will be heavier than the case above where the calculation of the screening was done for only two frequencies (zero and plasmonpole). As before, collect the valence and conduction bands as a function of the number of frequencies in a file called gw_ff.dat and plot the behaviour of the conduction and valence bands and the gap.
$ gnuplot gnuplot> p "gw_ff.dat" u 1:2 w lp t "Valence" gnuplot> p "gw_ff.dat" u 1:3 w lp t "Conduction" gnuplot> p "gw_ff.dat" u 1:($3$2) w lp t "Gap"
We can see that:
 Oscillations are still present, indicating that even more frequencies have to be considered. In general, a realaxis calculation is very demanding.
 The final result of the gap obtained up to now does not differ too much from the one obtained at the plasmonpole level (~50 meV)
Secant Solver
The real axis integration permits also to go beyond the linear expansion to solve the quasi particle equation. The QP equation is a nonlinear equation whose solution must be found using a suitable numerical algorithm.
The mostly used, based on the linearization of the selfenergy operator is the Newton method that is the one we have used up to now. Yambo can also perform a search of the QP energies using a nonlinear iterative method based on the Secant iterative Method.
In numerical analysis, the secant method is a rootfinding algorithm that uses a succession of roots of secant lines to better approximate a root of a function f. The secant method can be thought of as a finite difference approximation of Newton's method. The equation that defines the secant method is:
The first two iterations of the secant method are shown in the following picture. The red curve shows the function f and the blue lines are the secants.
To see if there is any nonlinear effect in the solution of the Dyson equation we compare the result of the calculation using the Newton solver as done before with the present case. In order to use the secant method you need to edit one of the the previous gw_ffN.in files e.g. gw_ff100.in and substitute:
DysSolver= "g"
to
DysSolver= "s"
Repeat the calculations in the same way as before:
$ yambo F gw_ff100.in J ff100
Note than now the screening will not be calculated again as it has been stored in the ffN directories as can be seen in the report file:
[05] Dynamical Dielectric Matrix ================================ [RD./ff10//ndb.em1d] ... [RD./ff10//ndb.em1d_fragment_1] ...
Comparing the output files, e.g. for the case with 100 freqs:
off100.qp Newton Solver:
# Kpoint Band Eo EEo ScEo # 7.00000 8.00000 0.41188 0.05933 2.94162 7.000000 9.000000 3.877976 1.393774 3.446402
off100.qp_01 Secant Solver:
# # Kpoint Band Eo EEo ScEo # 7.00000 8.00000 0.41188 0.05950 2.95728 7.000000 9.000000 3.877976 1.373645 3.759375
From the comparison we see that the effect is of the order of 20 meV, of the same order of magnitude of the accuracy of the GW calculations.
Step 3: Interpolating Band Structures
Up to now we have checked convergence for the gap. Now we want to calculate the quasiparticle corrections across the Brillouin zone in order to visualize the entire band structure along a path connecting high symmetry points.
To do that we start by calculating the QP correction in the plasmonpole approximation for all the k points of our sampling and for a number of bands around the gap. You can use a previous input file or generate a new one: yambo F gw_ppa_all_Bz.in x p p g n
and set the parameters found in the previous tests:
EXXRLvcs= 40 Ry % BndsRnXp 1  30  # [Xp] Polarization function bands % NGsBlkXp= 3 Ry # [Xp] Response block size % LongDrXp 1.000000  1.000000  1.000000  # [Xp] [cc] Electric Field % PPAPntXp= 27.21138 eV # [Xp] PPA imaginary energy % GbndRnge 1  40  # [GW] G[W] bands range %
and we calculate it for all the available kpoints by setting:
%QPkrange # [GW] QP generalized Kpoint/Band indices 1 14 611 %
Note that as we have already calculated the screening with these parameters you can save time and use that database either by running in the previous directory by using J 30b_3Ry
or if you prefer to put the new databases in the new all_Bz directory you can create a new directory and copy there the screening databases:
$ mkdir all_Bz $ cp ./30b_3Ry/ndb.pp* ./all_Bz/
and launch the calculation:
$ yambo F gw_ppa_all_Bz.in J all_Bz
Now we can inspect the output and see that it contains the correction for all the k points for the bands we asked:
# Kpoint Band Eo EEo ScEo # 1.000000 6.000000 1.299712 0.215759 3.788044 1.000000 7.000000 1.296430 0.238232 3.788091 1.000000 8.000000 1.296420 0.239959 3.785949 1.000000 9.000000 4.832399 0.951691 3.679265 1.00000 10.00000 10.76416 2.10457 4.38743 1.00000 11.00000 11.36167 2.48237 3.91021 ....
By plotting some of the 'oall_Bz.qp" columns it is possible to discuss some physical properties of the hBN QPs. Using columns 3 and (3+4), ie plotting the GW energies with respect to the LDA energies we can deduce the band gap renormalization and the stretching of the conduction/valence bands:
$ gnuplot gnuplot> p "oall_Bz.qp" u 3:($3+$4) w p t "Eqp vs Elda"
Essentially we can see that the effect of the GW self energy is the opening of the gap and a linear stretching of the conduction/valence bands that can be estimated by performing a linear fit of the positive and negative energies (the zero is set at top of the valence band).
In order to calculate the band structure, however, we need to interpolate the values we have calculated above on a given path. In Yambo the interpolation is done by the executable ypp
(Yambo Post Processing).
By typing:
$ ypp H
you will recognize that in order to interpolate the bands we need to build a ypp input file using
$ ypp s b
Before editing the ypp.in input file and running the interpolation, it is important to know that ypp
uses an algorithm ^{[2]} that cannot be used in presence of timereversal (TR) symmetry.
As a first step we therefore remove the TR symmetry by typing:
$ ypp y
and we uncomment the corresponding line to remove the TR.
fixsyms # [R] Reduce Symmetries #RmAllSymm # Remove all symmetries #RmTimeRev # Remove Time Reversal
and launch
$ ypp
This will create a new directory called FixSymm
where a SAVE
directory containing the electronic structure in the absence of TR is present.
We will calculate the band structure in the FixSymm
directory.
$ cd FixSymmm
After having performed the usual setup
$yambo
we can generate the input for the band interpolation:
$ypp s b F ypp_bands.in
and edit the ypp_bands.in file:
electrons # [R] Electrons (and holes) bnds # [R] Bands cooIn= "rlu" # Points coordinates (in) cc/rlu/iku/alat BANDS_steps=0 # Number of divisions % INTERPGrid 1 1 1  # Interpolation BZ Grid % ShellFac= 20.00000 # The bigger it is a higher number of shells is used %QPkrange # generalized Kpoint/Band indices 1 20 1100 % %BKpts # K points of the bands circuit 0.00000 0.00000 0.00000  %
We modify the following lines:
BANDS_steps=30 %QPkrange # generalized Kpoint/Band indices 1 20 611 % %BKpts # K points of the bands circuit 0.33300 .66667 0.00000  0.00000 0.00000 0.00000  0.50000 .50000 0.00000  0.33300 .66667 0.00000  0.33300 .66667 0.50000  0.00000 0.00000 0.50000  0.50000 .50000 0.50000  %
which means we assign 30 point in each segment, we ask to interpolate 3 occupied and 3 empty bands and we assign the following path passing from the high symmetry points: K Γ M K H A L. Launching:
$ ypp F ypp_bands.in
will produce the output file o.bands_interpolated containing:
# # k b6 b7 b8 b9 b10 b11 kx ky kz # # 0.00000 7.22092 0.13402 0.13395 4.67691 4.67694 10.08905 0.33300 0.66667 0.00000 0.03725 7.18857 0.17190 0.12684 4.66126 4.71050 10.12529 0.32190 0.64445 0.00000
...
and we can plot the bands using gnuplot:
$ gnuplot gnuplot> p "o.bands_interpolated" u 0:2 w l, "o.bands_interpolated" u 0:3 w l, ...
and you can recognize the index of the high symmetry point by inspecting the last three columns.
Note that up to now we have interpolated the LDA band structure. In order to plot the GW band structure we need to tell ypp
in the input file where the ndb.QP database is found. This is achieved by adding in the ypp_bands.in file the line:
GfnQPdb= "E < ../all_Bz/ndb.QP"
and relaunch
$ ypp F ypp_bands.in
Now the file o.bands_interpolated_01 contains the GW interpolated band structure. We can plot the LDA and GW band structure together by using the gnuplot script bands.gnu.
 As expected the effect of the GW correction is to open the gap.
 Comparing the obtained band structure with the one found in the literature by Arnaud and coworkers ^{[3]} we found a very nice qualitative agreement.
 Quantitatively we found a smaller gap: about 5.2 eV (indirect gap), 5.7 eV (direct gap) while in Ref.^{[3]} is found 5.95 eV for the indirect gap and a minimum direct band gap of 6.47 eV. Other values are also reported in the literature depending on the used pseudopotentials, starting functional and type of selfconsistency (see below).
 The present tutorial has been done with a small k point grid which is an important parameter to be checked, so convergence with respect the k point sampling has to be validated.
Step 4: Summary of the convergence parameters
We have calculated the band structure of hBN starting from a DFT calculation, here we summarize the main variable we have checked to achieve convergence:

EXXRLvcs
# [XX] Exchange RL components
Number of Gvectors in the exchange. This number should be checked carefully. Generally a large number is needed as the QP energies show a slow convergence. The calcualtion of the exchange part is rather fast.
BndsRnXp
#[Xp] Polarization function bands
Number of bands in the independent response function form which the dielectric matrix is calculated. Also this paramater has to be checked carefully,together with NGsBlkXp as the two variables are interconnected
NGsBlkXp
# [Xp] Response block size
Number of Gvectors block in the dielectric constant. Also this paramater has to be checked carefully, to be checked together with BndsRnXp. A large number of bands and block can make the calculation very demanding.
LongDrXp
# [Xp] [cc] Electric Field
Direction of the electric field for the calculation of the q=0 component of the dielectric constant e(q,w). In a bulk can be set to (1,1,1), attention must be paid for non 3D systems.
PPAPntXp
# [Xp] Plasmon pole imaginary energy: this is the second frequency used to fit the GodbyNeeds plasmonpole model (PPM). If results depend consistently by changing this frequency, the PPM is not adequate for your calculation and it is need to gp beyond that, e.g. Realaxis.
GbndRnge
# [GW] G[W] bands range
Number of bands used to expand the Green's function. This number is usually larger than the number of bands used to calculated the dielectricconstant. Single quasiparticle energies converge slowly with respect GbndRnge, energy difference behave better. You can use terminator technique to mitigate the slow dependence.
GDamping
# [GW] G[W] damping
Small damping in the Green's function definition, the delta parameter. The final result shouuld not depend on that, usually set at 0.1 eV
dScStep
# [GW]
Energy step to evaluate Z factors
DysSolver
# [GW] Dyson Equation solver ("n","s","g")
Parameters related to the solution of the Dyson equation, "n" Newton linearization, 's' non linear secant method
GTermKind
[GW] GW terminator
Terminator for the selfenergy^{[1]} . We have seen how this spped up the convergence with respect empty bands.
QPkrange
# [GW] QP generalized Kpoint/Band indices
Kpoints and band range where you want to calculate the GW correction. The syntax is first kpoint  last kpoint  first band  last band
Step 5: Eigenvalue only selfconsistent GW0
Here we want to see how we can compute an eigenvalue only GW0 correction in Yambo. By eigenvalue only we mean that the quasiparticle wavefunctions are assumed to be the LDA wavefunctions and are not changed during the loop, but only the quasiparticle energies are updated. Moreover, we will not update the screening, but only the Green function. It is also possible to update the screening, but we leave that as an exercise.
We will use the same parameters used before for the screening and number of bands. In contrast to the calculations done before now we need to calculate the correction for all the bands that will enter in the Green function so we can copy the screening from before:
$ mkdir SC $ cp ./all_Bz/ndb.p* ./SC/.
and copy the gw_ppa_all_Bz.in file to gw_evsc.in in which we set:
%QPkrange # [GW] QP generalized Kpoint/Band indices 1 14 611 %
and run the first run of the cycle, exactly as before:
$ yambo F gw_evsc.in J SC
Once the calculation is terminated we move the QP database in the current directory
$ mv ./SC/ndb.QP .
and we modify the input file in order to tell yambo to read the quasiparticle energies calculated in the previous run instead of the KS eigenvalue. This is achieved by adding in the input file:
GfnQPdb= "E < ndb.QP"
and launch yambo again. Next we move again the new quasi particle in the current directory and repeat the loop. Each run will last around 8 minutes in serial in the cecam cluster. We compare the output file "oSC.qp", "oSC.qp_01" and we stop when desired accuracy is reached. Looking for example at the last filled valence and conduction bands at Gamma point:
$ grep " 1.000000 8.000000 " oSC.qp* oSC.qp: 1.000000 8.000000 1.296420 0.239959 3.785949 oSC.qp_01: 1.000000 8.000000 1.296420 0.275293 3.745405 oSC.qp_02: 1.000000 8.000000 1.296420 0.287520 3.731213 oSC.qp_03: 1.000000 8.000000 1.296420 0.289822 3.728570
$ grep " 1.000000 9.000000 " oSC.qp* oSC.qp: 1.000000 9.000000 4.832399 0.951691 3.679265 oSC.qp_01: 1.000000 9.000000 4.832399 1.115555 3.503588 oSC.qp_02: 1.000000 9.000000 4.832399 1.136954 3.480806 oSC.qp_03: 1.000000 9.000000 4.832399 1.139693 3.477701
we can see that the calculation is converged after 4 iterations. We can now interpolate the band structure using the last ndb.QP:
$ cd FixSymm
Edit the "ypp_bands.in" setting:
GfnQPdb= "E < ../SC/ndb.QP"
and launch
$ ypp F ypp_bands.in
We can then compare the G0W0 with the new GW0 bands structure (you can use the gnuplot script bands.gnu substituting the new band file o.bands_interpolated_02
 The eigenvalue self consistency in G, i.e. GW0 gives a correction of about 0.2/0.25 eV giving an indirect gap of 5.45 eV and a direct gap of 5.9 eV in agreement with the calcualation in Wirtz et al. ^{[4]})
 Important effects on the gap are found in the literature when considering Self consistency also for the screening, e.g. N. Berseneva et al. ^{[5]}). Eigenvalue GW self consistency can be calculated with Yambo repeating the same procedure seen for the GW0 method and adding in the input file the keyword
XfnQPdb= "E < ./SC/ndb.QP"
and removing the ndb.PP* databases after each iteration, or alternatively redirecting the databases in different directories e.g. launching
$ yambo F gw_evsc.in J SC1 $ cp ./SC1/ndb.QP . $ yambo F gw_evsc.in J SC2 $ cp ./SC2/ndb.QP . ...
Step 6: Exercise: Convergence with respect K points
As an exercise now you can check the convergence with respect the K point sampling:
 perform a new nonscf calculation with a bigger k point grid: 9x9x3 and 12x12x4 ...
 convert wave functions and electronic structure to Yambo databases in a different directory as explained in the DFT and p2y module,
 Initialize the Yambo databases,
 Redo the steps explained in this section (exchange self energy, plasmon pole GW, band structure interpolation)
 The PPAGW calculation using 12x12x4 grid takes around 12 minutes in serial mode in the Cecam cluster. You can think to perform the exercise after having learned some basic on the parallelization strategy.
Prev: Tutorials Home  Now: Tutorials Home > GW  Next: If you did everything, choose another tutorial in the menu 
References
 ↑ ^{1.0} ^{1.1} ^{1.2} ^{1.3} F. Bruneval and X. Gonze, Physical Review B 78, 085125 (2008 )
 ↑ Warren E. Pickett, Henry Krakauer, and Philip B. Allen Phys. Rev. B 38, 2721
 ↑ ^{3.0} ^{3.1} ^{3.2} B. Arnaud, S. Lebegue,P. Rabiller, and M. Alouani Phys, Rev. Lett. 96, 026402 (2006)
 ↑ L. Wirtz, A. Marini, M. Gruning, A. Rubio in arXiv:condmat/0508421 [condmat.mtrlsci]
 ↑ N. Berseneva, A. Gulans, AV. Krasheninnikov and R. M. Nieminen1 Phys. Rev. B 87, 035404 (2013)