X-Ray Ab-Initio Correction Scheme

As mentioned in X-Ray Absorption from TDDFT, XAS LR-TDDFT results need to be rigidly shifted to match experiments. This is due to self-interaction error and the lack of orbital relaxation upon the creation of the core hole. An ab-initio correction scheme was developed to address these issues. Theory and benchmarks were published in Bussy2021b. Please cite this paper if you were to use this method.

Brief theory recap

XAS LR-TDDFT yield excitation energies as correction to ground state Kohn-Sham orbital energy differences, namely:

\[ \omega = \varepsilon_a - \varepsilon_I + \Delta_{xc}, \]

where \(\varepsilon_a\) is the orbital energy of a virtual MO and \(\varepsilon_I\) the energy of the donor core MO. Under Koopman’s condition, these energies are interpreted as the electron affinity and and the ionization potential (IP). However, DFT is notoriously bad at predicting accurate absolute orbital eigenvalues. Therefore, and because \(|\varepsilon_I| >> |\varepsilon_a|\), excitation energies are expected to be widely improved if the DFT energy \(\varepsilon_I\) were to be replace by an accurate value of the IP.

The IP can be accurately calculated using the second-order electron propagator equation:

\[ \text{IP}_I = -\varepsilon_I - \frac{1}{2} \sum_{ajk}\frac{|\langle Ia||jk\rangle|^2}{-\text{IP}_I + \varepsilon_a -\varepsilon_j -\varepsilon_k} - \frac{1}{2}\sum_{abj}\frac{|\langle Ij||ab\rangle|^2}{-\text{IP}_I + \varepsilon_j - \varepsilon_a - \varepsilon_b} \]

where \(a, b\) refer to virtual Hartree-Fock spin-orbitals and \(j,k\) to occupied HF spin-orbitals. The DFT generalization of this theory is known as GW2X (Shigeta2001). It involves calculating the Generalized Fock matrix and the rotation of the occupied and virtual DFT orbitals separately, such that they become pseudocanonical. Alternatively, the diagonal elements of the generalized Fock matrix can be used as approximations for the orbital energies (thus saving on the orbital rotation). This is known as the GW2X* method.

The GW2X input subsection

The parameters defining the GW2X correction to XAS LR-TDDFT are found in the GW2X subsection of XAS_TDP. GW2X will only work with hybrid functionals (or full Hartree-Fock), as the machinery necessary for the calculation of the generalized Fock matrix is not available otherwise.

There are not many parameters to set for the GW2X correction. Simply adding an empty GW2X subsection is usually enough. The electron propagator equation is solved iteratively with a Newton-Raphson scheme. EPS_GW2X controls the convergence threshold and MAX_GW2X_ITER the maximum number of iterations allowed. The PSEUDO_CANONICAL keyword controls whether the original GW2X scheme or its simplified GW2X* version is run (by default, the original GW2X is on). C_SS and C_OS allow to scale the same- and opposite-spin components (as in SOS- and SCS-MP2). Finally, if XPS_ONLY is set, only the core IP is calculated and the XAS LR-TDDFT calculation is skipped.

Simple examples

OCS molecule (L-edge + SOC)

This example covers GW2X corrected L-edge spectroscopy with spin-orbit coupling.

&GLOBAL
  PROJECT OCS
  PRINT_LEVEL MEDIUM
  RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
  METHOD Quickstep
  &DFT
    BASIS_SET_FILE_NAME BASIS_GW2X
    POTENTIAL_FILE_NAME POTENTIAL
    AUTO_BASIS RI_XAS MEDIUM

    &MGRID
      CUTOFF 800
      REL_CUTOFF 50
      NGRIDS 5
    &END MGRID
    &QS
      METHOD GAPW
    &END QS

    &POISSON
      PERIODIC NONE
      PSOLVER MT
    &END

    &SCF
      EPS_SCF 1.0E-8
      MAX_SCF 50
    &END SCF

    &XC
      &XC_FUNCTIONAL                ! The PBEh(45%) functional
         &LIBXC
            FUNCTIONAL GGA_C_PBE
         &END LIBXC
         &LIBXC
            FUNCTIONAL GGA_X_PBE
            SCALE 0.55
         &END LIBXC
      &END XC_FUNCTIONAL

      &HF
         FRACTION 0.45
      &END HF
    &END XC

    &XAS_TDP
      &DONOR_STATES
         DEFINE_EXCITED BY_KIND
         KIND_LIST S
         STATE_TYPES 2p          ! Need to look for the S 2p states within the 7 MOs with lowest energy;
         N_SEARCH 7              ! one S 1s, one S 2s, three S 2s , one C 1s and one O 1s
         LOCALIZE                ! Localization is required
      &END DONOR_STATES

      EXCITATIONS RCS_SINGLET
      EXCITATIONS RCS_TRIPLET
      SOC

      GRID S  300 500

      N_EXCITED 150
      TAMM_DANCOFF

      &GW2X                      ! This is the only difference in the input file with respect to a
      &END GW2X                  ! standard XAS_TDP calculation (defaults parameters are used)

      &KERNEL
         RI_REGION 3.0
         &XC_FUNCTIONAL
            &LIBXC
               FUNCTIONAL GGA_C_PBE
            &END LIBXC
            &LIBXC
               FUNCTIONAL GGA_X_PBE
               SCALE 0.55
            &END LIBXC
         &END XC_FUNCTIONAL
         &EXACT_EXCHANGE
            FRACTION 0.45
         &END EXACT_EXCHANGE
      &END KERNEL

    &END XAS_TDP
  &END DFT
  &SUBSYS
    &CELL
      ABC 10.0 10.0 10.0
      PERIODIC NONE
    &END CELL
    &COORD
      C         5.0000000209        4.9999999724        5.2021372095
      O         5.0000000094        5.0000000290        6.3579624316
      S         5.0000000207        5.0000000007        3.6399034216
    &END COORD
    &KIND C
      BASIS_SET aug-pcX-2
      POTENTIAL ALL
    &END KIND
    &KIND O
      BASIS_SET aug-pcX-2
      POTENTIAL ALL
    &END KIND
    &KIND S
      BASIS_SET aug-pcX-2
      POTENTIAL ALL
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

The only difference between the above input file and that of a standard XAS LR-TDDFT calculation is the addition of the GW2X subsection. In this case, only default parameters are used, which corresponds to the original GW2X scheme with a convergence threshold of 0.01 eV. Note that the core specific all-electron aug-pcX-2 basis set is used (triple zeta quality). This inputs corresponds to an entry of table II in Bussy2021b, although slacker parameters are used here (in order to make this tutorial cheap and easy to run, this particular calculations takes ~2 minutes on 4 cores).

In the output file, the correction for each S \(2p\) is displayed. Note that the correction amounts to a shift of 1.9 eV compared to standard XAS LR-TDDFT, leading to a first singlet excitation energy of 164.4 eV (at the L\(_3\) edge). This fits experimental results within 0.1 eV. thus clearly improving the XAS LR-TDDFT result. Note that the core IPs, including spin-orbit coupling effects, are also provided. These can be directly used to produce a XPS spectrum. The content of the OCS.spectrum file yields the corrected spectrum directly.

    - GW2X correction for donor MO with spin  1 and MO index    5:
                             iteration                convergence (eV)
                                     1                       10.047536
                                     2                        1.237503
                                     3                        0.014416
                                     4                       -0.000000

      Final GW2X shift for this donor MO (eV):   1.927146


    - GW2X correction for donor MO with spin  1 and MO index    6:
                             iteration                convergence (eV)
                                     1                        6.197650
                                     2                        4.963008
                                     3                        0.241838
                                     4                        0.000439

      Final GW2X shift for this donor MO (eV):   1.907648


    - GW2X correction for donor MO with spin  1 and MO index    7:
                             iteration                convergence (eV)
                                     1                        6.197650
                                     2                        4.963008
                                     3                        0.241838
                                     4                        0.000439

      Final GW2X shift for this donor MO (eV):   1.907648


    Calculations done:

    First singlet XAS excitation energy (eV):                165.014087
    First triplet XAS excitation energy (eV):                164.681850
    First SOC XAS excitation energy (eV):                    164.396537

    Ionization potentials for XPS (GW2X + SOC):              170.602279
                                                             169.457339
                                                             169.367465

Solid NH\(_3\) (K-edge, periodic)

This is a much larger example of a periodic system, namely solid ammonia. This example is much heavier to run (~45 minutes on 24 cores).


&GLOBAL
  PROJECT NH3
  RUN_TYPE ENERGY
  PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
  METHOD QS
  &DFT
    BASIS_SET_FILE_NAME BASIS_GW2X
    BASIS_SET_FILE_NAME BASIS_ADMM
    BASIS_SET_FILE_NAME BASIS_MOLOPT
    POTENTIAL_FILE_NAME POTENTIAL
    AUTO_BASIS RI_XAS MEDIUM

    &QS
      METHOD GAPW
    &END QS

    &MGRID
      CUTOFF 600
      REL_CUTOFF 50
      NGRIDS 5
    &END MGRID

    &SCF
      SCF_GUESS RESTART
      EPS_SCF 1.0E-8
      MAX_SCF 30

      &OT
         MINIMIZER CG
         PRECONDITIONER FULL_ALL
      &END OT

      &OUTER_SCF
         MAX_SCF 6
         EPS_SCF 1.0E-8
      &END OUTER_SCF

    &END SCF

    &AUXILIARY_DENSITY_MATRIX_METHOD
      ADMM_PURIFICATION_METHOD NONE
    &END AUXILIARY_DENSITY_MATRIX_METHOD

    &XC
      &XC_FUNCTIONAL
         &LIBXC
            FUNCTIONAL GGA_X_PBE
            SCALE 0.55
         &END
         &LIBXC
            FUNCTIONAL GGA_C_PBE
         &END
      &END XC_FUNCTIONAL
      &HF
         FRACTION 0.45
         &INTERACTION_POTENTIAL
            POTENTIAL_TYPE TRUNCATED
            CUTOFF_RADIUS 5.0
         &END INTERACTION_POTENTIAL
      &END HF
    &END XC

    &XAS_TDP
      &DONOR_STATES
         DEFINE_EXCITED BY_KIND
         KIND_LIST Nx
         STATE_TYPES 1s
         N_SEARCH 1
         LOCALIZE
      &END DONOR_STATES

      TAMM_DANCOFF
      GRID Nx 300 500
      E_RANGE 30.0

      &GW2X
      &END

      &KERNEL
         &XC_FUNCTIONAL
            &LIBXC
               FUNCTIONAL GGA_X_PBE
               SCALE 0.55
            &END
            &LIBXC
               FUNCTIONAL GGA_C_PBE
            &END
         &END XC_FUNCTIONAL
         &EXACT_EXCHANGE
            OPERATOR TRUNCATED
            CUTOFF_RADIUS 5.0
            FRACTION 0.45
         &END EXACT_EXCHANGE
      &END KERNEL

    &END XAS_TDP
  &END DFT
  &SUBSYS
    &CELL
      ABC   10.016118  10.016118  10.016118
    &END CELL
    &TOPOLOGY
      COORD_FILE_FORMAT XYZ
      COORD_FILE_NAME NH3.xyz
    &END TOPOLOGY
    &KIND H
      BASIS_SET DZVP-MOLOPT-SR-GTH
      BASIS_SET AUX_FIT FIT3
      POTENTIAL GTH-PBE
    &END KIND
    &KIND N
      BASIS_SET DZVP-MOLOPT-SR-GTH
      BASIS_SET AUX_FIT FIT3
      POTENTIAL GTH-PBE
    &END KIND
    &KIND Nx
      ELEMENT N
      BASIS_SET aug-pcseg-2
      BASIS_SET AUX_FIT aug-admm-2
      POTENTIAL ALL
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

Again, the only difference with respect to a standard XAS-LRTDDFT input file is the GW2X subsection. This input file corresponds exactly to figure 3 a) of Bussy2021b. In this case, the GW2X correction amounts to a blue shift of 3.7 eV, aligning the calculated spectrum to the experimental one remarkably well. The NH3 structure file as well as the necessary basis set file (also for the OCS example) are available here.

FAQ

How can I make the GW2X correction run faster ?

The GW2X correction scheme scales cubically with the number of MOs in the system. Therefore, the best way to improve performance is to reduce that number. Because an accurate description of the core region is only necessary for the exited atoms, all other atoms can be described with pseudopotentials. This drastically reduces the number of MOs since only valence states are kept. In the solid NH3 example above, all nitrogen atoms are equivalent under symmetry. Therefore, their individual contribution to the XAS spectrum is bound to be the same. This allows for the description of a single nitrogen atom at the all-electron level, while all others (and the hydrogens) use pseudopotentials. Note that the ADMM approximation is also utilized. This greatly reduces the cost of the underlying hybrid DFT calculation, as well as the evaluation of the generalized Fock matrix as required by GW2X.

Why don’t I get the absolute core IP in periodic systems ?

For molecules in non-periodic boundary conditions, the potential is such that it is zero far away. In the periodic case, the zero is ill defined. As a consequence, all Kohn-Sham eigenvalues end up shifted by some unknown, constant amount. Therefore, their absolute values and that of the calculated IP cannot be interpreted in a physical manner. However, the correction scheme depends on the difference \(|\varepsilon_a-\varepsilon_I|\), where the shift cancels out.

Why is the LOCALIZE keyword required ?

In order to efficiently evaluate the antisymmetric integrals of the type \(\langle Ia || jk \rangle\), the same local RI scheme as XAS_TDP is used. Therefore, the core state \(I\) needs to be local in space. However, the rotation required to get the pseudocanonical orbitals needed for the original GW2X scheme may break this localization, provided that there are other equivalent atoms in the system. To prevent that from happening, all core states localized on other atoms are ignored for the rotation and the subsequent IP calculation. This has negligible impact since core states belonging to different atoms only weakly interact. It is however important to keep the value of the LOCALIZE keyword to a minimum to insure that only core states are ignored.