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