Core ionization in the CO and N2 molecules¶
Introduction¶
We would like to study K-edge energies of CO and N2 molecules.
Carbon monoxide¶
A first approximation to the O1s and C1s ionization energies is provided by Koopmans’ theorem. Starting from the molecular input file CO.mol
INTGRL
Carbon monoxide
...
C 2 A
6. 1
C 0.0000000000000 0.0000000000000 0.0000000000000
LARGE BASIS dyall.3zp
8. 1
O 0.0000000000000 0.0000000000000 1.128323
LARGE BASIS dyall.3zp
FINISH
and the menu file CO.inp
**DIRAC
.WAVE FUNCTION
.ANALYZE
**WAVE FUNCTION
.SCF
*SCF
.CLOSED SHELL
14
**ANALYZE
.MULPOP
**END OF
we first run a Hartree-Fock calculation of the ground electronic state of carbon monoxide:
pam --inp=CO --mol=CO --get "DFCOEF=cf.CO"
We find that \(-\varepsilon (O1s)=20.6838\,E_h=562.8\) eV and \(-\varepsilon (C1s)=11.3672 \,E_h=309.3\) eV. This is significantly off the values of 542.1 eV and 295.9 eV, respectively, reported by Kai Siegbahn and co-workers [Siegbahn1969] .
Koopmans’ theorem work quite well for valence excitations because the two major sources, that is, i) lack of electron correlation and ii) lack of orbital relaxation, typically are of about the same magnitude, but of opposite sign and hence tend to cancel out. This is not the case for core ionization, where orbital relaxation is very much more important than electron correlation. To introduce orbital relaxation of the core-ionized state we carry out a \(\Delta SCF\) calculation, as pioneered by Paul Bagus and Henry F. Schaefer [Bagus_JCP1971].
In order to calculate the O1s ionized system we set up the menu file CO_O1s.inp
**DIRAC
.WAVE FUNCTION
.ANALYZE
**WAVE FUNCTION
.SCF
.REORDER
2..7,1
*SCF
.CLOSED SHELL
12
.OPEN SHELL
1
1/2
.OPENFAC
1.0
.OVLSEL
.NODYNSEL
**ANALYZE
.MULPOP
**END OF
and use the command:
pam --inp=CO_O1s --mol=CO --put "cf.CO=DFCOEF"
It is seen that we start from the coefficients of the neutral CO molecule. We then converge to the core-ionized system using reordering and overlap selection: In the input we have specified twelve electrons in six 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 non-dynamic overlap selection, which means that we order orbitals according to their overlap with the starting orbitals.
The analogous menu file for C1s ionization is CO_C1s.inp
**DIRAC
.WAVE FUNCTION
.ANALYZE
**WAVE FUNCTION
.SCF
.REORDER
1,3..7,2
*SCF
.CLOSED SHELL
12
.OPEN SHELL
1
1/2
.OPENFAC
1.0
.OVLSEL
.NODYNSEL
**ANALYZE
.MULPOP
**END OF
which puts the C1s on top of the occupied orbitals, and then use the command:
pam --inp=CO_C1s --mol=CO --put "cf.CO=DFCOEF"
The HF energies of the neutral, the O1s-ionized and the C1s-ionized systems are -112.857912 \(E_h\), -92.939954 \(E_h\) and -101.933919 \(E_h\), respectively, from which we calculate IP(O1s) = 19.9180 \(E_h\) = 542.0 eV and IP(O1s) = 10.9240 \(E_h\) = 297.3 eV, which agrees remarkably well with the experimental values, in particular for the oxygen K-edge.
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”).
Nitrogen molecule¶
We now switch to the isoelectronic nitrogen molecule. Starting from the molecular input file N2.mol
DIRAC
N2
cc-pVTZ basis
C 1 A
7. 2
N1 0.0 0.0 0.00000
N2 0.0 0.0 1.09768
LARGE BASIS dyall.3zp
FINISH
and the menu file N2.inp
**DIRAC
.WAVE FUNCTION
.ANALYZE
**GENERAL
.ACMOUT
**WAVE FUNCTION
.SCF
*SCF
.CLOSED SHELL
6 8
**ANALYZE
.MULPOP
**END OF
we run a Hartree-Fock calculation of the ground electronic state of the nitrogen molecule:
pam --inp=N2 --mol=N2 --get "DFCOEF=cf.N2 DFACMO=ac.N2"
(note that we are also getting the \(C_1\) coefficients; we shall come back to that). We find that \(-\varepsilon (1s\sigma_g)=15.6942\,E_h=427.1\) eV and \(-\varepsilon (1s\sigma_u)=15.6907 \,E_h=427.0\) eV. This is significantly off the value of 409.9 eV reported by Kai Siegbahn and co-workers [Siegbahn1969] ; the \(1s\sigma_g\) and \(1s\sigma_u\) are indistinguishable in the X-ray PES study (see also [Wight_JEPRP1972]).
We therefore proceed to \(\Delta SCF\) calculations: for the \(1s\sigma_g\) ionization we use the menu file N2_Kg.inp
**DIRAC
.WAVE FUNCTION
.ANALYZE
**GENERAL
.ACMOUT
**WAVE FUNCTION
.SCF
.REORDER
2,3,1
1,2,3,4
*SCF
.CLOSED SHELL
4 8
.OPEN SHELL
1
1/2,0
.OVLSEL
.NODYNSEL
**ANALYZE
.MULPOP
**END OF
and the command:
pam --inp=N2_Kg --mol=N2 --put "cf.N2=DFCOEF"
and equivalent menu file N2_Ku.inp for the \(1s\sigma_u\) ionization:
**DIRAC
.WAVE FUNCTION
.ANALYZE
**GENERAL
.ACMOUT
**WAVE FUNCTION
.SCF
.REORDER
1,2,3
2,3,4,1
*SCF
.CLOSED SHELL
6 6
.OPEN SHELL
1
1/0,2
.OVLSEL
.NODYNSEL
**ANALYZE
.MULPOP
**END OF
The HF energies of the neutral, the \(1s\sigma_g\) -ionized and the \(1s\sigma_u\) -ionized systems are -109.051015 \(E_h\), -93.626378 \(E_h\) and -93.630989 \(E_h\), respectively, from which we calculate IP(\(1s\sigma_g\)) = 15.4246 \(E_h\) = 419.7 eV and IP(\(1s\sigma_u\)) = 15.4200 \(E_h\) = 419.6 eV. However, this is still far from the experimental value of 409.9 eV !
Bagus and Schaefer found a similar discrepancy in their study of the \(O_2^+\) molecule [Bagus_JCP1972] and found that results improved dramatically by localization of the core hole, so let us try this. Localization breaks molecular symmetry, so we go down to \(C_1\) using the molecule input file N2_C1.mol
DIRAC
N2
cc-pVTZ basis
C 1 0 A
7. 2
N1 0.0 0.0 0.00000
N2 0.0 0.0 1.09768
LARGE BASIS dyall.3zp
FINISH
This is the reason why we saved the \(C_1\) coefficients in the initial calculation ! We next set up the menu file N2loc.inp
**DIRAC
#.WAVE FUNCTION
.ANALYZE
**WAVE FUNCTION
.SCF
*SCF
.CLOSED SHELL
14
**ANALYZE
.MULPOP
.LOCALIZE
*LOCALIZE
.THGRAD
1.0D-09
.SELECT
1,2
**END OF
for Pipek-Mezey localization. Note that we select only the a little trick: We select only the \(1s\sigma_g\) and \(1s\sigma_u\) orbitals for localization. We now localize using the command:
pam --inp=N2loc --mol=N2_C1 --put "ac.N2=DFCOEF" --get "DFCOEF=ac.N2loc"
We calculate the system with a localized core hole using the menu file N2_K.inp
**DIRAC
.WAVE FUNCTION
.ANALYZE
**GENERAL
.ACMOUT
**WAVE FUNCTION
.SCF
.REORDER
2..7,1
*SCF
.CLOSED SHELL
12
.OPEN SHELL
1
1/2
.OPENFAC
1.0
.OVLSEL
.NODYNSEL
**ANALYZE
.MULPOP
**END OF
and the command:
pam --inp=N2_K --mol=N2_C1 --put "ac.N2loc=DFCOEF"
In fact, we can do slightly better because we can go up to \(C_{2v}\) symmetry using that DIRAC allows the import of coefficients in \(C_1\) symmetry from the unformatted file DFACMO to current symmetry. This assumes that the current symmetry is lower than the symmetry used for obtaining the original coefficients, which is the case here. We now calculate the system with a localized core hole using the menu file N2_KC2v.inp
**DIRAC
.WAVE FUNCTION
.ANALYZE
**GENERAL
.ACMOIN
**WAVE FUNCTION
.SCF
.REORDER
2..7,1
*SCF
.CLOSED SHELL
12
.OPEN SHELL
1
1/2
.OPENFAC
1.0
.OVLSEL
.NODYNSEL
**ANALYZE
.MULPOP
**END OF
the molecule input file N2_C2v.mol
DIRAC
N2
cc-pVTZ basis
C 1 2 Y X A
7. 2
N1 0.0 0.0 0.00000
N2 0.0 0.0 1.09768
LARGE BASIS dyall.3zp
FINISH
and the command:
pam --inp=N2_KC2v --mol=N2_C2v --put "ac.N2loc=DFACMO"
The calculated energy of the core ionized system is -93.971190 \(E_h\) and we now obtain an ionization energy of 15.0978 \(E_h\) = 410.3 eV, in much better agreement with experiment. The price we had to pay was to break symmetry. It is possible to keep symmetry, but now at the price of going to more complex wave functions, as shown by Agren, Bagus and Roos [Agren_CPL1981].