Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

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., Norman & Dreuw (2018).

import copy
import adcc
import gator
import matplotlib.pyplot as plt
import multipsi as mtp
import numpy as np
import veloxchem as vlx
from matplotlib import gridspec
from scipy.interpolate import interp1d
import py3Dmol as p3d

# au to eV conversion factor
au2ev = 27.211386
Source
# These will be removed and included through VeloxChem

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

    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

class MockTask:
    def __init__(self, mol, basis, comm, ostream):
        self.molecule = mol
        self.ao_basis = basis
        self.mpi_comm = comm
        self.ostream = ostream

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 Vidal et al. (2019).

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 1s1s IE of:

water_xyz = """3

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
molecule = vlx.Molecule.read_xyz_string(water_xyz)
basis = vlx.MolecularBasis.read(molecule, "6-31G")
Loading...
# Perform SCF calculation
scf_gs = vlx.ScfRestrictedDriver()
scf_results = scf_gs.compute(molecule, basis)

# Extract orbital energies
orbital_energies = scf_results["E_alpha"]
print(f'1s IE from Koopmans\' theorem: {au2ev * orbital_energies[0]: 5.2f} eV')
1s IE from Koopmans' theorem: -559.50 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 ΔSCF\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 Lykos & Schmeising (1961)King et al. (1967)Bagus & Schaefer (1971)Gilbert et al. (2008). 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

O=(Cold)†SCnew\textbf{O} = \left( \textbf{C}^{\textrm{old}} \right)^{\dagger} \textbf{SC}^{\textrm{new}}
(1)

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

pj=∑inOij=∑aN(∑bN(∑inCibold)Sbc)Ccjoldp_j = \sum_i^n O_{ij} = \sum_a^N ( \sum_b^N ( \sum_i^n C_{ib}^{\textrm{old}} ) S_{bc} ) C_{cj}^{\textrm{old}}
(2)

The occupied orbitals are then selected by largest pjp_j values. The original formulation of MOM used the previous iteration as the old space, but more recent development (and what is used in VeloxChem) includes the use of the initial guess as the old space throughout - this is also referred to as initial maximum overlap method (IMOM).

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.

IEs from Δ\DeltaSCF

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

# Copy molecule object
mol_ion = copy.deepcopy(molecule)

# Change charge and multiplicity
mol_ion.set_charge(1)
mol_ion.set_multiplicity(2)

# Lists of occupied alpha and beta MOs
occa = [0, 1, 2, 3, 4]
occb = [1, 2, 3, 4]

# Calculate SCF of core-ion
scf_ion = vlx.ScfUnrestrictedDriver()
scf_ion.maximum_overlap(mol_ion, basis, scf_gs.mol_orbs, occa, occb)
scf_ion_results = scf_ion.compute(mol_ion, basis)
# IE from energy difference
hf_ie = au2ev * (scf_ion.get_scf_energy() - scf_gs.get_scf_energy())
print(f'HF ionization energy: {hf_ie : 5.2f} eV')
HF ionization energy:  541.73 eV

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

mp2_gs = vlx.Mp2Driver()
mp2_gs.compute_conventional(molecule, basis, scf_gs.mol_orbs)

mp2_ion = vlx.Mp2Driver()
mp2_ion.compute_conventional(molecule, basis, scf_ion.mol_orbs)

mp2_ie = au2ev * ((scf_ion.get_scf_energy() + mp2_ion.e_mp2) - (scf_gs.get_scf_energy() + mp2_gs.e_mp2))
print(f'MP2 ionization energy: {mp2_ie : 5.2f} eV')
MP2 ionization energy:  542.64 eV

Alternatively, we can consider this with KS-DFT:

# Perform SCF calculation
dft_gs = vlx.ScfRestrictedDriver()
method_settings = {"xcfun": "b3lyp"}
dft_gs.update_settings({}, method_settings)
dft_results = dft_gs.compute(molecule, basis)

# Calculate SCF of core-ion
dft_ion = vlx.ScfUnrestrictedDriver()
dft_ion.maximum_overlap(mol_ion, basis, scf_gs.mol_orbs, occa, occb)
dft_ion.update_settings({}, method_settings)
dft_ion_results = dft_ion.compute(mol_ion, basis)
# IE from energy difference
dft_ie = au2ev * (dft_ion.get_scf_energy() - dft_gs.get_scf_energy())
print(f'DFT ionization energy: {dft_ie : 5.2f} eV')
DFT ionization energy:  542.58 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:

methodIE [eV]
HF538.93
MP2540.26
B3LYP539.19

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

Spectra from ΔSCF\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( (CF3_3-CO-O-CH2_2-CH3_3), also known as the ‘ESCA’ molecule, which possess very illustrative shifts in XPS features of the carbon KK-edge due to interactions with the very electronegative oxygen and fluorine atoms. We calculate the four carbon 1s IEs with ΔSCF\Delta\textrm{SCF}:

esca_xyz = """14

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 veloxchem mol and basis objects
molecule = vlx.Molecule.read_xyz_string(esca_xyz)
basis = vlx.MolecularBasis.read(molecule, "6-31G")
Loading...
# Perform SCF calculation
scf_gs = vlx.ScfRestrictedDriver()
scf_gs.compute(molecule, basis)
nocc = molecule.number_of_alpha_electrons()

# Copy molecule object
mol_ion = copy.deepcopy(molecule)

# Change charge and multiplicity
mol_ion.set_charge(1)
mol_ion.set_multiplicity(2)

esca_ies = []
for n in [5, 6, 7, 8]:
    occa = list(np.arange(nocc))
    occb = copy.copy(occa)
    occb.remove(n)
    scf_ion = vlx.ScfUnrestrictedDriver()
    scf_ion.maximum_overlap(mol_ion, basis, scf_gs.mol_orbs, occa, occb)
    scf_ion_results = scf_ion.compute(mol_ion, basis)
    hf_ie = au2ev * (scf_ion.get_scf_energy() - scf_gs.get_scf_energy())
    esca_ies.append(hf_ie)

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

# Broaden with a Lorentzian 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()
<Figure size 600x300 with 1 Axes>

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). Theoretical spectrum shifted by 0.11 eV (to account for relativistic effects) and broadened by 0.3 eV.

The features originate in, from low to high energy: H3_3C-, -CH2_2-, -CO-, and -CF3_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 Wenzel et al. (2015)Besley & Asmuruf (2010)Norman & Dreuw (2018).

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:

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

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 below. Here we show the global spectrum of water, using the 6-311G basis set. Results obtained with ADC(1), CVS-ADC(1), ADC(2), CVS-ADC(2), and CPP-DFT. CPP-DFT results obtained using the B3LYP exchange-correlation functional and a large γ\gamma of 1 eV. ADC results are shown as bar spectra.

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 Herbst & Fransson (2020), 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 Rehn et al. (2021), so here CVS-CPP-ADC might instead be used - thus reintroducing this small error).

Spectra from CVS-ADC

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

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

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

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

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

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

# Print information on eigenstates
print(adc_result.describe())

plt.figure(figsize=(6, 5))
# Convolute using functionalities available in gator and adcc
plt.subplot(211)
plt.title("Built-in broadening")
adc_result.plot_spectrum()

# Manually convoluted using excitation energies and oscillator strengths
plt.subplot(212)
plt.title("Custom broadening (0.5 eV)")
x, y = au2ev * adc_result.excitation_energy, adc_result.oscillator_strength
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.0178       0.8       0.2  |
|  1       19.7967      538.6956   0.0373    0.8087    0.1913  |
|  2      20.49351      557.6567   0.0099    0.7858    0.2142  |
|  3      20.50482      557.9647   0.1016    0.8441    0.1559  |
+--------------------------------------------------------------+
<Figure size 600x500 with 2 Axes>

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 ∼\sim558 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. Here the theoretical spectra has been shifted by 0.37 eV (to account for relativistic effects) and broadened by 0.3 eV.

Spectra from TDDFT

RPA and TDA calculations:

# Prepare molecule and basis objects
molecule = vlx.Molecule.read_xyz_string(water_xyz)
basis = vlx.MolecularBasis.read(molecule, "6-31G")

scf_drv = vlx.ScfRestrictedDriver()
scf_drv.xcfun = "b3lyp"
scf_drv.grid_level = 4
scf_result = scf_drv.compute(molecule, basis)

rpa_drv = vlx.LinearResponseEigenSolver()
rpa_drv.core_excitation = True
rpa_drv.num_core_orbitals = 1
rpa_drv.update_settings({"nstates": 6}, method_settings)
rpa_result = rpa_drv.compute(molecule, basis, scf_result)

tda_drv = vlx.TdaEigenSolver()
tda_drv.core_excitation = True
tda_drv.num_core_orbitals = 1
tda_drv.update_settings({"nstates": 6}, method_settings)
tda_result = tda_drv.compute(molecule, basis, scf_result)

Plotting the results:

plt.figure(figsize=(6, 3))
# RPA results
x = 27.2114*rpa_result['eigenvalues']
y = rpa_result['oscillator_strengths']
xi,yi = lorentzian(x,y,min(x)-1,min(x)+5,0.01,0.3)
plt.plot(xi,yi)
# TDA results
x = 27.2114*tda_result['eigenvalues']
y = tda_result['oscillator_strengths']
xi,yi = lorentzian(x,y,min(x)-1,min(x)+5,0.01,0.3)
plt.plot(xi,yi)
plt.legend(('RPA','TDA'))
plt.tight_layout()
plt.show()
<Figure size 600x300 with 1 Axes>

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 Norman et al. (2001). 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 Δ\DeltaSCF 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
molecule = vlx.Molecule.read_xyz_string(water_xyz)
basis = vlx.MolecularBasis.read(molecule, "6-31G")

# SCF settings and SCF optimization
scf_drv = vlx.ScfRestrictedDriver()
scf_settings = {"conv_thresh": 1.0e-6}
method_settings = {"xcfun": "b3lyp"}
scf_drv.update_settings(scf_settings, method_settings)
scf_result = 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_drv = vlx.rsplinabscross.LinearAbsorptionCrossSection(
    {"frequencies": ",".join(freqs_str), "damping": 0.3 / au2ev}, method_settings
)
cpp_drv.init_driver()
cpp_result = cpp_drv.compute(molecule, basis, scf_result)

# Extract the imaginary part of the complex response function and convert to absorption cross section
sigma = []
for w in freqs:
    axx = -cpp_drv.rsp_property["response_functions"][("x", "x", w)].imag
    ayy = -cpp_drv.rsp_property["response_functions"][("y", "y", w)].imag
    azz = -cpp_drv.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()
<Figure size 900x500 with 4 Axes>

Using broadened eigenstates, increasing the basis set size and comparing to experiment we get the following, where theoretical spectrum has been shifted by 14.5 eV and broadened by 0.3 eV:

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 (∼\simfs) 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), yielding 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 follow the following trends, with figure taken from wikipedia.org:

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 Kowalska et al. (2016):

  1. The 2p→1s2p \rightarrow 1s (Kα1K\alpha_1 and Kα2K\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→1s3p \rightarrow 1s (Kβ′K\beta ' and Kβ1,3K\beta_{1,3}) lines are higher in energy and more sensitive to the local environment than KαK\alpha, but about an order of magnitude weaker.

  3. Transitions from the ligand MOs (Kβ′′K\beta '' and Kβ2,5K\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 Besley & Asmuruf (2010)Fransson & Dreuw (2019)Besley (2021).

Spectra from ADC

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 Besley & Asmuruf (2010)Fransson & Dreuw (2019). 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 veloxchem mol and basis objects
molecule = vlx.Molecule.read_xyz_string(water_xyz)
basis = vlx.MolecularBasis.read(molecule, "6-311G")

# Perform SCF calculation
scf_gs = vlx.ScfRestrictedDriver()
scf_result = scf_gs.compute(molecule, basis)

# Copy molecule object
mol_ion = copy.deepcopy(molecule)

# Change charge and multiplicity
mol_ion.set_charge(2)
mol_ion.set_multiplicity(1)

# Lists of occupied alpha and beta MOs
occa = [1, 2, 3, 4]
occb = [1, 2, 3, 4]

# Calculate SCF of core-ion
scf_ion = vlx.ScfRestrictedDriver()
scf_ion.maximum_overlap(mol_ion, basis, scf_gs.mol_orbs, occa, occb)
scf_ion_result = scf_ion.compute(mol_ion, basis)

# IE from energy difference
hf_ie = au2ev * (scf_ion.get_scf_energy() - scf_gs.get_scf_energy())

# Set up mock task for adcc
task = MockTask(mol_ion, basis, scf_ion.comm, scf_ion.ostream)
scf_ion.task = task

# Four lowest (negative) eigenstates
adc_xes = adcc.adc2(scf_ion, n_states=4)

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
print(adc_xes.describe())

plt.figure(figsize=(6, 5))
# Convolute using functionalities available in gator and adcc
plt.subplot(211)
plt.title("Built-in broadening [note negative energies]")
adc_xes.plot_spectrum()

# Convoluted using information on excitation energies and oscillator strengths
plt.subplot(212)
plt.title("Custom broadening")
x, y = -au2ev * adc_xes.excitation_energy, adc_xes.oscillator_strength
xi, yi = lorentzian(x, y, 598, 603, 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     -22.20632     -604.2647   0.0000    0.9644   0.03563  |
|  1     -22.15933      -602.986   0.0000    0.9622   0.03785  |
|  2     -22.09183     -601.1495   0.1612    0.9588   0.04123  |
|  3     -22.04624     -599.9086   0.1409    0.9561   0.04391  |
+--------------------------------------------------------------+
<Figure size 600x500 with 2 Axes>

Using a cc-pCVTZ/cc-pVTZ basis set and comparing to experiment, where theoretical spectrum shifted by 0.37 eV (to account for relativistic effects) and broadened by 0.3 eV:

While the broadening between the different lines is noticeably 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) Fransson & Dreuw (2019).

Spectra from TDDFT

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 Besley (2021). As such, a computationally simple approach for modeling X-ray emission spectra has been developed, in which transition energies are estimated from:

ΔE=ϵv−ϵc\Delta E = \epsilon_v - \epsilon_c
(3)

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

f∝∣⟨ϕi∣μ^∣ϕc⟩∣2f \propto | \langle \phi_i | \hat{\mu} | \phi_c \rangle | ^2
(4)

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
molecule = vlx.Molecule.read_xyz_string(water_xyz)
basis = vlx.MolecularBasis.read(molecule, "6-31G")

# SCF settings and calculation
scf_drv = vlx.ScfRestrictedDriver()
scf_settings = {"conv_thresh": 1.0e-6}
method_settings = {"xcfun": "b3lyp"}
scf_drv.update_settings(scf_settings, method_settings)
scf_result = scf_drv.compute(molecule, basis)
# Extract orbital energies and number of occupied MOs
orbital_energies = scf_result["E_alpha"]
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_result["C_alpha"]

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

# 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.65707757  -27.6878829   -14.39950787   -9.70379307   -7.94953867]

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
<Figure size 600x300 with 1 Axes>

Using a cc-pCVTZ/cc-pVTZ basis set and comparing to experiment Weinhardt et al. (2012), where theoretical spectrum shifted by 15.5 eV and broadened by 0.3 eV:

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.

Two-step approach with ADC

An approximate treatment of the RIXS plane can be achieved by considering it as a two-step process Fouda et al. (2018), 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
adc_xas = adcc.cvs_adc2(scf_gs, n_singlets=4, core_orbitals=1)
print(adc_xas.describe())

# 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
adc_rxes = []
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.addons.mom_occ(scf_ion, mo0, occ0)
    scf_ion.kernel()
    # Calculate resonant X-ray emission spectrum
    adc_rxes.append(adcc.adc2(scf_ion, n_states=4))

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
    E_xas = au2ev * adc_xas.excitation_energy[indx]
    I_xas = adc_xas.oscillator_strength[indx]
    # From XES calculation
    for i in np.arange(len(adc_rxes[indx])):
        E_xes = -au2ev * adc_rxes[indx].excitation_energy[i]
        I_xes = adc_rxes[indx].oscillator_strength[i]
        adc_rixs.append([E_xas, E_xes, I_xas * I_xes])

And the RIXS plane can then be plotted:

plt.figure(figsize=(8, 6))
# Plot absorption spectrum
plt.subplot(222)
x, y = au2ev * adc_xas.excitation_energy, adc_xas.oscillator_strength
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
x, y = -au2ev * adc_rxes[0].excitation_energy, adc_rxes[0].oscillator_strength
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
x, y = -au2ev * adc_rxes[0].excitation_energy, adc_rxes[1].oscillator_strength
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)
for i in np.arange(len(adc_rixs)):
    tmp = adc_rixs[i]
    plt.plot(tmp[1], tmp[0], "ok", markersize=10000 * tmp[2])

# Zooming in to region spanned by XAS and XES projections
plt.xlim((519, 530))
plt.ylim((536, 544))

plt.tight_layout()
plt.show()
References
  1. Norman, P., & Dreuw, A. (2018). Simulating X-ray Spectroscopies and Calculating Core-Excited States of Molecules. Chem. Rev., 118, 7208–7248. 10.1021/acs.chemrev.8b00156
  2. Vidal, M. L., Krylov, A. I., & Coriani, S. (2019). Dyson orbitals within the fc-CVS-EOM-CCSD framework: theory and application to X-ray photoelectron spectroscopy of ground and excited states. Phys. Chem. Chem. Phys., 22, 2693. 10.1039/C9CP03695D
  3. Lykos, P. G., & Schmeising, H. N. (1961). Maximum Overlap Atomic and Molecular Orbitals. J. Phys. Chem., 35, 288. 10.1063/1.1731901
  4. King, H. F., Stanton, R. E., Kim, H., Wyatt, R. E., & Parr, R. G. (1967). Corresponding Orbitals and the Nonorthogonality Problem in Molecular Quantum Mechanics. J. Chem. Phys., 47, 1936. 10.1063/1.1712221
  5. Bagus, P. S., & Schaefer, H. F. (1971). Direct Near-Hartree–Fock Calculations of the 1s Hole States of NO+. J. Chem. Phys., 55, 1474. 10.1063/1.1676248
  6. Gilbert, A. T. B., Besley, N. A., & Gill, P. M. W. (2008). Self-Consistent Field Calculations of Excited States Using the Maximum Overlap Method (MOM). J. Phys. Chem. A, 112, 13164–13171. 10.1021/jp801738f
  7. Travnikova, O., Børve, K. J., Patanen, M., Söderström, J., Miron, C., Sæthre, L. J., Mårtensson, N., & Svensson, S. (2012). The ESCA molecule—Historical remarks and new results. Journal of Electron Spectroscopy and Related Phenomena, 185(8–9), 191–197. 10.1016/j.elspec.2012.05.009
  8. Travnikova, O., Børve, K. J., Patanen, M., Söderström, J., Miron, C., Sæthre, L. J., Mårtensson, N., & Svensson, S. (2012). The ESCA molecule—Historical remarks and new results. Journal of Electron Spectroscopy and Related Phenomena, 185(8–9), 191–197. 10.1016/j.elspec.2012.05.009
  9. Wenzel, J., Holzer, A., Wormit, M., & Dreuw, A. (2015). Analysis and Comparison of CVS-ADC Approaches up to Third Order for the Calculation of Core-Excited States. J. Chem. Phys., 142, 214104. 10.1063/1.4921841
  10. Besley, N. A., & Asmuruf, F. A. (2010). Time-Dependent Density Functional Theory Calculations of the Spectroscopy of Core Electrons. Phys. Chem. Chem. Phys., 12, 12024–12039. 10.1039/C002207A
  11. Cederbaum, L. S., Domcke, W., & Schirmer, J. (1980). Many-body theory of core holes. Phys. Rev. A, 22, 206. 10.1103/PhysRevA.22.206
  12. Vidal, M. L., Feng, X., Epifanovsky, E., Krylov, A. I., & Coriani, S. (2019). New and efficient equation-of-motion coupled-cluster framework for core-excited and core-ionized states. J. Chem. Theory Comput., 15, 3117–3133. 10.1021/acs.jctc.9b00039
  13. Lopata, K., van Kuiken, B. E., Khalil, M., & Govind, N. (2012). Linear-Response and Real-Time Time-Dependent Density Functional Theory Studies of Core-Level Near-Edge X-Ray Absorption. J. Chem. Theory Comput., 8, 3284–3292. 10.1021/ct3005613
  14. Ekström, U., Norman, P., Carravetta, V., & \AAgren, H. (2006). Polarization Propagator for X-Ray Spectra. Phys. Rev. Lett., 97(143001). 10.1103/PhysRevLett.97.143001
  15. Herbst, M. F., & Fransson, T. (2020). Quantifying the error of the core-valence separation approximation. J. Chem. Phys., 153, 054114. 10.1063/5.0013538
Tutorials and workflows
X-ray spectroscopies
Tutorials and workflows
Spectrum analysis