Exercise Week 10: Photoemission

Introduction

This exercise is about photoemission topic. We will use the Delta SCF and GW methods to calculate electronic excitation energies. Both approaches are possbile to run with CP2K.

As it was discussed in the lecture, the ionization energy (or negative HOMO) can be evaluated with the delta-self-consistent-field (∆SCF) approach. Following the definition of the ionization potential, \(I = E_{tot}^{PBE}(N-1) - E_{tot}^{PBE}(N)\), the total energy difference between the neutral E_tot(N) and positively E_tot(N-1) charged species are computed from two separate DFT PBE (or HF, PBE0) total energy calculations, where N is the number of electrons of the neutral molecule.

An improved description of charged electronic excitations is obtained by the perturbative inclusion of many-body effects through the self-energy Σ via GW approximation. More specially, we will use G₀W₀ approximation in this exercise.

The purpose of this exercise is to simulate H₂O and C₂H₄'s HOMO, HOMO-1 and HOMO-2 energies.

delta SCF input file example

@SET LIBDIR /scratch/CP2K/cp2k/data
&FORCE_EVAL
  METHOD Quickstep
  &DFT
    LSD
    BASIS_SET_FILE_NAME ./BASIS_def2_QZVP_TZVP
    POTENTIAL_FILE_NAME ${LIBDIR}/POTENTIAL

    ! for neutral system, set CHARGE to 1 for charged system
    CHARGE 0
    &MGRID
      CUTOFF 400
     # NGRIDS 5
      REL_CUTOFF 50.
    &END MGRID
    &POISSON
       PERIODIC NONE
       POISSON_SOLVER MT
    &END
    &QS
      METHOD GAPW
    &END QS
    &SCF
      MAX_SCF 300
      SCF_GUESS ATOMIC
      EPS_SCF 1.0E-6
      &DIAGONALIZATION
      &END
    &END SCF
    &XC
      &XC_FUNCTIONAL
       &PBE
         ! set 75% GGA exchange for PBE0, set to 1.0 for PBE, set to 0.0 for HF
         SCALE_X 0.75
         ! set 100% GGA correlation for both PBE0 and PBE, set to 0.0 for HF
         SCALE_C 1.0
       &END PBE
      &END XC_FUNCTIONAL
      &HF
        ! set 25% HF exchange for PBE0, set to 0.0 for PBE, set to 1.0 for HF
        FRACTION 0.25
        &SCREENING
          EPS_SCHWARZ 1.0E-6
        &END
      &END
    &END XC
    &PRINT
     &MO
       EIGENVALUES
       MO_RANGE 1 50
       OCCUPATION_NUMBERS
     &END
    &END
  &END DFT

  &SUBSYS
    &CELL
      ABC [angstrom]  10 10 10
      PERIODIC NONE
    &END CELL

    &COORD
      O  0.00000000  -0.00000000  -0.00614048
      H  0.76443318  -0.00000000  0.58917024
      H  -0.76443318  0.00000000  0.58917024
    &END
!    &COORD
!      C  0.00000000  -0.00000000  -0.00156818
!      C  0.00000000  -0.00000000  1.33056818
!      H  0.92757992  0.00000000  -0.57540585
!      H  -0.92757992  0.00000000  -0.57540585
!      H  0.92757992  0.00000000  1.90440585
!      H  -0.92757992  0.00000000  1.90440585
!    &END COORD

    &TOPOLOGY
      &CENTER_COORDINATES
      &END
    &END TOPOLOGY

    &KIND H
      BASIS_SET  def2-QZVP
      POTENTIAL ALL
    &END KIND

    &KIND O
      BASIS_SET  def2-QZVP
      POTENTIAL ALL
    &END KIND

    &KIND C
      BASIS_SET  def2-QZVP
      POTENTIAL ALL
    &END KIND

  &END SUBSYS
&END FORCE_EVAL
&GLOBAL
  PROJECT cp2k_job
  RUN_TYPE energy
&END GLOBAL

A few comments about this input file.

Change CHARGE 0 to CHARGE 1 if you calculate charged molecule with N-1 electrons.
The coordinates of both H₂O and C₂H₄ are given in the above example. Comment out the water one when you deal with C₂H₄.
The basis set file we use is BASIS_def2_QZVP_TZVP. Please download it from GitHub. Put it in the same directory as your cp2k input file. You could find the def2-QZVP and def2-TZVP basis sets from Basis Set Exchange website. Here I download them and merge them into one file for convenience.
In this exercise we will test different exchange correlation functionals, i.e. PBE0, PBE and HF’s exchange contribution only. We set them in the &XC_FUNCTIONAL and &HF subsections. If you use PBE0, it means there are 75% PBE exchange + 25% HF exchange + 100% PBE correlation. If you use PBE, it means there are 100% PBE exchange + 100% PBE correlation. If you use HF, it means there are 100% HF exchange only.

GW input file example

@SET LIBDIR /scratch/CP2K/cp2k/data
&FORCE_EVAL
  METHOD Quickstep
  &DFT
    BASIS_SET_FILE_NAME ./BASIS_def2_QZVP_TZVP
    POTENTIAL_FILE_NAME ${LIBDIR}/POTENTIAL
    &MGRID
      CUTOFF 400
      REL_CUTOFF 50
    &END MGRID
    &QS
      ! all electron calculation since GW100 is all-electron test
      METHOD GAPW
    &END QS
    &POISSON
      PERIODIC NONE
      PSOLVER MT
    &END
    &SCF
      EPS_SCF 1.0E-6
      SCF_GUESS ATOMIC
      MAX_SCF 200
    &END SCF
    &XC
      &XC_FUNCTIONAL
      ! PBE section is for PBE and PBE0, HF doesn't need it
       &PBE
         ! 75% GGA exchange is for PBE0
         SCALE_X 0.75
         ! 100% GGA exchange is for PBE
         ! SCALE_X 1.0
         ! 100% GGA correlation
         SCALE_C 1.0
       &END PBE
      &END XC_FUNCTIONAL
      ! HF section is for PBE0 and HF, PBE doesn't need it
      &HF
        ! 25 % HFX exchange is for PBE0
        FRACTION 0.25
        ! 100 % HFX exchange is for HF
        ! FRACTION 1.0
        ! Important to improve scaling from O(N^4) to O(N)
        &SCREENING
          EPS_SCHWARZ 1.0E-6
        &END
      &END
      !&XC_FUNCTIONAL PBE
      !&END XC_FUNCTIONAL
      ! GW is part of the WF_CORRELATION section
      &WF_CORRELATION
        ! RPA is used to compute the density response function
        METHOD RI_RPA_GPW
        ! Use Obara-Saika integrals instead of GPW integrals
        ! since OS is much faster
        !&INTEGRALS
        ERI_METHOD OS
        !&END INTEGRALS
        &RI_RPA
          ! use 100 quadrature points to perform the
          ! frequency integration in GW
          RPA_NUM_QUAD_POINTS 100
          ! SIZE_FREQ_INTEG_GROUP is a group size for parallelization and
          ! should be increased for large calculations to prevent out of memory.
          ! maximum for SIZE_FREQ_INTEG_GROUP is the number of MPI tasks
          SIZE_FREQ_INTEG_GROUP 1
          GW
          &RI_G0W0
          !&GW
           ! compute the G0W0@PBE energy of HOMO-9,
           ! HOMO-8, ... , HOMO-1, HOMO
           CORR_OCC   10
           ! compute the G0W0@PBE energy of LUMO,
           ! LUMO+1, ... , LUMO+20
           CORR_VIRT  20
           ! fit a Pade approximant to the correlation self-energy
           ! as function of imaginary frequency. this has been done
           ! in the GW100 benchmark set and turned out to be reliable
           ANALYTIC_CONTINUATION PADE
           ! for solving the quasiparticle equation, the Newton method
           ! is used as in the GW100 benchmark
           CROSSING_SEARCH NEWTON
           ! use the RI approximation for the exchange part of the self-energy
           RI_SIGMA_X
          &END RI_G0W0
          !&END GW
        &END RI_RPA
        ! NUMBER_PROC is a group size for parallelization and should
        ! be increased for large calculations
        NUMBER_PROC 1
      &END
    &END XC
  &END DFT
  &SUBSYS
    &CELL
      ABC 10.0 10.0 10.0
      PERIODIC NONE
    &END CELL
    &COORD
      O  0.0000  0.0000  0.1173
      H  0.0000  0.7572 -0.4692
      H  0.0000 -0.7572 -0.4692
    &END COORD
!    &COORD
!      C  0.00000000  -0.00000000  -0.00156818
!      C  0.00000000  -0.00000000  1.33056818
!      H  0.92757992  0.00000000  -0.57540585
!      H  -0.92757992  0.00000000  -0.57540585
!      H  0.92757992  0.00000000  1.90440585
!      H  -0.92757992  0.00000000  1.90440585
!    &END COORD
    &TOPOLOGY
      &CENTER_COORDINATES
      &END
    &END TOPOLOGY
    &KIND H
      ! def2-QZVP is the basis which has been used in the GW100 paper
      BASIS_SET def2-TZVP
      !BASIS_SET def2-QZVP
      POTENTIAL ALL
    &END KIND
    &KIND O
      BASIS_SET def2-TZVP
      !BASIS_SET def2-QZVP
      POTENTIAL ALL
    &END KIND
    &KIND C
      BASIS_SET def2-TZVP
      !BASIS_SET def2-QZVP
      POTENTIAL ALL
    &END KIND
  &END SUBSYS
&END FORCE_EVAL
&GLOBAL
  RUN_TYPE     ENERGY
  PROJECT      ALL_ELEC
  PRINT_LEVEL  MEDIUM
&END GLOBAL

This file goes similar as delta SCF one, but with an additional &WF_CORRELATION part.

When you do a G₀W₀ calculation, please use basis set extrapolation to find the energies corresponding to infinite size of basis sets. In this exercise, you could use def2-TZVP and def2-QZVP to do the extrapolation. You could find the number of basis sets in the cp2k output file. The basis set extrapolation is like Figure 16 B in this review.
The reference experimental data for H₂O; HOMO: 12.615 eV, HOMO-1: 14.729 eV, HOMO-2: 18.550 eV.
The reference experimental data for C₂H₄; HOMO: 10.68 eV, HOMO-1: 12.80 eV, HOMO-2: 14.80 eV.

Calculate the HOMO energes of H₂O with delta SCF(PBE), G₀W₀@PBE0, G₀W₀@PBE and G₀W₀@HF methods, as well as purely DFT (PBE) and HF. You don’t need to run the last two (i.e. DFT and HF) independently, these values are already in the G₀W₀@PBE and G₀W₀@HF calculation output files.
Take the HOMO-1, HOMO-2 energies of H₂O with G₀W₀@PBE0, G₀W₀@PBE and G₀W₀@HF methods. We omit delta SCF method because it is not suitable for states other than HOMO. These values are in the output files of step 1.
Compare these data with experimental data. Which method gives the closest value to the experimental data? In brief, this exercise step 1 to 3 expects you to fill the table like this one. (10P)
(optional) Doing steps 1-3 with C₂H₄. (5P)