ABINIT, sixth lesson of the tutorial:

The quasi-particle band structure of Silicon, in the GW approximation

This lesson aims at showing how to get self-energy corrections to the DFT Kohn-Sham eigenvalues in the GW approximation.
The GW formalism will NOT be explained in this tutorial. The reader might consult the review

"Quasiparticle calculations in solids", by Aulbur WG, Jonsson L, Wilkins JW,
in Solid State Physics 54, 1-218 (2000),
also available at http://www.physics.ohio-state.edu/~wilkins/vita/gw_review.ps

The different formulas of the GW formalism, that have been implemented in ABINIT, have been written in a pdf document by Valerio Olevano (who also wrote the first version of this tutorial), see ~ABINIT/Infos/Theory/gwa.pdf .

This lesson should take about 2 hours to be done.

Copyright (C) 2002-2004 ABINIT group (VOlevano,XG)
This file is distributed under the terms of the GNU General Public License, see ~ABINIT/Infos/copyright or http://www.gnu.org/copyleft/gpl.txt .
For the initials of contributors, see ~ABINIT/Infos/contributors .

Help files : New user's guide | Abinis (main) | Abinis (respfn) | Mrgddb | Anaddb

Content of lesson 6

6.1 General example of well converged GW calculation.
6.2 Calculation of the Kohn-Sham structure (KSS file) and of the screening (SCR file).
6.3 Convergence on the number of planewaves in the wavefunctions to calculate the Self-Energy.
6.4 Convergence on the number of planewaves to calculate Sigma_x.
6.5 Convergence on the number of bands to calculate the Self-Energy.
6.6 Convergence on the number of planewaves in the wavefunctions to calculate the screening (epsilon^-1).
6.7 Convergence on the number of bands to calculate the screening.
6.8 Convergence on the dimension of the epsilon^-1 matrix.
6.9 Calculation of the GW corrections for the band gap in Gamma.

6.1 Computation of the Silicon band gap at Gamma, using a GW calculation.

Before beginning, you might consider to work in a different subdirectory as for the other lessons. Why not "Work6" ?

At the end of lesson 3, you computed the Kohn-Sham band structure of Silicon. In this approximation, the variation of eigenvalues inside each band is reasonable, as well as the band widths, but the band gaps are known to be qualitatively wrong. Now, we will compute the band gaps much more accurately, using the so-called GW approximation.

We start by an example, in which we show how to perform in one shot the calculation of the ground state, the Kohn Sham electronic structure, the screening, and the Self-Energy matrix elements, that is, the GW corrections, for one given k-point, for the highest occupied and the lowest empty bands. We provide some reasonable parameters without checking convergence. You will see that this procedure is MUCH MORE time-consuming than the corresponding calculation of the Kohn-Sham eigenvalues.

So, let us run immediately this calculation, and while it is running, we will explain what has been done.

In directory ~ABINIT/Tutorial/Work6, copy the files ~ABINIT/Tutorial/t6x.files and t61.in, and modify the t6x.files file as usual.
Then, issue :

../../abinis < t6x.files >& t61.log &

It is very important to run this job in background. Indeed, a PC Intel PIV/2.2 GHz will take about 6 minutes to complete it. In the meantime, you should read the following.

6.1.a The three steps of a GW calculation.

In order to complete a standard GW calculation, one has to :

Run a converged Ground State calculation (at fixed lattice parameters and atomic positions), to get self-consistent density and potential, and Kohn-Sham eigenvalues and eigenfunctions at the relevant k-point as well as on a regular grid of k-points ;
On the basis of these available Kohn-Sham data, compute the independent-particle susceptibility matrix ("chi0"), on a regular grid of wavevectors, for at least two frequencies (usually, zero frequency and a large frequency - on the order of the plasmon frequency, a dozen of eV), then compute the dielectric matrix ("epsilon") in the same conditions, its inverse, and the Random-Phase susceptibility matrix ("chi") ;
On this basis, compute the self-energy operator ("sigma") at a given k-point, and derive the GW eigenvalues for the target states at this k-point.

The input file t61.in has precisely that structure : there are three datasets. The first dataset starts a rather usual SCF calculation, then will construct a specialized file, t6xo_DS1_KSS (_KSS for "Kohn-Sham Structure), that contains the needed information to start step 2. The second dataset drives the computation of susceptibility and dielectric matrices, giving another specialized file, t6xo_DS2_SCR (_SCR for "Epsilon Minus 1" - the inverse dielectric matrix). Then, in the third dataset, one builds the eigenvalues of the 4th and 5th bands at the Gamma point.

So, you can edit this t61.in file.

In the dataset-independent part of this file (the last half of the file), there is the usual set of input variables, describing the cell, atom types, number, position, planewave cut-off energy, SCF convergence parameters, than in the t35.in file, driving the Kohn-Sham band structure calculation. Then, for the three datasets, you will find specialized additional input variables.

6.1.b Generating the Kohn-Sham band structure : the KSS file.

In dataset 1, apart from the usual input variables we are acquainted to through the previous tutorials, there is a new input variable:

nbandkss -1        # Number of bands in KSS file (-1 means the maximum possible)

This input variable tells the program to calculate the Kohn-Sham electronic structure by the (in this case) full diagonalization of the Kohn-Sham Hamiltonian evaluated at the converged density and calculated in each one of the k-points of the grid. Note that this diagonalization is performed in a routine (outkss.f) separated from the usual SCF cycle, so that there is additional control of the wavefunction actually stored, if needed. In particular, the number of bands to be computed in this routine is NOT determined by the usual input variable nband.

nbandkss is the key variable to create a _KSS file. If it is zero, no _KSS file will be created. -1 lead to the generation and storage of the maximum possible number of states (or bands) common to all points. This might lead to quite time-consuming calculations. One can reduce the load in the diagonalization by requiring less states.

Another way to reduce the load of the diagonalization is made possible through the use of npwkss. It governs the size of the plane wave basis set in which the Hamiltonian matrix will be expressed and diagonalized. The default value leaves the number of plane wave equal to the one of the SCF ground state calculation.

Another relevant input variable, related also to the specific set up of the _KSS file is kssform.

In this first dataset, we asked also the self-consistent cycle to be done for nine bands.

nband1      9         # Number of (occ and empty) bands to be computed

Only four bands would be needed for Si. The purpose of defining more bands in the ground-state determination is to verify that at least the first Kohn-Sham eigenvalues obtained through the diagonalization are sufficiently close to those determined in the self-consistent procedure. At present, the comparison is not done automatically, so please check (well, sometimes ...) that the Kohn-Sham eigenvalues given in the self-consistency part (with a residual) are close to those given after the diagonalization.

6.1.c Generating the screening : the SCR file.

In dataset 2 the calculation of the screening (susceptibility, dielectric matrix) is performed. We need to set optdriver=3 to do that :

optdriver2  3        # Screening calculation

The getkss input variable is similar to other "get" input variables of ABINIT :

getkss2     -1       # Obtain KSS file from previous dataset

In this case, it tells the code to use the KSS file calculated in the previous dataset.

Then, three input variables describe the computation :

nband2      25       # Bands to be used in the screening calculation
ecutwfn2    2.1      # Cut-off energy of the planewave set to represent the wavefunctions
ecuteps2    3.6      # Cut-off energy of the planewave set to represent the dielectric matrix

In this case, we use 25 bands to calculate the Kohn-Sham response function $\chi^{(0)}_{KS}$; we use a cut-off ecutwfn=2.1 Hartree, giving 89 planewaves to represent the wavefunctions in the calculation of $\chi^{(0)}_{KS}$. The dimension of $\chi^{(0)}_{KS}$, as well as all the other matrices ($\chi$, $\espilon$) is determined by ecuteps=3.6 Hartree, giving 169 planewaves.

Finally we define the frequencies at which the screening must be evaluated : $\omega = 0.0 eV$ and the imaginary frequency $\omega = i 16.7 eV$. The latter is determined by the input variable plasfrq

plasfrq2    16.7 eV  # Imaginary frequency where to calculate the screening

The two frequencies are used to calculate the plasmon-pole model parameters. For the non-zero frequency it is recommanded to use a value close to the plasmon frequency for the plasmon-pole model to work well. Plasmons frequencies are usually close to 0.5 Hartree. The parameters for the screening calculation are not far from the ones that give converged Energy Loss Function (-Im \epsilon^-1_00) spectra, So that one can start up by using indications from EELS calculations existing in literature.

6.1.d Computing the GW energies.

In dataset 3 the calculation of the Self-Energy matrix elements is performed. One needs to define the driver option, as well as the _KSS and _SCR files.

optdriver3  4        # Self-Energy calculation
getkss3     -2       # Obtain KSS file from dataset 1
getscr3     -1       # Obtain SCR file from previous dataset

The getscr input variable is also similar to other "get" input variables of ABINIT.

Then, comes the definition of parameters needed to compute the self-energy. As for the computation of the susceptibility and dielectric matrices, one must define the set of bands, and two sets of planewaves :

nband3      100      # Bands to be used in the Self-Energy calculation
ecutwfn3    5.0      # Planewaves to be used to represent the wavefunctions
ecutsigx3    6.0      # Dimension of the G sum in Sigma_x
                     # (the dimension in Sigma_c is controlled by npweps)

In this case, nband controls the number of bands used to calculate the Self-Energy. ecutwfn defines (as for optdriver=3) the number of planewaves used to represent the wavefunctions. ecutmat gives the dimension of the planewave sum needed to calculate Sigma_x (the exchange part of the self-energy, which is diagonal). The size of the planewave set needed to compute Sigma_c (the correlation part of the self-energy) is controlled by the dimension of the screening matrix read in the SCR file. However, it is taken equal to the number of planewave of Sigma_x if the latter is smaller than the one for Sigma_c.

Then, come the parameters defining the k-points and states for which the electronic energy must be computed :

nkptgw3      1               # number of k-point where to calculate the GW correction
kptgw3                       # k-points
  -0.125    0.000    0.000
bdgw3       4  5             # calculate GW corrections for bands from 4 to 5

nkptgw tells the number of k-points for which the GW corrections must be computed. The k-points coordinates are given in kptgw. At present, they must belong to the grid of k-points defined with the same repetition parameters (kptrlatt, or ngkpt) than the GS one, but WITHOUT any shift. bdgw gives the minimum/maximum band whose energies are calculated.

There is an additional parameter, called zcut, related to the self-energy computation. It is meant to avoid some divergencies that might occur in the calculation due to integrable poles along the integration path.

6.1.e Examination of the output file.

Let us hope now that your calculation has been completed, and that we can examine the output file. So, please edit the t61.out file.

The first departure from the usual information present in the output file for usual GS calculations appears after the SCF cycles of DATASET 1 :

======================================================================
 Calculating and writing out Kohn-Sham electronic Structure file  
 Using diagonalized wavefunctions and energies (kssform=1)
 number of Gamma centered plane waves    483
 number of Gamma centered shells     25
 number of bands    283

This section was issued when the Hamiltonian at the different k points were diagonalized, after the SCF cycles, in order to generate the KSS file. Then, comes the output of the numerous eigenvalues at the different k-points. Finally, the normalisation and orthogonalisation of the eigenvectors is tested. One should obtain perfect normalisation and orthogonalisation at that stage :

 Test on the normalization of the wavefunctions
  min sum_G |a(n,k,G)| =  1.000000
  max sum_G |a(n,k,G)| =  1.000000
 Test on the orthogonalization of the wavefunctions
  min sum_G a(n,k,G)* a(n',k,G) =  0.000000
  max sum_G a(n,k,G)* a(n',k,G) =  0.000000

Then, follows the usual information for the dataset 1. The dataset 2 drives the computation of the susceptibility and dielectric matrices, in preparation of the GW energy calculation of dataset 3. After some general information (origin of KSS file, header, description of unit cell), comes the echo of Kohn-Sham eigenenergies (in eV), and then the evaluation of the wavefunction normalisation and orthogonalisation USING ONLY THE PLANEWAVE SET DEFINED BY ecutwfn, npwwfn, or nshwfn. Thus, there is no surprise that these relations are not fulfilled :

 test on the normalization of the wavefunctions
 min sum_G |a(n,k,G)| =  0.497559
 max sum_G |a(n,k,G)| =  0.995840
 test on the orthogonalization of the wavefunctions
 min sum_G a(n,k,G)* a(n",k,G) =  0.000000
 max sum_G a(n,k,G)* a(n",k,G) =  0.179460

The squared norm of one of the wavefunctions is even as low as one half ! This should lead us to question the choice of ecutwfn that we have made : we will need a convergence study, see later.

The parameters of the FFT grid needed to represent the wavefunction and compute their convolution (so as to get the screening matrices) are then given.

Then, the grid of q-point (in the Irreducible part of the Brillouin Zone) on which the susceptibility and dielectric matrices will be computed is given. It is a grid of points with the same repetition parameters (kptrlatt, or ngkpt) than the GS one, but WITHOUT any shift.

On the basis of only the average density, one can obtain the plasmon frequency of metallic jellium (homogeneous electron gas, placed in a neutralizing background). The next lines start from the average density of the system, and evaluate the r_s parameter of the jellium, then compute the plasmon frequency. THIS IS A ROUGH ESTIMATE. In particular, it will be questionable for strongly inhomogeneous systems ! Also, the choice of pseudopotential (inclusion of core states) will have an effect on this estimate. So, take it cautiously. It is better to try a few values of plasfrq than to rely blindly on this value ...

At the end of the screening calculation, the macroscopic dielectric constant is printed:

 dielectric constant = 13.8476
 dielectric constant without local fields = 15.5520

Note that the convergence in the dielectric constant DOES NOT guarantee the convergence in the GW correction values at the end of the calculation. In fact, the dielectric constant is representative of only one element, the head, of the full dielectric matrix. Even if the convergence on the dielectric constant with local fields takes somehow into account also other non-diagonal elements. In a GW calculation all the \epsilon^-1 matrix is used to build the Self-Energy operator.
The dielectric constant here reported is the so-called RPA dielectric constant due to the electrons. Although evaluated at zero frequency, it is understood that the ionic response is not included. This is to be contrasted with the one computed in ANADDB). The RPA dielectric constant restricted to electronic effects is also not the same as the one computed in the RESPFN part of ABINIT, that includes exchange-correlation effects.

We enter now the third dataset. As for dataset 2, after some general information (origin of KSS file, header, description of unit cell), the echo of Kohn-Sham eigenenergies (in eV), the evaluation of the wavefunction normalisation, the description of the FFT grid and jellium parameters, there is the echo of parameters for the plasmon-pole model, and the inverse dielectric function (the screening). The self-energy operator has been constructed, and one can evaluate the GW energies, for each of the states.

The results follows :

k =   -0.125   0.000   0.000
 Band     E0  <VxcLDA>   SigX SigC(E0)    Z dSigC/dE  Sig(E)    E-E0     E
    4   5.616 -11.132 -12.334   1.257   0.775  -0.290 -11.089   0.043   5.659
    5   8.357 -10.157  -5.951  -3.336   0.779  -0.284  -9.480   0.677   9.034

 E^0_gap          2.741
 E^GW_gap         3.375
 DeltaE^GW_gap    0.634

For the desired k-point, state 4, then state 5, one finds different information:

E0 is the Kohn-Sham eigenenergy
VxcLDA gives the average Kohn-Sham exchange-correlation potential
SigX gives the exchange contribution to the self-energy
SigC(E0) gives the correlation contribution to the self-energy, evaluated at the Kohn-Sham eigenenergy
Z is the renormalisation factor
dSigC/dE is the energy derivative of SigC with respect to the energy
SigC(E) gives the correlation contribution to the self-energy, evaluated at the GW energy
E-E0 is the difference between GW energy and Kohn-Sham eigenenergy
E is the GW energy

In this case, the gap is also analyzed : E^0_gap is the Kohn-Sham one, E^GW_gap is the GW one, and DeltaE^GW_gap is the difference.

It is seen that the average Kohn-Sham exchange-correlation potential for the state 4 (a valence state) is very close to the exchange self-energy correction. For that state, the correlation correction is small, and the difference between Kohn-Sham and GW energies is also small (43 meV). By contrast, the exchange self-energy is much smaller than the average Kohn-Sham potential for the state 5 (a conduction state), but the correlation correction is much larger than for state 4. On the whole, the difference between Kohn-Sham and GW energies is not very large, but nevertheless, it is quite important when compared with the size of the gap.

6.2 Preparing convergence studies : Kohn-Sham structure (KSS file) and screening (SCR file).

In the following sections, we will perform different convergence analyses, because this is such an important task. In order to keep the CPU time at a reasonable level, we will use fake KSS and screening data, by limiting ourselves to the Gamma point only. In that way, we will verify convergence aspects that could be very cumbersome (at least in the framework of a tutorial) if more k-points were used.

In directory ~ABINIT/Tutorial/Work6, copy the file ../t62.in, and modify the t6x.files file as usual. Edit the t62.in file, and take the time to examine it. Note that the SCF cycles have been disconnected from the generation of the KSS file.
Then, issue :

../../abinis < t6x.files >& t62.log &

This small job lasts about 10 secs on a PC PIV Intel 2.2 GHz.

After that step you will need the KSS and SCR files produced in this run for the next runs (up to 6.8). Move t6xo_DS2_KSS to t6xo_DS1_KSS and t6xo_DS3_SCR to t6xo_DS1_SCR.

The next 6 sections are intended to show you how to find the converged parameters for a GW calculation.

6.3 Convergence on the number of planewaves in the wavefunctions to calculate the Self-Energy.

We begin by the convergence study on the three parameters needed in the self-energy calculation (optdriver=4). This is because for these, we will not need a double dataset loop to check this convergence, and we will rely on the previously determined SCR file.

First, we check the convergence on the number of planewaves to describe the wavefunctions, in the calculation of the Self-Energy. This will be done by defining five datasets, with increasing ecutwfn:

ndtset     5
ecutwfn:  3.0
ecutwfn+  1.0

In directory ~ABINIT/Tutorial/Work6, copy the file ../t63.in, and modify the t6x.files file as usual. Edit the t63.in file, and take the time to examine it.
Then, issue :

../../abinis < t6x.files >& t63.log &

This small job lasts about 10 secs on a PC PIV Intel 2.2 GHz.

Edit the output file. The number of plane waves used for the wavefunctions in the computation of the self-energy is mentioned in the fragments of output :

 SIGMA fundamental parameters:
 PLASMON POLE MODEL
 number of plane-waves for SigmaX                  169
 number of plane-waves for SigmaC and W            169
 number of plane-waves for wavefunctions            59

Gathering the GW energies for each planewave set, one gets :

 number of plane-waves for wavefunctions            59
 Band     E0  VxcLDA    SigX SigC(E0)      Z dSigC/dE  Sig(E)    E-E0       E
    4   5.915 -11.651 -15.237   3.897   0.806  -0.240 -11.401   0.251   6.166
    5   8.445  -9.669  -3.222  -5.460   0.819  -0.222  -8.861   0.808   9.253

 number of plane-waves for wavefunctions           113 
 Band     E0  VxcLDA    SigX SigC(E0)      Z dSigC/dE  Sig(E)    E-E0       E
    4   5.915 -11.654 -15.244   3.789   0.804  -0.244 -11.495   0.159   6.075
    5   8.445  -9.691  -3.213  -5.564   0.817  -0.224  -8.944   0.747   9.192

 number of plane-waves for wavefunctions           137
 Band     E0  VxcLDA    SigX SigC(E0)      Z dSigC/dE  Sig(E)    E-E0       E
    4   5.915 -11.654 -15.244   3.779   0.804  -0.244 -11.502   0.151   6.066
    5   8.445  -9.702  -3.216  -5.577   0.817  -0.225  -8.960   0.743   9.188

 number of plane-waves for wavefunctions           169
 Band     E0  VxcLDA    SigX SigC(E0)      Z dSigC/dE  Sig(E)    E-E0       E
    4   5.915 -11.651 -15.242   3.770   0.804  -0.245 -11.508   0.144   6.059
    5   8.445  -9.718  -3.221  -5.584   0.817  -0.225  -8.972   0.745   9.190

 number of plane-waves for wavefunctions           259
 Band     E0  VxcLDA    SigX SigC(E0)      Z dSigC/dE  Sig(E)    E-E0       E
    4   5.915 -11.667 -15.253   3.766   0.803  -0.245 -11.522   0.145   6.060
    5   8.445  -9.716  -3.219  -5.591   0.816  -0.225  -8.977   0.740   9.185

So that npwwfn=137 (ecutwfn=5.0) can be considered converged within 0.01eV.

6.4 Convergence on the number of planewaves to calculate Sigma_x.

Second, we check the convergence on the number of planewaves in the calculation of the Sigma_X. This will be done by defining five datasets, with increasing ecutmat:

ndtset     7
ecutsigx:  3.0
ecutsigx+  1.0

In directory ~ABINIT/Tutorial/Work6, copy the file ../t64.in, and modify the t6x.files file as usual. Edit the t64.in file, and take the time to examine it.
Then, issue :

../../abinis < t6x.files >& t64.log &

This small job lasts about 12 secs on a PC PIV Intel 2.2 GHz.

Edit the output file. The number of plane waves used for Sigma_X is mentioned in the fragments of output :

 SIGMA fundamental parameters:
 PLASMON POLE MODEL
 number of plane-waves for SigmaX                   59
 number of plane-waves for SigmaC and W             59

Gathering the GW energies for each planewave set, one gets :

 number of plane-waves for SigmaX                   59
 number of plane-waves for SigmaC and W             59
 Band     E0  VxcLDA    SigX SigC(E0)      Z dSigC/dE  Sig(E)    E-E0       E
    4   5.915 -11.654 -15.195   3.862   0.806  -0.241 -11.395   0.259   6.174
    5   8.445  -9.702  -3.177  -5.595   0.818  -0.223  -8.941   0.761   9.206

 number of plane-waves for SigmaX                  113
 number of plane-waves for SigmaC and W            113
 Band     E0  VxcLDA    SigX SigC(E0)      Z dSigC/dE  Sig(E)    E-E0       E
    4   5.915 -11.654 -15.235   3.795   0.804  -0.244 -11.482   0.172   6.087
    5   8.445  -9.702  -3.210  -5.581   0.817  -0.224  -8.958   0.744   9.189

 number of plane-waves for SigmaX                  137
 number of plane-waves for SigmaC and W            137
 Band     E0  VxcLDA    SigX SigC(E0)      Z dSigC/dE  Sig(E)    E-E0       E
    4   5.915 -11.654 -15.241   3.785   0.804  -0.244 -11.495   0.159   6.074
    5   8.445  -9.702  -3.213  -5.577   0.817  -0.224  -8.958   0.745   9.190

 number of plane-waves for SigmaX                  169
 number of plane-waves for SigmaC and W            169
 Band     E0  VxcLDA    SigX SigC(E0)      Z dSigC/dE  Sig(E)    E-E0       E
    4   5.915 -11.654 -15.244   3.779   0.804  -0.244 -11.502   0.151   6.066
    5   8.445  -9.702  -3.216  -5.577   0.817  -0.225  -8.960   0.743   9.188

 number of plane-waves for SigmaX                  259
 number of plane-waves for SigmaC and W            169
 Band     E0  VxcLDA    SigX SigC(E0)      Z dSigC/dE  Sig(E)    E-E0       E
    4   5.915 -11.654 -15.247   3.779   0.804  -0.244 -11.504   0.150   6.065
    5   8.445  -9.702  -3.218  -5.577   0.817  -0.225  -8.961   0.741   9.186

 number of plane-waves for SigmaX                  283
 number of plane-waves for SigmaC and W            169
 Band     E0  VxcLDA    SigX SigC(E0)      Z dSigC/dE  Sig(E)    E-E0       E
    4   5.915 -11.654 -15.247   3.779   0.804  -0.244 -11.504   0.150   6.065
    5   8.445  -9.702  -3.218  -5.577   0.817  -0.225  -8.961   0.741   9.186

 number of plane-waves for SigmaX                  283
 number of plane-waves for SigmaC and W            169
 Band     E0  VxcLDA    SigX SigC(E0)      Z dSigC/dE  Sig(E)    E-E0       E
    4   5.915 -11.654 -15.247   3.779   0.804  -0.244 -11.504   0.150   6.065
    5   8.445  -9.702  -3.218  -5.577   0.817  -0.225  -8.961   0.741   9.186

So that npwsigx=169 (ecutsigx=6.0) can be considered converged within 0.01eV.

6.5 Convergence on the number of bands to calculate the Self-Energy.

At last, as concerns the computation of the sel-energy, we check the convergence on the number of bands in the calculation of the Sigma_X. This will be done by defining five datasets, with increasing nband:

ndtset  5
nband:  50
nband+  50

In directory ~ABINIT/Tutorial/Work6, copy the file ../t65.in, and modify the t6x.files file as usual. Edit the t65.in file, and take the time to examine it.
Then, issue :

../../abinis < t6x.files >& t65.log &

This small job lasts about 12 secs on a PC PIV Intel 2.2 GHz.

Edit the output file. The number of bands used for the self-energy is mentioned in the fragments of output :

 SIGMA fundamental parameters:
 PLASMON POLE MODEL
 number of plane-waves for SigmaX                  169
 number of plane-waves for SigmaC and W            169
 number of plane-waves for wavefunctions           137
 number of bands                                    50

Gathering the GW energies for each number of bands, one gets :

 number of bands                                   50
    4   5.915 -11.654 -15.244   3.853   0.804  -0.243 -11.443   0.211   6.126
    5   8.445  -9.702  -3.216  -5.507   0.817  -0.224  -8.902   0.800   9.246

 number of bands                                  100
    4   5.915 -11.654 -15.244   3.779   0.804  -0.244 -11.502   0.151   6.066
    5   8.445  -9.702  -3.216  -5.577   0.817  -0.225  -8.960   0.743   9.188

 number of bands                                  150
    4   5.915 -11.654 -15.244   3.771   0.804  -0.244 -11.509   0.145   6.060
    5   8.445  -9.702  -3.216  -5.585   0.817  -0.225  -8.966   0.736   9.182

 number of bands                                  200
    4   5.915 -11.654 -15.244   3.769   0.804  -0.244 -11.510   0.143   6.059
    5   8.445  -9.702  -3.216  -5.587   0.817  -0.225  -8.967   0.735   9.180

 number of bands                                  250
    4   5.915 -11.654 -15.244   3.769   0.804  -0.244 -11.510   0.143   6.058
    5   8.445  -9.702  -3.216  -5.587   0.817  -0.225  -8.967   0.735   9.180

So that nband=100 can be considered converged within 0.01eV.

At this stage, we know that for the self-energy computation, we need ecutwfn=5.0 ecutmat=6.0, nband=100 .

6.6 Convergence on the number of planewaves in the wavefunctions to calculate the screening (epsilon^-1).

Now, we come back to the calculation of the screening. Adequate convergence studies will couple the change of parameters for optdriver=3 with a computation of the GW energy changes. One cannot rely on the convergence of the macroscopic dielectric constant to assess the convergence of the GW energies.

As a consequence, we will define a double loop over the datasets:

ndtset      10
udtset      5  2

The datasets 12,22,32,42 and 52, drive the computation of the GW energies :

# Calculation of the Self-Energy matrix elements (GW corrections)
optdriver?2  4
getscr?2     -1
ecutwfn?2    5.0
ecutsigx      6.0
nband?2      100

The datasets 11,21,31,41 and 51, drive the corresponding computation of the screening :

# Calculation of the screening (epsilon^-1 matrix)
optdriver?1  3

In this latter series, we will have to vary the three different parameters ecutwfn, ecuteps and nband.

First, we check the convergence on the number of planewaves to describe the wavefunctions, in the calculation of the screening. This will be done by defining five datasets, with increasing ecutwfn:

ecutwfn:?   3.0
ecutwfn+?   1.0

In directory ~ABINIT/Tutorial/Work6, copy the file ../t66.in, and modify the t6x.files file as usual. Edit the t66.in file, and take the time to examine it.
Then, issue :

../../abinis < t6x.files >& t66.log &

This small job lasts about 15 secs on a PC PIV Intel 2.2 GHz.

Edit the output file. The number of plane waves used for the wavefunctions in the computation of the screening is mentioned in the fragments of output :

 EPSILON^-1 parameters (SCR file):
 dimension of the eps^-1 matrix                    169
 number of plane-waves for wavefunctions            59

Gathering the macroscopic dielectric constant and GW energies for each planewave set, one gets :

 dielectric constant = 101.5301
 dielectric constant without local fields = 147.3095 
 number of plane-waves for wavefunctions            59
    4   5.915 -11.654 -15.244   3.799   0.806  -0.241 -11.486   0.168   6.083
    5   8.445  -9.702  -3.216  -5.555   0.816  -0.225  -8.942   0.761   9.206

 dielectric constant =  99.5265
 dielectric constant without local fields = 143.7208 
 number of plane-waves for wavefunctions           113
    4   5.915 -11.654 -15.244   3.769   0.804  -0.244 -11.510   0.143   6.059
    5   8.445  -9.702  -3.216  -5.582   0.815  -0.226  -8.964   0.738   9.183

 dielectric constant =  98.2598
 dielectric constant without local fields = 142.5982
 number of plane-waves for wavefunctions           137
    4   5.915 -11.654 -15.244   3.762   0.801  -0.248 -11.517   0.137   6.052
    5   8.445  -9.702  -3.216  -5.588   0.815  -0.227  -8.970   0.733   9.178

 dielectric constant =  97.6265
 dielectric constant without local fields = 142.1664 
 number of plane-waves for wavefunctions           169
    4   5.915 -11.654 -15.244   3.759   0.804  -0.244 -11.519   0.135   6.050
    5   8.445  -9.702  -3.216  -5.590   0.815  -0.227  -8.972   0.731   9.176

 dielectric constant =  96.4286
 dielectric constant without local fields = 140.5466 
 number of plane-waves for wavefunctions           259
    4   5.915 -11.654 -15.244   3.760   0.803  -0.245 -11.518   0.136   6.051
    5   8.445  -9.702  -3.216  -5.592   0.815  -0.227  -8.973   0.730   9.175

So that npwwfn=113 (ecutwfn=4.0) can be considered converged within 0.01eV.

6.7 Convergence on the number of bands to calculate the screening.

Second, we check the convergence on the number of bands in the calculation of the screening. This will be done by defining five datasets, with increasing nband:

   
nband11  25
nband21  50
nband31  100
nband41  150
nband51  200

In directory ~ABINIT/Tutorial/Work6, copy the file ../t67.in, and modify the t6x.files file as usual. Edit the t67.in file, and take the time to examine it.
Then, issue :

../../abinis < t6x.files >& t67.log &

This small job lasts about 22 secs on a PC PIV Intel 2.2 GHz.

Edit the output file. The number of bands used for the wavefunctions in the computation of the screening is mentioned in the fragments of output :

 EPSILON^-1 parameters (SCR file):
 dimension of the eps^-1 matrix                    169
 number of plane-waves for wavefunctions           113
 number of bands                                    25

Gathering the macroscopic dielectric constant and GW energies for each number of bands, one gets :

 dielectric constant =  99.5265
 dielectric constant without local fields = 143.7208
 number of bands                                    25
    4   5.915 -11.654 -15.244   3.769   0.804  -0.244 -11.510   0.143   6.059
    5   8.445  -9.702  -3.216  -5.582   0.815  -0.226  -8.964   0.738   9.183

 dielectric constant = 100.6436
 dielectric constant without local fields = 143.7240
 number of bands                                    50
    4   5.915 -11.654 -15.244   3.587   0.804  -0.244 -11.657  -0.003   5.912
    5   8.445  -9.702  -3.216  -5.764   0.815  -0.227  -9.113   0.589   9.034

 dielectric constant = 101.1764
 dielectric constant without local fields = 143.7244   
 number of bands                                   100
    4   5.915 -11.654 -15.244   3.516   0.804  -0.244 -11.714  -0.060   5.855
    5   8.445  -9.702  -3.216  -5.846   0.811  -0.233  -9.182   0.520   8.965

 dielectric constant = 101.2028  
 dielectric constant without local fields = 143.7244  
 number of bands                                   150
    4   5.915 -11.654 -15.244   3.510   0.804  -0.244 -11.718  -0.065   5.850
    5   8.445  -9.702  -3.216  -5.853   0.810  -0.234  -9.189   0.514   8.959

 dielectric constant = 101.2128
 dielectric constant without local fields = 143.7244
 number of bands                                   200
    4   5.915 -11.654 -15.244   3.509   0.803  -0.246 -11.719  -0.065   5.850
    5   8.445  -9.702  -3.216  -5.854   0.812  -0.231  -9.188   0.514   8.960

So that the computation using 100 bands can be considered converged within 0.01eV.

6.8 Convergence on the dimension of the epsilon^-1 matrix.

Third, we check the convergence on the number of plane waves in the calculation of the screening. This will be done by defining six datasets, with increasing ecuteps:

   
ecuteps:?     3.0
ecuteps+?     1.0

In directory ~ABINIT/Tutorial/Work6, copy the file ../t68.in, and modify the t6x.files file as usual. Edit the t68.in file, and take the time to examine it.
Then, issue :

../../abinis < t6x.files >& t68.log &

This small job lasts about 25 secs on a PC PIV Intel 2.2 GHz.

Edit the output file. The number of bands used for the wavefunctions in the computation of the screening is mentioned in the fragments of output :

 EPSILON^-1 parameters (SCR file):
 dimension of the eps^-1 matrix                     59

Gathering the macroscopic dielectric constant and GW energies for each number of bands, one gets :

 dielectric constant = 102.1281
 dielectric constant without local fields = 143.7244
 dimension of the eps^-1 matrix                     59
    4   5.915 -11.654 -15.244   3.684   0.806  -0.241 -11.579   0.075   5.990
    5   8.445  -9.702  -3.216  -5.847   0.811  -0.232  -9.183   0.519   8.964

 dielectric constant = 101.2712
 dielectric constant without local fields = 143.7244
 dimension of the eps^-1 matrix                    113
    4   5.915 -11.654 -15.244   3.559   0.804  -0.243 -11.680  -0.026   5.889
    5   8.445  -9.702  -3.216  -5.850   0.811  -0.233  -9.185   0.517   8.962

 dielectric constant = 101.2649
 dielectric constant without local fields = 143.7244   
 dimension of the eps^-1 matrix                    137
    4   5.915 -11.654 -15.244   3.535   0.804  -0.244 -11.699  -0.045   5.870
    5   8.445  -9.702  -3.216  -5.846   0.811  -0.232  -9.182   0.520   8.965

 dielectric constant = 101.1764
 dielectric constant without local fields = 143.7244   
 dimension of the eps^-1 matrix                    169
    4   5.915 -11.654 -15.244   3.516   0.804  -0.244 -11.714  -0.060   5.855
    5   8.445  -9.702  -3.216  -5.846   0.811  -0.233  -9.182   0.520   8.965

 dielectric constant = 101.1384
 dielectric constant without local fields = 143.7244   
 dimension of the eps^-1 matrix                    259
    4   5.915 -11.654 -15.244   3.517   0.804  -0.244 -11.713  -0.060   5.855
    5   8.445  -9.702  -3.216  -5.845   0.811  -0.232  -9.182   0.521   8.966

So that npweps=169 (ecuteps=6.0) can be considered converged within 0.01eV.

At this stage, we know that for the screening computation, we need ecutwfn=4.0 ecuteps=6.0, nband=100 .

Of course, until now, we have skipped the most difficult part of the convergence tests : the number of k-points. It is as important to check the convergence on this parameter, than on the other ones. However, this might be very time consuming, since the CPU time scales as the square of the number of k points (roughly), and the number of k-points can increase very rapidly from one possible grid to the next denser one. This is why we will leave this out of the present tutorial, and consider that we already know a sufficient k-point grid, for the last calculation.

6.9 Calculation of the GW corrections for the band gap in Gamma.

Now we try to perform a GW calculation for a real problem: the calculation of the GW corrections for the direct band gap of bulk Silicon in Gamma.

In directory ~ABINIT/Tutorial/Work6, copy the file ../t69.in, and modify the t6x.files file as usual. DO NOT EDIT IT NOW.
Issue :

../../abinis < t6x.files >& t69.log &

This job lasts about 20 minutes on a PC PIV Intel 2.2 GHz. Because it is so long, it was worth to run it before the examination of the input file.

Now, you can examine it.
We need the usual part of the input file to perform a ground state calculation. This is done in dataset 1 and at the end we print out the density. We use a 4x4x4 FCC grid (so, 256 k points in the full Brillouin Zone), shifted, because it is the most economical. It gives 10 k-points in the Irreducible part of the Brillouin Zone. However, this k-point grid does not contains the Gamma point, and, at present, one cannot perform calculations of the self-energy corrections for other k points than those present in the grid of k-points in the KSS file.

Then in dataset 2 we perform a non self-consistent calculation to calculate the Kohn-Sham structure in a set of 19 k-points in the Irreducible Brillouin Zone. This set of k-points is also derived from a 4x4x4 FCC grid, but a NON-SHIFTED one. It has the same density of points as the 10 k-point set, but the symmetries are not used in a very efficient way. However, this set contains the Gamma point, which allows us to tackle the computation of the band gap at this point.

In dataset 3 we calculate the screening. The screening calculation is very time-consuming. So, we have decided to weaken a bit the parameters found in the previous convergence studies. Indeed, ecutwfn has been decreased from 4.0 to 3.6 . This is rather innocuous. Also, nband has been decreased from 100 to 25. This is a drastic change. The CPU time of this part is linear with respect to this paramater (or more exacly, wrt the number of conduction bands). Thus, the CPU time has been decreased by a factor of 4. Referring to our previous convergence study, we see that the absolute accuracy on the GW energies is now on the order of 0.2 eV only. However, the gap energy (difference between valence and conduction states) is likely correct within 0.02 eV.

Finally in dataset 4 we calculate the self-energy matrix element in Gamma, using the previously determined parameters.

You should obtain the following results:

 k =    0.000   0.000   0.000
 Band     E0  VxcLDA    SigX SigC(E0)      Z dSigC/dE  Sig(E)    E-E0       E
    4   5.915 -11.255 -12.425   0.861   0.771  -0.296 -11.493  -0.238   5.677
    5   8.445 -10.067  -5.858  -3.690   0.772  -0.296  -9.666   0.401   8.846

 E^0_gap          2.530
 E^GW_gap         3.169
 DeltaE^GW_gap    0.639

So that the LDA energy gap in Gamma is about 2.53eV, while the GW correction is about 0.64eV, so that the GW band gap found is 3.17eV.

One can compare now what have been obtained to what one can get from the litterature.

 EXP         3.40 eV   Landolt-Boernstein	

 LDA         2.57 eV   L. Hedin, Phys. Rev. 139, A796 (1965)
 LDA         2.57 eV   M.S. Hybertsen and S. Louie, PRL 55, 1418 (1985)
 LDA (FLAPW) 2.55 eV   N. Hamada, M. Hwang and A.J. Freeman, PRB 41, 3620 (1990)
 LDA (PAW)   2.53 eV   B. Arnaud and M. Alouani, PRB 62, 4464 (2000)
 LDA         2.53 eV   present work

 GW          3.27 eV   M.S. Hybertsen and S. Louie, PRL 55, 1418 (1985)
 GW          3.35 eV   M.S. Hybertsen and S. Louie, PRB 34, 5390 (1986)
 GW          3.30 eV   R.W. Godby, M. Schlueter, L.J. Sham, PRB 37, 10159 (1988)
 GW  (FLAPW) 3.30 eV   N. Hamada, M. Hwang and A.J. Freeman, PRB 41, 3620 (1990)
 GW  (PAW)   3.15 eV   B. Arnaud and M. Alouani, PRB 62, 4464 (2000)
 GW  (FLAPW) 3.12 eV   W. Ku and A.G. Eguiluz, PRL 89, 126401 (2002)
 GW          3.17 eV   present work

The values are spread over an interval of 0.2eV. They depend on the details of the calculation. In the case of pseudopotential calculations, They depend of course on the pseudopotential used. However, a GW result is hardly meaningful beyond 0.1 eV, in the present state of the art.