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.211386Source
# 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 = ostreamX-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 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")# 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 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
with being the overlap between the :th old and :th new orbital, and and are the MO coefficient matrices and basis function overlap matrix, respectively. The projections of the :th new orbital onto the old space is
The occupied orbitals are then selected by largest 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:
Calculate the wave-function of the original system (e.g. no core-hole).
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).
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 SCF¶
The IE of water using 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:
| method | IE [eV] |
|---|---|
| HF | 538.93 |
| MP2 | 540.26 |
| B3LYP | 539.19 |
Thus closer to experiment for the correlated methods (note that scalar-relativistic effects will shift these results upwards by 0.37 eV).
Spectra from ¶
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-CO-O-CH-CH), also known as the ‘ESCA’ molecule, which possess very illustrative shifts in XPS features of the carbon -edge due to interactions with the very electronegative oxygen and fluorine atoms. We calculate the four carbon 1s IEs with :
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")# 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()
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: HC-, -CH-, -CO-, and -CF. 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 (, Rydberg, and ). Transitions to 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 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:
The core-valence separation (CVS) approximation Cederbaum et al. (1980)Wenzel et al. (2015)Vidal et al. (2019)
Restricted energy window (which is largely a different flavour of CVS) Lopata et al. (2012)
Damped linear response theory, or the complex polarization propagator (CPP), where the response at arbitrary frequencies is directly probed Ekström et al. (2006)
Real-time propagation of the electronic structure Lopata et al. (2012)
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 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 |
+--------------------------------------------------------------+
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 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. 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()
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 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
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()
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:
Lack of relaxation, which yields a final state too high in energy and thus shifts the theoretical spectra upwards in energy.
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 (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), 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):
The ( and ) lines are most intense (due to the large overlap) and of lowest energy, and relatively insensitive to the valence region
The ( and ) lines are higher in energy and more sensitive to the local environment than , but about an order of magnitude weaker.
Transitions from the ligand MOs ( and ), 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 |
+--------------------------------------------------------------+
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:
i.e. from the difference in ground state MO energies. The intensities are taken to be proportional to the transition dipole moment:
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
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:
A core-excitation as modeled like a normal X-ray absorption spectrum
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()
Spectra from TDDFT¶
As the name suggests, resonant inelastic X-ray scattering (RIXS) Gel'mukhanov & Ågren (1999) involves the inelastic scattering of an X-ray photon by a molecule, meaning that there is an exchange of energy, e.g., leaving the molecule in an excited state with the emitted photon redshifted. RIXS is also called resonant X-ray Raman scattering, due to the close analogy with Raman spectroscopy.
In a somewhat simplified view, the RIXS process can be thought of as the initial absorption of an X-ray, taking the system from the electronic ground state into a short-lived, core-excited state. The core-excited state quickly decays onto the manifold of valence-excited states, and naturally also back to the ground state, at which point we are at the elastic line since no energy is exchanged in this case. It is important to note that RIXS is a coherent process, and the previous description might suggest something more like two independent events, but in a proper description the absorption and subsequent emission are coupled and cannot be described independently, formally speaking (though this can serve as a crude approximation to computing RIXS spectra), hence the caveat that this is a simplified picture. In principle, one could have a blue-shift of the emitted photon as well if the initial state is not the ground state, under normal conditions however, molecules are in the electronic ground state, as we will assume, though one could imagine a pump--probe type setup as well.
Our starting point for the analysis is the Kramers--Heisenberg--Dirac (KHD) expression for the RIXS scattering amplitude, which at resonance reads
Here is the Cartesian component of the length gauge electric dipole operator denote the electronic ground, core/ intermediate, and valence/final state, respectively, is the vertical excitation energy from the ground state to , and are the incident and outgoing photon frequencies, respectively. is the natural linewidth, or damping factor, or full width at half maximum (FWHM) of the core-excited/intermediate state , which is related to the lifetime of the state as . Physically it reflects the total spontaneous decay rate of the core-excited state (radiative + non-radiative). One should be aware that the equation is already in a simplified form under what is commonly known as the rotating wave approximation, where the second term (obtained by letting in the above equation) is discarded, a term that represents the alternate time-ordering of emission followed by absorption instead---which at resonance can be safely neglected.
For a specific final state , the scattering cross-section is proportional to the squared modulus of the scattering amplitude
The Dirac delta distribution serves to ensure energy conservation of the overall process. Note that there is no such limiting factor on the intermediate states---why these states are sometimes known as virtual states. The total cross-section for a given frequency is then given by the sum over all final states . The cross-sections give Lorentzian line-shapes, where determines the width, and so also possible overlap between adjacent excited states.
Practical difficulties & implementation details¶
Based on the KHD expression for the scattering amplitudes, we can note two difficulties: the matrix elements involve transitions between excited states and hence are, formally, outside the reach of linear response theory, and secondly the sum goes over all possible intermediate states, including the continuum. The former means we need to go to at least quadratic response to obtain such terms, however, in addition to the increased computational cost there are also issues related to the use of an adiabatic -kernel. Alternatively, in what is known as the pseudo-wavefunction formalism of linear-response TDDFT (LR-TDDFT) one can approximate even these excited to excited state transitions. The sum over all intermediate states, can, at resonance, be truncated to those states close in energy. Finally, seeing that we will need both core-excited and valence-excited states, neither one CVS calculation nor one ordinary linear-response calculation is sufficient---unless the latter is fully diagonalized, at which point one obtains both sets of states from one calculation. Fully diagonalizing the response matrix is, except for small systems/and or basis sets, not feasible, and some more efficient scheme to obtaining them is needed.
In light of this, there are two linear-response-based approaches implemented in VeloxChem (apart from fully diagonalizing), these are known as the two-shot- (2S) and restricted-subspace approximation (RSA) approaches. Both of these approaches, naturally, attempt to obtain both the valence- and core-excited states as efficiently as possible, while keeping as much of the physics intact as possible. It should be noted that the pros and cons of each with respect to each other is not obvious, and benchmarking is required for investigating which option to use, and when. Nevertheless, generally, these two approaches provide similar results, with appropriate orbital selection. In the following we will only be concerned with vertical RIXS. Inclusion of vibrational degrees of freedom is generally a prerequisite for quantitative comparison to experiment but the vertical approximation nonetheless often yield qualitatively correct spectra.
Before going into further detail, we will first compute the RIXS spectra of water with all three approaches and see how they compare, starting with the RSA. The first step is to run the RSA-LR-TDDFT computation to get our eigenvectors and values, with which we can obtain all quantities we need for the scattering amplitudes. A final consideration regarding the parameter is that it is state-dependent. However, a common approximation is to take as constant across the core-excited state manifold, leaving it a function of the element only---the value of which can be obtained experimentally. For our purposes of studying the K-edge of oxygen, we will use eV. is further related to how many intermediate states should be kept in the computation of the scattering amplitudes, but more on this at the end.
Restricted-subspace approximation (RSA) approach¶
Let us start with the RSA Cruz et al. (2021) approach, as this lends itself as a natural restriction on the orbital manifold for our purposes. The idea is to include only those orbitals deemed adequate for describing both the core- and valence-excited states. Here, we directly might take issue with the use of the word adequate, but again, defining it in this context is beyond the scope of this chapter. One restricts not only which of the occupied orbitals are active, but also the virtual space. The limit of RSA, in taking all orbitals into account is of course the full diagonalization scheme---meaning that one has to compute all states in-between the valence- and core-excited states---still, this does mean that RSA is systematically improvable.
We start by running the RSA-linear-response TDDFT calculation, note the attributes determining the active orbitals.
h2o_xyz = """3
C2v
O 0.000000000000 0.000000000000 0.000000000000
H 0.000000000000 0.740848095288 0.582094932012
H 0.000000000000 -0.740848095288 0.582094932012
"""
molecule = vlx.Molecule.read_xyz_string(h2o_xyz)
basis = vlx.MolecularBasis.read(molecule, 'def2-svpd')
scf_drv = vlx.ScfRestrictedDriver()
scf_drv.ostream.mute()
scf_drv.xcfun = 'rcam-b3lyp'
scf_res = scf_drv.compute(molecule, basis)nocc = molecule.number_of_alpha_electrons()
norb = scf_res['C_alpha'].shape[0]
nvir = norb - nocc
print(nocc,nvir,norb)5 34 39
orb_drv = vlx.OrbitalViewer()
orb_drv.plot(molecule, basis, scf_drv.mol_orbs)Output
rsp_drv = vlx.LinearResponseEigenSolver()
rsp_drv.restricted_subspace = True
rsp_drv.num_core_orbitals = 1
rsp_drv.num_valence_orbitals = 2
# 34 total virtual
rsp_drv.num_virtual_orbitals = 20
# 5 core-excited/intermediate states
rsp_drv.nstates = 2*20 + 5
rsa_rsp = rsp_drv.compute(molecule, basis, scf_res)Output
To note here is that the core orbitals one chooses defines the intermediate states, i.e., which states are summed over in RixsDriver. The specific states that these correspond to can be seen in the output of the rsa_drv.compute, where one can note that we still need to compute all states between our valence- and core-excited states of interest, although these are fewer than in the full case due to our reduced orbital space.
rsa_rsp['excitation_details'][['HOMO -> LUMO -0.9838'],
['HOMO -> LUMO+1 0.9528', 'HOMO -> LUMO+4 0.2432'],
['HOMO-1 -> LUMO 0.9592'],
['HOMO -> LUMO+3 0.9790'],
['HOMO -> LUMO+2 -0.9075', 'HOMO-1 -> LUMO+3 -0.3725'],
['HOMO-1 -> LUMO+1 0.9758'],
['HOMO -> LUMO+4 -0.9513', 'HOMO -> LUMO+1 0.2572'],
['HOMO-1 -> LUMO+2 -0.7623', 'HOMO -> LUMO+5 0.6244'],
['HOMO -> LUMO+5 -0.7679', 'HOMO-1 -> LUMO+2 -0.6266'],
['HOMO-1 -> LUMO+3 -0.8766',
'HOMO -> LUMO+2 0.3468',
'HOMO -> LUMO+8 -0.2722'],
['HOMO -> LUMO+6 -0.8759', 'HOMO-1 -> LUMO+4 -0.4489'],
['HOMO -> LUMO+7 -0.9675'],
['HOMO-1 -> LUMO+4 -0.8800', 'HOMO -> LUMO+6 0.4684'],
['HOMO-1 -> LUMO+5 -0.9668'],
['HOMO -> LUMO+8 -0.8460',
'HOMO-1 -> LUMO+7 0.4025',
'HOMO-1 -> LUMO+3 0.2530'],
['HOMO-1 -> LUMO+6 -0.9943'],
['HOMO-1 -> LUMO+8 -0.8234', 'HOMO -> LUMO+10 0.5419'],
['HOMO -> LUMO+9 -0.8629', 'HOMO -> LUMO+11 -0.4765'],
['HOMO -> LUMO+10 -0.8279', 'HOMO-1 -> LUMO+8 -0.5361'],
['HOMO-1 -> LUMO+7 -0.8889', 'HOMO -> LUMO+8 -0.3999'],
['HOMO -> LUMO+11 0.8470', 'HOMO -> LUMO+9 -0.4962'],
['HOMO-1 -> LUMO+10 0.9844'],
['HOMO-1 -> LUMO+9 -0.9185', 'HOMO-1 -> LUMO+11 -0.3674'],
['HOMO-1 -> LUMO+11 -0.9102', 'HOMO-1 -> LUMO+9 0.3837'],
['HOMO -> LUMO+12 0.9848'],
['HOMO -> LUMO+13 0.9658', 'HOMO-1 -> LUMO+14 0.2186'],
['HOMO -> LUMO+14 -0.8922', 'HOMO-1 -> LUMO+13 0.3973'],
['HOMO-1 -> LUMO+12 0.7782', 'HOMO -> LUMO+16 -0.6158'],
['HOMO -> LUMO+15 -0.9915'],
['HOMO-1 -> LUMO+13 0.8979',
'HOMO -> LUMO+14 0.3571',
'HOMO-1 -> LUMO+15 0.2046'],
['HOMO -> LUMO+16 -0.7754', 'HOMO-1 -> LUMO+12 -0.6218'],
['HOMO-1 -> LUMO+14 0.9701', 'HOMO -> LUMO+13 -0.2259'],
['HOMO-1 -> LUMO+16 -0.9948'],
['HOMO-1 -> LUMO+15 0.9585', 'HOMO -> LUMO+14 -0.2510'],
['HOMO -> LUMO+17 0.9969'],
['HOMO-1 -> LUMO+17 -0.9779', 'HOMO -> LUMO+19 0.2009'],
['HOMO -> LUMO+18 0.9972'],
['HOMO-1 -> LUMO+18 -0.9981'],
['HOMO-1 -> LUMO+19 0.9983'],
['HOMO -> LUMO+19 0.9769', 'HOMO-1 -> LUMO+17 0.2029'],
['core_1 -> LUMO 0.9742'],
['core_1 -> LUMO+1 -0.9208',
'core_1 -> LUMO+4 -0.3067',
'core_1 -> LUMO+11 -0.2115'],
['core_1 -> LUMO+2 -0.9750', 'core_1 -> LUMO+8 -0.2107'],
['core_1 -> LUMO+3 0.9702', 'core_1 -> LUMO+7 -0.2317'],
['core_1 -> LUMO+4 0.9213', 'core_1 -> LUMO+1 -0.3329']]Now, we can construct the RIXS spectrum by using the RixsDriver class, where we set the incoming photon energies to the first core-excited resonance---this is what the driver computes if the incoming photon energy is not initially set. Alternatively, one can identify the core-excited states from the excitation_details of the return object from the response calculation.
rixs_drv = vlx.RixsDriver()
rixs_drv.restricted_subspace = True
# incoming photon energy
# if the rixs_drv.photon_energy is not set, the photon energy is set to match
# the transition to the first core-excited state with oscillator strength >1e-3
rixs_drv.photon_energy = rsa_rsp['eigenvalues'][40]
# one can truncate the final states at a certain energy
rixs_drv.final_state_cutoff = 14 / vlx.hartree_in_ev()
# angle between incoming polarization vector and outgoing propagation vector
rixs_drv.theta = 0
# intermediate/core state lifetime, 0.16 for O 1s
rixs_drv.gamma = 0.16 / vlx.hartree_in_ev()
rixs_res_rsa = rixs_drv.compute(molecule, basis, scf_res, rsa_rsp)
print(rixs_res_rsa.keys())Output
Resonant Inelastic X-ray Scattering (RIXS) Setup
================================================
Scattering angle (theta) [rad] : 0.00
Lifetime broadening (gamma) [eV] : 0.16
Incoming photon energies [a.u.] : 19.1148
Incoming photon energies [eV] : 520.14
Number of intermediate states : 5
Number of final states : 10
Final state energy cutoff [eV] : 14.00
State index sets
Intermediate/core states : [40 ... 44]
Final/valence states : [0 ... 9]
Orbital sets
Core orbitals : [Core1]
Valence orbitals : [HOMO-1, HOMO]
Virtual orbitals : [LUMO ... LUMO+19]
* Info * Running RIXS calculation in the restricted-subspace approach
J. Phys. Chem. A, 2025, 129, 8783-8797, DOI: 10.1021/acs.jpca.5c04528
Phys. Chem. Chem. Phys., 2021, 23, 1835-1848, DOI: 10.1039/D0CP04726K
* Info * Computed RIXS cross-sections for 10 final states at photon energy: 520.14 eV.
RIXS cross-sections at incident X-ray energy 520.14 eV, energy-loss mode
------------------------------------------------------------------------
Ground State S0: 0.00000000 a.u. 0.00000 eV Cross-section 1.47e-06
Excited State S1: 0.28664083 a.u. 7.79989 eV Cross-section 1.36e-05
Excited State S2: 0.35582839 a.u. 9.68258 eV Cross-section 6.19e-08
Excited State S3: 0.36592044 a.u. 9.95720 eV Cross-section 5.20e-06
Excited State S4: 0.40630716 a.u. 11.05618 eV Cross-section 7.43e-09
Excited State S5: 0.43225725 a.u. 11.76232 eV Cross-section 1.38e-07
Excited State S6: 0.43655005 a.u. 11.87913 eV Cross-section 4.52e-08
Excited State S7: 0.44702116 a.u. 12.16407 eV Cross-section 3.95e-09
Excited State S8: 0.47474958 a.u. 12.91859 eV Cross-section 7.31e-09
Excited State S9: 0.48136239 a.u. 13.09854 eV Cross-section 6.83e-09
Excited State S10: 0.50497418 a.u. 13.74105 eV Cross-section 4.45e-08
dict_keys(['cross_sections', 'elastic_cross_sections', 'photon_energies', 'scattering_amplitudes', 'emission_energies', 'energy_losses', 'core_eigenvalues', 'core_osc_strengths', 'gamma_fwhm_ev', 'rsp_type'])
Full diagonalization of response matrix¶
Due to the small system, and rather small basis set, we can also run the full diagonalization, also with the use of restricted_subspace variant, as this allows us to implicitly define the intermediate states as the core-excited ones. Without the tag on the states of core, the RixsDriver does not know which states should be the intermediate ones. We can also skip passing explicit response results and instead have the RixsDriver compute the response needed.
rixs_drv = vlx.RixsDriver()
rixs_drv.restricted_subspace = True
rixs_drv.num_core_orbitals = 1
rixs_drv.num_valence_orbitals = 4
rixs_drv.num_virtual_orbitals = nvir
rixs_drv.nstates = 34 * (4+1) - 29
# incoming photon energy
# if the rixs_drv.photon_energy is not set, the photon energy is set to match
# the transition to the first core-excited state with oscillator strength >1e-3
rixs_drv.photon_energy = None
# one can truncate the final states at a certain energy
rixs_drv.final_state_cutoff = 13.5 / vlx.hartree_in_ev()
# angle between incoming polarization vector and outgoing propagation vector
rixs_drv.theta = 0
# intermediate/core state lifetime, 0.16 for O 1s
rixs_drv.gamma = 0.16 / vlx.hartree_in_ev()
rixs_res_fd = rixs_drv.compute(molecule, basis, scf_res, fd_rsp)Output
Resonant Inelastic X-ray Scattering (RIXS) Setup
================================================
Scattering angle (theta) [rad] : 0.00
Lifetime broadening (gamma) [eV] : 0.16
Incoming photon energies [a.u.] : 19.1146
Incoming photon energies [eV] : 520.13
Number of intermediate states : 5
Number of final states : 10
Final state energy cutoff [eV] : 13.50
State index sets
Intermediate/core states : [136 ... 140]
Final/valence states : [0 ... 9]
Orbital sets
Core orbitals : [Core1]
Valence orbitals : [HOMO-3 ... HOMO]
Virtual orbitals : [LUMO ... LUMO+33]
* Info * Running RIXS calculation in the restricted-subspace approach
J. Phys. Chem. A, 2025, 129, 8783-8797, DOI: 10.1021/acs.jpca.5c04528
Phys. Chem. Chem. Phys., 2021, 23, 1835-1848, DOI: 10.1039/D0CP04726K
* Info * Computed RIXS cross-sections for 10 final states at photon energy: 520.13 eV.
RIXS cross-sections at incident X-ray energy 520.13 eV, energy-loss mode
------------------------------------------------------------------------
Ground State S0: 0.00000000 a.u. 0.00000 eV Cross-section 1.71e-06
Excited State S1: 0.28563706 a.u. 7.77258 eV Cross-section 1.47e-05
Excited State S2: 0.35546905 a.u. 9.67281 eV Cross-section 6.58e-08
Excited State S3: 0.36464290 a.u. 9.92244 eV Cross-section 5.59e-06
Excited State S4: 0.40509763 a.u. 11.02327 eV Cross-section 9.21e-09
Excited State S5: 0.42394828 a.u. 11.53622 eV Cross-section 1.02e-07
Excited State S6: 0.43539474 a.u. 11.84769 eV Cross-section 5.02e-08
Excited State S7: 0.44687617 a.u. 12.16012 eV Cross-section 4.15e-09
Excited State S8: 0.47439462 a.u. 12.90894 eV Cross-section 1.15e-08
Excited State S9: 0.48041384 a.u. 13.07273 eV Cross-section 4.52e-09
Excited State S10: 0.49399553 a.u. 13.44230 eV Cross-section 3.72e-08
Note the tiny but non-zero difference in energy of the first core-excited state obtained here, eV, and the one from the RSA, eV.
Two-shot approach¶
Next we turn to the two-shot (2S) approach, developed by Nascimento et al. (2021) [Nascimento et al. (2021), DOI: 10.1021/acs.jctc.1c00144], where two separate linear-response calculations are employed to obtain the eigenvalues and eigenvectors required to compute the scattering amplitudes. One CVS-linear-response computation is performed to determine the core-excited states, and one normal linear-response computation for the valence-excited states of interest. Strictly speaking, the CVS approximation means that the core-excited states are eigenvectors of a different Hamiltonian---why these eigenvectors must be padded with zeros---hence these are not necessarily orthogonal to eigenvectors from the calculation without the CVS approximation but this should be a small effect.
# normal LR-TDDFT
rsp_drv = vlx.LinearResponseEigenSolver()
rsp_drv.nstates = 25
valence_rsp = rsp_drv.compute(molecule, basis, scf_res)Output
rsp_drv = vlx.LinearResponseEigenSolver()
rsp_drv.core_excitation = True
rsp_drv.num_core_orbitals = 1
rsp_drv.nstates = 5
cvs_rsp = rsp_drv.compute(molecule, basis, scf_res)Output
Linear Response EigenSolver Setup
===================================
Number of States : 5
Max. Number of Iterations : 150
Convergence Threshold : 1.0e-04
ERI Screening Threshold : 1.0e-12
Exchange-Correlation Functional : RCAM-B3LYP
Molecular Grid Level : 4
* Info * Using the RCAM-B3LYP functional.
A. J. Cohen, P. Mori-Sánchez, and W. Yang, J. Chem. Phys. 126, 191109 (2007)
* Info * Using the Libxc library (v7.1.2).
S. Lehtola, C. Steigemann, M. J. T. Oliveira, and M. A. L. Marques, SoftwareX 7, 1–5 (2018)
* Info * Using the following algorithm for XC numerical integration.
J. Kussmann, H. Laqua and C. Ochsenfeld, J. Chem. Theory Comput. 2021, 17, 1512-1521
* Info * 15 gerade trial vectors in reduced space
* Info * 15 ungerade trial vectors in reduced space
*** Iteration: 1 * Residuals (Max,Min): 9.74e-02 and 5.00e-02
* Info * 20 gerade trial vectors in reduced space
* Info * 20 ungerade trial vectors in reduced space
*** Iteration: 2 * Residuals (Max,Min): 2.31e-02 and 1.59e-03
* Info * 25 gerade trial vectors in reduced space
* Info * 24 ungerade trial vectors in reduced space
*** Iteration: 3 * Residuals (Max,Min): 2.83e-03 and 1.11e-04
* Info * 29 gerade trial vectors in reduced space
* Info * 29 ungerade trial vectors in reduced space
*** Iteration: 4 * Residuals (Max,Min): 1.67e-04 and 7.54e-14
* Info * 30 gerade trial vectors in reduced space
* Info * 30 ungerade trial vectors in reduced space
*** Iteration: 5 * Residuals (Max,Min): 1.74e-05 and 6.64e-14
*** Linear response converged in 5 iterations. Time: 0.91 sec
Electric Transition Dipole Moments (dipole length, a.u.)
--------------------------------------------------------
X Y Z
Excited State S1: 0.000000 0.000000 0.027385
Excited State S2: -0.000000 0.043168 0.000000
Excited State S3: 0.031564 0.000000 -0.000000
Excited State S4: -0.000000 0.000000 -0.028115
Excited State S5: -0.000000 0.017823 -0.000000
Electric Transition Dipole Moments (dipole velocity, a.u.)
----------------------------------------------------------
X Y Z
Excited State S1: 0.000000 0.000000 0.025923
Excited State S2: -0.000000 0.042291 0.000000
Excited State S3: 0.033281 0.000000 -0.000000
Excited State S4: -0.000000 0.000000 -0.026255
Excited State S5: -0.000000 0.016801 -0.000000
Magnetic Transition Dipole Moments (a.u.)
-----------------------------------------
X Y Z
Excited State S1: 0.000000 -0.000000 -0.000000
Excited State S2: -0.000399 -0.000000 0.000000
Excited State S3: -0.000000 0.000135 -0.000000
Excited State S4: 0.000000 -0.000000 -0.000000
Excited State S5: -0.000245 -0.000000 -0.000000
One-Photon Absorption
---------------------
Excited State S1: 19.11457626 a.u. 520.13412 eV Osc.Str. 0.0096
Excited State S2: 19.18269122 a.u. 521.98762 eV Osc.Str. 0.0238
Excited State S3: 19.23278737 a.u. 523.35081 eV Osc.Str. 0.0128
Excited State S4: 19.23630445 a.u. 523.44651 eV Osc.Str. 0.0101
Excited State S5: 19.27780803 a.u. 524.57588 eV Osc.Str. 0.0041
Electronic Circular Dichroism
-----------------------------
Excited State S1: Rot.Str. -0.000000 a.u. -0.0000 [10**(-40) cgs]
Excited State S2: Rot.Str. 0.000000 a.u. 0.0000 [10**(-40) cgs]
Excited State S3: Rot.Str. -0.000000 a.u. -0.0000 [10**(-40) cgs]
Excited State S4: Rot.Str. 0.000000 a.u. 0.0000 [10**(-40) cgs]
Excited State S5: Rot.Str. 0.000000 a.u. 0.0000 [10**(-40) cgs]
Character of excitations:
Excited state 1
---------------
core_1 -> LUMO -0.9740
Excited state 2
---------------
core_1 -> LUMO+1 0.9204
core_1 -> LUMO+4 0.3067
core_1 -> LUMO+11 0.2127
Excited state 3
---------------
core_1 -> LUMO+2 -0.9750
core_1 -> LUMO+8 -0.2106
Excited state 4
---------------
core_1 -> LUMO+3 -0.9678
core_1 -> LUMO+7 0.2389
Excited state 5
---------------
core_1 -> LUMO+4 0.9211
core_1 -> LUMO+1 -0.3330
rixs_drv = vlx.RixsDriver()
# incoming photon energy
rixs_drv.photon_energy = [cvs_rsp['eigenvalues'][0], cvs_rsp['eigenvalues'][1]]
# angle between incoming polarization vector and outgoing propagation vector
rixs_drv.theta = 0
# intermediate/core state lifetime, 0.16 for O 1s
rixs_drv.gamma = 0.16 / vlx.hartree_in_ev()
rixs_res_2s = rixs_drv.compute(molecule, basis, scf_res, valence_rsp, cvs_rsp)Output
Resonant Inelastic X-ray Scattering (RIXS) Setup
================================================
Scattering angle (theta) [rad] : 0.00
Lifetime broadening (gamma) [eV] : 0.16
Incoming photon energies [a.u.] : 19.1146, 19.1827
Incoming photon energies [eV] : 520.13, 521.99
Number of intermediate states : 5
Number of final states : 25
State index sets
Intermediate/core states : [0 ... 4]
Final/valence states : [0 ... 24]
Orbital sets
Core orbitals : [Core1]
Valence orbitals : [HOMO-3 ... HOMO]
Virtual orbitals : [LUMO ... LUMO+33]
* Info * Running RIXS calculation in the two-shot approach
J. Phys. Chem. A, 2025, 129, 8783-8797, DOI: 10.1021/acs.jpca.5c04528
* Info * Computed RIXS cross-sections for 25 final states at photon energy: 520.13 eV.
RIXS cross-sections at incident X-ray energy 520.13 eV, energy-loss mode
------------------------------------------------------------------------
Ground State S0: 0.00000000 a.u. 0.00000 eV Cross-section 1.65e-06
Excited State S1: 0.28563706 a.u. 7.77258 eV Cross-section 1.43e-05
Excited State S2: 0.35546905 a.u. 9.67281 eV Cross-section 6.37e-08
Excited State S3: 0.36464290 a.u. 9.92244 eV Cross-section 5.38e-06
Excited State S4: 0.40509763 a.u. 11.02327 eV Cross-section 9.01e-09
Excited State S5: 0.42394828 a.u. 11.53622 eV Cross-section 1.01e-07
Excited State S6: 0.43539474 a.u. 11.84769 eV Cross-section 4.91e-08
Excited State S7: 0.44687617 a.u. 12.16012 eV Cross-section 4.04e-09
Excited State S8: 0.47439462 a.u. 12.90894 eV Cross-section 1.08e-08
Excited State S9: 0.48041384 a.u. 13.07273 eV Cross-section 5.09e-09
Excited State S10: 0.49399553 a.u. 13.44230 eV Cross-section 3.84e-08
Excited State S11: 0.50233182 a.u. 13.66915 eV Cross-section 9.53e-06
Excited State S12: 0.51836208 a.u. 14.10535 eV Cross-section 8.21e-08
Excited State S13: 0.53324470 a.u. 14.51033 eV Cross-section 1.15e-08
Excited State S14: 0.54128381 a.u. 14.72908 eV Cross-section 3.78e-07
Excited State S15: 0.55373243 a.u. 15.06783 eV Cross-section 1.53e-09
Excited State S16: 0.56161556 a.u. 15.28234 eV Cross-section 6.44e-09
Excited State S17: 0.59232073 a.u. 16.11787 eV Cross-section 1.46e-11
Excited State S18: 0.59376072 a.u. 16.15705 eV Cross-section 9.77e-09
Excited State S19: 0.62237950 a.u. 16.93581 eV Cross-section 8.14e-09
Excited State S20: 0.62952796 a.u. 17.13033 eV Cross-section 2.03e-08
Excited State S21: 0.63498655 a.u. 17.27886 eV Cross-section 2.55e-09
Excited State S22: 0.64402344 a.u. 17.52477 eV Cross-section 2.63e-08
Excited State S23: 0.64619582 a.u. 17.58388 eV Cross-section 1.92e-10
Excited State S24: 0.64834114 a.u. 17.64226 eV Cross-section 1.10e-08
Excited State S25: 0.68067851 a.u. 18.52221 eV Cross-section 6.56e-11
* Info * Computed RIXS cross-sections for 25 final states at photon energy: 521.99 eV.
RIXS cross-sections at incident X-ray energy 521.99 eV, energy-loss mode
------------------------------------------------------------------------
Ground State S0: 0.00000000 a.u. 0.00000 eV Cross-section 1.03e-05
Excited State S1: 0.28563706 a.u. 7.77258 eV Cross-section 2.69e-08
Excited State S2: 0.35546905 a.u. 9.67281 eV Cross-section 3.51e-05
Excited State S3: 0.36464290 a.u. 9.92244 eV Cross-section 1.72e-08
Excited State S4: 0.40509763 a.u. 11.02327 eV Cross-section 4.55e-08
Excited State S5: 0.42394828 a.u. 11.53622 eV Cross-section 2.16e-07
Excited State S6: 0.43539474 a.u. 11.84769 eV Cross-section 2.67e-05
Excited State S7: 0.44687617 a.u. 12.16012 eV Cross-section 2.04e-07
Excited State S8: 0.47439462 a.u. 12.90894 eV Cross-section 2.88e-08
Excited State S9: 0.48041384 a.u. 13.07273 eV Cross-section 2.21e-08
Excited State S10: 0.49399553 a.u. 13.44230 eV Cross-section 5.90e-07
Excited State S11: 0.50233182 a.u. 13.66915 eV Cross-section 3.82e-08
Excited State S12: 0.51836208 a.u. 14.10535 eV Cross-section 1.37e-08
Excited State S13: 0.53324470 a.u. 14.51033 eV Cross-section 1.55e-10
Excited State S14: 0.54128381 a.u. 14.72908 eV Cross-section 1.01e-06
Excited State S15: 0.55373243 a.u. 15.06783 eV Cross-section 1.13e-07
Excited State S16: 0.56161556 a.u. 15.28234 eV Cross-section 1.81e-06
Excited State S17: 0.59232073 a.u. 16.11787 eV Cross-section 3.73e-09
Excited State S18: 0.59376072 a.u. 16.15705 eV Cross-section 6.09e-06
Excited State S19: 0.62237950 a.u. 16.93581 eV Cross-section 4.53e-08
Excited State S20: 0.62952796 a.u. 17.13033 eV Cross-section 3.45e-08
Excited State S21: 0.63498655 a.u. 17.27886 eV Cross-section 2.00e-10
Excited State S22: 0.64402344 a.u. 17.52477 eV Cross-section 1.75e-06
Excited State S23: 0.64619582 a.u. 17.58388 eV Cross-section 7.55e-08
Excited State S24: 0.64834114 a.u. 17.64226 eV Cross-section 1.06e-10
Excited State S25: 0.68067851 a.u. 18.52221 eV Cross-section 2.89e-08
Note how the CVS approximation yields a core-excited energy closer to the full-matrix one than the RSA did.
Comparison¶
Conceptually, in the 2S approach, via the use of the core-valence separation (CVS) approximation, the coupling between the core- and valence-excited states are neglected in the calculation of the core-excited states, impacting both energies and properties. However, this approach is motivated by that the CVS approximation has been thoroughly investigated and found to perturb these only slightly, supported by the small couplings between valence and core orbitals. The RSA, captures some of this coupling---at the expense of neglecting the couplings to parts of the virtual-state manifold instead, this can make the error somewhat harder to predict.
Source
fig, ax = plt.subplots(1, figsize=(10, 6))
valence_enes = valence_rsp['eigenvalues'] * vlx.hartree_in_ev()
valence_osc = valence_rsp['oscillator_strengths']
# broaden RIXS stick spectrum
rsa_x, rsa_y = add_broadening(rixs_res_rsa['energy_losses'][:,0]*vlx.hartree_in_ev(), rixs_res_rsa['cross_sections'][:,0], interval=(1, 15), step=.01, fwhm=.3)
# broaden RIXS at first core resonance
fd_x, fd_y = add_broadening(rixs_res_fd['energy_losses'][:,0]*vlx.hartree_in_ev(), rixs_res_fd['cross_sections'][:,0], interval=(1, 15), step=.01, fwhm=.3)
# broaden RIXS at the first core resonance
twos_x, twos_y = add_broadening(rixs_res_2s['energy_losses'][:,0]*vlx.hartree_in_ev(), rixs_res_2s['cross_sections'][:,0], interval=(1, 15), step=.01, fwhm=.3)
ax.plot(fd_x, fd_y, label='FD', linewidth=2.5)
ax.plot(rsa_x, rsa_y, label='RSA', linestyle='dashed', linewidth=2.5)
ax.plot(twos_x, twos_y, label='Two-shot', linestyle='dotted', linewidth=2.5)
ax.set_xlim(12,6)
ax.set_xlabel('Energy loss (eV)', fontsize=15)
ax.set_yticks([])
ax.set_ylabel('Intensity (arbitrary units)', fontsize=15)
ax.legend(fontsize=13)
plt.show()
The full diagonalization approach serves as the reference in this case. One can see how the two-shot performs slightly better than the RSA---for this particular choice of orbital space.
Targeting dipole-forbidden excited states with RIXS¶
As a second-order process, RIXS is not constrained by the usual one-photon selection rules. This is because of the indirect path to the final state . Even though the transition from might be dipole forbidden, transitioning to an intermediate core-excited state can open up channels to final states that are otherwise dark in the optical spectra. To see this more concretely, we can take a specific example from our computed water UV/Vis, and RIXS spectra. First, note that the transition to the second excited singlet state, , is dipole forbidden, that is
Source
print(f"S2: E = {valence_rsp['eigenvalues'][1]*vlx.hartree_in_ev():.2f} eV, Oscillator strength = {valence_rsp['oscillator_strengths'][1]:.3f}")S2: E = 9.67 eV, Oscillator strength = 0.000
This could have already been deduced by group theory, as the transition is essentially an excitation from HOMO LUMO+1, with the respective orbitals transforming as the and irreps. The direct product is , and since none of the Cartesian components transform as A we can deduce that this transition is dipole-forbidden, as that means that does not contain the totally symmetric irrep, , for any component , .
However, for the 1s oxygen orbital, belonging to the totally symmetric irrep, we can note that the 1s to LUMO+1 transition is not forbidden---with the hole in the 1s orbital and electron in LUMO+1 comprising the second core-excited state. If the now core-excited state relaxes back via a transition from HOMO 1s, the system winds up in the state. Via a similar analysis, we have, for the absorption step, (corresponding to ) we see that must belong to , that is, polarization along , i.e., . For the emission step, we instead have the HOMO (belonging to ) and this leads to required polarization along . If this is correct, we should have non-zero amplitudes of , for on the second core-excited resonance frequency, here indicated by (). The frequency dependence means one can target specific core-excited states and deduce symmetries and polarization dependence of the states involved. This has further implications: for instance, for solvated molecules it might be that specific absorption lines in a regular UV/Vis spectrum can be heavily mixed by absorption lines of the solvent. With RIXS one can still probe the solute, specifically, by going the indirect path to that same valence spectrum by tuning the frequency to a specific core-excitation instead---allowing for a “purification” of the spectrum (the chemical shift of the core orbitals further allows for sensitivity to the local chemical environment as well).
Source
# F_fg^xy.shape: (nr_final_states, nr_intermediate_states, 3, 3)
print(f'F_fg = \n{rixs_res_2s['scattering_amplitudes'][1,1]}')
print(f"\nF^xy: {rixs_res_2s['scattering_amplitudes'][1,1,0,1]:.2f}")F_fg =
[[ 0. +0.j -0.3024-307.6483j -0. -0.j ]
[-0.1055 +0.0062j -0. -0.j 0. -0.j ]
[ 0. -0.j -0. -0.j -0. -0.j ]]
F^xy: -0.30-307.65j
Indeed, only the -component has substantial amplitude! Non-resonant scattering contributes to the permuted entry , but the resonant contribution dominates. However, such a clean separation is not what is generally observed experimentally. Among the difficulties with these experiments and their comparison to theory is that samples in the gas- or liquid phase are randomly oriented, and additionally, that the emitted photon polarization is seldom experimentally measured. Hence, one typically averages over molecular orientations as well as over outgoing photon polarization, yielding a cross-section that depends on the angle ---between the incoming photon polarization vector and outgoing photon propagation vector, . This gives the averaged cross-section, for a particular final state
We can now plot the spectra to get a better view of what we have discussed so far
Source
fig, ax = plt.subplots(1, figsize=(10, 6))
enes = rixs_res_2s['energy_losses'] * vlx.hartree_in_ev()
xsection = rixs_res_2s['cross_sections']
valence_enes = valence_rsp['eigenvalues'] * vlx.hartree_in_ev()
valence_osc = valence_rsp['oscillator_strengths']
# broaden UV/vis stick spectrum
val_x, val_y = add_broadening(valence_enes, valence_osc, interval=(1, 15), step=.01, fwhm=.3)
# broaden RIXS at first core resonance
x1, y1 = add_broadening(enes[:,0], xsection[:,0], interval=(1, 15), step=.01, fwhm=.3)
# broaden RIXS at second core resonance
x2, y2 = add_broadening(enes[:,1], xsection[:,1], interval=(1, 15), step=.01, fwhm=.3)
# area-normalize the spectra for ease of comparison
area1 = scipy.integrate.simpson(y1)
area2 = scipy.integrate.simpson(y2)
area_val = scipy.integrate.simpson(val_y)
ax.plot(val_x, val_y/area_val, label='One-photon',linewidth=2.5)
ax.fill_between(val_x, 0, val_y/area_val, alpha=.3)
ax.plot(x1, y1/area1, label='1st core res.', color='orange',linewidth=2.5)
ax.plot(x2, y2/area2, label='2nd core res.', color='darkred',linewidth=2.5)
ax.set_xlim(14,6)
ax.set_xlabel('Energy loss (eV)', fontsize=15)
ax.set_yticks([])
ax.set_ylabel('Intensity (area normalized)', fontsize=15)
ax.text(valence_rsp['eigenvalues'][1] * vlx.hartree_in_ev() - .125, .0115, 'S2', fontsize=15)
ax.axvline(valence_rsp['eigenvalues'][1] * vlx.hartree_in_ev(), color='black', linewidth=1.2, ls='--')
ax.legend(fontsize=13)
plt.show()
We can observe the second excited state when the incoming photon energy is tuned to the second core-excited (CE) state resonance, in contrast to the valence- and RIXS (at the first core-excited resonance) spectra for which is dark.
RIXS map¶
The RIXS spectra we have looked at so far has been for one, or two, particular incoming photon frequencies. A common way of representing RIXS is by plotting the so-called RIXS-map, . The previous spectra, the ones at the first- and second CE resonances, represent one-dimensional slices of this map, that is, the map describes spectra for a range of possible incoming frequencies. We can construct such a map, starting with giving a range of input frequencies to the RIXS driver like in the following, where we take as range the values between our two CE resonances (and some below and above).
# step size between incoming frequencies [eV]
step = 0.1/2
core_eigvals_ev = cvs_rsp['eigenvalues'] * vlx.hartree_in_ev()
incoming_photon_enes = np.arange(core_eigvals_ev[0] - 1, core_eigvals_ev[2] + 1, step)
# give photon energies in a.u.
rixs_drv.photon_energy = list(incoming_photon_enes / vlx.hartree_in_ev())
rixs_drv.ostream.mute()
rixs_map_res = rixs_drv.compute(molecule, basis, scf_res, valence_rsp, cvs_rsp)Source
n_inc = incoming_photon_enes.shape[0]
rixs_map = np.empty((1400, n_inc))
for i in range(n_inc):
x, rixs_map[:, i] = add_broadening(
rixs_map_res["energy_losses"][:, i] * vlx.hartree_in_ev(),
rixs_map_res["cross_sections"][:, i],
interval=(1, 15),
step=0.01,
fwhm=0.3,
)
Z = rixs_map.T
y = incoming_photon_enes
xas_x, xas_y = add_broadening(
core_eigvals_ev,
cvs_rsp["oscillator_strengths"],
interval=(core_eigvals_ev[0] - 1, core_eigvals_ev[1] + 1),
step=0.01,
fwhm=0.3,
)
fig, (ax1, ax2) = plt.subplots(
1, 2,
sharey=True,
figsize=(10, 6),
gridspec_kw={"wspace": 0, "width_ratios": [4, 1.25]},
)
ax1.imshow(
Z,
extent=[x.min(), x.max(), y.min(), y.max()],
origin="lower",
aspect="auto",
cmap="cividis",
)
ax1.set_xlabel("Energy loss (eV)", fontsize=14)
ax1.set_ylabel("Photon energy (eV)", fontsize=14)
ax2.plot(xas_y, xas_x, lw=2.5)
ax2.fill_betweenx(xas_x, 0, xas_y, alpha=0.3)
ax1.set_xlim(x.max(), 5)
ax1.set_ylim(y.min(), core_eigvals_ev[1] + 0.75)
ax2.set_xlim(0)
ax2.set_xticks([])
ax1.set_title("RIXS map", fontsize=14)
ax2.set_title("XAS", fontsize=14)
ax2.set_xlabel("Intensity", fontsize=14)
ax2.yaxis.set_visible(False)
for idx in (0, 1):
E = core_eigvals_ev[idx]
ax1.axhline(E, ls="--", lw=1.2, color="w")
I_E = np.interp(E, xas_x, xas_y)
ax2.hlines(E, 0, I_E, linestyles="--", linewidth=1.2, colors="w")
plt.show()
The dashed lines indicate the 1-dimensional cuts we saw in the previous figure.
Effect of intermediate state lifetime¶
As a final curiosity we will have a short look at the role played by the intermediate state lifetime , in relation to the effect of detuning, i.e., of going off resonance. Detuning is usually written and is related to the notion of the scattering duration. As we saw initially, sets the effective width of the Lorentzians fo the cross-section contributions. This is closely related to interference between scattering channels. Note that the cross-section can be written
with the partial channel . Interference is contained in the off-diagonal terms. These can be important to include, depending on the detuning and core-hole lifetime. Discarding the off-diagonal terms serves as a valid approximation provided the intermediate states contributing non-negligibly to a given final state are well-separated, that is: , with the energy difference between adjacent core-excited states. For light elements, in particular, at or near resonance, the condition often holds, where the cross-terms become negligible in comparison to the diagonal, as each Lorentzian is effectively isolated. Neglecting the interference means we assume that coherence is lost and that the process can be described as absorption + emission---not as a coherent scattering process---so the amplitudes factorize.
We will not go into too much more detail here but we can have a look at what happens when becomes very (unphysically) large---on the order of ---for three frequencies, two on resonance and the final one the average of the two.
rixs_drv = vlx.RixsDriver()
rixs_drv.gamma = 0.16 / vlx.hartree_in_ev()
rixs_drv.restricted_subspace = False
rixs_drv.photon_energy = [cvs_rsp['eigenvalues'][0], cvs_rsp['eigenvalues'][1], (cvs_rsp['eigenvalues'][0] + cvs_rsp['eigenvalues'][1])/2]
rixs_res = rixs_drv.compute(molecule, basis, scf_res, valence_rsp, cvs_rsp)
rixs_drv.gamma = 1.5 / vlx.hartree_in_ev()
rixs_res_large_gamma = rixs_drv.compute(molecule, basis, scf_res, valence_rsp, cvs_rsp)Output
Source
fig, (ax1,ax2) = plt.subplots(1,2, figsize=(12, 6))
enes = rixs_res['energy_losses'] * vlx.hartree_in_ev()
xsection = rixs_res['cross_sections']
xsection_wide = rixs_res_large_gamma['cross_sections']
# broaden RIXS at first core resonance
x1, y1 = add_broadening(enes[:,0], xsection[:,0], interval=(1, 15), step=.01, fwhm=.3)
# broaden RIXS at second core resonance
x2, y2 = add_broadening(enes[:,1], xsection[:,1], interval=(1, 15), step=.01, fwhm=.3)
# broaden RIXS at the mean photon energy
x3, y3 = add_broadening(enes[:,1], xsection[:,2], interval=(1, 15), step=.01, fwhm=.3)
x1w, y1w = add_broadening(enes[:,0], xsection_wide[:,0], interval=(1, 15), step=.01, fwhm=.3)
# broaden RIXS at second core resonance
x2w, y2w = add_broadening(enes[:,1], xsection_wide[:,1], interval=(1, 15), step=.01, fwhm=.3)
# broaden RIXS at the mean photon energy
x3w, y3w = add_broadening(enes[:,1], xsection_wide[:,2], interval=(1, 15), step=.01, fwhm=.3)
ax1.plot(x1, y1, label='1st CE res.', lw=2.5)
ax1.plot(x2, y2, label='2nd CE res.', lw=2.5)
ax1.plot(x3, y3, label='Avg. CE res.', lw=2.5)
ax2.plot(x1w, y1w, label='1st CE res.', lw=2.5)
ax2.plot(x2w, y2w, label='2nd CE res.', lw=2.5)
ax2.plot(x3w, y3w, label='Avg. CE res.', lw=2.5)
ax1.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
ax2.ticklabel_format(axis='y', style='sci', scilimits=(0,0))
ax1.set_xlim(14,6)
ax2.set_xlim(14,6)
ax1.set_xlabel('Energy loss (eV)', fontsize=14)
ax2.set_xlabel('Energy loss (eV)', fontsize=14)
ax1.set_ylabel('Intensity (arbitrary units)', fontsize=14)
ax1.set_title('Lifetime: 0.16 eV', fontsize=14)
ax2.set_title('Lifetime: 1.5 eV', fontsize=14)
ax1.legend(fontsize=12)
plt.show()
One can clearly see that for the case of the large the spectra start mixing---approaching a common average spectral shape---even the detuned spectrum has gained considerable intensity, or rather the other two have lost intensity. In such a situation interference cannot be neglected, as the spectrum departs from a simple sum of isolated Lorentzians.
- 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
- 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
- Lykos, P. G., & Schmeising, H. N. (1961). Maximum Overlap Atomic and Molecular Orbitals. J. Phys. Chem., 35, 288. 10.1063/1.1731901
- 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
- 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
- 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
- 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
- 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
- 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
- Cederbaum, L. S., Domcke, W., & Schirmer, J. (1980). Many-body theory of core holes. Phys. Rev. A, 22, 206. 10.1103/PhysRevA.22.206
- 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
- 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
- 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
- Herbst, M. F., & Fransson, T. (2020). Quantifying the error of the core-valence separation approximation. J. Chem. Phys., 153, 054114. 10.1063/5.0013538
- Rehn, D. R., Rinkevicius, Z., Herbst, M., Li, M., X. Scheurer, Brand, M., Dempwolff, A. L., Brumboiu, I. E., Fransson, T., Dreuw, A., & Norman, P. (2021). Gator: a Python-driven program for spectroscopy simulations using correlated wave functions. WIREs Comput. Mol. Sci., e1528. 10.1002/wcms.1528