GW DOS calculations using QE and YAMBO

In this short tutorial I will show how to first perform a QE calculation of the molecule benzene in a vacuum box at the GGA level. Then, we will use YAMBO to calculate a single shot G0W0 correction to the Kohn-Sham eigenvalues. As a last step we calculate the density of states at the GW level as well as a projected density of states on the C:s orbitals.

The main purpose for this tutorial is to demonstrate the work-flow based on instructions by Daniele Varsano (here).

I will assume, for the sake of keeping this tutorial short, that one already knows the basic steps of creating the base folder for a YAMBO calculation by first performing an SCF and NSCF calculation. The crucial step for the GW-PDOS calculation is to perform a PDOS calculation using projwfc.x before creating the YAMBO base folder in p2y. That means the proper sequence is: scf (pw.x) > nscf (pw.x) >pdos (projwfc.x) > p2y > yambo initialization > yambo G0W0> ypp. After the p2y step there should be a new file in the base folder called ns.atom_proj_pwscf – this basically indicates that p2y found atom projections and created a database.

As always: this is a basic example and not a converged calculation (this will become obvious right away)

Let’s first see what the G0W0 (yambo -x -p p -g n) calculation does to the eigenvalues:

As you can see in above plot of the o.qp file, the gap between the highest occupied molecular orbital (HOMO) and lowest unoccupied MO (LUMO) at the GGA level is around 5 eV. This increases to around 8 eV at the G0W0 level (shown in red). A converged GW gap for this molecule is around 10 eV. [1]

Ignoring this discrepancy for now we move to the generation of a DOS input using ypp -s s -V all, increasing the verbosity to all helps with defining the input. A few lines deserve special attention:

To tell ypp that we attempt to create atom projected DOS we change the project mode line to: PROJECT_mode= “ATOM”. By itself, the GW corrected database will not be read, therefore the line related to that reads: GfnQPdb= “none”
last but not least you will see a larger block at the bottom of the file related to selecting atomic wave functions – I will come to that in a second.

To find out which atom types and wave functions are in our calculation we can take a look at the report generated during the initialization of the calculation:

``````  [02.01] Unit cells
==================

Cell kind             :  CUB
Atoms in the cell     :  C  H
number of C  atoms    :   6
number of H  atoms    :   6

List of atomic masses
Mass of C     :  12.01070
Mass of H     :  1.000000``````

As we can see, there are two types (C, H) and 6 atoms each (C6H6). The H pseudo potential will contain the 1s wave function with angular quantum number l=0, C will contain the 2s (l=0) and 2p (l=1). If we don’t care about the type of atom and wave function we can just leave them set to -1 and run ypp. This will generate the total DOS at the GGA level.

Advancing one step, we can insert the QP corrections by replacing GfnQPdb= “none” with

``GfnQPdb= "E < ./SAVE/ndb.QP"``

Note that the ndb.QP can also lie in a different folder if the GW step was called using a job identifier (-J)! The default folder is however SAVE. Running ypp one more time will now tell us in the output that it will apply corrections to the eigenvalues:

``[QP_apply] Action to be applied: E<./G0W0/ndb.QP[ PPA@E  27.21138 * XG 949 * Xb 1-100 * Scb 1-100]``

The last but very interesting step is to not calculate only the total DOS but select a specific atom type and angular character, for example the H:s or C:s or C:p. This can be done by changing the part at the bottom of the ypp input file:

``````% PDOS_kinds
1 |1 |``````

Allows us to project onto the first type of atom, which is carbon. The next setting

``````% PDOS_l
0 |0 |``````

Allows us to project onto the s (l=0) orbital of that atom. More fine-grained selections can be defined. The resulting data can be summarized in the following DOS plot

Shown are the GGA DOS (blue), which exhibits a significantly smaller gap than expected, the G0W0 DOS in orange as well as the C:s projected G0W0 DOS in green. Note that the DOS is aligned so that 0 eV corresponds to the highest occupied orbital.

A script to run the calculations and produce the plot (using gnuplot) is attached – have fun reproducing the results yourself!

[1] Neaton et al., https://arxiv.org/abs/cond-mat/0606640

Hybrid band structure calculations (QE 7.0 update)

Dear all,

I have received several emails and comments reporting issues with the pw.x + open_grid.x + wannier90 approach for band structure calculations with hybrid XC functionals. I decided to quickly test the latest version of QE and found that all is well – at least for my simple test case.

The attached run file will let you calculate the bands of bulk MgO with hybrid functionals (I tested HSE, B3LYP and PBE0 but that choice was entirely arbitrary and does not reflect any meaningful selection). I also tested that both, norm-conserving (ONCV) and ultrasoft/PAW pseudopotentials work at this stage and found no problem.

Below you will see the bands obtained for the ONCV case – note that this is just a test of the code’s capabilities and not converged (wiggles in the bands indicate that the wannierization failed to localize the states sufficiently which might require adjustment of the wannierization window)

And we can also quickly check if the hybrid functionals do what we expect them to do – fix the “band-gap problem”. Indeed, we find that the band-gap is increased from the PZ case (~4.1 eV) to a maximum value of 7.2 eV (using PBE0) – with literature values usually falling between 5.7 and 7.7 this is not a bad match at all 🙂

As always, use these with caution and make sure your calculations are converged! A run script that tests above-mentioned cases as well as an example of the output for HSE using ONCV pseudopotentials are available for download below.

Have fun!

Five easy ways to visualize xsf volume grids

Post-processing of the total charge density allows the visualization of many quantities such as charge density, spin-polarization, electron localization and many more (see the option plot_num in the manual). In many cases the resulting data will be a 3D grid representation (voxel grid) in either the CUBE format or XSF – I prefer the latter because I am more familiar with using it.

For this example I generated an xsf file of the spin-polarization (plot_num=6) of a molecule (I used nickelocene as an example as I had it readily at hand). By default, pp.x will generate the xsf on the entire FFT grid, which might be bigger than necessary. To avoid this, you can use “bspline” interpolation together with output_format = 3 which seems to work without any issue in QE 6.5.

The format of the xsf file is pretty self-explanatory and also detailed here. The header contains information about the cell dimensions, the coordinates of the atoms in the cell, followed by the grid specifications and the grid data “one pixel at a time”, which are written as explained in the link as follows:

` write(*,*) \$ (((value(ix,iy,iz),ix=1,nx),iy=1,ny),iz=1,nz)`

Now, let’s move on to how to visualize this data. The most “natural” way would probably be using xcrysden, which is excellent when dealing with xsf files, but I wanted to use MATLAB first to get a better feeling of how to read the xsf raw data so let’s look at a simple MATLAB script to read and visualize xsf voxel data first.

1) MATLAB (https://www.mathworks.com/

Reading in the voxel grid is rather easy; I opted to look for the keyword DATAGRID_3D_UNKNOWN first and then simply read in all following lines starting 6 lines after the beginning of the datagrid. There might be much better ways of doing this, but it was all I could come up with in 5 minutes or less. In the end we have to reshape the data which can be easily achieved by MATLAB’s own reshape routine. Isosurfaces can be constructed (and simultaneously plotted) by the isosurface() function. The script to read the xsf file and create the following plot is attached at the end of this post.

I have centered the molecule at (0/0/0) for plotting purposes and I would say the result isn’t bad for such a simple approach! If someone wants to manipulate or use the data to calculate volume overlaps or similar then MATLAB is definitely convenient, but when it comes to creating pretty plots I am not sure how much this can be improved. So let’s move to option

2) XCRYSDEN (http://www.xcrysden.org/

It does not get much easier to read and visualize xsf files than with xcrysden. Simply call xcrysden with the — xsf  switch (I am using xcrysden 1.6.2 on Ubuntu), press a to activate anti-aliasing to make it  a bit prettier and click on Tools -> Data Grid -> OK; input the isovalue, check to render positive and negative isovalues and you are ready to go; Unfortunately exporting the resulting file to an image is a bit limited and the pixel based images (png, TIFF etc.) are not the prettiest and their resolution is limited to the screen resolution. However, this is definitely a fast and convenient way.

3) VESTA (http://jp-minerals.org/vesta/en/

Slightly nicer plots can be made with next to no effort in VESTA as well. Just drag and drop the xsf file into a recent version of VESTA (I use version 3.4.4 dated March 2018) and adjust the settings under “Properties -> Isosurfaces” to your liking. What I particularly like in VESTA is the easy slicing using the “Boundary -> Cutoff planes -> new) function. To get a better resolution, export the raster image (File -> Export raster image) at larger sizes, then shrink it down.

4) OVITO (https://www.ovito.org/

I have only used this one for a short period of time compared to the others but OVITO (Open Visualization Tool) is hands down my favorite. The following plots were made with version 3.0.0 and there is also a “Pro” version available that will come with license costs in the near future. Simply drag-and-drop the xsf file into the main window, then add a “Create isosurface” modifier (Add modification -> Visualization -> Create Isosurface); do this twice, once for positive and once for negative values of the isosurface and you are ready to export the raster image with a render engine of your choice.  Here I have used Tachyon which already gives quite nice results.

The GUI offers a huge variety of customizations for atom color, radius, bond customization, slicing (Add Modification -> Modification -> Slice) making this tool really fun to work with. But if you want to go one step further you can use…

5) OVITOS

OK, I admit, this is not really another tool but rather the python scripting interface to OVITO. Nevertheless, this is an amazing addition that makes repeated tasks much less annoying. When using the script interface you can further take advance of numpy which makes “doing math” on the xsf file very easy. This approach really shines when you have to create a lot of similar plots for a publication and you want to make sure that plotting is consistent across these plots (for example atomic radii or colors). I have attached the python script used to create this plot which I hope is basically self-explanatory.

If you have any questions drop them in the comments and I will try to help you out!

Happy Plotting!

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. 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!