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


Fatbands with spin-orbit coupling

A recent question touched on the case of fatband representation when spin-orbit coupling is included in the calculation (lspinorb = .true.). I assumed that this would be straightforward and indeed there is no problem up to the last step, where plotband.x should read the band energies and projections. It seems that in QE 6.5 the projection file is not properly read

 Input file > bands.dat
Reading 30 bands at 401 k-points
file with projections not compatible with bands
Range: -5.6460 34.2420eV Emin, Emax >

However the projection file is correctly written as we can identify the (j,mj) blocks such as

1    1  Pb  5D     1    2  2.5  -2.5

To circumvent this incompatibility of plotband.x one can simply read the energy eigenvalues from the bands.x outfile and relative weight of the wave-function from the projection file, where every block looks something likes this:

    9    1  Pb  5D     2    2  1.5   0.5
       1       1        0.5475316611
       1       2        0.2162196288
       1       3        0.1468317944
       1       4        0.0885236565
       1       5        0.0000000000
       1       6        0.0000000000
     401      29        0.0000000000
     401      30        0.0000000000

The first column is the k-point index, followed by the band and the weight for each mj state. If we want to plot the fatbands for all the d-states we have to identify the “block” number (i.e. the preceding 9 in the example above), collect all the d-states projections and plot them on top of the d-bands. I found a script in the usergroup and the origin seems to trace back to but unfortunately I can’t read Chinese and cannot verify that information.

And there we have it: fatband representation for the d-bands of bulk lead with spin-orbit coupling.

A script that does all of this automatically is as always included, the pseudopotential used is the fully relativistic version of the ONCV pseudopotential which seems to be giving reliable results according to SSSP, however I did not converge the calculation as this seemed precise enough for a quick demonstration.

I hope that helps!