# Spectrum calculation#

This section focuses on the calculation of XPS, XAS, XES, and RIXS, mainly by considering these spectroscopies at the ADC and TDDFT levels of theory. For a more general overview of the large number of methods available to determine these properties, see, e.g., [ND18].

import copy
import gator
import matplotlib.pyplot as plt
import numpy as np
import veloxchem as vlx
from matplotlib import gridspec

# QC software suites and functionalities
from pyscf import cc, gto, mp, scf
from scipy.interpolate import interp1d
from veloxchem.rsplinabscross import LinearAbsorptionCrossSection

# for vlx
silent_ostream = vlx.OutputStream(None)
from mpi4py import MPI
comm = MPI.COMM_WORLD

# au to eV conversion factor
au2ev = 27.211386

def lorentzian(x, y, xmin, xmax, xstep, gamma):
"""

Call: xi,yi = lorentzian(energies, intensities, start energy,
end energy, energy step, gamma)
"""
xi = np.arange(xmin, xmax, xstep)
yi = np.zeros(len(xi))
for i in range(len(xi)):
for k in range(len(x)):
yi[i] = (
yi[i]
+ y[k] * gamma / ((xi[i] - x[k]) ** 2 + (gamma / 2.0) ** 2) / np.pi
)
return xi, yi

def gaussian(x, y, xmin, xmax, xstep, sigma):
"""

Call: xi,yi = gaussian(energies, intensities, start energy, end energy, energy step, gamma)
"""
xi = np.arange(xmin, xmax, xstep)
yi = np.zeros(len(xi))
for i in range(len(xi)):
for k in range(len(y)):
yi[i] = yi[i] + y[k] * np.e ** (-((xi[i] - x[k]) ** 2) / (2 * sigma ** 2))
return xi, yi


## X-ray photoemission spectroscopy#

In XPS the photoemission of electrons is measured, yielding information on core-electron binding energies. Here we focus on the calculation of ionization energies (IEs), from which the isotropic photoemission spectrum is formed by giving each IE the same intensity and then broadening with an appropriate convolution function. It is possible to go beyond such a simple approach using Dyson orbitals, but it typically does not make much of a difference [VKC19].

### IEs from Koopmans’ theorem#

A simple model for calculating ionization energies and electron affinities is by using Koopmans’ theorem, where they are estimated from the MO energies. For water, we then obtain a $$1s$$ IE of:

water_xyz = """
O       0.0000000000     0.0000000000     0.1178336003
H      -0.7595754146    -0.0000000000    -0.4713344012
H       0.7595754146     0.0000000000    -0.4713344012
"""

# Create veloxchem mol and basis objects

# Perform SCF calculation
scf_gs = vlx.ScfRestrictedDriver(comm, ostream=silent_ostream)
scf_gs.compute(mol_vlx, bas_vlx)

# Extract orbital energies
orbital_energies = scf_gs.scf_tensors["E"]
print("1s E from Koopmans' theorem:", np.around(au2ev * orbital_energies[0], 2), "eV")

1s E from Koopmans' theorem: -559.5 eV


This can be compared to the experimental value, which is 539.9 eV. The overestimation of almost 20 eV is due to lack of relaxation, which results from the reorganization of the valence electrons in response to a reduction in screening of the nuclei. By comparison, in valence spectroscopies the change in shielding is significantly smaller, yielding far smaller errors due to the lack of relaxation.

An approach for including relaxation effects in IE calculations is to use the $$\Delta\textrm{SCF}$$ approach, in which the difference in total energy between the neutral and explicitly optimized core-hole states yields estimates of core-electron binding energies.

### The maximum overlap method#

In a normal SCF cycle without constraints, the core-ionized state is not possible to converge since it collapses into a lower-energy valence-ionized state. To obtain a core-hole state, an approach which avoids the convergence to the valence-hole state is needed. This can be done using several approaches, but the focus here will be on the maximum overlap method (MOM). The idea behind MOM is to run an SCF optimization with the constraint of maximizing the overlap of the orbitals to some initial guess [BS71, GBG08, KSK+67, LS61]. The guess is typically taken as the converged wave function of the neutral (no core-hole) system, and new occupied orbitals at each iteration are chosen such that they overlap most with the span of old occupied orbitals. Defining the overlap matrix as

$\textbf{O} = \left( \textbf{C}^{\textrm{old}} \right)^{\dagger} \textbf{SC}^{\textrm{new}},$

with $$O_{ij}$$ being the overlap between the $$i$$:th old and $$j$$:th new orbital, and $$\textbf{C}$$ and $$\textbf{S}$$ are the MO coefficient matrices and basis function overlap matrix, respectively. The projections of the $$j$$:th new orbital onto the old space is

$p_j = \sum_i^n O_{ij} = \sum_a^N ( \sum_b^N ( \sum_i^n C_{ib}^{old} ) S_{bc} ) C_{cj}^{old}.$

The occupied orbitals are then selected by largest $$p_j$$ values.

With this, a converged core-hole wave-function is obtained with the following scheme:

1. Calculate the wave-function of the original system (e.g. no core-hole).

2. Use this as a starting guess for a new wave-function optimization, removing an electron from the desired orbital (and potentially occupying some other virtual orbital).

3. Optimize this state using the MOM restriction.

This scheme can consider many unstable states, including core-ionized and certain core-excited states. Some care should be observed, though:

• If excited states are considered, check that the final state is reasonably close to what is to be investigated.

• Delocalized core orbitals need to be localized, for example by employing ECPs.

• If post-HF methods are used, the ground state calculation can attempt to (partially) fill the core-hole of the non-Aufbau reference state, leading to a corrupt wave-function which represents an intermediate between a core-hole/valence-hole state.

Note

We currently use pyscf for the calculations involving explicit core-holes, as veloxchem presently lacks MOM.

### IEs from $$\Delta$$SCF#

The IE of water using $$\Delta\textrm{SCF}$$ is then computed as:

# Create pyscf mol object
mol = gto.Mole()
mol.atom = water_xyz
mol.basis = "6-31G"
mol.build()

# Perform unrestricted SCF calculation
scf_gs = scf.UHF(mol)
scf_gs.kernel()

# Copy molecular orbitals and occupations
mo0 = copy.deepcopy(scf_gs.mo_coeff)
occ0 = copy.deepcopy(scf_gs.mo_occ)

# Create 1s core-hole by setting alpha_0 population to zero
occ0[0][0] = 0.0

# Perform unrestricted SCF calculation with MOM constraint
scf_ion = scf.UHF(mol)
scf_ion.kernel()

# IE from energy difference
print(
"Ionization energy:",
np.around(au2ev * (scf_ion.energy_tot() - scf_gs.energy_tot()), 2),
"eV",
)

converged SCF energy = -75.9838703827193  <S^2> = 6.3353767e-12  2S+1 = 1

converged SCF energy = -56.0754789470865  <S^2> = 0.76257805  2S+1 = 2.0125387

Ionization energy: 541.73 eV

Overwritten attributes  get_occ  of <class 'pyscf.scf.uhf.UHF'>


This estimate is within 2 eV from experiment, and we can include electron correlation using MP2:

# Run MP2 on neutral and core-hole state
mp_res = mp.MP2(scf_gs).run()
mp_ion = mp.MP2(scf_ion).run()

# IE from energy difference
print("Ionization energy:", np.around(au2ev * (mp_ion.e_tot - mp_res.e_tot), 2), "eV")

E(UMP2) = -76.113048395549  E_corr = -0.129178012829714

E(UMP2) = -56.1523709631035  E_corr = -0.0768920160169647

Ionization energy: 543.16 eV


Alternatively, we can consider this with KS-DFT:

# Perform unrestricted SCF calculation
scf_gs = scf.UKS(mol)
scf_gs.xc = "b3lyp"
scf_gs.kernel()

# Copy molecular orbitals
mo0 = copy.deepcopy(scf_gs.mo_coeff)
occ0 = copy.deepcopy(scf_gs.mo_occ)

# Create 1s core-hole by setting alpha_0 population to zero
occ0[0][0] = 0.0

# Perform unrestricted SCF calculation with MOM constraint
scf_ion = scf.UKS(mol)
scf_ion.xc = "b3lyp"
scf_ion.kernel()

# IE from energy difference
print(
"Ionization energy:",
np.around(au2ev * (scf_ion.energy_tot() - scf_gs.energy_tot()), 2),
"eV",
)

converged SCF energy = -76.3480474320838  <S^2> = 2.5509372e-11  2S+1 = 1

Overwritten attributes  get_occ  of <class 'pyscf.dft.uks.UKS'>

converged SCF energy = -56.4125421194362  <S^2> = 0.75441021  2S+1 = 2.0044054

Ionization energy: 542.47 eV


MP2 and DFT currently bring us away from experiment when comparing to HF, but we note that the basis set is small and ill-suited for core properties. Using a cc-pCVTZ basis we obtain:

• HF: 538.93 eV

• MP2: 540.26 eV

• B3LYP: 539.19 eV

Thus closer to experiment for the correlated methods (note that scalar-relativistic effects will shift these results upwards by $$\sim$$0.37 eV).

### Spectra from $$\Delta\textrm{SCF}$$#

The X-ray photoemission spectrum can be formed by computing all relavant core-electron binding energies (IEs), giving each unique ionization equal intensity, and broadening over the result with an appropriate convolution function.

We here consider ethyl trifluoroacetate( (CF$$_3$$-CO-O-CH$$_2$$-CH$$_3$$), also known as the ‘ESCA’ molecule, which possess very illustrative shifts in XPS features of the carbon $$K$$-edge due to interactions with the very electronegative oxygen and fluorine atoms. We calculate the four carbon 1s IEs with $$\Delta\textrm{SCF}$$:

esca_xyz = """
C         1.1183780041   -0.3449927563    0.0140172822
C         0.2599020293    0.9430366372   -0.0248385056
F         2.4019839124   -0.0408200726   -0.0368443870
F         0.8399255955   -1.1560630574   -1.0169777359
F         0.8945563482   -1.0267696528    1.1485539196
O         0.7954196415    2.0187110930    0.0011424145
O        -1.0675509573    0.8144007049   -0.0614767265
C        -1.7486063992   -0.4602333713   -0.1739210991
C        -3.1986940727   -0.1850294900    0.1302033614
H        -1.3268717284   -1.1720377773    0.5289059043
H        -1.6185520565   -0.8344576449   -1.1859838788
H        -3.3116641442    0.1959612359    1.1419001202
H        -3.7716619472   -1.1057672384    0.0380096706
H        -3.6012985602    0.5468300237   -0.5651778685
"""

# Create pyscf mol object
mol = gto.Mole()
mol.atom = esca_xyz
mol.basis = "6-31G"
mol.build()

# Perform unrestricted SCF calculation
scf_gs = scf.UHF(mol)
scf_gs.kernel()

# Perform core-hole calculation for each carbon atom (MOs 5-8)
esca_ies = []
for n in [5, 6, 7, 8]:
# Copy molecular orbitals
mo0 = copy.deepcopy(scf_gs.mo_coeff)
occ0 = copy.deepcopy(scf_gs.mo_occ)
# Create 1s core-hole by setting alpha_n population to zero
occ0[0][n] = 0.0
# Perform unrestricted SCF calculation with MOM constraint
scf_ion = scf.UHF(mol)
scf_ion.kernel()
# Append resulting IEs to an array
esca_ies.append(au2ev * (scf_ion.energy_tot() - scf_gs.energy_tot()))

# IEs from above calculation
esca_ies = [302.187, 299.362, 295.405, 293.277]

# Broaden with a Lorenztian and plot
plt.figure(figsize=(6, 3))
x, y = esca_ies, np.ones((len(esca_ies)))
# Lorentzian function takes energies, intensities, energy range, step-size, and a broadening parameter
xi, yi = lorentzian(x, y, min(x) - 2, max(x) + 2, 0.01, 0.5)
plt.plot(xi, yi)
plt.xlabel("Photon energy [eV]")
plt.ylabel("Intensity [arb. units]")
plt.tight_layout()
plt.show()


For the comparison to experimental measurements, a cc-pVTZ basis set is selected, as augmented with additional core-polarizing functions for the carbon atoms (i.e. cc-pCVTZ):

The features originate in, from low to high energy: H$$_3$$C-, -CH$$_2$$-, -CO-, and -CF$$_3$$. We note a good agreement for the first two features, but the comparison deteriorates for higher energies. This is largely an effect of the more significant change in electronic structure imposed by the very electronegative fluorine and oxygen atoms, which the current level of theory is not capable of describing as well.

## X-ray absorption spectroscopy#

In XAS the energy-dependent photoabsorption of the sample is measured, probing excitation of core electrons to bound or continuum states. These two different final states then provide information on unoccupied states and local chemical environment, respectively, with a prototypical spectrum of a conjugated system:

We see two disctinct regions with different features and physical origins:

• NEXAFS (near-edge X-ray absorption fine structure): the region in which transitions to bound states occurs (this region is also known as XANES, X-ray absorption near-edge structure).

• EXAFS (extended X-ray absorption fine structure): which consists of transitions to continuum states, as modulated by constructive or destructive interference between the emitted electron and the environment. This thus provide information on local chemical bonds.

Included is also the position of the IE, as well assignment of typical features ($$\pi^*$$, Rydberg, and $$\sigma^*$$). Transitions to $$\pi^*$$ states are typically low in energy and very intense, and thus easy to identify. By comparison, Rydberg states are typically weak and converge to the ionization energy, and for condensed phase they tend to be atttenuated. Finally, transitions to $$\sigma^*$$ states can be found above the ionization energy.

The total spectrum is referred to as the XAFS (X-ray absorption fine structure) region. For the remainder of this tutorial we will focus on NEXAFS, for which good starting points for modeling with ADC and TDDFT can be found in [BA10, ND18, WHWD15].

### Decoupling from the valence continuum#

A complication that arises when considering core excitations is the embedding of the core-excited states in the continuum of valence-ionized states. This means that the targetted core-excited states cannot easily be resolved with, e.g., a Davidson approach, as a large number of valence-excited and -ionized states need to be converged first. A number of approaches have been developed to remove or circumvent this issue, including:

• The core-valence separation (CVS) approximation

• Restricted energy window (which is largely a different flavour of CVS) [LvKKG12]

• Damped linear response theory, or the complex polarization propagator (CPP), where the response at arbitrary frequencies is directly probed

• Real-time propagation of the electronic structure [LvKKG12]

Alternatively, when using methods constructing explicit excited states this can be avoided by some constraint scheme, such as MOM.

Warning

Within a damped response framework and when using real-time propagation there can be spurious mixing with valene-ionized states, which needs to be avoided.

The focus in this tutorial is the use of the CVS approximation and the CPP approach, with which interior eigenstates and response at arbitrary frequencies can be addressed, as illustrated in Fig. 8:

Here we see the global spectrum of water, using a rather limited basis set. Even so, ADC(1) requires more than 50 states to resolve the core-excitation region (situated at about 535 eV for water), and for ADC(2) we could only converge abot 50 states, which did not bring us to the core-excitations. With the CVS approximation we directly target the core-excitation region, albeit at the cost of introducing an error. However, this error has been shown to be relatively small (provided a reasonable basis set) and stable for ADC [HF20], and for TDHF and TDDFT it is all but negligible. Using CPP-DFT, this error is avoided all together (although CPP-ADC has been noted to be hard to converge for core-excited state [RRH+21], so here CVS-CPP-ADC might instead be used - thus reintroducing this small error).

With CVS-ADC, the core-excitation spectrum is calculate by specifying the index of the probed core orbitals(s):

# Construct structure and basis objects
struct = gator.get_molecule(water_xyz)
basis = gator.get_molecular_basis(struct, "6-31G")

# Perform SCF calculation
scf_gs = gator.run_scf(struct, basis)

# Calculate the 6 lowest eigenstates with CVS restriction to MO #1 (oxygen 1s)
struct, basis, scf_gs, method="cvs-adc2x", singlets=4, core_orbitals=1
)

* Info * Reading basis set from file: /home/emi/miniconda3/envs/echem/lib/python3.9/site-packages/veloxchem/basis/6-31G



Self Consistent Field Driver Setup
====================================

Wave Function Model             : Spin-Restricted Hartree-Fock
Initial Guess Model             : Superposition of Atomic Densities
Convergence Accelerator         : Two Level Direct Inversion of Iterative Subspace
Max. Number of Iterations       : 50
Max. Number of Error Vectors    : 10
Convergence Threshold           : 1.0e-06
ERI Screening Scheme            : Cauchy Schwarz + Density
ERI Screening Mode              : Dynamic
ERI Screening Threshold         : 1.0e-12
Linear Dependence Threshold     : 1.0e-06

* Info * Nuclear repulsion energy: 9.1561447194 a.u.

* Info * Overlap matrix computed in 0.00 sec.

* Info * Kinetic energy matrix computed in 0.00 sec.

* Info * Nuclear potential matrix computed in 0.00 sec.


* Info * Orthogonalization matrix computed in 0.00 sec.


* Info * SAD initial guess computed in 0.00 sec.


* Info * Starting Reduced Basis SCF calculation...
* Info * ...done. SCF energy in reduced basis set: -75.983870205310 a.u. Time: 0.03 sec.


* Info * Overlap matrix computed in 0.00 sec.

* Info * Kinetic energy matrix computed in 0.00 sec.

* Info * Nuclear potential matrix computed in 0.00 sec.


* Info * Orthogonalization matrix computed in 0.00 sec.



Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change
--------------------------------------------------------------------------------------------
1       -75.983870373952    0.0000000000      0.00006826      0.00001638      0.00000000

                  2       -75.983870375702   -0.0000000017      0.00002636      0.00000475      0.00006304

                  3       -75.983870375765   -0.0000000001      0.00000396      0.00000061      0.00000524

                  4       -75.983870375769   -0.0000000000      0.00000029      0.00000008      0.00000296


*** SCF converged in 4 iterations. Time: 0.02 sec.


               Spin-Restricted Hartree-Fock:
-----------------------------
Total Energy                       :      -75.9838703758 a.u.
Electronic Energy                  :      -85.1400150952 a.u.
Nuclear Repulsion Energy           :        9.1561447194 a.u.
------------------------------------

Ground State Information
------------------------
Charge of Molecule            :  0.0
Multiplicity (2S+1)           :  1.0
Magnetic Quantum Number (M_S) :  0.0

Spin Restricted Orbitals
------------------------

Molecular Orbital No.   1:
--------------------------
Occupation: 2.000 Energy:  -20.56128 a.u.
(   1 O   1s  :    -1.00)

Molecular Orbital No.   2:
--------------------------
Occupation: 2.000 Energy:   -1.35434 a.u.
(   1 O   1s  :     0.21) (   1 O   2s  :    -0.47) (   1 O   3s  :    -0.48)

Molecular Orbital No.   3:
--------------------------
Occupation: 2.000 Energy:   -0.70774 a.u.
(   1 O   1p+1:     0.51) (   1 O   2p+1:     0.27) (   2 H   1s  :    -0.26)
(   3 H   1s  :     0.26)

Molecular Orbital No.   4:
--------------------------
Occupation: 2.000 Energy:   -0.56024 a.u.
(   1 O   2s  :    -0.18) (   1 O   3s  :    -0.31) (   1 O   1p0 :    -0.55)
(   1 O   2p0 :    -0.40)

Molecular Orbital No.   5:
--------------------------
Occupation: 2.000 Energy:   -0.50122 a.u.
(   1 O   1p-1:     0.64) (   1 O   2p-1:     0.51)

Molecular Orbital No.   6:
--------------------------
Occupation: 0.000 Energy:    0.20280 a.u.
(   1 O   3s  :     1.17) (   1 O   1p0 :    -0.23) (   1 O   2p0 :    -0.48)
(   2 H   2s  :    -0.99) (   3 H   2s  :    -0.99)

Molecular Orbital No.   7:
--------------------------
Occupation: 0.000 Energy:    0.29874 a.u.
(   1 O   1p+1:    -0.34) (   1 O   2p+1:    -0.82) (   2 H   2s  :    -1.38)
(   3 H   2s  :     1.38)

Molecular Orbital No.   8:
--------------------------
Occupation: 0.000 Energy:    1.05548 a.u.
(   1 O   2p+1:     0.74) (   2 H   1s  :     0.97) (   2 H   2s  :    -0.50)
(   3 H   1s  :    -0.97) (   3 H   2s  :     0.50)

Molecular Orbital No.   9:
--------------------------
Occupation: 0.000 Energy:    1.16442 a.u.
(   1 O   1p-1:    -0.96) (   1 O   2p-1:     1.04)

Molecular Orbital No.  10:
--------------------------
Occupation: 0.000 Energy:    1.18237 a.u.
(   1 O   2s  :     0.26) (   1 O   3s  :    -0.18) (   1 O   1p0 :    -0.67)
(   1 O   2p0 :     0.22) (   2 H   1s  :    -0.84) (   2 H   2s  :     0.55)
(   3 H   1s  :    -0.84) (   3 H   2s  :     0.55)


SCF converged in 4 iterations.
Total Energy: -75.9838703758 au

Starting cvs-adc2x singlet Jacobi-Davidson ...
Niter n_ss  max_residual  time  Ritz values

  1     8        1.2621   95ms  [20.27035896 20.30081829 20.86977777 20.91121528]

  2    16      0.038976  180ms  [19.72495586 19.80691963 20.50881    20.51569984]

  3    24     0.0035541  234ms  [19.71673642 19.79707758 20.49577433 20.5050343 ]

  4    32    0.00029894  269ms  [19.71639059 19.79671156 20.49373717 20.50484016]

  5    40    3.1991e-05  310ms  [19.7163757  19.79669977 20.49353402 20.5048266 ]
=== Restart ===

  6    16    5.9779e-06  267ms  [19.7163752  19.79669907 20.49351316 20.50482466]

  7    24    1.1003e-06  236ms  [19.71637516 19.79669898 20.49350798 20.50482383]

  8    32    5.7208e-07  268ms  [19.71637516 19.79669897 20.49350754 20.50482373]
=== Converged ===
Number of matrix applies:    64
Total solver time:             1s 881ms


The results can be printed as a table, and convoluted and plotted using built-in functionalities or a custom broadening:

# Print information on eigenstates

plt.figure(figsize=(6, 5))
# Convolute using functionalities available in gator and adcc
plt.subplot(211)

# Manually convoluted using excitation energies and oscillator strengths
plt.subplot(212)
xi, yi = lorentzian(x, y, 533, 543, 0.01, 0.5)
plt.plot(xi, yi)
plt.xlabel("Photon energy [eV]")
plt.ylabel("Cross section [au]")
plt.tight_layout()
plt.show()

+--------------------------------------------------------------+
| cvs-adc2x                               singlet ,  converged |
+--------------------------------------------------------------+
|  #        excitation energy     osc str    |v1|^2    |v2|^2  |
|          (au)           (eV)                                 |
|  0      19.71638      536.5099   0.0175       0.8       0.2  |
|  1       19.7967      538.6956   0.0368    0.8087    0.1913  |
|  2      20.49351      557.6567   0.0098    0.7858    0.2142  |
|  3      20.50482      557.9647   0.1007    0.8441    0.1559  |
+--------------------------------------------------------------+


The built-in functionality uses a Lorentzian broadening with a small broadening parameter, and plots the full region which is resolved, while we might want to focus on lower-energy regions and use custom broadening schemes (as done in the lower panel). We note that high-energy features can often be a result of a discretized continuum region, as seen by the unphysically intense feature at $$\sim$$558 eV. A larger basis set will flatten out this region, but care should be taken for any analysis above the ionization energy.

Increasing the basis set to aug-cc-pCVTZ/cc-pVTZ (for oxygen and hydrogen, respectively), we obtain results in good agreement with experiment:

### Spectra from TDDFT#

Note

To be added. For now, see the CVS example in non-standard calculations.

### Spectra from CPP-DFT#

With the complex polarization propagator (CPP), also called the damped linear response function, the absorption and dispersion spectrum can be evaluated at arbitrary energies [NBJO01]. As such, the X-ray absorption spectrum can be directly calculated by calculating the response at X-ray photon energies. For this, we need to decide on which energy region to use, which can be estimated by, e.g.:

• Screening extended energy regions with large step size. This could potentially be done using a smaller basis set.

• With a method lacking relaxation: starting the region a few eV below the IE from Koopmans’ theorem. This will not work for methods capable of including relaxation, where IEs from $$\Delta$$SCF could instead be considered as a first rough estimate (noting that it will work better or worse for different levels of, e.g., ADC theory).

Choosing an energy region of 514-525 eV, using a step length of 0.1 eV and a damping factor of 0.3 eV, we calculate the CPP-DFT (with the B3LYP functional) water spectrum as:

# Prepare molecule and basis objects

# SCF settings and SCF optimization
scf_drv = vlx.ScfRestrictedDriver(comm, ostream=silent_ostream)
scf_settings = {"conv_thresh": 1.0e-6}
method_settings = {"xcfun": "b3lyp"}
scf_drv.update_settings(scf_settings, method_settings)
scf_drv.compute(molecule, basis)

# Define spectrum region to be resolved
freqs = np.arange(515.0, 525.0, 0.1) / au2ev
freqs_str = [str(x) for x in freqs]

# Calculate the response
cpp_prop = LinearAbsorptionCrossSection(
{"frequencies": ",".join(freqs_str), "damping": 0.3 / au2ev}, method_settings
)
cpp_prop.init_driver(comm, ostream=silent_ostream)
cpp_prop.compute(molecule, basis, scf_drv.scf_tensors)

# Extract the imaginary part of the complex response function and convert to absorption cross section
sigma = []
for w in freqs:
axx = -cpp_prop.rsp_property["response_functions"][("x", "x", w)].imag
ayy = -cpp_prop.rsp_property["response_functions"][("y", "y", w)].imag
azz = -cpp_prop.rsp_property["response_functions"][("z", "z", w)].imag
alpha_bar = (axx + ayy + azz) / 3.0
sigma.append(4.0 * np.pi * w * alpha_bar / 137.035999)


Resulting spectra can be plotted as raw, or splined for smoother figures:

# Make figure with panels of 3:1 width
plt.figure(figsize=(9, 5))
gs = gridspec.GridSpec(2, 2, width_ratios=[3, 1])

# Raw results for the full region
plt.subplot(gs[0])
plt.plot(au2ev * freqs, sigma, "bx-")

# Raw results for a zoomed-in region
plt.subplot(gs[1])
plt.plot(au2ev * freqs, sigma, "bx-")
plt.xlim((519.5, 520.5))
plt.ylim((0, 0.10))

# Raw and splined spectra for the full region
plt.subplot(gs[2])
plt.plot(au2ev * freqs, sigma, "bx")
x = np.arange(min(au2ev * freqs), max(au2ev * freqs), 0.01)
y = interp1d(au2ev * freqs, sigma, kind="cubic")
plt.plot(x, y(x), "r")
plt.legend(("Raw", "Splined"))

# Raw and splined spectra for a zoomed-in region
plt.subplot(gs[3])
plt.plot(au2ev * freqs, sigma, "bx")
plt.plot(x, y(x), "r")
plt.xlim((519.5, 520.5))
plt.ylim((0, 0.10))
plt.show()


Using broadened eigenstates, increasing the basis set size and comparing to experiment we get:

As can be seen, the absolute excitation energies are significantly off, with a shift of 14.5 eV required to yield an approximate alignment to experiment. Part of this shift is due to lack of relativistic effects (~0.37 eV), and a small discrepancy can be expected due to the use of a limited basis set, but the main part of the error is due to two effects:

1. Lack of relaxation, which yields a final state too high in energy and thus shifts the theoretical spectra upwards in energy.

2. The self-interaction error (SIE), which yields an erroneous self-repulsion of the dense core orbitals, and results in too low core-excitation energies.

For B3LYP, the combined effect of these two errors results in too low transition energies, meaning that the SIE dominates. By changing the amount of HF exchange, this balance will shift and the error eventually becomes positive.

Despite this large absolute error, the TDDFT spectrum is in reasonable agreement with experiment once it has been shifted. In general, different xc-functionals can exhibit good precision (relative errors), despite poor accuracy (absolute errors), such that a constant shift in energy can be applied.

## X-ray emission spectroscopy#

Following the creation of a core-hole, the system will rapidly ($$\sim$$fs) decay to a valence-hole or valence-excited state, depending on whether the system is core-ionized or core-excited. These two different intermediate states are probed when the X-ray pulse is in resonance (excited regime) or not (ionized regime), yieldig resonant/non-resonant X-ray emission spectroscopy. This section will focus on the latter, with the former following in the next section.

In the decay process the system gains large amounts of energy which need to be discarded, as can occur by emission of high-energy photons (fluorescence) or high-lying electrons (Auger), with relative probability depending on edge and element. For the K-edge the relative probabilities are:

Auger decay dominates for light elements, which has impact on both our ability of obtaining well-resolved spectra for these elements, as well as the impact of X-ray damage on biological matter.

In non-resonant XES we measure the fluorescent decay of core-ionized states, through which we probe the (element-specific) nature of the local chemical environment. For light elements, the valence space is probed directly, but for heavier elements we obtain a number of different decay channels from the valence region and the outer core region. These lines vary considerably in energy, intensity, and sensitivity to local environment [KLP+16]:

1. The $$2p \rightarrow 1s$$ ($$K\alpha_1$$ and $$K\alpha_2$$) lines are most intense (due to the large overlap) and of lowest energy, and relatively insensitive to the valence region

2. The $$3p \rightarrow 1s$$ ($$K\beta '$$ and $$K\beta_{1,3}$$) lines are higher in energy and more sensitive to the local environment than $$K\alpha$$, but about an order of magnitude weaker.

3. Transitions from the ligand MOs ($$K\beta ''$$ and $$K\beta_{2,5}$$), or valence-to-core (VtC), yield the weakest, highest-energy features, which are the most sensitive to the local environment

This tutorial focuses on spectra of light elements, as modelled by ground state DFT, TDDFT, and ADC. Good starting points for learning about such modeling techniques can be found in [Bes21, BA10, FD19].

One approach for considering X-ray emission spectra is to converge a core-hole state and use this core-ionized system as the reference state for a normal excited state calculation [BA10, FD19]. The core-decay transitions are then the first negative eigenstates of the system, such that a standard Davidson calculation can be run on top of this reference state.

Considering water with ADC(2), the oxygen X-ray emission spectrum is calculated as:

# Create pyscf mol object
mol = gto.Mole()
mol.atom = water_xyz
mol.basis = "6-31G"
mol.build()

# Perform unrestricted SCF calculation
scf_res = scf.UHF(mol)
scf_res.kernel()

# Copy molecular orbitals
mo0 = copy.deepcopy(scf_res.mo_coeff)
occ0 = copy.deepcopy(scf_res.mo_occ)

# Create 1s core-hole by setting alpha_0 population to zero
occ0[0][0] = 0.0

# Perform unrestricted SCF calculation with MOM constraint
scf_ion = scf.UHF(mol)
scf_ion.kernel()


converged SCF energy = -75.9838703827193  <S^2> = 6.3327121e-12  2S+1 = 1

converged SCF energy = -56.0754789470865  <S^2> = 0.76257805  2S+1 = 2.0125387

Starting adc2  Jacobi-Davidson ...
Niter n_ss  max_residual  time  Ritz values

  1     8       0.77768  851ms  [-19.29082067 -19.19940205 -18.96181907 -18.3637012 ]
2    16     0.0010764   89ms  [-19.44427093 -19.36689642 -19.16826546 -18.57947501]
3    24    4.4301e-07  126ms  [-19.44429764 -19.3669284  -19.16831914 -18.57954987]

=== Converged ===
Number of matrix applies:    24
Total solver time:             1s  72ms


Tip

There may be convergence issues if you include transitions of very different types, e.g. when resolving transitions into the core and out to the unoccupied valence region simultaneously. As such, the number of required states should not be too high.

The results can be printed as a table, and convoluted and plotted using built-in functionalities or a custom broadening. Note that the eigenvalues are negative, and a sign change is thus appropriate when constructing the spectra.

# Print information on eigenstates

plt.figure(figsize=(6, 5))
# Convolute using functionalities available in gator and adcc
plt.subplot(211)

# Convoluted using information on excitation energies and oscillator strengths
plt.subplot(212)
xi, yi = lorentzian(x, y, 519, 531, 0.01, 0.5)
plt.plot(xi, yi)
plt.xlabel("Photon energy [eV]")
plt.ylabel("Cross section [au]")
plt.tight_layout()
plt.show()

+--------------------------------------------------------------+
| adc2                                        any ,  converged |
+--------------------------------------------------------------+
|  #        excitation energy     osc str    |v1|^2    |v2|^2  |
|          (au)           (eV)                                 |
|  0      -19.4443     -529.1063   0.0568    0.9548    0.0452  |
|  1     -19.36693      -527.001   0.0458    0.9503   0.04966  |
|  2     -19.16832     -521.5965   0.0419    0.9376   0.06245  |
|  3     -18.57955     -505.5753   0.0031    0.9335   0.06655  |
+--------------------------------------------------------------+


Using a cc-pCVTZ/cc-pVTZ basis set and comparing to experiment:

While the broadening between the different lines is noticably different, it is clear that the energy spacing and integrated intensities of the features are in good agreement with experiment. The theoretical spectrum is obtained at the ADC(2) level of theory, which has been seen to yield better absolute energies than ADC(2)-x and ADC(3/2) (for XES) [FD19].

Note

### Spectra from ground state MOs#

While we have previously seen that Koopmans’ theorem is in significant disagreement with experimental ionization energies (due to lack of relaxation), it has been noted that the relative position of the (valence) MO energies can still be quite reasonable [Bes21]. As such, a computationally simple approach for modeling X-ray emission spectra has been developed, in which transition energies are estimated from:

$\Delta E = \epsilon_v - \epsilon_c$

i.e. from the difference in ground state MO energies. The intensities are taken to be proportional to the transition dipole moment:

$f \propto | \langle \phi_i | \hat{\mu} | \phi_c \rangle | ^2$

This independent-particle picture method is very attractive in its simplicity, requiring only a single ground state calculation and the construction of a number of transition dipole moments in order to construct the full spectrum. However, the absolute energies are significantly off, and any differences in relaxation between the different channels will not be accounted for, and multichannel and many-body effects are neglected. Nevertheless, this approach has been noted to yield good (relative) agreement to experimental measurements, and it is capable of considering even massive molecular systems.

The calculation of an emission spectra thus requires the calculation of the ground state and of all relevant transition dipole moments:

# Prepare molecule and basis objects

# SCF settings and calculation
scf_drv = vlx.ScfRestrictedDriver(comm, ostream=silent_ostream)
scf_settings = {"conv_thresh": 1.0e-6}
method_settings = {"xcfun": "b3lyp"}
scf_drv.update_settings(scf_settings, method_settings)
scf_drv.compute(molecule, basis)

# Extract orbital energies and number of occupied MOs
orbital_energies = scf_drv.scf_tensors["E"]
nocc = molecule.number_of_alpha_electrons()
print("Occupied MO energies:", au2ev * orbital_energies[:nocc])

# Define probed core MO and calculate energy differences
n_core = 0
energy_diff = []
for n_mo in np.arange(nocc):
if n_mo != n_core:
energy_diff.append(au2ev * (orbital_energies[n_mo] - orbital_energies[n_core]))

# Extract MO coefficients
mo_C = scf_drv.scf_tensors["C"]

# Load drivers for transition dipole moment
dipole_drv = vlx.ElectricDipoleIntegralsDriver(comm)

# Calculated transition dipole moments and convert to MO basis
dipole_matrices = dipole_drv.compute(molecule, basis)
x_ao = dipole_matrices.x_to_numpy()
x_mo = np.matmul(mo_C.T, np.matmul(x_ao, mo_C))
y_ao = dipole_matrices.y_to_numpy()
y_mo = np.matmul(mo_C.T, np.matmul(y_ao, mo_C))
z_ao = dipole_matrices.z_to_numpy()
z_mo = np.matmul(mo_C.T, np.matmul(z_ao, mo_C))

# Construct oscillator strengths
xx, yy, zz, sigma = [], [], [], []
for n_mo in np.arange(nocc):
if n_mo != n_core:
x_tmp, y_tmp, z_tmp = x_mo[n_core, n_mo], y_mo[n_core, n_mo], z_mo[n_core, n_mo]
xx.append(x_tmp)
yy.append(y_tmp)
zz.append(z_tmp)
sigma.append(x_tmp ** 2 + y_tmp ** 2 + z_tmp ** 2)

Occupied MO energies: [-520.65709503  -27.6878928   -14.39952882   -9.70377531   -7.94952958]


The results can be printed as a table, and convoluted and plotted using a custom broadening. We here focus on the high-energy features:

# Print results as a table
print("Energy   Osc. str.   x-component   y-component   z-component")
for i in np.arange(len(energy_diff)):
e, os, x, y, z = energy_diff[i], sigma[i], xx[i], yy[i], zz[i]
print(
"{:.3f}  {:8.5f}     {:8.5f}      {:8.5f}      {:8.5f}".format(e, os, x, y, z)
)

# Plot X-ray emission spectrum
plt.figure(figsize=(6, 3))
x, y = energy_diff, sigma
xi, yi = lorentzian(x, y, min(x) - 5, max(x) + 5, 0.01, 0.4)
plt.plot(xi, yi)
plt.xlabel("Photon energy [eV]")
plt.ylabel("Intensity [arb. units]")
plt.xlim((504, 514))
plt.show()

Energy   Osc. str.   x-component   y-component   z-component
492.969   0.00022      0.00000      -0.00000      -0.01478
506.258   0.00238     -0.04878       0.00000      -0.00000
510.953   0.00255      0.00000      -0.00000       0.05050
512.708   0.00342     -0.00000      -0.05848      -0.00000


Using a cc-pCVTZ/cc-pVTZ basis set and comparing to experiment [WBM+12]:

The resulting spectrum has poor absolute energies, but the relative agreement to experiment is quite impressive.

## Resonant inelastic X-ray scattering#

We will here focus on calculating the electronic part of the RIXS plane, and skip cases where one is more interested in probing vibrational modes.

An approximate treatment of the RIXS plane can be achieved by considering it as a two-step process [FPB18], with:

1. A core-excitation as modeled like a normal X-ray absorption spectrum

2. A core-decay process from the core-hole system with an electron put in each relevant final state probed in (1), modeled as a (non-resonant) X-ray emission spectrum

The RIXS plane is then constructed from the above absorption and emission energies, with the RIXS intensity from the multiplication of the two respective oscillator strengths.

Note that this approach is an approximation, as it disregards potentially important detuning effects. In the example given below we also ignore the coupling between the polarization of the incoming and outgoing photons. Furthermore, the calculation requires the convergence of core-excited states, which is much trickier than core-ionized, and will likely only work for the first few core-excited states. Finally, depending on how the intermediate core-excited state is constructed, issues with spin-contamination can become highly influential.

Nevertheless, the scheme is quite easy to understand, and introducing core-excited dynamics is relatively straightforward. As such, we here illustrate how the RIXS plane of water can be constructed. First, calculate the X-ray absorption spectrum, as well as the emission spectra with a core-excited reference state:

# Create pyscf mol object
mol = gto.Mole()
mol.atom = water_xyz
mol.basis = "6-31G"
mol.build()

# Perform restricted SCF calculation
scf_gs = scf.RHF(mol)
scf_gs.kernel()

# Calculate X-ray absorption spectrum

# Perform unrestricted SCF calculation
scf_gs = scf.UHF(mol)
scf_gs.kernel()

# Calculate resonant X-ray emission spectrum
#   for states where the core-electron is moved to LUMO or LUMO+1
for n in [5, 6]:
# Copy molecular orbitals
mo0 = copy.deepcopy(scf_gs.mo_coeff)
occ0 = copy.deepcopy(scf_gs.mo_occ)
# Move 1s electron to LUMO / LUMO+1
occ0[0][0] = 0.0
occ0[0][n] = 1.0
# Perform unrestricted SCF calculation with MOM constraint
scf_ion = scf.UHF(mol)
scf_ion.kernel()
# Calculate resonant X-ray emission spectrum


Now the correct combinations of transition energies and intensities need to be calculated:

adc_rixs = []

# Probing LUMO and LUMO+1
for indx in [0, 1]:
# From XAS calculation
# From XES calculation


And the RIXS plane can then be plotted:

plt.figure(figsize=(8, 6))
# Plot absorption spectrum
plt.subplot(222)
xi, yi = lorentzian(x, y, min(x) - 2, min(x) + 5, 0.01, 0.3)
plt.plot(yi, xi)

# Plot emission specta, with vertical shift from absorbed energy
plt.subplot(223)
# LUMO
xi, yi = lorentzian(x, 5.0 * y, x[1] - 9, x[1] + 2, 0.01, 0.3)
plt.plot(xi, yi + au2ev * adc_xas.excitation_energy[0])
# LUMO+1
xi, yi = lorentzian(x, 5.0 * y, x[1] - 9, x[1] + 2, 0.01, 0.3)
plt.plot(xi, yi + au2ev * adc_xas.excitation_energy[1])

# Plot RIXS plane, here with marker size scaling with intensity
plt.subplot(221)