Band structure calculations in QE using hybrid functionals

In recent years the interest in hybrid functionals (that is the incorporation of parts of Hartree-Fock exchange in calculations based on common approximations of exchange-correlations such as LDA, GGA and so on) has steadily increased owing to its improvement over most common functionals, especially when it comes to band-gap calculation of extended solids – see for example this work for a comparison. Quantum Espresso now offers a variety of hybrid functionals (for a complete list see the header of funct.f90) but currently the code can only use hybrid functionals for self-consistent calculations. The reason being that:
The problem is quite fundamental, because in order to get the Fock operator at a certain k-point you need the wavefunctions on a grid that is commensurate with it, this can only be done self-consistently. (L. Paulatto)
But there is a quite elegant solution via maximally localized Wannier functions. If you worked with wannier functions before you know they required the “full grid” (conveniently generated using the tool kmesh.pl) but a recently introduced tool allows to “unfold” the SCF calculation of the reduced grid to the full grid using the executable open_grid.x. An example is provided with recent versions of QE as well. In the following I will use this approach to improve the band-gap of MgO. Using PBE pseudopotentials the band-gap is underestimated at about 4.3 eV but using the HSE hybrid the band-gap can be improved to about 7 eV, which is within acceptable range of the experimental value of ~7.7 eV. The workflow is fairly straightforward:
  1. run a (converged) SCF calculation with input_dft=’HSE’ and a number of empty bands. You have to ensure convergence with respect to the usual parameters (k-points, cutoff, …) AND the mesh for the Fock operator (nqx)
  2.  unfold the reduced grid onto the full grid using open_grid.x
  3.  wannierize the obtained wave-functions using wannier90.x and plot the band-structure along a desired high symmetry path
The wannierization is by far the most tricky part in this particular example but by projecting on O:p and Mg:s one can accurately describe  the valence band and an additional single conduction band. For the resulting band-structure see the figure below.
Figure 1. Comparison of the band-structure of MgO on PBE level (purple) and the corresponding wannierzation (grey) as well as the final HSE hybrid calculation with vastly improved band-gap of 7.4 eV (blue)
So this is actually quite straight-forward. please note that in above calculation the Fock operator was calculated on a very coarse 1x1x1 grid. The calculation on such a coarse grid actually seems to over-estimate the band-gap and a converged energy can be obtained on a 6x6x6 nqx grid with ~6.7 eV. The entire calculation can be run using the attached script – I hope you find it helpful!

From “Fatbands” to Wannier Functions of Silver

Electronic bands – the energy solutions εi(k) of the Kohn-Sham equation of a non-interacting electron in an effective potential, are, despite of approximations made in practical DFT calculations, very successful in describing the physics of condensed matter. It is sometimes necessary (or simply practical) to describe the band-structure (“Spaghetti”) in terms of local orbitals, that is in a language that is more close to the “chemist’s picture” of atomic s,p,d orbitals. Generally, the s,p,d nomenclature is linked to the angular momentum l of a radially symmetrical potential (“Hydrogen model”) – and not a periodic crystal potential. It can be shown, however, that the resulting electronic bands can often be assigned to a certain l, at least for the low-energy core bands.

Most programs available allow to obtain this information either as a function of energy alone (the  density of states, DOS and projected density of states, PDOS) or resolved in energy and crystal momentum which is often called a “fat-band” representation as the weight of the orbital s,p,d characteristic is shown as the line thickness for practical reasons.  As a simple example I have calculate the PDOS of bulk silver using the ELK code, a full-potential (all-electrons) linearised augmented-plane wave (FP-LAPW) DFT code. The relevant task numbers in this case are 0 (the ground-state), 10 (the total and projected density of states) and 21 which the manual lists conveniently to give us the Band structure plot which includes angular momentum characters for every atom. Details can be found in the attached elk.in file.

fatbands_ELK_Ag

The figure shows the band structure on the left-hand side and the corresponding integrated information with the same color-code on the right-hand side. I included arrows to highlight how certain features in the PDOS can be interpreted based on the band structure and vice versa. The first thing to notice is that most bands above the Fermi level (set as 0 in the ELK code) are a mix of various l-states whereas just below the Fermi level the d-states form a fairly flat (localized) set of bands.

Now I usually work with the plane-wave pseudopotential method as implemented in Quantum-Espresso (QE 6.3) and not with full-potentials but I chose above calculation as an entry point to show what can be a problem with pseudopotentials – bear with me. Turning to QE I generate the same input structure as previously used in ELK using cif2cell with the program switch -p pwscf. I use my favorite set of norm-conserving pseudopotentials from the ONCV library and generate the wave-functions using the latest version (3.3.1) of the ONCV code. For silver I find that a regular k-point grid of 11x11x11 and a cutoff of the wave-functions of about 80 Ry gives converged results. The calculation of fatbands is a bit more difficult in QE than ELK and involves the following steps:

  1. calculation of the band-structure using pw.x and bands.x
  2. calculation of the “projection on wave-functions” using projwfc.x
  3. generation of the weights by calculating the overlap between projections and wave functions using plotband.x

All steps are automated in the bands.sh script attached but I want to highlight that the setting of  lsym=.true./.false. in the input for bands.x and projwfc.x is important. When the projection file is found by plotband.x it asks for a list of wave functions to calculate the overlaps. The proper numbers can be found in the output of projwfc.x and looks something like

state # 1: atom 1 (Ag ), wfc 1 (l=0 m= 1)
state # 2: atom 1 (Ag ), wfc 2 (l=1 m= 1)
state # 3: atom 1 (Ag ), wfc 2 (l=1 m= 2)
state # 4: atom 1 (Ag ), wfc 2 (l=1 m= 3)
state # 5: atom 1 (Ag ), wfc 3 (l=0 m= 1)
state # 6: atom 1 (Ag ), wfc 4 (l=2 m= 1)
state # 7: atom 1 (Ag ), wfc 4 (l=2 m= 2)
state # 8: atom 1 (Ag ), wfc 4 (l=2 m= 3)
state # 9: atom 1 (Ag ), wfc 4 (l=2 m= 4)
state # 10: atom 1 (Ag ), wfc 4 (l=2 m= 5)

The ordering here is the same as in the pseudopotential file, so in case there is any doubt consult the header of the pseudo and confirm the valence configuration:

#n l occ
4 0 2.00
4 1 6.00
4 2 10.00
5 0 1.00

The following figure shows the fatbands as obtained by QE and the band structure obtained from maximally localized wannier functions (MLWF) via wannier90.

bands-wannier_ONCV_QE

If you are familiar with wannier90 you know that the definition of an initial “guess” of wavefunctions is usually critical in order to obtained well localized wannier functions. Looking at the fatbands of QE the character of the bands just above the Fermi level remains, unfortunately a mystery. Overall the band-structure agrees very well with the results obtained from ELK which tells us that the pseudo has no problem but the wave functions generated with ONCV are less than perfect. If I had read the release notes carefully I would have been able to expect this:

The upf output now includes pseudo wave funtions for the occupied states.
This will allow LDA+U calculations in some codes. Please note that
ONCVPSP potentials have not been tested or benchmarked for this use.
Also be aware that LDA+U is a semi-empirical mean field theory, and
does NOT represent many-body physics.

Ultimately I had no problem defining proper projections (with the guidance of the ELK fatbands) but above should be kept in mind when performing, for example, LDA+U calculations in QE (note that this does not really pose a problem here as the highest-l manifold, the d-bands, have proper wave-functions!).

For wannier90 initial projections of

begin projections
Ag: s;p;d
end projections

result in an average spread of Ω~1.1 per atom which is fairly well localized and the resulting tight-binding bands agree well with the DFT bands, as can be seen on the right-hand side above.

Hope this was helpful – if there is any problem running attached scripts please leave it in the comments!

 

PS: as correctly mentioned in the comments, a “better” pseudopotential would be able to produce a more accurate fat-band representation. Therefore I attached a small comparison of 3 easily available PPs (from the PSLibrary, the GBRV library and ONCV library) – an updated “run script” is attached as well. HTH!

comparison_Ag_pseudos

fatbands_Ag

Ag_bands_comparison.sh

 

Work function calculation of MgO in QE 6.2.1

Using a slab-vacuum system and the dipole correction as previously described we aim to calculate the work-function (WF) of MgO. The experimental value for the WF of a low-index surface such as (200) is about 5 eV according to Choi et al.

pwi2xsf

By constructing a 6 layer MgO slab padded by 6 layers of vacuum we create an artificial surface that is infinitely repeated by periodic boundary conditions. The dipole correction is applied by setting tefield and dipfield to true. We know that the position of the dipole is critical as the discontinuity has to fall into the vacuum region. By moving the dipole across the vacuum region it can be shown that if the dipole is too close to the slab the dipole correction will fail.

pot-1

If placed correctly the electrostatic potential will be flat in the vacuum region. Two failing cases are also shown by dashed lines. The slab is repeated to the negative region to emphasize the region in which the dipole correction is physical which extends from emaxpos+eopreg-1 (in units of entire cell extending from [0-1]). The difference between vacuum level (here set to 0) and the flat region of the electrostatic potential is about 5 eV which is in acceptable agreement with experiment. Find the script to run this job here.

 

Dipole correction in Quantum-Espresso

For surfaces that have a non-vanishing dipole density quantum-espresso (QE) offers a dipole correction in combination with the sawtooth potential. The relevant variables in the pw.x input are dipfield and tefield, which both have to be set to true. Unfortunately this is not all there is to do – we also have to define where the dipole is located and the location (and size) of the dipole play an important role. Since this confuses me regularly I decided to write it down here in the hope it helps me and others.

Lets assume we have centered the slab in the center of our cell and the cell extends from 0 (the bottom) to 1 (the top). As QE uses periodic boundary conditions everything that exceeds the length of the cell will re-enter the cell at the bottom, i.e. (1+x)=x in relative z-coordinates. As the manual tells us one has to take care of 3 parameters: emaxpos, eopreg and eamp. Since we only want to apply a dipole correction the amplitude eamp=0. The field linearly increases until emaxpos and then decreases to reach 0 at emaxpos+eopreg, as such they represent the positive and negative charges that create the dipole. It is important that the change of slope, that is emaxpos, lies in the vacuum region and so must emaxpos+eopreg (the opposite charge) – otherwise the correction will not work. A very useful setup is to center the slab in the cell (~0.5 in relative z-coordinates) and the dipole around 0 of the slab. One can then slowly increase the length of the dipole and see if the resulting electrostatic potential (plot_num=11 in pp.x) becomes flat in the vacuum region (thanks for Dr. Thomas Brumme for the hint!). As long as the dipole and the charge density of the electrons do not overlap the correction works and the necessary amplitude of the dipole field decreases linearly with the dipole length – as to be expected since p=q*z (see the right side of Figure 1). When the dipole gets too close the correction fails and the potential is no longer flat in the vacuum region. Please note that this does not mean that SCF convergence cannot be reached – in fact, the last case converges faster than all the cases were the results are physically meaningful in this particular case.pot-1

Figure 1. Potential as obtained by pp.x (plot_num=11) and averaged in xy-plane (average.x) as a function of the dipole length with the dipole centered at 0 of the cell. On the right side the amplitude of the dipole field is shown as function of the dipole cell.

In recent versions of QE a different correction has been implemented by Sohier et al. as well. This is activated with the input-flag assume_isolated =’2d’ and gives basically the same result. As an advantage the new correction seems to be working better with phonons of 2-dimensional systems but not all codes that build on QE output are able to read it yet (for example the resulting electron-phonon matrix elements in EPW are too large by a factor of at least 100).

Attached shell script will run a calculation for dipoles of varying length – by collecting the results one can reproduce above figure.

Hope this was helpful – if not, let me know in the comments!

run