First tutorial on GW¶
The quasiparticle band structure of Silicon in the GW approximation.¶
This tutorial aims at showing how to calculate selfenergy corrections to the DFT KohnSham (KS) eigenvalues in the GW approximation.
A brief description of the formalism and of the equations implemented in the code can be found in the GW_notes. The different formulas of the GW formalism have been written in a pdf document by Valerio Olevano who also wrote the first version of this tutorial. For a much more consistent discussion of the theoretical aspects of the GW method we refer the reader to the review article Quasiparticle calculations in solids by W.G Aulbur et al also available here.
It is suggested to acknowledge the efforts of developers of the GW part of ABINIT, by citing the 2005 ABINIT publication.
The user should be familiarized with the four basic tutorials of ABINIT, see the tutorial home page. After this first tutorial on GW, you should read the second GW tutorial.
This tutorial should take about 2 hours.
Note
Supposing you made your own installation of ABINIT, the input files to run the examples are in the ~abinit/tests/ directory where ~abinit is the absolute path of the abinit toplevel directory. If you have NOT made your own install, ask your system administrator where to find the package, especially the executable and test files.
In case you work on your own PC or workstation, to make things easier, we suggest you define some handy environment variables by executing the following lines in the terminal:
export ABI_HOME=Replace_with_absolute_path_to_abinit_top_level_dir # Change this line
export PATH=$ABI_HOME/src/98_main/:$PATH # Do not change this line: path to executable
export ABI_TESTS=$ABI_HOME/tests/ # Do not change this line: path to tests dir
export ABI_PSPDIR=$ABI_TESTS/Psps_for_tests/ # Do not change this line: path to pseudos dir
Examples in this tutorial use these shell variables: copy and paste
the code snippets into the terminal (remember to set ABI_HOME first!) or, alternatively,
source the set_abienv.sh
script located in the ~abinit directory:
source ~abinit/set_abienv.sh
The ‘export PATH’ line adds the directory containing the executables to your PATH so that you can invoke the code by simply typing abinit in the terminal instead of providing the absolute path.
To execute the tutorials, create a working directory (Work*
) and
copy there the input files of the lesson.
Most of the tutorials do not rely on parallelism (except specific tutorials on parallelism). However you can run most of the tutorial examples in parallel with MPI, see the topic on parallelism.
1 General example of an almost converged GW calculation¶
Before beginning, you might consider to work in a different subdirectory as for the other tutorials. Why not Work_gw1?
At the end of tutorial 3, we computed the KS band structure of silicon. In this approximation, the band dispersion as well as the band widths are reasonable but the band gaps are qualitatively wrong. Now we will compute the band gaps much more accurately, using the socalled GW approximation.
We start by an example, in which we show how to perform in a single input file the calculation of the ground state density, the Kohn Sham band structure, the screening, and the GW corrections. We use reasonable values for the parameters of the calculation. The discussion on the convergence tests is postponed to the next paragraphs. We will see that GW calculations are much more timeconsuming than the computation of the KS eigenvalues.
So, let us run immediately this calculation, and while it is running, we will explain what has been done.
mkdir Work_gw1
cd Work_gw1
cp $ABI_TESTS/tutorial/Input/tgw1_1.abi .
Then, issue:
abinit tgw1_1.abi > log 2> err &
Please run this job in background because it takes about 1 minute. In the meantime, you should read the following.
1.a The four steps of a GW calculation.¶
In order to perform a standard oneshot GW calculation one has to:

Run a converged Ground State calculation to obtain the selfconsistent density.

Perform a non selfconsistent run to compute the KS eigenvalues and the eigenfunctions including several empty states. Note that, unlike standard band structure calculations, here the KS states must be computed on a regular grid of kpoints. (This limitation is also present with hybrid functional calculations).

Use optdriver = 3 to compute the independentparticle susceptibility \chi^0 on a regular grid of qpoints, for at least two frequencies (usually, \omega=0 and a purely imaginary frequency  of the order of the plasmon frequency, a dozen of eV). The inverse dielectric matrix \epsilon^{1} is then obtained via matrix inversion and stored in an external file (SCR). The list of qpoints is automatically defined by the kmesh used to generate the KS states in the previous step.

Use optdriver = 4 to compute the selfenergy \Sigma matrix elements for a given set of kpoints in order to obtain the GW quasiparticle energies. Note that the kpoint must belong to the kmesh used to generate the WFK file in step 2.
The flowchart diagram of a standard oneshot run is depicted in the figure below.
The input file tgw1_1.abi has precisely that structure: there are four datasets.
The first dataset performs the SCF calculation to get the density. The second dataset reads the previous density file and performs a NSCF run including several empty states. The third dataset reads the WFK file produced in the previous step and drives the computation of susceptibility and dielectric matrices, producing another specialized file, tgw1_xo_DS2_SCR (_SCR for “Screening”, actually the inverse dielectric matrix \epsilon^{1}). Then, in the fourth dataset, the code calculates the quasiparticle energies for the 4^{th} and 5^{th} bands at the \Gamma point.
So, you can edit this tgw1_1.abi file.
# Crystalline silicon # Calculation of the GW corrections # Dataset 1: ground state calculation to get the density # Dataset 2: NSCF run to produce the WFK file for 10 kpoints in IBZ # Dataset 3: calculation of the screening (epsilon^1 matrix for W) # Dataset 4: calculation of the SelfEnergy matrix elements (GW corrections) ndtset 4 ############ # Dataset 1 ############ # SCFGS run nband1 6 tolvrs1 1.0e10 ############ # Dataset 2 ############ # Definition of parameters for the calculation of the WFK file nband2 100 # Number of (occ and empty) bands to be computed nbdbuf2 20 # Do not apply the convergence criterium to the last 20 bands (faster) iscf2 2 getden2 1 tolwfr2 1.0d12 # Will stop when this tolerance is achieved ############ # Dataset 3 ############ # Calculation of the screening (epsilon^1 matrix) optdriver3 3 # Screening calculation getwfk3 1 # Obtain WFK file from previous dataset nband3 60 # Bands to be used in the screening calculation ecuteps3 3.6 # Cutoff energy of the planewave set to represent the dielectric matrix. # It is important to adjust this parameter, that is usually smaller than ecut, # and between 5 and 10 Ha.. ppmfrq3 16.7 eV # Imaginary frequency where to calculate the screening. # It is easier (and safer) to let ABINIT compute and use the Drude plasma frequency, # instead of selecting a value by hand. This would be done thanks to the default value ppmfrq 0.0 . ############ # Dataset 4 ############ # Calculation of the SelfEnergy matrix elements (GW corrections) optdriver4 4 # SelfEnergy calculation getwfk4 2 # Obtain WFK file from dataset 1 getscr4 1 # Obtain SCR file from previous dataset nband4 80 # Bands to be used in the SelfEnergy calculation ecutsigx4 8.0 # Dimension of the G sum in Sigma_x. # ecutsigx = ecut is usually a wise choice # (the dimension in Sigma_c is controlled by ecuteps) nkptgw4 1 # number of kpoint where to calculate the GW correction kptgw4 # kpoints in reduced coordinates 0.000 0.000 0.000 bdgw4 4 5 # calculate GW corrections for bands from 4 to 5 # Data common to the three different datasets # Definition of the unit cell: fcc acell 3*10.26 # Experimental lattice constants rprim 0.0 0.5 0.5 # FCC primitive vectors (to be scaled by acell) 0.5 0.0 0.5 0.5 0.5 0.0 # Definition of the atom types ntypat 1 # There is only one type of atom znucl 14 # The keyword "znucl" refers to the atomic number of the # possible type(s) of atom. The pseudopotential(s) # mentioned in the "files" file must correspond # to the type(s) of atom. Here, the only type is Silicon. # Definition of the atoms natom 2 # There are two atoms typat 1 1 # They both are of type 1, that is, Silicon. xred # Reduced coordinate of atoms 0.0 0.0 0.0 0.25 0.25 0.25 # Definition of the kpoint grid ngkpt 2 2 2 nshiftk 4 shiftk 0.0 0.0 0.0 # These shifts will be the same for all grids 0.0 0.5 0.5 0.5 0.0 0.5 0.5 0.5 0.0 istwfk *1 # This is mandatory in all the GW steps. # Definition of the planewave basis set (at convergence 16 Rydberg 8 Hartree) ecut 8.0 # Maximal kinetic energy cutoff, in Hartree # Definition of the SCF procedure nstep 20 # Maximal number of SCF cycles diemac 12.0 # Although this is not mandatory, it is worth to # precondition the SCF cycle. The model dielectric # function used as the standard preconditioner # is described in the "dielng" input variable section. # Here, we follow the prescription for bulk silicon. pp_dirpath "$ABI_PSPDIR" pseudos "Pseudodojo_nc_sr_04_pbe_standard_psp8/Si.psp8" ############################################################## # This section is used only for regression testing of ABINIT # ############################################################## #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% [files] #%% files_to_test = #%% tgw1_1.abo, tolnlines= 10, tolabs= 0.03, tolrel= 1.500e01 #%% [paral_info] #%% max_nprocs = 4 #%% [extra_info] #%% authors = V. Olevano, F. Bruneval, M. Giantomassi #%% keywords = GW #%% description = #%% Crystalline silicon #%% Calculation of the GW corrections #%% Dataset 1: ground state calculation to get the density #%% Dataset 2: NSCF run to produce the WFK file for 10 kpoints in IBZ #%% Dataset 3: calculation of the screening (epsilon^1 matrix for W) #%% Dataset 4: calculation of the SelfEnergy matrix elements (GW corrections) #%%<END TEST_INFO>
In the first half of the file, you will find specialized input variables for the datasets 1 to 4.
In the second half of the file, one find the datasetindependent information, namely, input variables describing the cell, atom types, number, position, planewave cutoff energy, SCF convergence parameters driving the KS band structure calculation.
1.b Generating the KohnSham band structure: the WFK file.¶
Dataset 1 drives a rather standard SCF calculation. It is worth noticing that we use tolvrs to stop the SCF cycle because we want a wellconverged KS potential to be used in the subsequent NSCF calculation. Dataset 2 computes 100 bands and we set nbdbuf to 20 so that only the first 80 states must be converged within tolwfr. The 20 highest energy states are simply not considered when checking the convergence.
###########
# Dataset 1
############
# SCFGS run
nband1 6
tolvrs1 1.0e10
############
# Dataset 2
############
# Definition of parameters for the calculation of the WFK file
nband2 100 # Number of (occ and empty) bands to be computed
nbdbuf2 20 # Do not apply the convergence criterium to the last 20 bands (faster)
iscf2 2
getden2 1
tolwfr2 1.0d12 # Will stop when this tolerance is achieved
Important
The nbdbuf trick allows us to save several minimization steps because the last bands usually require more iterations to converge in the iterative diagonalization algorithms. Also note that it is a very good idea to increase significantly the value of nbdbuf when computing many empty states. As a rule of thumb, use 10% of nband or even more in complicated systems. This can really make a huge difference at the level of the wall time.
1.c Generating the screening: the SCR file.¶
In dataset 3, the calculation of the screening (KS susceptibility \chi^0 and then inverse dielectric matrix \epsilon^{1}) is performed. We need to set optdriver=3 to do that:
optdriver3 3 # Screening calculation
The getwfk input variable is similar to other “get” input variables of ABINIT:
getwfk3 1 # Obtain WFK file from previous dataset
In this case, it tells the code to use the WFK file calculated in the previous dataset.
Then, three input variables describe the computation:
nband3 60 # Bands used in the screening calculation
ecuteps3 3.6 # Cutoff energy of the planewave set to represent the dielectric matrix
In this case, we use 60 bands to calculate the KS response function \chi^{0}. The dimension of \chi^{0}, as well as all the other matrices (\chi, \epsilon^{1}) is determined by the cutoff energy ecuteps = 3.6 Hartree, which yields 169 planewaves in our case.
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 ppmfrq
ppmfrq3 16.7 eV # Imaginary frequency where to calculate the screening
The two frequencies are used to calculate the plasmonpole model parameters. For the nonzero frequency, it is recommended to use a value close to the plasmon frequency for the plasmonpole 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 Electron Energy Loss Function (\mathrm{Im} \epsilon^{1}_{00}) spectra, so that one can start up by using indications from EELS calculations existing in literature. Alternatively, ABINIT can compute an approximate plasmon frequency using the Drude formula. This is activated by letting ppmfrq to its default value. It is actually safer to use the Drude value than to use blindly a value like 16.7 eV for other materials than silicon.
1.d Computing the GW energies.¶
In dataset 4 the calculation of the SelfEnergy matrix elements is performed. One needs to define the driver option as well as the _WFK and _SCR files.
optdriver4 4 # SelfEnergy calculation
getwfk4 2 # Obtain WFK file from dataset 2
getscr4 1 # Obtain SCR file from previous dataset
The getscr input variable is similar to other “get” input variables of ABINIT.
Then, comes the definition of parameters needed to compute the selfenergy. As for the computation of the susceptibility and dielectric matrices, one must define the set of bands and two sets of planewaves:
nband4 80 # Bands to be used in the SelfEnergy calculation
ecutsigx4 8.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 correlation part of the SelfEnergy while ecutsigx gives the number of planewaves used to calculate \Sigma_x (the exchange part of the selfenergy). The size of the planewave set used to compute \Sigma_c (the correlation part of the selfenergy) is controlled by ecuteps and cannot be larger than the value used to generate the SCR file. For the initial convergence studies, it is advised to set ecutsigx to a value as high as ecut since, anyway, this parameter is not much influential on the total computational time. Note that the exact treatment of the exchange part requires, in principle, ecutsigx = 4 * ecut.
Then, come the parameters defining the kpoints and the band indices for which the quasiparticle energies will be computed:
nkptgw4 1 # number of kpoint where to calculate the GW correction
kptgw4 0.00 0.00 0.00 # kpoints
bdgw4 4 5 # calculate GW corrections for bands from 4 to 5
nkptgw defines the number of kpoints for which the GW corrections will be computed. The kpoint reduced coordinates are specified in kptgw while bdgw gives the minimum/maximum band whose energies are calculated for each selected kpoint.
Important
These kpoints must belong to the kmesh used to generate the WFK file. Hence if you wish the GW correction in a particular kpoint, you should choose a grid containing it. Usually this is done by taking the kpoint grid where the convergence is achieved and shifting it such as at least one kpoint is placed on the wished position in the Brillouin zone.
There is an additional parameter, called zcut, (not studied here) related to the selfenergy computation. It is meant to avoid some divergences that might occur in the calculation due to integrable poles along the integration path.
1.e Examination of the output file.¶
Your calculation should have ended now. Let’s examine the output file. Open tgw1_1.abo in your preferred editor and find the section corresponding to DATASET 3.
.Version 9.7.2 of ABINIT .(MPI version, prepared for a x86_64_linux_gnu9.3 computer) .Copyright (C) 19982024 ABINIT group . ABINIT comes with ABSOLUTELY NO WARRANTY. It is free software, and you are welcome to redistribute it under certain conditions (GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt). ABINIT is a project of the Universite Catholique de Louvain, Corning Inc. and other collaborators, see ~abinit/doc/developers/contributors.txt . Please read https://docs.abinit.org/theory/acknowledgments for suggested acknowledgments of the ABINIT effort. For more information, see https://www.abinit.org . .Starting date : Thu 10 Feb 2022.  ( at 14h13 )  input file > /home/buildbot/ABINIT/alps_gnu_9.3_openmpi/trunk__beuken/tests/TestBot_MPI1/tutorial_tgw1_1/tgw1_1.abi  output file > tgw1_1.abo  root for input files > tgw1_1i  root for output files > tgw1_1o DATASET 1 : space group Fd 3 m (#227); Bravais cF (facecenter cubic) ================================================================================ Values of the parameters that define the memory need for DATASET 1. intxc = 0 ionmov = 0 iscf = 7 lmnmax = 6 lnmax = 6 mgfft = 20 mpssoang = 3 mqgrid = 3001 natom = 2 nloc_mem = 1 nspden = 1 nspinor = 1 nsppol = 1 nsym = 48 n1xccc = 2501 ntypat = 1 occopt = 1 xclevel = 2  mband = 6 mffmem = 1 mkmem = 6 mpw = 303 nfft = 8000 nkpt = 6 ================================================================================ P This job should need less than 3.487 Mbytes of memory. Rough estimation (10% accuracy) of disk space for files : _ WF disk file : 0.168 Mbytes ; DEN or POT disk file : 0.063 Mbytes. ================================================================================ DATASET 2 : space group Fd 3 m (#227); Bravais cF (facecenter cubic) ================================================================================ Values of the parameters that define the memory need for DATASET 2. intxc = 0 ionmov = 0 iscf = 2 lmnmax = 6 lnmax = 6 mgfft = 20 mpssoang = 3 mqgrid = 3001 natom = 2 nloc_mem = 1 nspden = 1 nspinor = 1 nsppol = 1 nsym = 48 n1xccc = 2501 ntypat = 1 occopt = 1 xclevel = 2  mband = 100 mffmem = 1 mkmem = 6 mpw = 303 nfft = 8000 nkpt = 6 ================================================================================ P This job should need less than 5.406 Mbytes of memory. Rough estimation (10% accuracy) of disk space for files : _ WF disk file : 2.776 Mbytes ; DEN or POT disk file : 0.063 Mbytes. ================================================================================ DATASET 3 : space group Fd 3 m (#227); Bravais cF (facecenter cubic) ================================================================================ Values of the parameters that define the memory need for DATASET 3. intxc = 0 ionmov = 0 iscf = 7 lmnmax = 6 lnmax = 6 mgfft = 20 mpssoang = 3 mqgrid = 3001 natom = 2 nloc_mem = 1 nspden = 1 nspinor = 1 nsppol = 1 nsym = 48 n1xccc = 2501 ntypat = 1 occopt = 1 xclevel = 2  mband = 60 mffmem = 1 mkmem = 6 mpw = 303 nfft = 8000 nkpt = 6 ================================================================================ P This job should need less than 5.083 Mbytes of memory. Rough estimation (10% accuracy) of disk space for files : _ WF disk file : 1.666 Mbytes ; DEN or POT disk file : 0.063 Mbytes. ================================================================================ DATASET 4 : space group Fd 3 m (#227); Bravais cF (facecenter cubic) ================================================================================ Values of the parameters that define the memory need for DATASET 4. intxc = 0 ionmov = 0 iscf = 7 lmnmax = 6 lnmax = 6 mgfft = 20 mpssoang = 3 mqgrid = 3001 natom = 2 nloc_mem = 1 nspden = 1 nspinor = 1 nsppol = 1 nsym = 48 n1xccc = 2501 ntypat = 1 occopt = 1 xclevel = 2  mband = 80 mffmem = 1 mkmem = 6 mpw = 303 nfft = 8000 nkpt = 6 ================================================================================ P This job should need less than 5.724 Mbytes of memory. Rough estimation (10% accuracy) of disk space for files : _ WF disk file : 2.221 Mbytes ; DEN or POT disk file : 0.063 Mbytes. ================================================================================   Echo of variables that govern the present computation     outvars: echo of selected default values  iomode0 = 0 , fftalg0 =312 , wfoptalg0 = 0   outvars: echo of global parameters not present in the input file  max_nthreads = 0  outvars: echo values of preprocessed input variables  acell 1.0260000000E+01 1.0260000000E+01 1.0260000000E+01 Bohr amu 2.80855000E+01 bdgw4 4 5 diemac 1.20000000E+01 ecut 8.00000000E+00 Hartree ecuteps1 0.00000000E+00 Hartree ecuteps2 0.00000000E+00 Hartree ecuteps3 3.60000000E+00 Hartree ecuteps4 0.00000000E+00 Hartree ecutsigx1 0.00000000E+00 Hartree ecutsigx2 0.00000000E+00 Hartree ecutsigx3 0.00000000E+00 Hartree ecutsigx4 8.00000000E+00 Hartree ecutwfn1 0.00000000E+00 Hartree ecutwfn2 0.00000000E+00 Hartree ecutwfn3 8.00000000E+00 Hartree ecutwfn4 8.00000000E+00 Hartree  fftalg 312 getden1 0 getden2 1 getden3 0 getden4 0 getscr1 0 getscr2 0 getscr3 0 getscr4 1 getwfk1 0 getwfk2 0 getwfk3 1 getwfk4 2 iscf1 7 iscf2 2 iscf3 7 iscf4 7 istwfk 0 0 1 0 1 1 ixc 11 jdtset 1 2 3 4 kpt 2.50000000E01 2.50000000E01 0.00000000E+00 2.50000000E01 2.50000000E01 0.00000000E+00 5.00000000E01 5.00000000E01 0.00000000E+00 2.50000000E01 5.00000000E01 2.50000000E01 5.00000000E01 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 kptgw4 0.00000000E+00 0.00000000E+00 0.00000000E+00 kptrlatt 2 2 2 2 2 2 2 2 2 kptrlen 2.05200000E+01 P mkmem 6 natom 2 nband1 6 nband2 100 nband3 60 nband4 80 nbdbuf1 0 nbdbuf2 20 nbdbuf3 0 nbdbuf4 0 ndtset 4 ngfft 20 20 20 nkpt 6 nkptgw1 0 nkptgw2 0 nkptgw3 0 nkptgw4 1 npweps1 0 npweps2 0 npweps3 89 npweps4 0 npwsigx1 0 npwsigx2 0 npwsigx3 0 npwsigx4 283 npwwfn1 0 npwwfn2 0 npwwfn3 283 npwwfn4 283 nstep 20 nsym 48 ntypat 1 occ1 2.000000 2.000000 2.000000 2.000000 0.000000 0.000000 occ3 2.000000 2.000000 2.000000 2.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 occ4 2.000000 2.000000 2.000000 2.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 optdriver1 0 optdriver2 0 optdriver3 3 optdriver4 4 ppmfrq1 0.00000000E+00 Hartree ppmfrq2 0.00000000E+00 Hartree ppmfrq3 6.13713734E01 Hartree ppmfrq4 0.00000000E+00 Hartree rprim 0.0000000000E+00 5.0000000000E01 5.0000000000E01 5.0000000000E01 0.0000000000E+00 5.0000000000E01 5.0000000000E01 5.0000000000E01 0.0000000000E+00 spgroup 227 symrel 1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 1 0 1 1 0 1 0 1 1 0 0 1 1 0 1 0 1 1 0 1 0 0 1 0 1 1 1 0 1 0 0 1 0 1 1 1 0 0 1 1 1 0 1 0 0 1 0 1 1 1 0 1 0 0 1 1 0 0 1 1 0 1 0 1 1 0 0 1 1 0 1 0 1 0 1 1 1 1 0 0 1 0 0 1 1 1 1 0 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 0 0 1 0 1 0 0 1 1 0 0 1 1 0 1 0 1 1 0 0 1 1 0 1 1 0 1 1 1 0 1 0 0 1 0 1 1 1 0 1 0 0 0 1 0 1 1 0 0 1 1 0 1 0 1 1 0 0 1 1 1 0 1 0 0 1 0 1 1 1 0 1 0 0 1 0 1 1 0 1 0 0 0 1 1 0 0 0 1 0 0 0 1 1 0 0 1 0 1 0 1 1 0 0 1 1 0 1 0 1 1 0 0 1 0 1 0 0 1 1 1 1 0 0 1 0 0 1 1 1 1 0 1 0 1 1 0 0 1 1 0 1 0 1 1 0 0 1 1 0 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 0 0 1 0 0 1 0 1 1 1 0 1 0 0 1 0 1 1 1 0 1 1 1 0 0 1 1 0 1 0 1 1 0 0 1 1 0 1 0 0 0 1 1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 1 1 0 1 0 0 1 0 1 1 1 0 1 0 0 1 0 1 0 0 1 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 1 1 0 0 1 0 0 1 1 1 1 0 0 1 0 0 1 1 0 0 1 1 0 1 0 1 1 0 0 1 1 0 1 0 1 1 1 1 0 1 0 1 1 0 0 1 1 0 1 0 1 1 0 0 tnons 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 tolvrs1 1.00000000E10 tolvrs2 0.00000000E+00 tolvrs3 0.00000000E+00 tolvrs4 0.00000000E+00 tolwfr1 0.00000000E+00 tolwfr2 1.00000000E12 tolwfr3 0.00000000E+00 tolwfr4 0.00000000E+00 typat 1 1 wtk 0.18750 0.37500 0.09375 0.18750 0.12500 0.03125 xangst 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 1.3573395400E+00 1.3573395400E+00 1.3573395400E+00 xcart 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 2.5650000000E+00 2.5650000000E+00 2.5650000000E+00 xred 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 2.5000000000E01 2.5000000000E01 2.5000000000E01 znucl 14.00000 ================================================================================ chkinp: Checking input parameters for consistency, jdtset= 1. chkinp: Checking input parameters for consistency, jdtset= 2. chkinp: Checking input parameters for consistency, jdtset= 3. chkinp: Checking input parameters for consistency, jdtset= 4. ================================================================================ == DATASET 1 ==================================================================  mpi_nproc: 1, omp_nthreads: 1 (1 if OMP is not activated)  !DatasetInfo iteration_state: {dtset: 1, } dimensions: {natom: 2, nkpt: 6, mband: 6, nsppol: 1, nspinor: 1, nspden: 1, mpw: 303, } cutoff_energies: {ecut: 8.0, pawecutdg: 1.0, } electrons: {nelect: 8.00000000E+00, charge: 0.00000000E+00, occopt: 1.00000000E+00, tsmear: 1.00000000E02, } meta: {optdriver: 0, ionmov: 0, optcell: 0, iscf: 7, paral_kgb: 0, } ... Exchangecorrelation functional for the present dataset will be: GGA: PerdewBurkeErnzerhof functional  ixc=11 Citation for XC functional: J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996) Real(R)+Recip(G) space primitive vectors, cartesian coordinates (Bohr,Bohr^1): R(1)= 0.0000000 5.1300000 5.1300000 G(1)= 0.0974659 0.0974659 0.0974659 R(2)= 5.1300000 0.0000000 5.1300000 G(2)= 0.0974659 0.0974659 0.0974659 R(3)= 5.1300000 5.1300000 0.0000000 G(3)= 0.0974659 0.0974659 0.0974659 Unit cell volume ucvol= 2.7001139E+02 bohr^3 Angles (23,13,12)= 6.00000000E+01 6.00000000E+01 6.00000000E+01 degrees getcut: wavevector= 0.0000 0.0000 0.0000 ngfft= 20 20 20 ecut(hartree)= 8.000 => boxcut(ratio)= 2.16515  Pseudopotential description   pspini: atom type 1 psp file is /home/buildbot/ABINIT/alps_gnu_9.3_openmpi/trunk__beuken/tests/Psps_for_tests/Pseudodojo_nc_sr_04_pbe_standard_psp8/Si.psp8  pspatm: opening atomic psp file /home/buildbot/ABINIT/alps_gnu_9.3_openmpi/trunk__beuken/tests/Psps_for_tests/Pseudodojo_nc_sr_04_pbe_standard_psp8/Si.psp8  Si ONCVPSP3.2.3.1 r_core= 1.60303 1.72197 1.91712  14.00000 4.00000 170510 znucl, zion, pspdat 8 11 2 4 600 0.00000 pspcod,pspxc,lmax,lloc,mmax,r2well 5.99000000000000 4.00000000000000 0.00000000000000 rchrg,fchrg,qchrg nproj 2 2 2 extension_switch 1 pspatm : epsatm= 9.34321699  l ekb(1:nproj) > 0 5.168965 0.829883 1 2.571282 0.578307 2 2.427311 0.488097 pspatm: atomic psp has been read and splines computed 1.49491472E+02 ecore*ucvol(ha*bohr**3)  _setup2: Arith. and geom. avg. npw (full set) are 291.094 290.895 ================================================================================  !BeginCycle iteration_state: {dtset: 1, } solver: {iscf: 7, nstep: 20, nline: 4, wfoptalg: 0, } tolerances: {tolvrs: 1.00E10, } ... iter Etot(hartree) deltaE(h) residm vres2 ETOT 1 8.4477770574543 8.448E+00 6.039E03 2.180E+00 ETOT 2 8.4508649767972 3.088E03 7.993E04 3.375E02 ETOT 3 8.4508905681654 2.559E05 4.495E04 3.880E04 ETOT 4 8.4508907628105 1.946E07 3.877E04 1.445E05 ETOT 5 8.4508907674908 4.680E09 6.981E05 4.700E07 ETOT 6 8.4508907676238 1.330E10 8.334E05 2.407E09 ETOT 7 8.4508907676242 3.588E13 1.045E05 1.231E11 At SCF step 7 vres2 = 1.23E11 < tolvrs= 1.00E10 =>converged. Cartesian components of stress tensor (hartree/bohr^3) sigma(1 1)= 4.48073872E05 sigma(3 2)= 0.00000000E+00 sigma(2 2)= 4.48073872E05 sigma(3 1)= 0.00000000E+00 sigma(3 3)= 4.48073872E05 sigma(2 1)= 0.00000000E+00  !ResultsGS iteration_state: {dtset: 1, } comment : Summary of ground state results lattice_vectors:  [ 0.0000000, 5.1300000, 5.1300000, ]  [ 5.1300000, 0.0000000, 5.1300000, ]  [ 5.1300000, 5.1300000, 0.0000000, ] lattice_lengths: [ 7.25492, 7.25492, 7.25492, ] lattice_angles: [ 60.000, 60.000, 60.000, ] # degrees, (23, 13, 12) lattice_volume: 2.7001139E+02 convergence: {deltae: 3.588E13, res2: 1.231E11, residm: 1.045E05, diffor: null, } etotal : 8.45089077E+00 entropy : 0.00000000E+00 fermie : 1.62353054E01 cartesian_stress_tensor: # hartree/bohr^3  [ 4.48073872E05, 0.00000000E+00, 0.00000000E+00, ]  [ 0.00000000E+00, 4.48073872E05, 0.00000000E+00, ]  [ 0.00000000E+00, 0.00000000E+00, 4.48073872E05, ] pressure_GPa: 1.3183E+00 xred :  [ 0.0000E+00, 0.0000E+00, 0.0000E+00, Si]  [ 2.5000E01, 2.5000E01, 2.5000E01, Si] cartesian_forces: # hartree/bohr  [ 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, ]  [ 0.00000000E+00, 0.00000000E+00, 0.00000000E+00, ] force_length_stats: {min: 0.00000000E+00, max: 0.00000000E+00, mean: 0.00000000E+00, } ... Integrated electronic density in atomic spheres:  Atom Sphere_radius Integrated_density 1 2.00000 1.72529250 2 2.00000 1.72529250 ================================================================================ iterations are completed or convergence reached Mean square residual over all n,k,spin= 29.053E08; max= 10.447E06 reduced coordinates (array xred) for 2 atoms 0.000000000000 0.000000000000 0.000000000000 0.250000000000 0.250000000000 0.250000000000 rms dE/dt= 2.1985E28; max dE/dt= 2.6926E28; dE/dt below (all hartree) 1 0.000000000000 0.000000000000 0.000000000000 2 0.000000000000 0.000000000000 0.000000000000 cartesian coordinates (angstrom) at end: 1 0.00000000000000 0.00000000000000 0.00000000000000 2 1.35733954003335 1.35733954003335 1.35733954003335 cartesian forces (hartree/bohr) at end: 1 0.00000000000000 0.00000000000000 0.00000000000000 2 0.00000000000000 0.00000000000000 0.00000000000000 frms,max,avg= 0.0000000E+00 0.0000000E+00 0.000E+00 0.000E+00 0.000E+00 h/b cartesian forces (eV/Angstrom) at end: 1 0.00000000000000 0.00000000000000 0.00000000000000 2 0.00000000000000 0.00000000000000 0.00000000000000 frms,max,avg= 0.0000000E+00 0.0000000E+00 0.000E+00 0.000E+00 0.000E+00 e/A length scales= 10.260000000000 10.260000000000 10.260000000000 bohr = 5.429358160133 5.429358160133 5.429358160133 angstroms prteigrs : about to open file tgw1_1o_DS1_EIG Fermi (or HOMO) energy (hartree) = 0.16235 Average Vxc (hartree)= 0.34044 Eigenvalues (hartree) for nkpt= 6 k points: kpt# 1, nband= 6, wtk= 0.18750, kpt= 0.2500 0.2500 0.0000 (reduced coord) 0.23833 0.03292 0.09196 0.09196 0.20352 0.27776 prteigrs : prtvol=0 or 1, do not print more kpoints.  !EnergyTerms iteration_state : {dtset: 1, } comment : Components of total free energy in Hartree kinetic : 3.08753412466338E+00 hartree : 5.59850380435293E01 xc : 3.09781415724427E+00 Ewald energy : 8.40046478618609E+00 psp_core : 5.53648753925927E01 local_psp : 2.30656115212577E+00 non_local_psp : 1.15291606890733E+00 total_energy : 8.45089076762421E+00 total_energy_eV : 2.29960432636752E+02 band_energy : 1.88175119850873E01 ... Cartesian components of stress tensor (hartree/bohr^3) sigma(1 1)= 4.48073872E05 sigma(3 2)= 0.00000000E+00 sigma(2 2)= 4.48073872E05 sigma(3 1)= 0.00000000E+00 sigma(3 3)= 4.48073872E05 sigma(2 1)= 0.00000000E+00 Cartesian components of stress tensor (GPa) [Pressure= 1.3183E+00 GPa]  sigma(1 1)= 1.31827862E+00 sigma(3 2)= 0.00000000E+00  sigma(2 2)= 1.31827862E+00 sigma(3 1)= 0.00000000E+00  sigma(3 3)= 1.31827862E+00 sigma(2 1)= 0.00000000E+00 ================================================================================ == DATASET 2 ==================================================================  mpi_nproc: 1, omp_nthreads: 1 (1 if OMP is not activated)  !DatasetInfo iteration_state: {dtset: 2, } dimensions: {natom: 2, nkpt: 6, mband: 100, nsppol: 1, nspinor: 1, nspden: 1, mpw: 303, } cutoff_energies: {ecut: 8.0, pawecutdg: 1.0, } electrons: {nelect: 8.00000000E+00, charge: 0.00000000E+00, occopt: 1.00000000E+00, tsmear: 1.00000000E02, } meta: {optdriver: 0, ionmov: 0, optcell: 0, iscf: 2, paral_kgb: 0, } ... mkfilename : getden/=0, take file _DEN from output of DATASET 1. Exchangecorrelation functional for the present dataset will be: GGA: PerdewBurkeErnzerhof functional  ixc=11 Citation for XC functional: J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996) Real(R)+Recip(G) space primitive vectors, cartesian coordinates (Bohr,Bohr^1): R(1)= 0.0000000 5.1300000 5.1300000 G(1)= 0.0974659 0.0974659 0.0974659 R(2)= 5.1300000 0.0000000 5.1300000 G(2)= 0.0974659 0.0974659 0.0974659 R(3)= 5.1300000 5.1300000 0.0000000 G(3)= 0.0974659 0.0974659 0.0974659 Unit cell volume ucvol= 2.7001139E+02 bohr^3 Angles (23,13,12)= 6.00000000E+01 6.00000000E+01 6.00000000E+01 degrees getcut: wavevector= 0.0000 0.0000 0.0000 ngfft= 20 20 20 ecut(hartree)= 8.000 => boxcut(ratio)= 2.16515  ================================================================================ prteigrs : about to open file tgw1_1o_DS2_EIG NonSCF case, kpt 1 ( 0.25000 0.25000 0.00000), residuals and eigenvalues= 2.31E13 4.27E13 6.12E14 6.35E14 2.22E13 8.99E13 1.42E13 3.53E13 2.37E14 4.21E13 3.96E13 3.90E13 5.13E13 1.83E13 3.55E13 6.54E13 5.91E13 1.42E14 2.52E13 2.43E14 4.70E14 5.95E13 3.93E13 3.30E13 4.65E13 1.33E13 4.29E14 2.82E14 2.91E13 4.11E13 3.53E13 7.85E14 6.45E14 9.71E14 6.27E13 6.93E14 3.75E14 1.66E13 1.02E13 3.35E13 7.76E14 1.70E13 2.21E13 9.01E14 1.13E13 2.02E13 1.62E13 2.19E13 3.65E13 1.41E13 2.24E13 6.83E13 8.06E13 8.98E13 5.94E13 4.08E13 4.94E13 7.40E13 8.48E13 5.05E13 7.21E13 7.08E13 6.87E14 1.56E13 4.29E13 3.83E13 2.76E13 5.63E13 7.96E13 4.79E13 5.38E13 9.97E13 1.36E13 8.37E14 9.49E14 3.45E13 2.70E13 1.89E13 7.06E13 3.75E13 2.20E13 8.36E13 1.65E12 2.63E12 3.06E11 3.75E11 1.29E10 6.99E10 2.69E09 7.79E09 4.48E08 9.65E09 1.88E07 2.10E06 7.78E07 1.20E06 2.58E06 1.89E05 2.14E04 2.52E04 2.3833E01 3.2916E02 9.1963E02 9.1963E02 2.0352E01 2.7776E01 3.7458E01 3.7458E01 4.5966E01 4.9437E01 5.8246E01 6.4334E01 6.4334E01 6.6963E01 8.0649E01 8.0649E01 8.4299E01 9.1578E01 9.3904E01 1.0878E+00 1.1459E+00 1.1772E+00 1.1772E+00 1.2524E+00 1.2524E+00 1.2656E+00 1.2843E+00 1.4509E+00 1.4596E+00 1.4770E+00 1.4770E+00 1.5394E+00 1.5426E+00 1.5426E+00 1.6688E+00 1.6688E+00 1.6755E+00 1.6845E+00 1.6928E+00 1.7847E+00 1.7847E+00 1.7911E+00 1.8781E+00 1.9163E+00 1.9163E+00 1.9725E+00 1.9784E+00 2.0346E+00 2.0466E+00 2.0466E+00 2.0795E+00 2.1289E+00 2.2433E+00 2.2787E+00 2.2787E+00 2.3791E+00 2.4811E+00 2.4811E+00 2.4898E+00 2.4932E+00 2.4953E+00 2.4953E+00 2.5176E+00 2.6118E+00 2.6504E+00 2.6703E+00 2.6703E+00 2.7489E+00 2.7600E+00 2.7881E+00 2.8100E+00 2.8100E+00 2.8579E+00 2.8844E+00 2.8844E+00 2.9090E+00 2.9736E+00 2.9915E+00 3.0730E+00 3.0730E+00 3.0827E+00 3.1663E+00 3.1663E+00 3.2167E+00 3.2686E+00 3.3076E+00 3.3076E+00 3.3803E+00 3.4345E+00 3.4345E+00 3.4498E+00 3.4534E+00 3.5229E+00 3.5229E+00 3.5259E+00 3.5404E+00 3.5494E+00 3.5495E+00 3.5903E+00 3.6047E+00 prteigrs : prtvol=0 or 1, do not print more kpoints.  !ResultsGS iteration_state: {dtset: 2, } comment : Summary of ground state results lattice_vectors:  [ 0.0000000, 5.1300000, 5.1300000, ]  [ 5.1300000, 0.0000000, 5.1300000, ]  [ 5.1300000, 5.1300000, 0.0000000, ] lattice_lengths: [ 7.25492, 7.25492, 7.25492, ] lattice_angles: [ 60.000, 60.000, 60.000, ] # degrees, (23, 13, 12) lattice_volume: 2.7001139E+02 convergence: {deltae: 0.000E+00, res2: 0.000E+00, residm: 9.975E13, diffor: 0.000E+00, } etotal : 8.45089077E+00 entropy : 0.00000000E+00 fermie : 1.62353054E01 cartesian_stress_tensor: null pressure_GPa: null xred :  [ 0.0000E+00, 0.0000E+00, 0.0000E+00, Si]  [ 2.5000E01, 2.5000E01, 2.5000E01, Si] cartesian_forces: null force_length_stats: {min: null, max: null, mean: null, } ... Integrated electronic density in atomic spheres:  Atom Sphere_radius Integrated_density 1 2.00000 1.72529250 2 2.00000 1.72529250 ================================================================================ iterations are completed or convergence reached Mean square residual over all n,k,spin= 31.196E14; max= 99.749E14 reduced coordinates (array xred) for 2 atoms 0.000000000000 0.000000000000 0.000000000000 0.250000000000 0.250000000000 0.250000000000 cartesian coordinates (angstrom) at end: 1 0.00000000000000 0.00000000000000 0.00000000000000 2 1.35733954003335 1.35733954003335 1.35733954003335 length scales= 10.260000000000 10.260000000000 10.260000000000 bohr = 5.429358160133 5.429358160133 5.429358160133 angstroms prteigrs : about to open file tgw1_1o_DS2_EIG Eigenvalues (hartree) for nkpt= 6 k points: kpt# 1, nband=100, wtk= 0.18750, kpt= 0.2500 0.2500 0.0000 (reduced coord) 0.23833 0.03292 0.09196 0.09196 0.20352 0.27776 0.37458 0.37458 0.45966 0.49437 0.58246 0.64334 0.64334 0.66963 0.80649 0.80649 0.84299 0.91578 0.93904 1.08776 1.14592 1.17716 1.17716 1.25239 1.25239 1.26564 1.28431 1.45093 1.45964 1.47700 1.47700 1.53936 1.54263 1.54263 1.66881 1.66881 1.67555 1.68455 1.69275 1.78466 1.78466 1.79107 1.87812 1.91631 1.91631 1.97250 1.97836 2.03457 2.04662 2.04662 2.07950 2.12889 2.24331 2.27866 2.27866 2.37907 2.48107 2.48107 2.48980 2.49323 2.49531 2.49531 2.51758 2.61177 2.65043 2.67032 2.67032 2.74895 2.76000 2.78805 2.80998 2.80998 2.85793 2.88444 2.88444 2.90905 2.97358 2.99145 3.07297 3.07297 3.08272 3.16627 3.16627 3.21675 3.26856 3.30759 3.30759 3.38032 3.43452 3.43452 3.44983 3.45344 3.52290 3.52290 3.52587 3.54037 3.54944 3.54946 3.59033 3.60475 prteigrs : prtvol=0 or 1, do not print more kpoints. ================================================================================ == DATASET 3 ==================================================================  mpi_nproc: 1, omp_nthreads: 1 (1 if OMP is not activated)  !DatasetInfo iteration_state: {dtset: 3, } dimensions: {natom: 2, nkpt: 6, mband: 60, nsppol: 1, nspinor: 1, nspden: 1, mpw: 303, } cutoff_energies: {ecut: 8.0, pawecutdg: 1.0, } electrons: {nelect: 8.00000000E+00, charge: 0.00000000E+00, occopt: 1.00000000E+00, tsmear: 1.00000000E02, } meta: {optdriver: 3, gwcalctyp: 0, } ... mkfilename : getwfk/=0, take file _WFK from output of DATASET 2. Exchangecorrelation functional for the present dataset will be: GGA: PerdewBurkeErnzerhof functional  ixc=11 Citation for XC functional: J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996) SCREENING: Calculation of the susceptibility and dielectric matrices Based on a program developped by R.W. Godby, V. Olevano, G. Onida, and L. Reining. Incorporated in ABINIT by V. Olevano, G.M. Rignanese, and M. Torrent. .Using double precision arithmetic ; gwpc = 8 Real(R)+Recip(G) space primitive vectors, cartesian coordinates (Bohr,Bohr^1): R(1)= 0.0000000 5.1300000 5.1300000 G(1)= 0.0974659 0.0974659 0.0974659 R(2)= 5.1300000 0.0000000 5.1300000 G(2)= 0.0974659 0.0974659 0.0974659 R(3)= 5.1300000 5.1300000 0.0000000 G(3)= 0.0974659 0.0974659 0.0974659 Unit cell volume ucvol= 2.7001139E+02 bohr^3 Angles (23,13,12)= 6.00000000E+01 6.00000000E+01 6.00000000E+01 degrees  ==== Kmesh for the wavefunctions ==== Number of points in the irreducible wedge : 6 Reduced coordinates and weights : 1) 2.50000000E01 2.50000000E01 0.00000000E+00 0.18750 2) 2.50000000E01 2.50000000E01 0.00000000E+00 0.37500 3) 5.00000000E01 5.00000000E01 0.00000000E+00 0.09375 4) 2.50000000E01 5.00000000E01 2.50000000E01 0.18750 5) 5.00000000E01 0.00000000E+00 0.00000000E+00 0.12500 6) 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.03125 Together with 48 symmetry operations and timereversal symmetry yields 32 points in the full Brillouin Zone. ==== Qmesh for the screening function ==== Number of points in the irreducible wedge : 6 Reduced coordinates and weights : 1) 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.03125 2) 5.00000000E01 5.00000000E01 0.00000000E+00 0.09375 3) 5.00000000E01 2.50000000E01 2.50000000E01 0.37500 4) 0.00000000E+00 5.00000000E01 0.00000000E+00 0.12500 5) 5.00000000E01 2.50000000E01 2.50000000E01 0.18750 6) 0.00000000E+00 2.50000000E01 2.50000000E01 0.18750 Together with 48 symmetry operations and timereversal symmetry yields 32 points in the full Brillouin Zone. setmesh: FFT mesh size selected = 20x 20x 20 total number of points = 8000  screening: taking advantage of timereversal symmetry  Maximum band index for partially occupied states nbvw = 4  Remaining bands to be divided among processors nbcw = 56  Number of bands treated by each node ~56 Number of electrons calculated from density = 8.0000; Expected = 8.0000 average of density, n = 0.029628 r_s = 2.0048 omega_plasma = 16.6039 [eV] calculating chi0 at frequencies [eV] : 1 0.000000E+00 0.000000E+00 2 0.000000E+00 1.670000E+01  qpoint number 1 q = ( 0.000000, 0.000000, 0.000000) [r.l.u.]  chi0(G,G') at the 1 th omega 0.0000 0.0000 [eV] 1 2 3 4 5 6 7 8 9 1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 2 0.000 18.770 0.000 0.072 0.000 0.072 0.000 0.072 0.000 0.000 0.000 5.064 0.000 0.279 0.000 0.279 0.000 0.279 chi0(G,G') at the 2 th omega 0.0000 16.7000 [eV] 1 2 3 4 5 6 7 8 9 1 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 2 0.000 7.757 0.000 0.039 0.000 0.039 0.000 0.039 0.000 0.000 0.000 1.029 0.000 0.122 0.000 0.122 0.000 0.122 For qpoint: 0.000010 0.000020 0.000030 dielectric constant = 22.4176 dielectric constant without local fields = 24.7005 Average fulfillment of the sum rule on Im[epsilon] for qpoint 1 : 82.54 [%] Heads and wings of the symmetrical epsilon^1(G,G') Upper and lower wings at the 1 th omega 0.0000 0.0000 [eV] 1 2 3 4 5 6 7 8 9 0.045 0.004 0.004 0.012 0.012 0.012 0.012 0.004 0.004 0.000 0.004 0.004 0.012 0.012 0.012 0.012 0.004 0.004 1 2 3 4 5 6 7 8 9 0.045 0.004 0.004 0.012 0.012 0.012 0.012 0.004 0.004 0.000 0.004 0.004 0.012 0.012 0.012 0.012 0.004 0.004 Upper and lower wings at the 2 th omega 0.0000 16.7000 [eV] 1 2 3 4 5 6 7 8 9 0.492 0.007 0.007 0.022 0.022 0.022 0.022 0.007 0.007 0.000 0.007 0.007 0.022 0.022 0.022 0.022 0.007 0.007 1 2 3 4 5 6 7 8 9 0.492 0.007 0.007 0.022 0.022 0.022 0.022 0.007 0.007 0.000 0.007 0.007 0.022 0.022 0.022 0.022 0.007 0.007  qpoint number 2 q = ( 0.500000, 0.500000, 0.000000) [r.l.u.]  chi0(G,G') at the 1 th omega 0.0000 0.0000 [eV] 1 2 3 4 5 6 7 8 9 1 18.001 3.139 1.157 1.157 3.139 1.157 3.139 3.139 1.157 0.000 3.139 1.157 1.157 3.139 1.157 3.139 3.139 1.157 2 3.139 17.803 0.000 0.018 0.000 0.018 0.000 0.152 0.000 3.139 0.000 2.461 0.000 0.491 0.000 0.491 0.000 0.263 chi0(G,G') at the 2 th omega 0.0000 16.7000 [eV] 1 2 3 4 5 6 7 8 9 1 4.372 0.918 0.223 0.223 0.918 0.223 0.918 0.918 0.223 0.000 0.918 0.223 0.223 0.918 0.223 0.918 0.918 0.223 2 0.918 8.743 0.000 0.031 0.000 0.031 0.000 0.138 0.000 0.918 0.000 0.578 0.000 0.033 0.000 0.034 0.000 0.061 Average fulfillment of the sum rule on Im[epsilon] for qpoint 2 : 86.77 [%]  qpoint number 3 q = ( 0.500000, 0.250000, 0.250000) [r.l.u.]  chi0(G,G') at the 1 th omega 0.0000 0.0000 [eV] 1 2 3 4 5 6 7 8 9 1 14.492 2.533 0.521 2.302 2.302 0.521 2.533 2.302 2.302 0.000 2.533 0.521 2.302 2.302 0.521 2.533 2.302 2.302 2 2.533 18.009 0.000 0.299 0.000 0.041 0.000 0.299 0.000 2.533 0.000 3.373 0.000 0.106 0.000 0.334 0.000 0.106 chi0(G,G') at the 2 th omega 0.0000 16.7000 [eV] 1 2 3 4 5 6 7 8 9 1 2.571 0.727 0.255 0.414 0.414 0.255 0.727 0.414 0.414 0.000 0.727 0.255 0.414 0.414 0.255 0.727 0.414 0.414 2 0.727 8.786 0.000 0.050 0.000 0.031 0.000 0.050 0.000 0.727 0.000 0.808 0.000 0.055 0.000 0.047 0.000 0.055 Average fulfillment of the sum rule on Im[epsilon] for qpoint 3 : 87.33 [%]  qpoint number 4 q = ( 0.000000, 0.500000, 0.000000) [r.l.u.]  chi0(G,G') at the 1 th omega 0.0000 0.0000 [eV] 1 2 3 4 5 6 7 8 9 1 14.428 1.961 2.660 1.961 2.660 1.961 2.660 2.274 3.152 0.000 1.961 2.660 1.961 2.660 1.961 2.660 2.274 3.152 2 1.961 20.236 0.000 0.280 0.000 0.280 0.000 0.340 0.000 1.961 0.000 2.681 0.000 0.211 0.000 0.211 0.000 0.897 chi0(G,G') at the 2 th omega 0.0000 16.7000 [eV] 1 2 3 4 5 6 7 8 9 1 3.518 0.343 0.716 0.343 0.716 0.343 0.716 0.910 0.537 0.000 0.343 0.716 0.343 0.716 0.343 0.716 0.910 0.537 2 0.343 7.557 0.000 0.135 0.000 0.135 0.000 0.066 0.000 0.343 0.000 0.648 0.000 0.044 0.000 0.044 0.000 0.102 Average fulfillment of the sum rule on Im[epsilon] for qpoint 4 : 86.99 [%]  qpoint number 5 q = ( 0.500000,0.250000, 0.250000) [r.l.u.]  chi0(G,G') at the 1 th omega 0.0000 0.0000 [eV] 1 2 3 4 5 6 7 8 9 1 19.327 3.001 1.292 3.006 2.417 2.417 3.006 1.292 3.001 0.000 3.001 1.292 3.006 2.417 2.417 3.006 1.292 3.001 2 3.001 16.162 0.000 0.366 0.000 0.549 0.000 0.242 0.000 3.001 0.000 2.479 0.000 0.207 0.000 0.320 0.000 0.063 chi0(G,G') at the 2 th omega 0.0000 16.7000 [eV] 1 2 3 4 5 6 7 8 9 1 5.045 1.050 0.043 0.904 0.625 0.625 0.904 0.043 1.050 0.000 1.050 0.043 0.904 0.625 0.625 0.904 0.043 1.050 2 1.050 8.424 0.000 0.139 0.000 0.077 0.000 0.008 0.000 1.050 0.000 0.510 0.000 0.040 0.000 0.034 0.000 0.063 Average fulfillment of the sum rule on Im[epsilon] for qpoint 5 : 86.85 [%]  qpoint number 6 q = ( 0.000000,0.250000,0.250000) [r.l.u.]  chi0(G,G') at the 1 th omega 0.0000 0.0000 [eV] 1 2 3 4 5 6 7 8 9 1 10.610 2.292 0.014 0.014 2.292 2.292 0.014 0.014 2.292 0.000 2.292 0.014 0.014 2.292 2.292 0.014 0.014 2.292 2 2.292 18.972 0.000 0.413 0.000 0.354 0.000 0.413 0.000 2.292 0.000 3.503 0.000 0.184 0.000 0.374 0.000 0.184 chi0(G,G') at the 2 th omega 0.0000 16.7000 [eV] 1 2 3 4 5 6 7 8 9 1 1.418 0.449 0.089 0.089 0.449 0.449 0.089 0.089 0.449 0.000 0.449 0.089 0.089 0.449 0.449 0.089 0.089 0.449 2 0.449 8.551 0.000 0.048 0.000 0.044 0.000 0.048 0.000 0.449 0.000 0.896 0.000 0.063 0.000 0.100 0.000 0.063 Average fulfillment of the sum rule on Im[epsilon] for qpoint 6 : 88.68 [%] ================================================================================ == DATASET 4 ==================================================================  mpi_nproc: 1, omp_nthreads: 1 (1 if OMP is not activated)  !DatasetInfo iteration_state: {dtset: 4, } dimensions: {natom: 2, nkpt: 6, mband: 80, nsppol: 1, nspinor: 1, nspden: 1, mpw: 303, } cutoff_energies: {ecut: 8.0, pawecutdg: 1.0, } electrons: {nelect: 8.00000000E+00, charge: 0.00000000E+00, occopt: 1.00000000E+00, tsmear: 1.00000000E02, } meta: {optdriver: 4, gwcalctyp: 0, } ... mkfilename : getwfk/=0, take file _WFK from output of DATASET 2. mkfilename : getscr/=0, take file _SCR from output of DATASET 3. Exchangecorrelation functional for the present dataset will be: GGA: PerdewBurkeErnzerhof functional  ixc=11 Citation for XC functional: J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996) SIGMA: Calculation of the GW corrections Based on a program developped by R.W. Godby, V. Olevano, G. Onida, and L. Reining. Incorporated in ABINIT by V. Olevano, G.M. Rignanese, and M. Torrent. .Using double precision arithmetic ; gwpc = 8 Real(R)+Recip(G) space primitive vectors, cartesian coordinates (Bohr,Bohr^1): R(1)= 0.0000000 5.1300000 5.1300000 G(1)= 0.0974659 0.0974659 0.0974659 R(2)= 5.1300000 0.0000000 5.1300000 G(2)= 0.0974659 0.0974659 0.0974659 R(3)= 5.1300000 5.1300000 0.0000000 G(3)= 0.0974659 0.0974659 0.0974659 Unit cell volume ucvol= 2.7001139E+02 bohr^3 Angles (23,13,12)= 6.00000000E+01 6.00000000E+01 6.00000000E+01 degrees  ==== Kmesh for the wavefunctions ==== Number of points in the irreducible wedge : 6 Reduced coordinates and weights : 1) 2.50000000E01 2.50000000E01 0.00000000E+00 0.18750 2) 2.50000000E01 2.50000000E01 0.00000000E+00 0.37500 3) 5.00000000E01 5.00000000E01 0.00000000E+00 0.09375 4) 2.50000000E01 5.00000000E01 2.50000000E01 0.18750 5) 5.00000000E01 0.00000000E+00 0.00000000E+00 0.12500 6) 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.03125 Together with 48 symmetry operations and timereversal symmetry yields 32 points in the full Brillouin Zone. ==== Qmesh for screening function ==== Number of points in the irreducible wedge : 6 Reduced coordinates and weights : 1) 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.03125 2) 5.00000000E01 5.00000000E01 0.00000000E+00 0.09375 3) 5.00000000E01 2.50000000E01 2.50000000E01 0.37500 4) 0.00000000E+00 5.00000000E01 0.00000000E+00 0.12500 5) 5.00000000E01 2.50000000E01 2.50000000E01 0.18750 6) 0.00000000E+00 2.50000000E01 2.50000000E01 0.18750 Together with 48 symmetry operations and timereversal symmetry yields 32 points in the full Brillouin Zone. setmesh: FFT mesh size selected = 20x 20x 20 total number of points = 8000 Number of electrons calculated from density = 8.0000; Expected = 8.0000 average of density, n = 0.029628 r_s = 2.0048 omega_plasma = 16.6039 [eV] === KS Band Gaps === >>>> For spin 1 Minimum direct gap = 2.5435 [eV], located at kpoint : 0.0000 0.0000 0.0000 Fundamental gap = 0.7097 [eV], Top of valence bands at : 0.0000 0.0000 0.0000 Bottom of conduction at : 0.5000 0.5000 0.0000 SIGMA fundamental parameters: PLASMON POLE MODEL 1 number of planewaves for SigmaX 283 number of planewaves for SigmaC and W 89 number of planewaves for wavefunctions 283 number of bands 80 number of independent spin polarizations 1 number of spinorial components 1 number of kpoints in IBZ 6 number of qpoints in IBZ 6 number of symmetry operations 48 number of kpoints in BZ 32 number of qpoints in BZ 32 number of frequencies for dSigma/dE 9 frequency step for dSigma/dE [eV] 0.25 number of omega for Sigma on real axis 0 max omega for Sigma on real axis [eV] 0.00 zcut for avoiding poles [eV] 0.10 EPSILON^1 parameters (SCR file): dimension of the eps^1 matrix on file 89 dimension of the eps^1 matrix used 89 number of planewaves for wavefunctions 283 number of bands 60 number of qpoints in IBZ 6 number of frequencies 2 number of real frequencies 1 number of imag frequencies 1 matrix elements of selfenergy operator (all in [eV]) Perturbative Calculation  !SelfEnergy_ee iteration_state: {dtset: 4, } kpoint : [ 0.000, 0.000, 0.000, ] spin : 1 KS_gap : 2.544 QP_gap : 3.110 Delta_QP_KS: 0.567 data: !SigmaeeData  Band E0 <VxcDFT> SigX SigC(E0) Z dSigC/dE Sig(E) EE0 E 2 4.418 11.332 13.262 1.352 0.766 0.305 11.775 0.443 3.975 3 4.418 11.332 13.262 1.352 0.766 0.305 11.775 0.443 3.975 4 4.418 11.332 13.262 1.352 0.766 0.305 11.775 0.443 3.975 5 6.961 10.028 5.550 4.316 0.766 0.305 9.904 0.124 7.085 6 6.961 10.028 5.550 4.316 0.766 0.305 9.904 0.124 7.085 7 6.961 10.028 5.550 4.316 0.766 0.305 9.904 0.124 7.085 ... == END DATASET(S) ============================================================== ================================================================================ outvars: echo values of variables after computation  acell 1.0260000000E+01 1.0260000000E+01 1.0260000000E+01 Bohr amu 2.80855000E+01 bdgw4 4 5 diemac 1.20000000E+01 ecut 8.00000000E+00 Hartree ecuteps1 0.00000000E+00 Hartree ecuteps2 0.00000000E+00 Hartree ecuteps3 3.60000000E+00 Hartree ecuteps4 0.00000000E+00 Hartree ecutsigx1 0.00000000E+00 Hartree ecutsigx2 0.00000000E+00 Hartree ecutsigx3 0.00000000E+00 Hartree ecutsigx4 8.00000000E+00 Hartree ecutwfn1 0.00000000E+00 Hartree ecutwfn2 0.00000000E+00 Hartree ecutwfn3 8.00000000E+00 Hartree ecutwfn4 8.00000000E+00 Hartree etotal1 8.4508907676E+00 etotal3 0.0000000000E+00 etotal4 0.0000000000E+00 fcart1 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 fcart3 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 fcart4 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00  fftalg 312 getden1 0 getden2 1 getden3 0 getden4 0 getscr1 0 getscr2 0 getscr3 0 getscr4 1 getwfk1 0 getwfk2 0 getwfk3 1 getwfk4 2 iscf1 7 iscf2 2 iscf3 7 iscf4 7 istwfk 0 0 1 0 1 1 ixc 11 jdtset 1 2 3 4 kpt 2.50000000E01 2.50000000E01 0.00000000E+00 2.50000000E01 2.50000000E01 0.00000000E+00 5.00000000E01 5.00000000E01 0.00000000E+00 2.50000000E01 5.00000000E01 2.50000000E01 5.00000000E01 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 kptgw4 0.00000000E+00 0.00000000E+00 0.00000000E+00 kptrlatt 2 2 2 2 2 2 2 2 2 kptrlen 2.05200000E+01 P mkmem 6 natom 2 nband1 6 nband2 100 nband3 60 nband4 80 nbdbuf1 0 nbdbuf2 20 nbdbuf3 0 nbdbuf4 0 ndtset 4 ngfft 20 20 20 nkpt 6 nkptgw1 0 nkptgw2 0 nkptgw3 0 nkptgw4 1 npweps1 0 npweps2 0 npweps3 89 npweps4 0 npwsigx1 0 npwsigx2 0 npwsigx3 0 npwsigx4 283 npwwfn1 0 npwwfn2 0 npwwfn3 283 npwwfn4 283 nstep 20 nsym 48 ntypat 1 occ1 2.000000 2.000000 2.000000 2.000000 0.000000 0.000000 occ3 2.000000 2.000000 2.000000 2.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 occ4 2.000000 2.000000 2.000000 2.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 optdriver1 0 optdriver2 0 optdriver3 3 optdriver4 4 ppmfrq1 0.00000000E+00 Hartree ppmfrq2 0.00000000E+00 Hartree ppmfrq3 6.13713734E01 Hartree ppmfrq4 0.00000000E+00 Hartree rprim 0.0000000000E+00 5.0000000000E01 5.0000000000E01 5.0000000000E01 0.0000000000E+00 5.0000000000E01 5.0000000000E01 5.0000000000E01 0.0000000000E+00 spgroup 227 strten1 4.4807387228E05 4.4807387228E05 4.4807387228E05 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 strten3 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 strten4 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 symrel 1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 1 0 1 1 0 1 0 1 1 0 0 1 1 0 1 0 1 1 0 1 0 0 1 0 1 1 1 0 1 0 0 1 0 1 1 1 0 0 1 1 1 0 1 0 0 1 0 1 1 1 0 1 0 0 1 1 0 0 1 1 0 1 0 1 1 0 0 1 1 0 1 0 1 0 1 1 1 1 0 0 1 0 0 1 1 1 1 0 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 0 0 1 0 1 0 0 1 1 0 0 1 1 0 1 0 1 1 0 0 1 1 0 1 1 0 1 1 1 0 1 0 0 1 0 1 1 1 0 1 0 0 0 1 0 1 1 0 0 1 1 0 1 0 1 1 0 0 1 1 1 0 1 0 0 1 0 1 1 1 0 1 0 0 1 0 1 1 0 1 0 0 0 1 1 0 0 0 1 0 0 0 1 1 0 0 1 0 1 0 1 1 0 0 1 1 0 1 0 1 1 0 0 1 0 1 0 0 1 1 1 1 0 0 1 0 0 1 1 1 1 0 1 0 1 1 0 0 1 1 0 1 0 1 1 0 0 1 1 0 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 0 0 1 0 0 1 0 1 1 1 0 1 0 0 1 0 1 1 1 0 1 1 1 0 0 1 1 0 1 0 1 1 0 0 1 1 0 1 0 0 0 1 1 0 0 0 1 0 0 0 1 1 0 0 0 1 0 1 1 0 1 0 0 1 0 1 1 1 0 1 0 0 1 0 1 0 0 1 0 1 0 1 0 0 0 0 1 0 1 0 1 0 0 1 1 0 0 1 0 0 1 1 1 1 0 0 1 0 0 1 1 0 0 1 1 0 1 0 1 1 0 0 1 1 0 1 0 1 1 1 1 0 1 0 1 1 0 0 1 1 0 1 0 1 1 0 0 tnons 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 0.0000000 0.0000000 0.0000000 0.2500000 0.2500000 0.2500000 tolvrs1 1.00000000E10 tolvrs2 0.00000000E+00 tolvrs3 0.00000000E+00 tolvrs4 0.00000000E+00 tolwfr1 0.00000000E+00 tolwfr2 1.00000000E12 tolwfr3 0.00000000E+00 tolwfr4 0.00000000E+00 typat 1 1 wtk 0.18750 0.37500 0.09375 0.18750 0.12500 0.03125 xangst 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 1.3573395400E+00 1.3573395400E+00 1.3573395400E+00 xcart 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 2.5650000000E+00 2.5650000000E+00 2.5650000000E+00 xred 0.0000000000E+00 0.0000000000E+00 0.0000000000E+00 2.5000000000E01 2.5000000000E01 2.5000000000E01 znucl 14.00000 ================================================================================  Timing analysis has been suppressed with timopt=0 ================================================================================ Suggested references for the acknowledgment of ABINIT usage. The users of ABINIT have little formal obligations with respect to the ABINIT group (those specified in the GNU General Public License, http://www.gnu.org/copyleft/gpl.txt). However, it is common practice in the scientific literature, to acknowledge the efforts of people that have made the research possible. In this spirit, please find below suggested citations of work written by ABINIT developers, corresponding to implementations inside of ABINIT that you have used in the present run. Note also that it will be of great value to readers of publications presenting these results, to read papers enabling them to understand the theoretical formalism and details of the ABINIT implementation. For information on why they are suggested, see also https://docs.abinit.org/theory/acknowledgments.   [1] The Abinit project: Impact, environment and recent developments.  Computer Phys. Comm. 248, 107042 (2020).  X.Gonze, B. Amadon, G. Antonius, F.Arnardi, L.Baguet, J.M.Beuken,  J.Bieder, F.Bottin, J.Bouchet, E.Bousquet, N.Brouwer, F.Bruneval,  G.Brunin, T.Cavignac, J.B. Charraud, Wei Chen, M.Cote, S.Cottenier,  J.Denier, G.Geneste, Ph.Ghosez, M.Giantomassi, Y.Gillet, O.Gingras,  D.R.Hamann, G.Hautier, Xu He, N.Helbig, N.Holzwarth, Y.Jia, F.Jollet,  W.LafargueDitHauret, K.Lejaeghere, M.A.L.Marques, A.Martin, C.Martins,  H.P.C. Miranda, F.Naccarato, K. Persson, G.Petretto, V.Planes, Y.Pouillon,  S.Prokhorenko, F.Ricci, G.M.Rignanese, A.H.Romero, M.M.Schmitt, M.Torrent,  M.J.van Setten, B.Van Troeye, M.J.Verstraete, G.Zerah and J.W.Zwanzig  Comment: the fifth generic paper describing the ABINIT project.  Note that a version of this paper, that is not formatted for Computer Phys. Comm.  is available at https://www.abinit.org/sites/default/files/ABINIT20.pdf .  The licence allows the authors to put it on the Web.  DOI and bibtex: see https://docs.abinit.org/theory/bibliography/#gonze2020   [2] Optimized normconserving Vanderbilt pseudopotentials.  D.R. Hamann, Phys. Rev. B 88, 085117 (2013).  Comment: Some pseudopotential generated using the ONCVPSP code were used.  DOI and bibtex: see https://docs.abinit.org/theory/bibliography/#hamann2013   [3] ABINIT: Overview, and focus on selected capabilities  J. Chem. Phys. 152, 124102 (2020).  A. Romero, D.C. Allan, B. Amadon, G. Antonius, T. Applencourt, L.Baguet,  J.Bieder, F.Bottin, J.Bouchet, E.Bousquet, F.Bruneval,  G.Brunin, D.Caliste, M.Cote,  J.Denier, C. Dreyer, Ph.Ghosez, M.Giantomassi, Y.Gillet, O.Gingras,  D.R.Hamann, G.Hautier, F.Jollet, G. Jomard,  A.Martin,  H.P.C. Miranda, F.Naccarato, G.Petretto, N.A. Pike, V.Planes,  S.Prokhorenko, T. Rangel, F.Ricci, G.M.Rignanese, M.Royo, M.Stengel, M.Torrent,  M.J.van Setten, B.Van Troeye, M.J.Verstraete, J.Wiktor, J.W.Zwanziger, and X.Gonze.  Comment: a global overview of ABINIT, with focus on selected capabilities .  Note that a version of this paper, that is not formatted for J. Chem. Phys  is available at https://www.abinit.org/sites/default/files/ABINIT20_JPC.pdf .  The licence allows the authors to put it on the Web.  DOI and bibtex: see https://docs.abinit.org/theory/bibliography/#romero2020   [4] Recent developments in the ABINIT software package.  Computer Phys. Comm. 205, 106 (2016).  X.Gonze, F.Jollet, F.Abreu Araujo, D.Adams, B.Amadon, T.Applencourt,  C.Audouze, J.M.Beuken, J.Bieder, A.Bokhanchuk, E.Bousquet, F.Bruneval  D.Caliste, M.Cote, F.Dahm, F.Da Pieve, M.Delaveau, M.Di Gennaro,  B.Dorado, C.Espejo, G.Geneste, L.Genovese, A.Gerossier, M.Giantomassi,  Y.Gillet, D.R.Hamann, L.He, G.Jomard, J.Laflamme Janssen, S.Le Roux,  A.Levitt, A.Lherbier, F.Liu, I.Lukacevic, A.Martin, C.Martins,  M.J.T.Oliveira, S.Ponce, Y.Pouillon, T.Rangel, G.M.Rignanese,  A.H.Romero, B.Rousseau, O.Rubel, A.A.Shukri, M.Stankovski, M.Torrent,  M.J.Van Setten, B.Van Troeye, M.J.Verstraete, D.Waroquier, J.Wiktor,  B.Xu, A.Zhou, J.W.Zwanziger.  Comment: the fourth generic paper describing the ABINIT project.  Note that a version of this paper, that is not formatted for Computer Phys. Comm.  is available at https://www.abinit.org/sites/default/files/ABINIT16.pdf .  The licence allows the authors to put it on the Web.  DOI and bibtex: see https://docs.abinit.org/theory/bibliography/#gonze2016   And optionally:   [5] ABINIT: Firstprinciples approach of materials and nanosystem properties.  Computer Phys. Comm. 180, 25822615 (2009).  X. Gonze, B. Amadon, P.M. Anglade, J.M. Beuken, F. Bottin, P. Boulanger, F. Bruneval,  D. Caliste, R. Caracas, M. Cote, T. Deutsch, L. Genovese, Ph. Ghosez, M. Giantomassi  S. Goedecker, D.R. Hamann, P. Hermet, F. Jollet, G. Jomard, S. Leroux, M. Mancini, S. Mazevet,  M.J.T. Oliveira, G. Onida, Y. Pouillon, T. Rangel, G.M. Rignanese, D. Sangalli, R. Shaltaf,  M. Torrent, M.J. Verstraete, G. Zerah, J.W. Zwanziger  Comment: the third generic paper describing the ABINIT project.  Note that a version of this paper, that is not formatted for Computer Phys. Comm.  is available at https://www.abinit.org/sites/default/files/ABINIT_CPC_v10.pdf .  The licence allows the authors to put it on the Web.  DOI and bibtex: see https://docs.abinit.org/theory/bibliography/#gonze2009   Proc. 0 individual time (sec): cpu= 11.6 wall= 11.7 ================================================================================ Calculation completed. .Delivered 1 WARNINGs and 8 COMMENTs to log file. +Overall time at end (sec) : cpu= 11.6 wall= 11.7
After the description of the unit cell and of the pseudopotentials, you will find the list of kpoints used for the electrons and the grid of qpoints (in the Irreducible part of the Brillouin Zone) on which the susceptibility and dielectric matrices will be computed.
==== Kmesh for the wavefunctions ====
Number of points in the irreducible wedge : 6
Reduced coordinates and weights :
1) 2.50000000E01 2.50000000E01 0.00000000E+00 0.18750
2) 2.50000000E01 2.50000000E01 0.00000000E+00 0.37500
3) 5.00000000E01 5.00000000E01 0.00000000E+00 0.09375
4) 2.50000000E01 5.00000000E01 2.50000000E01 0.18750
5) 5.00000000E01 0.00000000E+00 0.00000000E+00 0.12500
6) 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.03125
Together with 48 symmetry operations and timereversal symmetry
yields 32 points in the full Brillouin Zone.
==== Qmesh for the screening function ====
Number of points in the irreducible wedge : 6
Reduced coordinates and weights :
1) 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.03125
2) 5.00000000E01 5.00000000E01 0.00000000E+00 0.09375
3) 5.00000000E01 2.50000000E01 2.50000000E01 0.37500
4) 0.00000000E+00 5.00000000E01 0.00000000E+00 0.12500
5) 5.00000000E01 2.50000000E01 2.50000000E01 0.18750
6) 0.00000000E+00 2.50000000E01 2.50000000E01 0.18750
Together with 48 symmetry operations and timereversal symmetry
yields 32 points in the full Brillouin Zone.
The qmesh is the set of all the possible momentum transfers. These points are obtained as all the possible differences among the kpoints ( \mathbf{q} =\mathbf{k}\mathbf{k}' ) of the grid chosen to generate the WFK file. From the last statement it is clear the importance of choosing homogeneous kpoint grids in order to minimize the number of qpoints is clear.
After this section, the code prints the parameters of the FFT grid needed to represent the wavefunctions and to compute their convolution (required for the screening matrices). Then we have some information about the MPI distribution of the bands and the total number of valence electrons computed by integrating the density in the unit cell.
setmesh: FFT mesh size selected = 20x 20x 20
total number of points = 8000
 screening: taking advantage of timereversal symmetry
 Maximum band index for partially occupied states nbvw = 4
 Remaining bands to be divided among processors nbcw = 56
 Number of bands treated by each node ~56
With the valence density, one can obtain the classical Drude plasma frequency. The next lines calculate the average density of the system, and evaluate the Wigner radius r_s, then compute the Drude plasma frequency, reported as omega_plasma.
Number of electrons calculated from density = 7.9999; Expected = 8.0000
average of density, n = 0.029628
r_s = 2.0048
omega_plasma = 16.6039 [eV]
This omega_plasma is the value used when the default for ppmfrq, namely 0.0, is specified. It is in fact the second frequency where the code calculates the dielectric matrix to adjust the plasmonpole model parameters.
It has been found that Drude plasma frequency is a reasonable value where to adjust the model. The control over this parameter is however left to the user in order to check that the result does not change when changing ppmfrq. One has to be careful with finite systems or with systems having semicore electrons. If the result depends much on ppmfrq, then the plasmonpole model is not appropriate and one should go beyond it by taking into account a full dynamical dependence in the screening (see later, the contourdeformation method). However, the plasmonpole model has been found to work well for a very large range of solidstate systems when focusing only on the real part of the GW corrections in the band gap region.
At the end of the screening calculation, the macroscopic dielectric constant is printed:
dielectric constant = 22.4176
dielectric constant without local fields = 24.7005
Note
Note that the convergence in the dielectric constant does not guarantee the convergence in the GW corrections and viceversa. In fact, the dielectric constant is representative of only one element i.e. the head of the full dielectric matrix. Even if the convergence on the dielectric constant with local fields takes somehow into account also other nondiagonal elements. In a GW calculation the whole \epsilon^{1} matrix is used to build the SelfEnergy operator.
The dielectric constant reported here is the socalled RPA dielectric constant due to the electrons. Although evaluated at zero frequency, it is understood that the ionic response is not included (this term can be computed with DFPT and ANADDB). The RPA dielectric constant restricted to electronic effects is also not the same as the one computed in the DFPT part of ABINIT, that includes exchangecorrelation effects.
We now enter the fourth dataset. As for dataset 3, after some general information (origin of WFK file, header, description of unit cell, kpoints, qpoints), the description of the FFT grid and jellium parameters, there is the echo of parameters for the plasmonpole model, and the inverse dielectric function (the screening). The selfenergy operator has been constructed, and one can evaluate the GW energies for each state.
The final results are:
 !SelfEnergy_ee
kpoint : [ 0.000, 0.000, 0.000, ]
spin : 1
KS_gap : 2.544
QP_gap : 3.110
Delta_QP_KS: 0.567
data: !SigmaeeData 
Band E0 <VxcDFT> SigX SigC(E0) Z dSigC/dE Sig(E) EE0 E
2 4.418 11.332 13.262 1.352 0.766 0.305 11.775 0.443 3.975
3 4.418 11.332 13.262 1.352 0.766 0.305 11.775 0.443 3.975
4 4.418 11.332 13.262 1.352 0.766 0.305 11.775 0.443 3.975
5 6.961 10.028 5.550 4.316 0.766 0.305 9.904 0.124 7.085
6 6.961 10.028 5.550 4.316 0.766 0.305 9.904 0.124 7.085
7 6.961 10.028 5.550 4.316 0.766 0.305 9.904 0.124 7.085
...
For the desired kpoint (\Gamma point), for state 4, then state 5, one finds different information in the SigmaeeData section:
 E0 is the KS eigenenergy
 VxcDFT gives the KS exchangecorrelation potential expectation value
 SigX gives the exchange contribution to the selfenergy
 SigC(E0) gives the correlation contribution to the selfenergy, evaluated at the KS 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 selfenergy, evaluated at the GW energy
 EE0 is the difference between GW energy and KS eigenenergy
 E is the final GW quasiparticle energy
In this case, prior to the SigmaeeData section, the direct band gap was also analyzed: KS_gap is the direct KS gap at that particular kpoint (and spin, in the case of spinpolarized calculations), QP_gap is the GW one, and Delta_QP_KS is the difference. This direct gap is always computed between the band whose number is equal to the number of electrons in the cell divided by two (integer part, in case of spinpolarized calculation), and the next one. This means that the value reported by the code may be wrong if the final QP energies obtained in the perturbative approach are not ordered by increasing energy anymore. So it’s always a good idea to check that the “gap” reported by the code corresponds to the real QP direct gap.
Warning
For a metal, these two bands do not systematically lie below and above the KS Fermi energy  but the concept of a direct gap is not relevant in that case. Moreover one should compute the Fermi energy of the QP system.
It is seen that the KS exchangecorrelation potential expectation value for the state 4 (a valence state) is rather close to the exchange selfenergy correction. For that state, the correlation correction is small, and the difference between KS and GW energies is also small (0.128 eV). By contrast, the exchange selfenergy is much smaller than the average KohnSham 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 KohnSham and GW energies is not very large, but nevertheless, it is quite important when compared with the size of the gap.
If AbiPy is installed on your machine, you can use the abiopen.py script
with the print
option to extract the results from the SIGRES.nc file
and print them to terminal:
abiopen.py tgw1_1o_DS4_SIGRES.nc p
================================= Structure =================================
Full Formula (Si2)
Reduced Formula: Si
abc : 3.839136 3.839136 3.839136
angles: 60.000000 60.000000 60.000000
Sites (2)
# SP a b c
    
0 Si 0 0 0
1 Si 0.25 0.25 0.25
Abinit Spacegroup: spgid: 0, num_spatial_symmetries: 48, has_timerev: True, symmorphic: True
============================== KohnSham bands ==============================
Number of electrons: 8.0, Fermi level: 4.773 (eV)
nsppol: 1, nkpt: 6, mband: 80, nspinor: 1, nspden: 1
smearing scheme: none, tsmear_eV: 0.272, occopt: 1
Direct gap:
Energy: 2.544 (eV)
Initial state: spin: 0, kpt: [+0.000, +0.000, +0.000], weight: 0.031, band: 3, eig: 4.418, occ: 2.000
Final state: spin: 0, kpt: [+0.000, +0.000, +0.000], weight: 0.031, band: 4, eig: 6.961, occ: 0.000
Fundamental gap:
Energy: 0.710 (eV)
Initial state: spin: 0, kpt: [+0.000, +0.000, +0.000], weight: 0.031, band: 3, eig: 4.418, occ: 2.000
Final state: spin: 0, kpt: [+0.500, +0.500, +0.000], weight: 0.094, band: 4, eig: 5.128, occ: 0.000
Bandwidth: 11.985 (eV)
Valence maximum located at:
spin: 0, kpt: [+0.000, +0.000, +0.000], weight: 0.031, band: 3, eig: 4.418, occ: 2.000
Conduction minimum located at:
spin: 0, kpt: [+0.500, +0.500, +0.000], weight: 0.094, band: 4, eig: 5.128, occ: 0.000
TIP: Call set_fermie_to_vbm() to set the Fermi level to the VBM if this is a nonmagnetic semiconductor
=============================== QP direct gaps ===============================
QP_dirgap: 3.110 (eV) for Kpoint: [+0.000, +0.000, +0.000] $\Gamma$, spin: 0
============== QP results for each kpoint and spin (all in eV) ==============
Kpoint: [+0.000, +0.000, +0.000] $\Gamma$, spin: 0
band e0 qpe qpe_diago vxcme sigxme sigcmee0 vUme ze0
1 1 4.418 3.975 3.840 11.332 13.262 1.352 0.0 0.766
2 2 4.418 3.975 3.840 11.332 13.262 1.352 0.0 0.766
3 3 4.418 3.975 3.840 11.332 13.262 1.352 0.0 0.766
4 4 6.961 7.085 7.123 10.028 5.550 4.316 0.0 0.766
5 5 6.961 7.085 7.123 10.028 5.550 4.316 0.0 0.766
6 6 6.961 7.085 7.123 10.028 5.550 4.316 0.0 0.766
For further details about the SIGRES.nc file and the AbiPy API see the Sigres notebook .
2 Preparing convergence studies: KohnSham structure (WFK file) and screening (SCR file)¶
In the following sections, we will perform different convergence studies. In order to keep the CPU time at a reasonable level, we will use fake WFK and SCR data. We will focus on the GW correction for \Gamma point to determine the values of the GW parameters needed to reach the convergence. Indeed, we will use a coarse kpoint grid with one shift only, and we will not vary ecut. This is a common strategy to find the adequate specific GW parameters before the final calculations, that should be done with a sufficiently fine kpoint grid, and an adequate ecut, in addition to adequate specific GW parameters.
In directory Work_gw1, copy the file tgw1_2.abi:
cp $ABI_TESTS/tutorial/Input/tgw1_2.abi .
Edit the tgw1_2.abi file, and take the time to examine it. Then, issue:
abinit tgw1_2.abi > tgw1_2.log 2> err &
After this step you will need the WFK and SCR files produced in this run for the next runs. Move tgw1o_DS2_WFK to tgw1o_DS1_WFK and tgw1o_DS3_SCR to tgw1o_DS1_SCR.
The next sections are intended to show you how to find the converged values of parameters that are specific of a GW calculation. The following parameters might be used to decrease the CPU time and/or the memory requirements, in addition to the wellknown k point sampling and ecut. For optdriver = 3, one needs to study the convergence with respect to ecuteps and nband simultaneously, while for optdriver = 4, only the behaviour with respect to nband should be monitored. As mentioned above, the global convergence with respect to ecut and to the number of k points has to be monitored as well, but the determination of the adequate parameters can be done independently from the determination of the adequate values for ecuteps and nband. Altogether, one has to determine the adequate values of four parameters in GW calculations, instead of only two in groundstate calculations (ecut and the number of k points). The adequate values of ecut and the number of k points for converged results might perhaps be the same as for groundstate calculations, but this is not always the case !
We will test the convergence with respect to nband and ecuteps, simultaneously for optdriver=3 and =4. As a side note, there are actually other technical parameters like ecutwfn or ecutsigx. However, for them, one can use the default value of ecut. For PAW, pawecutdg can be tuned as well.
We begin by the convergence study with respect to nband, the most important parameter needed in the selfenergy calculation, optdriver = 4. This is because for the selfenergy calculation, we will not need a double dataset loop to check this convergence (as ecuteps is not a parameter of the optdriver = 4 calculation), and we will rely on the previously determined SCR file.
3 Convergence of the selfenergy with respect to the number of bands¶
Let us check the convergence of the band gap with respect to the number of bands in the calculation of \Sigma_c with a fixed screening file. This convergence study is very important. However most of the time, the converged nband is similar for \Sigma_c and for \chi_0 so that the same value is taken for both. Here we will proceed carefully and converge the two occurences of nband independently.
The convergence on the number of bands to calculate \Sigma_c will be done by defining five datasets, with increasing nband:
ndtset 5
nband: 50
nband+ 50
In directory Work_gw1, copy the file tgw1_3.abi:
cp $ABI_TESTS/tutorial/Input/tgw1_3.abi .
Edit the tgw1_3.abi file, and take the time to examine it. Then, issue:
cp tgw1_2o_DS2_WFK tgw1_3o_DS2_WFK
cp tgw1_2o_DS3_SCR tgw1_3o_DS3_SCR
abinit tgw1_3.abi > tgw1_3.log 2> err &
# Crystalline silicon # Calculation of the GW corrections # Convergence with respect to the number of bands in Sigmac ndtset 5 nband: 50 nband+ 50 # Calculation of the SelfEnergy matrix elements (GW corrections) optdriver 4 getwfk 2 getscr 3 ecutsigx 8.0 nkptgw 1 kptgw 0.000 0.000 0.000 bdgw 4 5 # Definition of the unit cell: fcc acell 3*10.26 # Experimental lattice constants rprim 0.0 0.5 0.5 # FCC primitive vectors (to be scaled by acell) 0.5 0.0 0.5 0.5 0.5 0.0 # Definition of the atom types ntypat 1 # There is only one type of atom znucl 14 # The keyword "znucl" refers to the atomic number of the # possible type(s) of atom. The pseudopotential(s) # mentioned in the "files" file must correspond # to the type(s) of atom. Here, the only type is Silicon. # Definition of the atoms natom 2 # There are two atoms typat 1 1 # They both are of type 1, that is, Silicon. xred # Reduced coordinate of atoms 0.0 0.0 0.0 0.25 0.25 0.25 # Definition of the planewave basis set (at convergence 16 Rydberg 8 Hartree) ecut 8.0 # Maximal kinetic energy cutoff, in Hartree # Definition of the kpoint grid ngkpt 2 2 2 nshiftk 1 shiftk 0.0 0.0 0.0 istwfk *1 # This is mandatory in all the GW steps. pp_dirpath "$ABI_PSPDIR" pseudos "Pseudodojo_nc_sr_04_pbe_standard_psp8/Si.psp8" ############################################################## # This section is used only for regression testing of ABINIT # ############################################################## #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% test_chain = tgw1_2.abi, tgw1_3.abi, tgw1_4.abi, tgw1_5.abi #%% [files] #%% files_to_test = #%% tgw1_3.abo, tolnlines= 70, tolabs= 9.000e03, tolrel= 3.000e02 #%% [paral_info] #%% max_nprocs = 4 #%% [extra_info] #%% authors = Unknown #%% keywords = GW #%% description = #%% Crystalline silicon #%% Calculation of the GW corrections #%%<END TEST_INFO>
Edit the output file. The number of bands used for the selfenergy is mentioned in the fragments of output:
SIGMA fundamental parameters:
PLASMON POLE MODEL
number of planewaves for SigmaX 283
number of planewaves for SigmaC and W 169
number of planewaves for wavefunctions 283
number of bands 50
Gathering the GW energies for each number of bands, one gets:
number of bands 50
4 4.665 11.412 13.527 1.904 0.785 0.273 11.578 0.165 4.500
5 7.108 9.962 4.945 4.344 0.793 0.261 9.428 0.534 7.643
number of bands 100
4 4.665 11.412 13.527 1.768 0.784 0.275 11.684 0.271 4.394
5 7.108 9.962 4.945 4.470 0.792 0.263 9.528 0.434 7.542
number of bands 150
4 4.665 11.412 13.527 1.741 0.784 0.275 11.705 0.293 4.372
5 7.108 9.962 4.945 4.494 0.792 0.263 9.547 0.415 7.523
number of bands 200
4 4.665 11.412 13.527 1.733 0.784 0.275 11.711 0.299 4.366
5 7.108 9.962 4.945 4.500 0.792 0.263 9.553 0.410 7.518
number of bands 250
4 4.665 11.412 13.527 1.731 0.784 0.275 11.713 0.300 4.365
5 7.108 9.962 4.945 4.502 0.792 0.263 9.554 0.408 7.516
So that nband = 100 can be considered converged within 30 meV, which is fair to compare with experimental accuracy.
With AbiPy , one can use the abicomp.py script provides to compare multiple SIGRES.nc files
Use the expose
option to visualize of the QP gaps extracted from the different netcdf files:
$ abicomp.py sigres tgw1_3o_*_SIGRES.nc e sns
Output of robot.get_dataframe():
nsppol qpgap nspinor nspden nband nkpt \
tgw1_3o_DS1_SIGRES.nc 1 3.142871 1 1 50 3
tgw1_3o_DS2_SIGRES.nc 1 3.148588 1 1 100 3
tgw1_3o_DS3_SIGRES.nc 1 3.151012 1 1 150 3
tgw1_3o_DS4_SIGRES.nc 1 3.151603 1 1 200 3
tgw1_3o_DS5_SIGRES.nc 1 3.151485 1 1 250 3
ecutwfn ecuteps ecutsigx scr_nband sigma_nband \
tgw1_3o_DS1_SIGRES.nc 8.0 5.062893 8.0 60 50
tgw1_3o_DS2_SIGRES.nc 8.0 5.062893 8.0 60 100
tgw1_3o_DS3_SIGRES.nc 8.0 5.062893 8.0 60 150
tgw1_3o_DS4_SIGRES.nc 8.0 5.062893 8.0 60 200
tgw1_3o_DS5_SIGRES.nc 8.0 5.062893 8.0 60 250
gwcalctyp scissor_ene
tgw1_3o_DS1_SIGRES.nc 0 0.0
tgw1_3o_DS2_SIGRES.nc 0 0.0
tgw1_3o_DS3_SIGRES.nc 0 0.0
tgw1_3o_DS4_SIGRES.nc 0 0.0
tgw1_3o_DS5_SIGRES.nc 0 0.0
Invoking the script without options will open an ipython terminal to interact with the AbiPy robot.
Use the nb
option to automatically generate a jupyter notebook that will open in your browser.
For further details about the API provided by SigRes Robots see the Sigres notebook
and the notebook with the GW lesson for GW calculations powered by AbiPy.
4 Convergence of the screening with respect to the number of bands¶
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 SelfEnergy matrix elements (GW corrections)
optdriver?2 4
getscr?2 1
ecutsigx 8.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 two different parameters ecuteps and nband.
Let us begin with nband. This convergence study is rather important. It can be done at the same time as the convergence study for the number of bands for the selfenergy. Note that the number of bands used to calculate both the screening and the selfenergy can be lowered by a large amount by resorting to the extrapolar technique (see the input variable gwcomp).
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 Work_gw1, copy the file tgw1_4.abi:
cp $ABI_TESTS/tutorial/Input/tgw1_4.abi .
Edit the tgw1_4.abi file, and take the time to examine it. Then, issue:
cp tgw1_2o_DS2_WFK tgw1_4o_DS2_WFK
abinit tgw1_4.abi > tgw1_4.log 2> err &
# Crystalline silicon # Calculation of the GW corrections ndtset 10 udtset 5 2 # Convergence with respect to the number of bands to calculate epsilon^1 nband11 25 nband21 50 nband31 100 nband41 150 nband51 200 # Calculation of the screening (epsilon^1 matrix) optdriver?1 3 ecuteps 6.0 # Calculation of the SelfEnergy matrix elements (GW corrections) optdriver?2 4 getscr?2 1 ecutsigx 8.0 nband?2 100 nkptgw?2 1 kptgw?2 0.000 0.000 0.000 bdgw?2 4 5 getwfk 2 ppmfrq 16.7 eV # It is easier (and safer) to let ABINIT compute and use the Drude plasma frequency, # instead of selecting a value by hand. This would be done thanks to the default value ppmfrq 0.0 . # Definition of the unit cell: fcc acell 3*10.26 # Experimental lattice constants rprim 0.0 0.5 0.5 # FCC primitive vectors (to be scaled by acell) 0.5 0.0 0.5 0.5 0.5 0.0 # Definition of the atom types ntypat 1 # There is only one type of atom znucl 14 # The keyword "znucl" refers to the atomic number of the # possible type(s) of atom. The pseudopotential(s) # mentioned in the "files" file must correspond # to the type(s) of atom. Here, the only type is Silicon. # Definition of the atoms natom 2 # There are two atoms typat 1 1 # They both are of type 1, that is, Silicon. xred # Reduced coordinate of atoms 0.0 0.0 0.0 0.25 0.25 0.25 # Definition of the planewave basis set (at convergence 16 Rydberg 8 Hartree) ecut 8.0 # Maximal kinetic energy cutoff, in Hartree # Sampling of the BZ ngkpt 2 2 2 nshiftk 1 shiftk 0.0 0.0 0.0 istwfk *1 # This is mandatory in all the GW steps. pp_dirpath "$ABI_PSPDIR" pseudos "Pseudodojo_nc_sr_04_pbe_standard_psp8/Si.psp8" ############################################################## # This section is used only for regression testing of ABINIT # ############################################################## #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% test_chain = tgw1_2.abi, tgw1_3.abi, tgw1_4.abi, tgw1_5.abi #%% [files] #%% files_to_test = #%% tgw1_4.abo, tolnlines= 70, tolabs= 7.000e02, tolrel= 3.000e02 #%% [paral_info] #%% max_nprocs = 4 #%% [extra_info] #%% authors = Unknown #%% keywords = GW #%% description = #%% Crystalline silicon #%% Calculation of the GW corrections #%%<END TEST_INFO>
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 on file 169
dimension of the eps^1 matrix used 169
number of planewaves for wavefunctions 283
number of bands 25
Gathering the GW energies for each number of bands, one gets:
number of bands 25
4 4.665 11.412 13.527 1.968 0.786 0.273 11.527 0.115 4.550
5 7.108 9.962 4.945 4.279 0.796 0.257 9.375 0.587 7.696
number of bands 50
4 4.665 11.412 13.527 1.798 0.784 0.275 11.661 0.248 4.417
5 7.108 9.962 4.945 4.446 0.789 0.268 9.512 0.451 7.559
number of bands 100
4 4.665 11.412 13.527 1.708 0.784 0.276 11.731 0.318 4.347
5 7.108 9.962 4.945 4.523 0.791 0.264 9.571 0.391 7.499
number of bands 150
4 4.665 11.412 13.527 1.685 0.783 0.276 11.749 0.336 4.329
5 7.108 9.962 4.945 4.542 0.790 0.265 9.586 0.376 7.485
number of bands 200
4 4.665 11.412 13.527 1.678 0.783 0.277 11.754 0.341 4.324
5 7.108 9.962 4.945 4.549 0.790 0.265 9.592 0.370 7.479
So that the computation using 100 bands can be considered converged within 30 meV. Note that the value of nband that gives a converged dielectric matrix is usually of the same order of magnitude than the one that gives a converged \Sigma_c.
5 Convergence of the screening matrix with respect to the number of planewaves¶
Then, 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 Work_gw1, get the file tgw1_5.abi:
cp $ABI_TESTS/tutorial/Input/tgw1_5.abi .
Edit the tgw1_5.abi file, and take the time to examine it. Then, issue:
cp tgw1_2o_DS2_WFK tgw1_5o_DS2_WFK
abinit tgw1_5.abi > tgw1_5.log 2> err &
# Crystalline silicon # Calculation of the GW corrections ndtset 12 udtset 6 2 # Calculation of the screening (epsilon^1 matrix) optdriver?1 3 nband?1 100 # Convergence with respect to the dimension of epsilon^1 matrix ecuteps:? 3.0 ecuteps+? 1.0 # Calculation of the SelfEnergy matrix elements (GW corrections) optdriver?2 4 getscr?2 1 ecutsigx 8.0 nband?2 100 nkptgw 1 kptgw 0.000 0.000 0.000 bdgw 4 5 # GW calculation general parameters getwfk 2 ppmfrq 16.7 eV # It is easier (and safer) to let ABINIT compute and use the Drude plasma frequency, # instead of selecting a value by hand. This would be done thanks to the default value ppmfrq 0.0 . # Definition of the unit cell: fcc acell 3*10.26 # Experimental lattice constants rprim 0.0 0.5 0.5 # FCC primitive vectors (to be scaled by acell) 0.5 0.0 0.5 0.5 0.5 0.0 # Definition of the atom types ntypat 1 # There is only one type of atom znucl 14 # The keyword "znucl" refers to the atomic number of the # possible type(s) of atom. The pseudopotential(s) # mentioned in the "files" file must correspond # to the type(s) of atom. Here, the only type is Silicon. # Definition of the atoms natom 2 # There are two atoms typat 1 1 # They both are of type 1, that is, Silicon. xred # Reduced coordinate of atoms 0.0 0.0 0.0 0.25 0.25 0.25 # Definition of the planewave basis set (at convergence 16 Rydberg 8 Hartree) ecut 8.0 # Maximal kinetic energy cutoff, in Hartree # Sampling of the BZ ngkpt 2 2 2 nshiftk 1 shiftk 0.0 0.0 0.0 istwfk *1 # This is mandatory in all the GW steps. pp_dirpath "$ABI_PSPDIR" pseudos "Pseudodojo_nc_sr_04_pbe_standard_psp8/Si.psp8" ############################################################## # This section is used only for regression testing of ABINIT # ############################################################## #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% test_chain = tgw1_2.abi, tgw1_3.abi, tgw1_4.abi, tgw1_5.abi #%% [files] #%% files_to_test = #%% tgw1_5.abo, tolnlines= 70, tolabs= 7.000e02, tolrel= 7.000e02 #%% [paral_info] #%% max_nprocs = 4 #%% [extra_info] #%% authors = Unknown #%% keywords = GW #%% description = #%% Crystalline silicon #%% Calculation of the GW corrections #%%<END TEST_INFO>
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 GW energies for each number of bands, one gets:
dimension of the eps^1 matrix 59
4 4.665 11.412 13.527 1.899 0.786 0.272 11.582 0.169 4.496
5 7.108 9.962 4.945 4.462 0.791 0.264 9.523 0.440 7.548
dimension of the eps^1 matrix 113
4 4.665 11.412 13.527 1.755 0.784 0.275 11.694 0.282 4.383
5 7.108 9.962 4.945 4.513 0.791 0.264 9.563 0.400 7.508
dimension of the eps^1 matrix 137
4 4.665 11.412 13.527 1.728 0.784 0.276 11.715 0.303 4.362
5 7.108 9.962 4.945 4.518 0.791 0.264 9.567 0.395 7.504
dimension of the eps^1 matrix 169
4 4.665 11.412 13.527 1.708 0.784 0.276 11.731 0.318 4.347
5 7.108 9.962 4.945 4.523 0.791 0.264 9.571 0.391 7.499
dimension of the eps^1 matrix 259
4 4.665 11.412 13.527 1.696 0.784 0.276 11.740 0.328 4.338
5 7.108 9.962 4.945 4.527 0.791 0.264 9.574 0.388 7.496
dimension of the eps^1 matrix 283
4 4.665 11.412 13.527 1.695 0.784 0.276 11.741 0.329 4.337
5 7.108 9.962 4.945 4.527 0.791 0.264 9.575 0.388 7.496
So that ecuteps = 6.0 (%npweps = 169) can be considered converged within 10 meV.
At this stage, we know that for the screening computation, we need ecuteps = 6.0 Ha and nband = 100.
Of course, until now, we have skipped the most difficult part of the convergence tests: the convergence in the number of kpoints. 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 kpoints (roughly), and the number of kpoints 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 kpoint grid, for the last calculation.
As discussed in [Setten2017], the convergence study for kpoints the number of bands and the cutoff energies can be decoupled in the sense that one can start from a reasonaby coarse kmesh to find the converged values of nband, ecuteps, ecutsigx and then fix these values and look at the convergence with respect to the BZ mesh.
6 Calculation of the GW corrections for the band gap at the zone center¶
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 at the \Gamma point.
In directory Work_gw1, get the file tgw1_6.abi:
cp $ABI_TESTS/tutorial/Input/tgw1_6.abi .
Then, edit the tgw1_6.abi file, and, without examining it, comment the line
ngkpt 2 2 2 # Density of k points used for the automatic tests of the tutorial
and uncomment the line
#ngkpt 4 4 4 # Density of k points needed for a converged calculation
Then, issue:
abinit tgw1_6.abi > tgw1_6.log 2> err &
This job lasts a couple of minutes or so. It is worth to run it before the examination of the input file. Now, you can examine it.
# Crystalline silicon # Calculation of the GW correction to the direct band gap in Gamma # Dataset 1: ground state calculation # Dataset 2: calculation of the WFK file # Dataset 3: calculation of the screening (epsilon^1 matrix for W) # Dataset 4: calculation of the SelfEnergy matrix elements (GW corrections) ndtset 4 ngkpt 2 2 2 # Density of k points used for the automatic tests of the tutorial #ngkpt 4 4 4 # Density of k points needed for a converged calculation nshiftk 4 shiftk 0.0 0.0 0.0 # This grid contains the Gamma point, which is the point at which 0.0 0.5 0.5 # we will compute the (direct) band gap. There are 19 k points 0.5 0.0 0.5 # in the irreducible Brillouin zone, if ngkpt 4 4 4 is used. 0.5 0.5 0.0 istwfk *1 # For the GW computations, do not take advantage of the # specificities of k points to reduce the number of components of the # wavefunction. # Dataset1: usual selfconsistent groundstate calculation # Definition of the kpoint grid nband1 6 tolvrs1 1e10 # Dataset2: calculation of WFK file # Definition of kpoints iscf2 2 # Non selfconsistent calculation getden2 1 # Read previous density file nband2 120 nbdbuf2 20 tolwfr2 1.0d12 # Dataset3: Calculation of the screening (epsilon^1 matrix) optdriver3 3 # Screening calculation getwfk3 1 # Obtain WFK file from previous dataset nband3 50 # Bands to be used in the screening calculation ecuteps3 6.0 # Dimension of the screening matrix ppmfrq3 16.7 eV # Imaginary frequency where to calculate the screening # It is easier (and safer) to let ABINIT compute and use the Drude plasma frequency, # instead of selecting a value by hand. This would be done thanks to the default value ppmfrq 0.0 . # Dataset4: Calculation of the SelfEnergy matrix elements (GW corrections) optdriver4 4 # SelfEnergy calculation getwfk4 2 # Obtain WFK file from dataset 1 getscr4 1 # Obtain SCR file from previous dataset nband4 100 # Bands to be used in the SelfEnergy calculation ecutsigx4 8.0 # Dimension of the G sum in Sigma_x nkptgw4 1 # number of kpoint where to calculate the GW correction kptgw4 # kpoints 0.000 0.000 0.000 # (Gamma) bdgw4 4 5 # calculate GW corrections for bands from 4 to 5 # Definition of the unit cell: fcc acell 3*10.26 # Experimental lattice constants rprim 0.0 0.5 0.5 # FCC primitive vectors (to be scaled by acell) 0.5 0.0 0.5 0.5 0.5 0.0 # Definition of the atom types ntypat 1 # There is only one type of atom znucl 14 # The keyword "znucl" refers to the atomic number of the # possible type(s) of atom. The pseudopotential(s) # mentioned in the "files" file must correspond # to the type(s) of atom. Here, the only type is Silicon. # Definition of the atoms natom 2 # There are two atoms typat 1 1 # They both are of type 1, that is, Silicon. xred # Reduced coordinate of atoms 0.0 0.0 0.0 0.25 0.25 0.25 # Definition of the planewave basis set (at convergence 16 Rydberg 8 Hartree) ecut 8.0 # Maximal kinetic energy cutoff, in Hartree # Definition of the SCF procedure nstep 10 # Maximal number of SCF cycles diemac 12.0 # Although this is not mandatory, it is worth to # precondition the SCF cycle. The model dielectric # function used as the standard preconditioner # is described in the "dielng" input variable section. # Here, we follow the prescription for bulk silicon. pp_dirpath "$ABI_PSPDIR" pseudos "Pseudodojo_nc_sr_04_pbe_standard_psp8/Si.psp8" ############################################################## # This section is used only for regression testing of ABINIT # ############################################################## #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% [files] #%% files_to_test = #%% tgw1_6.abo, tolnlines= 70, tolabs= 7.000e02, tolrel= 3.000e02, fld_options= ridiculous #%% [paral_info] #%% max_nprocs = 4 #%% [extra_info] #%% authors = V. Olevano, F. Bruneval, M. Giantomassi #%% keywords = GW #%% description = #%% Crystalline silicon #%% Calculation of the GW correction to the direct band gap in Gamma #%% Dataset 1: ground state calculation #%% Dataset 2: calculation of the WFK file #%% Dataset 3: calculation of the screening (epsilon^1 matrix for W) #%% Dataset 4: calculation of the SelfEnergy matrix elements (GW corrections) #%%<END TEST_INFO>
We need the usual part of the input file to perform a ground state calculation. This is done in datasets 1 and 2. At the end of dataset 2, we print out the density and wavefunction files. We use a set of 19 kpoints in the Irreducible Brillouin Zone. This set of kpoints is not shifted so it contains the \Gamma point.
In dataset 3 we calculate the screening. The screening calculation is very timeconsuming. So, we have decided to decrease a bit the parameters found in the previous convergence studies. Indeed, nband has been decreased from 100 to 50. The CPU time of this part is linear with respect to this parameter (or more exactly, with the number of conduction bands). Thus, the CPU time has been decreased by a factor of 2. 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. This would be annoying for the absolute positioning of the band energy as required for bandoffset or ionization potential of finite systems. However, as long as we are only interested in the gap energy that is fine enough.
Finally, in dataset 4, we calculate the selfenergy matrix element at \Gamma, using the previously determined parameters.
You should obtain the following results:
 !SelfEnergy_ee
iteration_state: {dtset: 4, }
kpoint : [ 0.000, 0.000, 0.000, ]
spin : 1
KS_gap : 2.564
QP_gap : 3.196
Delta_QP_KS: 0.632
data: !SigmaeeData 
Band E0 <VxcDFT> SigX SigC(E0) Z dSigC/dE Sig(E) EE0 E
2 4.369 11.316 12.769 0.817 0.765 0.308 11.803 0.487 3.882
3 4.369 11.316 12.769 0.817 0.765 0.308 11.803 0.487 3.882
4 4.369 11.316 12.769 0.817 0.765 0.308 11.803 0.487 3.882
5 6.933 10.039 5.840 4.010 0.765 0.307 9.894 0.144 7.078
6 6.933 10.039 5.840 4.010 0.765 0.307 9.894 0.144 7.078
7 6.933 10.039 5.840 4.010 0.765 0.307 9.894 0.144 7.078
...
So that the DFT energy gap in \Gamma is about 2.564 eV, while the GW correction is about 0.632 eV, so that the GW band gap found is 3.196 eV.
One can compare now what have been obtained to what one can get from the literature.
EXP 3.40 eV LandoltBoernstein
DFT (LDA)
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 (FLAPW) 3.12 eV W. Ku and A.G. Eguiluz, PRL 89, 126401 (2002)
GW 3.20 eV present work
The values are spread over an interval of 0.2 eV. They depend on the details of the calculations. In the case of pseudopotential calculations, they depend of course on the pseudopotential used. However, a GW result is hardly more accurate than 0.1 eV, in the present state of the art. But this goes also with the other source of inaccuracy, the choice of the pseudopotential, that can arrive up to even 0.2 eV. This can also be taken into account when choosing the level of accuracy for the convergence parameters in the GW calculation. As a reasonable target, the numerical sources of errors, due to insufficient ecuteps, nband, k point grid, should be kept lower than 0.02 or 0.03 eV.
7 How to compute GW band structures¶
Finally, it is possible to calculate a full GW band plot of a system via interpolation. There are three possible techniques.
The first one is based on the use of Wannier functions to interpolate a few selected points in the IBZ obtained using the direct GW approach [Hamann2009]. You need to have the Wannier90 plugin installed. See the directory tests/wannier90, test case 03, for an example of a file where a GW calculation is followed by the use of Wannier90.
# Selfconsistent GW test including wannier90 interface for GW quasiparticles # This test is poorly converged (see GW and wannier90 tutorials). # Silicon structure acell 10.263 10.263 10.263 rprim 0.00 0.50 0.50 0.50 0.00 0.50 0.50 0.50 0.00 natom 2 xred 0.00 0.00 0.00 0.25 0.25 0.25 ntypat 1 typat 1 1 znucl 14.00 symmorphi 0 symsigma 0 # Parameters common to all runs ecut 6.00 ecutsigx 5.00000000 ecuteps 1.49923969 istwfk 8*1 ngkpt 4 4 4 nstep 100 nshiftk 1 shiftk 0.00 0.00 0.00 enunit 2 gw_icutcoul 6 # To preserve the results of older tests: current default gw_icutcoul=3 ndtset 7 gwpara 1 # Selfconsistent run to get the density toldfe1 1.00d6 # Nonselfconsistent run to get all cg wavefunctions getden2 1 getwfk2 1 iscf2 2 tolwfr2 1.0d10 nband2 30 # Calculation of the dielectric matrix  iteration 1 optdriver3 3 gwcalctyp3 28 getwfk3 2 nband3 10 ecutwfn3 5.00 # Increased value to stabilize the test awtr3 0 # Note : the default awtr 1 is better # Calculation of the model GW corrections  iteration 1 optdriver4 4 gwcalctyp4 28 getwfk4 2 getscr4 3 nband4 10 ecutwfn4 5.00 # Increased value to stabilize the test gw_icutcoul4 3 # old deprecated value of icutcoul, only used for legacy # Calculation of the dielectric matrix  iteration 2 optdriver5 3 gwcalctyp5 28 getwfk5 2 getqps5 4 nband5 10 ecutwfn5 5.00 # Increased value to stabilize the test awtr5 0 # Note : the default awtr 1 is better # Calculation of the model GW corrections  iteration 2 optdriver6 4 gwcalctyp6 28 getwfk6 2 getqps6 4 getscr6 5 nband6 10 ecutwfn6 5.00 # Increased value to stabilize the test gw_icutcoul6 3 # old deprecated value of icutcoul, only used for legacy # Calculation of the quasiparticle Wannier functions getden7 1 getwfk7 2 # Must point to lda wavefunction file completely # consistent with _WFK file generated for GW getqps7 6 kptopt7 3 # Must use fullzone k mesh for wannier90 istwfk7 64*1 iscf7 2 nstep7 0 # Irreduciblezone wavefunctions will be transformed # using symmetry operations to fill the full zone, # and must not be modified (for consistency with GW) tolwfr7 1.0d10 # Dummy, but necessary nband7 10 # Must be consistent with nband in quasiparticle GW above prtwant7 3 # Flag for wannier90 interface with quaisparticles w90iniprj7 2 # Flag to use hydrogenic or gaussian orbital initial # projectors (to be specified in *.win file) w90prtunk7 0 # Flag for producing UNK files necessary for plotting # (suppressed here by 0 value) #Common to all model GW calculations rhoqpmix 0.5 nkptgw 8 kptgw 0.00000000E+00 0.00000000E+00 0.00000000E+00 2.50000000E01 0.00000000E+00 0.00000000E+00 5.00000000E01 0.00000000E+00 0.00000000E+00 2.50000000E01 2.50000000E01 0.00000000E+00 5.00000000E01 2.50000000E01 0.00000000E+00 2.50000000E01 2.50000000E01 0.00000000E+00 5.00000000E01 5.00000000E01 0.00000000E+00 2.50000000E01 5.00000000E01 2.50000000E01 bdgw 1 8 # Only bands 18 are quasiparticle. LDA will be used for # bands 9 and 10 in the wannier90 calculation. 1 8 1 8 1 8 1 8 1 8 1 8 1 8 pp_dirpath "$ABI_PSPDIR" pseudos "PseudosTM_pwteter/14si.pspnc" #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% [files] #%% files_to_test = #%% t03.abo, tolnlines = 30, tolabs = 7.0e03, tolrel = 1.05e00 #%% extra_inputs = t03o_DS7_w90.win #%% [paral_info] #%% max_nprocs = 1 #%% [extra_info] #%% authors = D. Hamann #%% keywords = GW #%% description = #%% Cannot be executed in parallel (mlwfovl_qp) #%% Si fcc, in primitive cell (2 atoms/cell) #%% Test of selfconsistent model GW (2 iterations) following Faleev et al., #%% [PRL 93, 126406 (2004)] followed by construction of quasiparticle #%% maximallylocalized wannier functions [Hamann & Vanderbilt, #%% arXiv:0810.3616v1 (condmat.mtrlsci)]. Cutoffs are set for test #%% acceptable speed, and the results are poorly converged. The input #%% file is sufficiently annotated to serve as a model. Note that well #%% converged GW calculations are extremely time consuming, and in general #%% it is advisable to run the SCGW part separately on a parallel system, #%% and then run a separate serial job modeled on the last dataset, #%% substituting "irdwfk" and "irdqps" for "getwfk" and "getqps," with #%% appropriate links to the files produced in the serial run. Note that #%% the _DEN file from the first dataset is also needed as input, although #%% the discontinued "irdden" input variable is not needed or supported. #%% Note that acceptable names for the secondary input file needed by the #%% wannier90 library are wannier90.win, w90.win, abo_DSn_w90.win (ndtset #%% >0) and abo_w90.win (ndtset=0), where abo is the 4th line of the .files #%% file and n is the wannier dataset. #%% topics = Wannier #%%<END TEST_INFO>
The wannier interpolation is a very accurate method, can handle band crossings but it may require additional work to obtain well localized wannier functions. Another practical way follows from the fact that the QP energies, similarly to the KS eigenvalues, must fulfill the symmetry properties:
and
where \GG is a reciprocal lattice vector and S is a rotation of the point group of the crystal. Therefore it’s possible to employ the starfunction interpolation by Shankland, Koelling and Wood [Euwema1969], [Koelling1986] in the improved version proposed by [Pickett1988] to fit the abinitio results. This interpolation technique, by construction, passes through the initial points and satisfies the basic symmetry property of the band energies. It should be stressed, however, that this Fourierbased method can have problems in the presence of band crossings that may cause unphysical oscillations between the abinitio points. To reduce this spurious effect, we suggest to interpolate the GW corrections instead of the GW energies. The corrections, indeed, are usually smoother in kspace and the resulting fit is more stable. A python example showing how to construct an interpolated scissor operator with AbiPy is available here
The third method uses the fact that the GW corrections are usually linear with the energy, for each group of bands. This is evident when reporting on a plot the GW correction with respect to the 0order KS energy for each state. One can then simply correct the KS band structure at any point, by using a GW correction for the kpoints where it has not been calculated explicitly, using a fit of the GW correction at a sparse set of points. A python example showing how to construct an energydependent scissor operator with AbiPy is available here.
8 Advanced features of GW calculations¶
The user might switch to the second GW tutorial before coming back to the present section.
Calculations without using the PlasmonPole model¶
In order to circumvent the plasmonpole model, the GW frequency convolution has to be calculated explicitly along the real axis. This is a tough job, since G and W have poles along the real axis. Therefore it is more convenient to use another path of integration along the imaginary axis plus the residues enclosed in the path.
Consequently, it is better to evaluate the screening for imaginary frequencies (to perform the integration) and also for real frequencies (to evaluate the contributions of the residues that may enter into the path of integration). The number of imaginary frequencies is set by the input variable nfreqim. The regular grid of real frequencies is determined by the input variables nfreqre, which sets the number of real frequencies, and freqremax, which indicates the maximum real frequency used.
The method is particularly suited to output the spectral function (contained in file out.sig). The grid of real frequencies used to calculate the spectral function is set by the number of frequencies (input variable nfreqsp) and by the maximum frequency calculated (input variable freqspmax).
Selfconsistent calculations¶
The details in the implementation and the justification for the approximations retained can be found in [Bruneval2006]. The only added input variables are getqps and irdqps. These variables concerns the reading of the _QPS file, that contains the eigenvalues and the unitary transform matrices of a previous quasiparticle calculation. QPS stands for “QuasiParticle Structure”.
The only modified input variables for selfconsistent calculations are gwcalctyp and bdgw. When the variable gwcalctyp is in between 0 and 9, The code calculates the quasiparticle energies only and does not output any QPS file (as in a standard GW run). When the variable gwcalctyp is in between 10 and 19, the code calculates the quasiparticle energies only and outputs them in a QPS file. When the variable gwcalctyp is in between 20 and 29, the code calculates the quasiparticle energies and wavefunctions and outputs them in a QPS file.
For a full selfconsistency calculation, the quasiparticle wavefunctions are expanded in the basis set of the KS wavefunctions. The variable bdgw now indicates the size of all matrices to be calculated and diagonalized. The quasiparticle wavefunctions are consequently linear combinations of the KS wavefunctions in between the min and max values of bdgw.
A correct selfconsistent calculation should consist of the following runs:
 1) Selfconsistent KS calculation: outputs a WFK file
 2) Screening calculation (with KS inputs): outputs a SCR file
 3) Sigma calculation (with KS inputs): outputs a QPS file
 4) Screening calculation (with the WFK, and QPS file as an input): outputs a new SCR file
 5) Sigma calculation (with the WFK, QPS and the new SCR files): outputs a new QPS file
 6) Screening calculation (with the WFK, the new QPS file): outputs a newer SCR file
 7) Sigma calculation (with the WFK, the newer QPS and SCR files): outputs a newer QPS
 ............ and so on, until the desired accuracy is reached
Note that for HartreeFock calculations a dummy screening is required for initialization reasons. Therefore, a correct HF calculations should look like
 1) Selfconsistent KS calculation: outputs a WFK file
 2) Screening calculation using very low convergence parameters (with KS inputs): output a dummy SCR file
 3) Sigma calculation (with KS inputs): outputs a QPS file
 4) Sigma calculation (with the WFK and QPS files): outputs a new QPS file
 5) Sigma calculation (with the WFK and the new QPS file): outputs a newer QPS file
 ............ and so on, until the desired accuracy is reached
In the case of a selfconsistent calculation, the output is slightly more complex: For instance, at iteration 2
 !SelfEnergy_ee
iteration_state: {dtset: 3, }
kpoint : [ 0.500, 0.250, 0.000, ]
spin : 1
KS_gap : 3.684
QP_gap : 5.764
Delta_QP_KS: 2.080
data: !SigmaeeData 
Band E_DFT <VxcDFT> E(N1) <Hhartree> SigX SigC[E(N1)] Z dSigC/dE Sig[E(N)] DeltaE E(N)_pert E(N)_diago
1 3.422 10.273 3.761 6.847 15.232 4.034 1.000 0.000 11.198 0.590 4.351 4.351
2 0.574 10.245 0.850 9.666 13.806 2.998 1.000 0.000 10.807 0.291 1.141 1.141
3 2.242 9.606 2.513 11.841 11.452 1.931 1.000 0.000 9.521 0.193 2.320 2.320
4 3.595 10.267 4.151 13.866 11.775 1.842 1.000 0.000 9.933 0.217 3.934 3.934
5 7.279 8.804 9.916 16.078 4.452 1.592 1.000 0.000 6.044 0.119 10.034 10.035
6 10.247 9.143 13.462 19.395 4.063 1.775 1.000 0.000 5.838 0.095 13.557 13.557
7 11.488 9.704 15.159 21.197 4.061 1.863 1.000 0.000 5.924 0.113 15.273 15.273
8 11.780 9.180 15.225 20.958 3.705 1.893 1.000 0.000 5.598 0.135 15.360 15.360
...
The columns are
 Band: Index of the band
 E_DFT: DFT eigenvalue
 VxcDFT: Diagonal expectation value of the xc potential in between DFT bra and ket
 E(N1): Quasiparticle energy of the previous iteration (equal to DFT for the first iteration)
 Hhartree: Diagonal expectation value of the Hartree Hamiltonian (equal to E_DFT  VxcDFT for the first iteration only)
 SigX: Diagonal expectation value of the exchange selfenergy
 SigC[E(N1)]: Diagonal expectation value of the correlation selfenergy (evaluated for the energy of the preceeding iteration)
 Z: Quasiparticle renormalization factor Z (taken equal to 1 in methods HF, SEX, COHSEX and model GW)
 dSigC/dE: Derivative of the correlation selfenergy with respect to the energy
 Sig[E(N)]: Total selfenergy for the new quasiparticle energy
 DeltaE: Energy difference with respect to the previous step
 E(N)_pert: QP energy as obtained by the usual perturbative method
 E(N)_diago: QP energy as obtained by the full diagonalization