Skip to content

Core Levels with GW

In this tutorial, we use the \(GW\) method to calculate the core levels of a methanol molecule in the gas phase. The implementation of GW for core levels in FHI-aims is described in Golze et al. 1 2. We employ the relativistic corrections described by Keller et al. 3

Background

Core levels are typically observed in x-ray photoemission spectroscopy, in which the absorbed photon excites a core electron whose kinetic energy is observed. The absolute binding energy of the emitted quasiparticle is given as a correction to the DFT orbital energy by replacing the exchange-correlation energy \(v^{XC}\) with the correlation part of the \(G_0W_0\) approximation to the self-energy operator \(\Sigma^{G_0W_0}\). The quasiparticle energy for the \(n\)th level is given by the quasiparticle equation

\[ \epsilon_n^{G_0W_0} = \epsilon_n + Re\Sigma_{n}^{G_0W_0}(\epsilon_n^{G_0W_0}) - v^{XC}_n \]

where \(\epsilon_n\) is the Kohn-Sham eigenvalue. This is a non-linear equation, and is solved iteratively.

There are a number of core-level specific concerns that must be accounted for:

  • More accurate integration techniques are required in the core region. We use the contour-deformation approach for the self-energy integration 1

  • A DFT Generalized Gradient Approximation starting point yields spurious multisolution behavior. This can be circumvented either by using a hybrid starting point with a high fraction of exact exchange, or by iterating the Green's function calculation to eigenvalue self-consistency (eigenvalue self-consistent \(GW_0\), henceforth \(evGW_0\)). The former approach leaves the fraction of exact exchange as a tuning parameter, which for the 1s cores of the second period elements is optimized at 0.45. The latter approach is entirely ab-initio, but computationally more cumbersome. We leave the latter as an optional additional exercise. 12

  • As is generally the case with correlated methods, the correlation energy slowly converges with the size of the basis set. We perform a two-point extrapolation utilizing the cc-pVTZ and cc-pVQZ Dunning correlation consistent basis sets; see 2,4 and refs. therein.

  • Relativistic effects are more pronounced in the core-region, and must be accurately accounted for. The most common scalar relativistic approach, ZORA, is insufficiently accurate for the deep core. Fortunately, as these effects are fairly independent of the valence region, a simple atomic correction dependent only on the atomic species is sufficient, and this correction does not increase the computational effort of the calculation. 3

Exercise

In this exercise, we will perform an \(G_0W_0\) calculation based on PBEh exchange-correlation functional DFT calculation (\(G_0W_0\)@PBEh(\(\alpha =0.45\))) employing Dunning correlation-consistent Gaussian basis sets.

The quantities of interest in this case are the core-level quasiparticle energies. You can proceed as follows:

As a first step, generate the geometry.in file for methanol. Here are the fully relaxed positions from reference 2:

    atom        -0.72621163       -0.01396083        0.00081861  C
    atom         0.70098012        0.00264395       -0.00142482  O
    atom        -1.02758845       -1.06791437        0.00399808  H
    atom        -1.15279322        0.46446619       -0.89700965  H
    atom        -1.14999624        0.46912430        0.89744859  H
    atom         0.99167641        0.92622076       -0.00390182  H

For the control.in file please use the following settings:

#########################################
#            SCF Settings               #
#########################################

  xc                      pbe0 
  hybrid_xc_coeff         0.45
  relativistic            atomic_zora scalar

#########################################
#  Quasiparticle calculation Settings   #
#########################################

  qpe_calc                 gw 

#########################################
#         GW calculation Settings       #
#########################################

  anacon_type                 1
  n_anacon_par                16
  frequency_points            200
  contour_def_gw              1 2
  state_lower_limit           1
  post_adjust_qp_relativistic  1  0.3593 2 0.1001

Let us have a look at the meaning of this control.in file:

  • anacon_type 1: This chooses the Pade approximation for the analytical continuation.
  • n_anacon_par 16: This sets the order of Pade approximation
  • frequency_points 200: Number of frequency points on the imaginary axis for analytical continuation, as well as for the integration over the imaginary axis in the contour-deformation routines. The default value is 100. This parameter should be in principle converged. However, 200 frequency points is typically a very safe setting. Convergence is in most cases achieved earlier.
  • contour_def_gw 1 2: states to calculate with contour deformation (inner core states). Here, for the states 1 and 2, ordered by increasing binding energy.
  • state_lower_limit 1: Include all states in eigenvalue calculation
  • post_adjust_qp_relativistic: Post-GW relativistic correction to inner core states. Using atomic ZORA, the DFT eigenvalues and GW quasiparticle energies will be too negative compared to the fully relativistic reference, while they will be too positive when performing a non-relativistic calculation, see Figure 1 in ref 3. Since in our example we use atomic ZORA, a positive correction of 0.3593 eV is added to the quasiparticle energy of state 1 (O1s state) and positive correction of 0.1001 eV to state 2 (C1s state). The values for the relativistic corrections are given in ref. 3

Now, run the calculation using the Dunning cc-pVnZ basis sets with n=T,Q. The basis sets can be found in AIMS_DIR/species_defaults/non-standard/gaussian_tight_770/cc-pVnZ, where AIMS_DIR refers to the path of the FHI-aims root directory.

After the calculation has finished, look for the final listing of 1s quasiparticle energies for the states 1 and 2 (the 1s states of carbon and oxygen) in the output file. Only the states specified in the `control.in' file are computed with contour deformation. The other states are computed with analytic continuation. Below you see the output for the cc-pVTZ basis set.

QP_eps_output

After extracting the 1s quasiparticle energies for both calculations, extract also the number of basis functions. You will have to look for the following output block during the initialization of the calculation:

| Total number of basis functions :      116

We have used valence correlation consistent basis functions, which allows us to extrapolate to an infinite basis set and, thus, obtain results independent of the underlying basis set (also confirm with the previous chapter). For the infinite-basis-set extrapolation, use the formula (see 2,4 and references therein):

\[ IP_n = IP_n^{\infty} - \frac{a}{N_\text{basis}} \]

where \(N_{basis}\) is the number of basis functions and \(IP\) the ionization potential, which corresponds to the negative of the quasiparticle energy, i.e., \(IP_n=-\epsilon_n^{G_0W_0}\). For each core-level, make a separate plot with \(1/N_{basis}\) on the x-axis, and the IP on the y-axis. The two data points (cc-pVTZ, cc-pVQZ) determine a line. Where does this line intercept the y-axis? Shown below is an example of such an extrapolation for the highest occupied molecular orbital of a Benzene molecule. Compare to the experimental binding energies of 538.88 eV and 292.3 eV 5.

QP_BS_extrapolation

Reproduced from Golze et. al. 4

Optional additional exercise.

Perform the same calculation and extrapolation using the \(evGW_0\) approach with a PBE starting point. The settings for the control.in with comments are given here.

##########################################
#     Control.in for evGW_0@PBE          #
##########################################

#########################################
#            SCF Settings               #
#########################################

  xc                      pbe
  spin                    none
  relativistic            atomic_zora scalar

#########################################
#  Quasiparticle calculation Settings   #
#########################################

  qpe_calc           ev_scgw0 

#########################################
#         GW calculation Settings       #
#########################################


  anacon_type  1              #  Pade approximation
  n_anacon_par 16             #  order of Pade approximant
  frequency_points 200        #  Frequency integration grid for Self-Energy calculation
  contour_def_gw   1 2        #  states to calculate with contour deformation (inner core states)
  state_lower_limit 1         #  Include all states in eigenvalue calculation
  pre_adjust_qp_relativistic 1  0.3593 2 0.1001

  nocc_sc_cd 1000             # include all occupied states
  contour_eta 0.05            # broadening parameter for contour deformation.
                              # this is now raised above the default value for quicker convergence
                  # but qp energies should be converged as eta -> 0