Core electron excitations and ionization in H2O at the HF and DFT levels¶
Introduction¶
We want to study excitation of the oxygen 1s electron of water.
Restricted excitation window TDHF (REWTDHF)¶
Starting from the molecular input file H2O.mol
INTGRL
Water
ccpV5Z basis (note that a tight p has NOT been added)
C 2
8. 1
O .0000000000 0.0000000000 .2249058930
LARGE BASIS dyall.3zp
1. 2
H1 1.4523499293 .0000000000 .8996235720
H2 1.4523499293 .0000000000 .8996235720
LARGE BASIS dyall.3zp
FINISH
and the menu file H2O.inp
**DIRAC
.TITLE
H2O
.WAVE FUNCTIONS
.ANALYZE
**HAMILTONIAN
.X2C
**INTEGRALS
*TWOINT
.SCREEN
1.0D12
*READIN
.UNCONT
**WAVE FUNCTIONS
.SCF
*SCF
.CLOSED SHELL
10
**ANALYZE
.MULPOP
*MULPOP
.VECPOP
1..oo
*END OF
we first run a HartreeFock calculation of the ground electronic state:
pam mol=H2O inp=H2O outcmo
From Mulliken population analysis we confirm that the oxygen 1s orbital is the first orbital:
* Electronic eigenvalue no. 1: 20.581089118212 (Occupation : f = 1.0000)
==========================================================================================
* Gross populations greater than 0.00010
Gross Total  L A1 O s

alpha 1.0000  1.0000
beta 0.0000  0.0000
We now focus on electric dipole allowed transitions. The molecular symmetry is \(C_{2v}\) where the components of the electric dipole operator \(e(x,y,z)\) span irreps \((B_1,B_2,A_1)\). This leads us to set up the following input
**DIRAC
.TITLE
H2O
.PROPERTIES
**HAMILTONIAN
.X2C
**INTEGRALS
*TWOINT
.SCREEN
1.0D12
*READIN
.UNCONT
**WAVE FUNCTIONS
.SCF
*SCF
.CLOSED SHELL
10
**PROPERTIES
*EXCITATION ENERGIES
.OCCUP
1
.EXCITATION
1 10
.EXCITATION
2 10
.EXCITATION
3 10
.INTENS
0
.ANALYZE
*END OF
We run our calculation:
pam inp=H2O_O1s mol=H2O incmo
and find the isotropically averaged oscillator strenghts
ISOTROPIC DLDL OSCILLATOR STRENGTHS (f)
==========================================
DL = dipole length
Rate = Dipole radiation rate (s1)
Lifetime = corresponding radiation lifetime (s)
Level Frequency f Rate Lifetime

2 20.1880843814 0.00000002 0.20483E+06 0.48822E05 6(E1 ) 45.42%; 11(E1 ) 30.51%
3 20.1880853618 0.00000012 0.16192E+07 0.61758E06 6(E1 ) 45.42%; 11(E1 ) 30.51%
4 20.2106058884 0.00000020 0.26305E+07 0.38016E06 12(E1 ) 32.69%; 7(E1 ) 30.96%; 9(E1 ) 19.02%
6 20.2637822952 0.04094463 0.54019E+12 0.18512E11 6(E1 ) 56.01%; 11(E1 ) 27.09%
7 20.2840331058 0.07566515 0.10003E+13 0.99973E12 7(E1 ) 38.54%; 12(E1 ) 28.96%; 9(E1 ) 20.02%
9 20.4453126149 0.00000002 0.28737E+06 0.34799E05 10(E1 ) 87.89%
10 20.4725521853 0.00000026 0.35482E+07 0.28184E06 11(E1 ) 49.36%; 6(E1 ) 22.09%
12 20.4935960507 0.04949243 0.66786E+12 0.14973E11 10(E1 ) 91.75%
13 20.5086683295 0.03037128 0.41044E+12 0.24364E11 11(E1 ) 58.52%; 6(E1 ) 23.96%
14 20.5624885103 0.00000003 0.44279E+06 0.22584E05 7(E1 ) 55.82%; 9(E1 ) 37.98%
15 20.5624885107 0.00000002 0.33427E+06 0.29916E05 7(E1 ) 55.82%; 9(E1 ) 37.98%
16 20.5662290804 0.00436872 0.59371E+11 0.16843E10 7(E1 ) 48.37%; 9(E1 ) 45.36%
17 20.5701122256 0.00000001 0.13954E+06 0.71663E05 8(E1 ) 63.89%; 6(E1 ) 28.82%
18 20.5701122415 0.00000005 0.64527E+06 0.15497E05 8(E1 ) 63.89%; 6(E1 ) 28.82%
19 20.5851294613 0.00000477 0.64894E+08 0.15410E07 8(E1 ) 77.46%; 6(E1 ) 17.30%
22 20.6853090865 0.01353373 0.18606E+12 0.53746E11 12(E1 ) 56.08%; 9(E1 ) 27.59%; 7(E1 ) 12.87%
25 20.8084949905 0.00056602 0.78746E+10 0.12699E09 13(E1 ) 95.28%
31 20.8954525590 0.00057906 0.81233E+10 0.12310E09 15(E1 ) 92.81%

Sum of oscillator strenghts: 0.21553
where we have added by hand the most important virtual orbitals. For reference, orbital densities of the fifteen first canonical HF orbitals of water are given below
We may simulating the spectrum using a Lorentzian lineshape with halfwidth at halfmaximum (HWHM) equal to 0.005 a.u.
Complex response¶
We can alternatively obtain the above spectrum from complex response. The isotropically averaged ocillator strength associated within the electric dipole approximation is related to the isotropic electric dipole polarizability through the relation
where \(\gamma\) is a damping parameter corresponding to the halfwidth at halfmaximum (HWHM) of the Lorentzian lineshape.
We can calculate the real and imaginary parts of the polarizability by complex response. By choosing a frequency window we can directly access the region of the spectrum of interest.
We use the input file
**DIRAC
.TITLE
H2O
.PROPERTIES
**HAMILTONIAN
.X2C
**INTEGRALS
*TWOINT
.SCREEN
1.0D12
*READIN
.UNCONT
**WAVE FUNCTIONS
.SCF
*SCF
.CLOSED SHELL
10
**PROPERTIES
.POLARI
*LINEAR RESPONSE
.FREQ INTERVAL
20.0 21.0 0.02
.DAMPING
0.005
*END OF
and the command:
pam incmo mol=H2O inp=cpp
To get enough data points we do a second run using:
.FREQ INTERVAL
20.01 21.0 0.02
Extracting the imaginary part of the frequencydependent electric dipole polarizatibility and converting to oscillator strength we can directly plot the spectrum as below
In the above graph we have also included the results from REWTDDFT and they completely overlap, except that a keen eyemay note that REWTDDFT is missing the peak around 570 eV, simply because not enough excitations were specified in the input.
Localizing the K edge¶
An interesting question is how to localize the K edge in the XANES spectrum. A first approximation to the 1s ionization energy is provided by Koopmans’ theorem. We find \(IP_{1s}\approx \varepsilon_{1s} =\) 20.581089 \(E_h\) = 560.04 eV. This is singificantly off the value 539.7 eV reported by Kai Siegbahn and coworkers [ESCA applied to free molecules (1969)] [Siegbahn1969] ( see also here , but here the work function of the reference metal must be subtracted). reported by Kai Siegbahn and coworkers in 1977(?). We know that Koopman’s theorem ignores correlation and orbital relaxation. For valence ionization Koopman’s theorem often provides a reasonable approximation since the errors tend to cancel each other. For core excitations orbital relaxation dominates such that Koopman’s theorem greatly overestimates ionization energies.
To see this we carry our a averageofconfiguration (AOC) calculation of the 1s coreionized system. Starting from the coefficients from the neutral system and the input
**DIRAC
.TITLE
H2O
.WAVE FUNCTIONS
.ANALYZE
**HAMILTONIAN
.X2C
**INTEGRALS
*TWOINT
.SCREEN
1.0D12
*READIN
.UNCONT
**WAVE FUNCTIONS
.SCF
.REORDER
2..5,1
*SCF
.CLOSED SHELL
8
.OPEN SHELL
1
1/2
.OVLSEL
.NODYNSEL
**ANALYZE
.MULPOP
*MULPOP
.VECPOP
1..oo
*END OF
and the command:
pam incmo mol=H2O inp=H2O_1s
we find a total energy of 56.287847 \(E_h\) for the core ionized system compared to 76.115149 \(E_h\) for the neutral system. This corresponds to a \(\Delta\) SCF value of 539.52 eV for the ionization energy, tantalizingly close to experiment.
In passing we that note that in the first iteration of this calculation we obtain a total energy of 55.534060 \(E_h\). This is the energy of the coreionized system obtained using the orbitals of the neutral system. Koopman’s theorem is obtained by subtracting this energy from the energy of the neutral system. For the neutral system we obtained 76.115149 \(E_h\), so by taking the difference we obtain 20.581089 \(E_h\), which is exactly the 1s orbital energy in the neutral system.
You should also note that we easily converge to the coreionized system using reordering and overlap selection: In the input we have specified eight electrons in four inactive (closed) orbitals, followed by a single electron in an active (open) orbital. We want the O1s to be the active orbital, but this is not achieved automatically since DIRAC will normally order orbitals according to their energy. We therefore start be reordering the orbitals such that the O1s orbital from the previous calculation on the neutral system comes out on top of the occupied orbitals. However, this is not enough to converge to the desired state since after the first diagonalization DIRAC will again by default order orbitals according to their energy. This is why we use overlap selection, that is, we ask DIRAC to rather order orbitals according to their overlap with some reference orbitals. By default (dynamic overlap selection) this will be the orbitals from the previous iteration. However, in this case we activate nondynamic overlap selection, which means that we order orbitals according to their overlap with the starting orbitals.
Note
Overlap selection is nowadays marketed hard as MOM (Maximum Orbital Method, see [Gilbert_JPCA2008]), but this method has been included in DIRAC for at least two decades and goes back to the pioneering work of Paul Bagus It was used in [Bagus_JCP1971], but not reported explicitly. However, it is for instance documented in the 1975 manual of the ALCHEMY program (On pdf page 15 you find a description of keyword MOORDR using a “maximum overlap criterion”).
The K edge obtained by \(\Delta SCF\) does not correspond to that of our TDHF or complex reponse calculations since they only allow linear reponse (orbital relaxation).
In the figure below we have plotted oscillator strength per atom for 1s core excitation spectrum for water, taken from the Gas Phase Core Excitation Database and recorded at 0.7 eV fwhm.
The vertical orange line corresponds to the O1s binding energy reported by Siegbahn and coworkers and seems to be too early. We have also plotted the spectrum obtained by REWTDHF with a Lorentzian wideshape corresponding to the fwhm of the experiment. In green we plot the original spectrum, whereas in red it has been shifted so that our linear response estimate for the ionization energy has been aligned with the experimental O1s binding energy. We can see that the shift improves agreement, but the spacing and relative intensities of peaks do not agree with experiment.
Static Exchange Approximation¶
In order to incorporate orbital relaxation we carry out a STEX calculation. At the moment symmetry is not implemented, so we turn off symmetry in the nolecular input file
INTGRL
Water
dyall.3zp basis
C 2 0
8. 1
O .0000000000 0.0000000000 .2249058930
LARGE BASIS dyall.3zp
1. 2
H1 1.4523499293 .0000000000 .8996235720
H2 1.4523499293 .0000000000 .8996235720
LARGE BASIS dyall.3zp
FINISH
We then first run the ground state:
pam inp=H2O mol=H2O_C1 outcmo
followed by the 1s coreionized state:
pam mol=H2O_C1 inp=H2O_1s incmo get "DFCOEF=DFCOEF.ION"
We now run STEX using the input
**DIRAC
.TITLE
H2O
.WAVE FUNCTIONS
.PROPERTIES
**HAMILTONIAN
.X2C
**INTEGRALS
*TWOINT
.SCREEN
1.0D12
*READIN
.UNCONT
**WAVE FUNCTIONS
.SCF
*SCF
.CLOSED SHELL
10
**PROPERTIES
*STEX
.HOLES
1
5
*END OF
and the command:
pam inp=stex mol=H2O_C1 put "DFCOEF DFCOEF.ION"
We have plotted the STEX spectrum below together with the experimental one
Switching to DFT¶
Let us now look at what we can do with DFT. The first thing to note is that Koopman’s theorem does not hold:
HF(AOC) 
HF(focc) 
LDA 
PBE 
PBE0 

Neutral 
Energy 
76.115149 
76.115149 
75.962033 
76.438943 
76.437412 
\(\varepsilon_{1s}\) 
20.581089 
20.581089 
18.620721 
18.766629 
19.223080 

\(1s^{1}\) 
Energy 
56.287847 
55.070646 
55.980781 
56.300366 
56.069079 
Energy(0) 
55.534060 
54.347068 
55.231846 
55.540835 
55.319556 

\(\Delta E\) 
relax 
19.827303 
21.044503 
19.981252 
20.138577 
20.368333 
norelax 
20.581089 
21.768082 
20.730186 
20.898108 
21.117856 
BP86 
BLYP 
B3LYP 
CAMB3LYP 

Neutral 
Energy 
76.525093 
76.508410 
76.487092 
76.496078 
\(\varepsilon_{1s}\) 
18.786866 
18.791935 
19.145210 
19.219699 

\(1s^{1}\) 
Energy 
56.366800 
56.340103 
56.148003 
56.115449 
Energy(0) 
55.606191 
55.582410 
55.398934 
55.366886 

\(\Delta E\) 
relax 
20.158293 
20.168307 
20.339089 
20.380629 
norelax 
20.918902 
20.926000 
21.088158 
21.129191 
In the above table the entry Energy(0) is the total energy of the coreionized system calculated using the orbitals of the neutral system. When we calculate the energy difference between the neutral and coreionized system using the orbitals of the neutral system we reproduce the 1s orbital energy to the cited decimals, but this is not the case of any other method, including HF using fractional occupation.
Below we give the \(\Delta SCF\) numbers is eV. Interestingly AOCHF at539.5 eV easily comes closest to the experimental value 539.7 eV. It can furthermore be seen that \(\Delta SCF\) values obtained with DFT functionals show the trend LDA < GGA < hybrid with LDA closest, but still far from experiment. Switching from averageofconfiguration to fractional occupation at the HF level leads to a dramatic deterioration of the agreement with experiment.
HF(AOC) 
HF(focc) 
LDA 
PBE 
PBE0 
BP86 
BLYP 
B3LYP 
CAMB3LYP 

\(\varepsilon_{1s}\) 
560.0 
560.0 
506.7 
510.7 
523.1 
511.2 
511.4 
521.0 
523.0 
\(\Delta E\) (relax) 
539.5 
572.7 
564.1 
568.7 
574.6 
569.2 
569.4 
573.8 
575.0 
\(\Delta E\) (no relax) 
560.0 
592.3 
543.7 
548.0 
554.3 
548.5 
548.8 
553.5 
554.6 
At the DFT level we have employed an energy expression based on fractional occupation, which is the model that leads to Janak’s theorem, namely that the derivative of the energy with respect to occupation number \(n_i\) gives the energy of the corresponding orbital, that is
We may investigate Janak’s theorem numerically. We set up a script to do DFT calculations with fractional occupation from 1.0 to 1.9 of the oxygen 1s orbital
#!/bin/bash
for dft in LDA BLYP B3LYP CAMB3LYP PBE PBE0 BP86
do
cat <<EOF > H2O_dft.inp
**DIRAC
.TITLE
H2O
.WAVE FUNCTIONS
.ANALYZE
**HAMILTONIAN
.X2C
.DFT
${dft}
**INTEGRALS
*TWOINT
.SCREEN
1.0D12
*READIN
.UNCONT
**WAVE FUNCTIONS
.SCF
*SCF
.CLOSED SHELL
10
**ANALYZE
.MULPOP
*MULPOP
.VECPOP
1..oo
*END OF
EOF
pam inp=H2O_dft mol=H2O outcmo suffix=${dft}
for occ in 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9
do
cat <<EOF > H2O_1s_dft.inp
**DIRAC
.TITLE
H2O
.WAVE FUNCTIONS
.ANALYZE
**HAMILTONIAN
.X2C
.DFT
${dft}
**INTEGRALS
*TWOINT
.SCREEN
1.0D12
*READIN
.UNCONT
**WAVE FUNCTIONS
.SCF
.REORDER
2..5,1
*SCF
.CLOSED SHELL
8
.OPEN SHELL
1
${occ}/2
.OVLSEL
.NODYNSEL
.FOCC
**ANALYZE
.MULPOP
*MULPOP
.VECPOP
1..oo
*END OF
EOF
pam inp=H2O_1s mol=H2O incmo suffix=${dft}_occ${occ}
done
done
exit 0
and also add the energy of the neutral system to our data set. We then carry out polynomial fits to various orders and calculate the derivative of the energy at occupation 2.0 with respect to 1s occupatio number. We then obtain
HF(AOC) 
HF(focc) 
LDA 
PBE 
PBE0 
BP86 
BLYP 
B3LYP 
CAMB3LYP 

1 
19.824314 
21.044193 
19.981185 
20.138613 
20.368425 
20.158286 
20.168248 
20.339084 
20.380650 
2 
18.196672 
20.579045 
18.618588 
18.765178 
19.222619 
18.785047 
18.789831 
19.143977 
19.218655 
3 
18.220409 
20.581507 
18.619168 
18.764949 
19.221929 
18.785159 
18.790354 
19.144055 
19.218529 
4 
18.220112 
20.581073 
18.620943 
18.766866 
19.223251 
18.787109 
18.792168 
19.145388 
19.219883 
5 
18.220135 
20.581093 
18.620699 
18.766601 
19.223061 
18.786846 
18.791909 
19.145191 
19.219680 
6 
18.220131 
20.581090 
18.620724 
18.766633 
19.223084 
18.786872 
18.791940 
19.145213 
19.219703 
7 
18.220130 
20.581089 
18.620720 
18.766628 
19.223080 
18.786869 
18.791934 
19.145209 
19.219699 
8 
18.220130 
20.581089 
18.620721 
18.766629 
19.223080 
18.786857 
18.791935 
19.145210 
19.219699 
9 
18.220130 
20.581089 
18.620721 
18.766629 
19.223080 
18.786879 
18.791935 
19.145210 
19.219699 
\(\varepsilon_{1s}\) 
20.581089 
20.581089 
18.620721 
18.766629 
19.223080 
18.786866 
18.791935 
19.145210 
19.219699 
We see that a linear fit is clearly is insufficient, whereas a quadratic fit is reasonable. However, a 8th order fit is in general needed to converge the energy derivative to the 1s orbital energy with the cited number of decimals.