{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Special aspects\n",
"\n",
"In this section we will discuss a number of issues and questions relating to core spectroscopies, ranging from the physical effect of a core-transition, to technical aspects relating to lowering computational costs. They are ordered by approximate relevance, with some topics being pertient only for certain cases.\n",
"\n",
"\n",
"(sec:xray_topics_relax)=\n",
"## Relaxation\n",
"\n",
"Core spectroscopies exhibit strong relaxation effects, originating in the significant change in shielding of the nuclear charge arising from the creation/annihilation of a core-hole. These effects need to be properly accounted for in order to get good absolute energies, as the lack of this will mean that, *e.g.* an unrelaxed core-hole system is too high in energy. This explains the previously [observed](sec:xray_calc_iekoop) large error obtained using Koopmans' theorem.\n",
"\n",
"### Relaxation of core-ionized systems\n",
"\n",
"In {numref}`relax_form` we see the effects of relaxation for creating a core-hole in the oxygen/carbon 1s of formaldehyde:\n",
"\n",
"```{figure} ../../img/xray/relax_form.png\n",
"---\n",
"name: relax_form\n",
"---\n",
"The total electron density of formaldehyde, featuring the neutral system and core-ionized at the carbon and oxygen 1s.\n",
"```\n",
"\n",
"We see that there is a noticable attraction of electron density toward the core-holes, as marked with a **+**. These effects can be modeled with the $Z+1$ approximation, within which the change in screening is modeled by substituting the ionized atom with the next element in the periodic table. This can be illustrated with the following radial distribution plots (using uncontracted 6-31G basis sets, in order to increase flexibility in the core region):\n",
"\n",
"```{note}\n",
"To be added.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"remove-cell"
]
},
"source": [
"```python\n",
"water_xyz = \"\"\"\n",
"O 0.0000000000 0.0000000000 0.1178336003\n",
"H -0.7595754146 -0.0000000000 -0.4713344012\n",
"H 0.7595754146 0.0000000000 -0.4713344012\n",
"\"\"\"\n",
"water_fluorine_xyz = \"\"\"\n",
"F 0.0000000000 0.0000000000 0.1178336003\n",
"H -0.7595754146 -0.0000000000 -0.4713344012\n",
"H 0.7595754146 0.0000000000 -0.4713344012\n",
"\"\"\"\n",
"\n",
"### SCF of neutral and core-hole water\n",
"# Create pyscf mol object\n",
"mol = gto.Mole()\n",
"mol.atom = water_xyz\n",
"mol.basis = \"unc-6-31G\"\n",
"mol.build()\n",
"# Perform unrestricted SCF calculation\n",
"scf_gs = scf.UHF(mol)\n",
"scf_gs.kernel()\n",
"# Copy molecular orbitals\n",
"mo0 = copy.deepcopy(scf_gs.mo_coeff)\n",
"occ0 = copy.deepcopy(scf_gs.mo_occ)\n",
"# Create 1s core-hole by setting beta_0 population to zero\n",
"occ0[1][0] = 0.0\n",
"# Perform unrestricted SCF calculation with MOM constraint\n",
"scf_ion = scf.UHF(mol)\n",
"scf.addons.mom_occ(scf_ion, mo0, occ0)\n",
"scf_ion.kernel()\n",
"\n",
"### SCF of FO2 (Z+1 model)\n",
"# Create pyscf mol object\n",
"mol_fl = gto.Mole()\n",
"mol_fl.atom = water_fluorine_xyz\n",
"mol_fl.basis = \"unc-6-31G\"\n",
"mol_fl.charge = 1 # Closed-shell\n",
"mol_fl.build()\n",
"# Perform unrestricted SCF calculation\n",
"scf_gs_fl = scf.UHF(mol_fl)\n",
"scf_gs_fl.kernel()\n",
"\n",
"Ay = 0.0 # oxygen y-position (in au)\n",
"Az = 0.1178336003 / 0.529177249 # oxygen z-position (in au)\n",
"\n",
"# Create coordinate object along x-axis (intersecting with O/F atom)\n",
"coords = []\n",
"for ix in np.arange(-4, 4, 0.001):\n",
" coords.append((ix, Ay, Az))\n",
"coords = np.array(coords)\n",
"\n",
"# AO values and MO values on given grids\n",
"# Neutral water\n",
"ao_gs = mol.eval_gto(\"GTOval_sph\", coords)\n",
"mo_gs = ao_gs.dot(scf_gs.mo_coeff)\n",
"\n",
"# Water core-hole\n",
"ao_ion = mol.eval_gto(\"GTOval_sph\", coords)\n",
"mo_ion = ao_ion.dot(scf_ion.mo_coeff)\n",
"\n",
"# FO2\n",
"ao_Z1 = mol_fl.eval_gto(\"GTOval_sph\", coords)\n",
"mo_Z1 = ao_Z1.dot(scf_gs_fl.mo_coeff)\n",
"\n",
"plt.figure(figsize=(10, 6))\n",
"plt.subplot(221)\n",
"plt.title(r\"MO$_{\\alpha 1}$ ($1a_1$)\")\n",
"mo_nr = 0\n",
"r, rho1, rho2, rho3 = [], [], [], []\n",
"for i in np.arange(len(coords)):\n",
" r.append(coords[i][0])\n",
" rho1.append(coords[i][0] ** 2 * mo_gs[i][0][mo_nr] ** 2)\n",
" rho2.append(coords[i][0] ** 2 * mo_ion[i][0][mo_nr] ** 2)\n",
" rho3.append(coords[i][0] ** 2 * mo_Z1[i][0][mo_nr] ** 2)\n",
"plt.plot(r, rho1)\n",
"plt.plot(r, rho2)\n",
"plt.plot(r, rho3)\n",
"plt.legend((\"Ground-state\", \"Core-hole\", \"Ground-state [Z+1]\"))\n",
"plt.xlim((0, 0.75))\n",
"plt.xlabel(\"Radius [a.u.]\")\n",
"plt.ylabel(r\"Density [$r^2 \\rho$]\")\n",
"\n",
"plt.subplot(222)\n",
"plt.title(r\"MO$_{\\alpha 2}$ ($2a_1$)\")\n",
"mo_nr = 1\n",
"r, rho1, rho2, rho3 = [], [], [], []\n",
"for i in np.arange(len(coords)):\n",
" r.append(coords[i][0])\n",
" rho1.append(coords[i][0] ** 2 * mo_gs[i][0][mo_nr] ** 2)\n",
" rho2.append(coords[i][0] ** 2 * mo_ion[i][0][mo_nr] ** 2)\n",
" rho3.append(coords[i][0] ** 2 * mo_Z1[i][0][mo_nr] ** 2)\n",
"plt.plot(r, rho1)\n",
"plt.plot(r, rho2)\n",
"plt.plot(r, rho3)\n",
"plt.xlim((0, 2.5))\n",
"plt.xlabel(\"Radius [a.u.]\")\n",
"plt.ylabel(r\"Density [$r^2 \\rho$]\")\n",
"\n",
"plt.subplot(224)\n",
"plt.title(r\"MO$_{\\beta 2}$ ($2a_1$)\")\n",
"mo_nr = 1\n",
"r, rho1, rho2, rho3 = [], [], [], []\n",
"for i in np.arange(len(coords)):\n",
" r.append(coords[i][0])\n",
" rho1.append(coords[i][0] ** 2 * mo_gs[i][1][mo_nr] ** 2)\n",
" rho2.append(coords[i][0] ** 2 * mo_ion[i][1][mo_nr] ** 2)\n",
" rho3.append(coords[i][0] ** 2 * mo_Z1[i][1][mo_nr] ** 2)\n",
"plt.plot(r, rho1)\n",
"plt.plot(r, rho2)\n",
"plt.plot(r, rho3)\n",
"plt.xlim((0, 2.5))\n",
"plt.xlabel(\"Radius [a.u.]\")\n",
"plt.ylabel(r\"Density [$r^2 \\rho$]\")\n",
"\n",
"plt.tight_layout()\n",
"plt.show()\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"import copy\n",
"import adcc\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"au2ev = 27.211386\n",
"\n",
"\n",
"def lorentzian(x, y, xmin, xmax, xstep, gamma):\n",
" \"\"\"\n",
" Lorentzian broadening function\n",
"\n",
" Call: xi,yi = lorentzian(energies, intensities, start energy,\n",
" end energy, energy step, gamma)\n",
" \"\"\"\n",
" xi = np.arange(xmin, xmax, xstep)\n",
" yi = np.zeros(len(xi))\n",
" for i in range(len(xi)):\n",
" for k in range(len(x)):\n",
" yi[i] = (\n",
" yi[i]\n",
" + y[k] * gamma / ((xi[i] - x[k]) ** 2 + (gamma / 2.0) ** 2) / np.pi\n",
" )\n",
" return xi, yi"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As can be seen, the radial distribution of particularly $\\alpha_2$ for the core-ionized system is closer to $Z+1$ than to the neutral state, while for $\\alpha_1$ and $\\beta_2$ it is closer to that of the ground state. Also note the differences in the x-scales, with $\\alpha_1$ being more compressed and thus having a larger peak max.\n",
"\n",
"\n",
"### Relaxation of core-(de)excitation\n",
"\n",
"For X-ray absorption and emission spectroscopy the influence of the relaxation is more complex, on account of the interaction with the created singly-occupied state. Comparing the two spectroscopies, the relaxation processes are partially inverse to each other, with\n",
"\n",
"\n",
"- Relaxation in XAS leading to a more contracted final state, which lowers absorption energies\n",
"\n",
"- Relaxation in XES leading to a more decontracted final state, which increases emission energies\n",
"\n",
"In addition to these strong changes in shielding, there are also additional couplings due to the population of a previous unoccupied state (XAS) or the creation of a valence-hole (XES). These effects are smaller and partially counteract the stronger change in shielding from the occupied/unoccupied core orbital, leading to a total picture as [follows](https://pubs.acs.org/doi/10.1021/acs.jctc.8b01046):\n",
"\n",
"```{figure} ../../img/xray/relax_processes.svg\n",
"---\n",
"name: relax_processes\n",
"---\n",
"Relaxation effects involved in X-ray absorption (left) and emission (right) spectroscopies. Showing which orbitals primarily contribute to the respective effects, as well as the result in MO energies (or electron densities).\n",
"```\n",
"\n",
"We here see the four main effects:\n",
"\n",
"- X-ray absorption spectroscopy\n",
"\n",
" 1. A strong relaxation effect due to the creation of a core-hole, which leads to the contraction of electron density and scales with element\n",
" \n",
" 2. A weak polarization effect due to coupling with the newly occupied valence state, which partially decontracts the outer valence density and does not scale appreciably with the element\n",
" \n",
"- X-ray emission spectroscopy\n",
"\n",
" 1. A strong repulsion effect due to the annihilation of a core-hole, which leads to decontraction of electron density and scales with element\n",
" \n",
" 2. A weaker outer relaxation effect due to the creation of a valence-hole, which leads to contraction of the (outer) valence density and scales with element\n",
"\n",
"These effects balances in different manners, and the exact strength depends heavily on the element and which edge is being probed. As such, the accuracy of a particular method depends on which spectroscopy, element, and edge is considered.\n",
"\n",
"\n",
"(sec:xray_delta_rel)=\n",
"### Including relaxation with $\\Delta$SCF methods\n",
"\n",
"$\\Delta$SCF methods provide the most direct way of including relaxation, as the core-hole state is explicitly optimized. However, these approaches need to consider each transition explicitly, and will have difficulties in calculating intensities and in properly converging higher excited states. As such, they are mainly useful for considering ionization energies and for providing a reference state for XES calculations. Comparing the performance of calculating IEs of the ten-electron series using Koopmans' theorem and $\\Delta$SCF:\n",
"\n",
"```{note}\n",
"pyscf version - to be updated.\n",
"```\n",
"\n",
"\n",
"```python\n",
"methane_xyz = \"\"\"\n",
"C 0.0000000000 0.0000000000 0.0000000000\n",
"H 0.6265918120 0.6265918120 0.6265918120\n",
"H -0.6265918120 -0.6265918120 0.6265918120\n",
"H -0.6265918120 0.6265918120 -0.6265918120\n",
"H 0.6265918120 -0.6265918120 -0.6265918120\n",
"\"\"\"\n",
"ammonia_xyz = \"\"\"\n",
"N 0.0000000000 0.0000000000 0.1175868000\n",
"H 0.9323800000 0.0000000000 -0.2743692000\n",
"H -0.4661900000 -0.8074647660 -0.2743692000\n",
"H -0.4661900000 0.8074647660 -0.2743692000\n",
"\"\"\"\n",
"water_xyz = \"\"\"\n",
"O 0.0000000000 0.0000000000 0.1178336003\n",
"H -0.7595754146 -0.0000000000 -0.4713344012\n",
"H 0.7595754146 0.0000000000 -0.4713344012\n",
"\"\"\"\n",
"hydrofluoride_xyz = \"\"\"\n",
"H 0.0000000000 0.0000000000 -0.8261856000\n",
"F 0.0000000000 0.0000000000 0.0917984000\n",
"\"\"\"\n",
"neon_xyz = \"\"\"\n",
"Ne 0.0000000000 0.0000000000 0.0000000000\n",
"\"\"\"\n",
"\n",
"\n",
"def Koopman_vs_delta(molecule, basis):\n",
" # Create pyscf mol object\n",
" mol = gto.Mole()\n",
" mol.atom = molecule\n",
" mol.basis = basis\n",
" mol.build()\n",
" # Perform unrestricted SCF calculation\n",
" scf_gs = scf.UHF(mol)\n",
" scf_gs.kernel()\n",
" # Copy molecular orbitals\n",
" mo0 = copy.deepcopy(scf_gs.mo_coeff)\n",
" occ0 = copy.deepcopy(scf_gs.mo_occ)\n",
" # Create 1s core-hole by setting alpha_0 population to zero\n",
" occ0[0][0] = 0.0\n",
" # Perform unrestricted SCF calculation with MOM constraint\n",
" scf_ion = scf.UHF(mol)\n",
" scf.addons.mom_occ(scf_ion, mo0, occ0)\n",
" scf_ion.kernel()\n",
" # Return MO energy and deltaSCF\n",
" return -au2ev * scf_gs.mo_energy[0][0], au2ev * (\n",
" scf_ion.energy_tot() - scf_gs.energy_tot()\n",
" )\n",
"\n",
"\n",
"# Basis set and container for errors\n",
"basis = \"unc-6-31G\"\n",
"err_koopman, err_delta = [], []\n",
"\n",
"expt, rel = 290.76, 0.11 # experimental IE and scalar-relativistic effects\n",
"E_Koopman, E_delta = Koopman_vs_delta(methane_xyz, basis)\n",
"err_koopman.append(E_Koopman + rel - expt)\n",
"err_delta.append(E_delta + rel - expt)\n",
"\n",
"expt, rel = 405.52, 0.21 # experimental IE and scalar-relativistic effects\n",
"E_Koopman, E_delta = Koopman_vs_delta(ammonia_xyz, basis)\n",
"err_koopman.append(E_Koopman + rel - expt)\n",
"err_delta.append(E_delta + rel - expt)\n",
"\n",
"expt, rel = 539.90, 0.37 # experimental IE and scalar-relativistic effects\n",
"E_Koopman, E_delta = Koopman_vs_delta(water_xyz, basis)\n",
"err_koopman.append(E_Koopman + rel - expt)\n",
"err_delta.append(E_delta + rel - expt)\n",
"\n",
"expt, rel = 694.10, 0.61 # experimental IE and scalar-relativistic effects\n",
"E_Koopman, E_delta = Koopman_vs_delta(hydrofluoride_xyz, basis)\n",
"err_koopman.append(E_Koopman + rel - expt)\n",
"err_delta.append(E_delta + rel - expt)\n",
"\n",
"expt, rel = 870.09, 0.94 # experimental IE and scalar-relativistic effects\n",
"E_Koopman, E_delta = Koopman_vs_delta(neon_xyz, basis)\n",
"err_koopman.append(E_Koopman + rel - expt)\n",
"err_delta.append(E_delta + rel - expt)\n",
"```\n",
"\n",
"```python\n",
"plt.figure(figsize=(7, 4))\n",
"# Plot error from MO energies and delta-SCF\n",
"plt.plot([1, 2, 3, 4, 5], err_koopman, \"bs\")\n",
"plt.plot([1, 2, 3, 4, 5], err_delta, \"ro\")\n",
"plt.plot([0.5, 5.5], [0, 0], \"k--\", zorder=0)\n",
"\n",
"plt.legend((\"Koopmans' theorem\", r\"$\\Delta$SCF\"))\n",
"plt.ylabel(\"Error wrt experiment [eV]\")\n",
"plt.xticks([1, 2, 3, 4, 5], (r\"CH$_4$\", r\"NH$_3$\", r\"H$_2$O\", \"HF\", \"Ne\"))\n",
"plt.xlim((0.5, 5.5))\n",
"plt.tight_layout()\n",
"plt.show()\n",
"```\n",
"\n",
"```{figure} ../../img/xray/relax_scf.svg\n",
"---\n",
"name: xray_relax_scf\n",
"scale: 6%\n",
"---\n",
"```\n",
"\n",
"We see that the error from orbital energies is large (14.4-22.3 eV) and clearly scales with atomic number, while the error of $\\Delta$SCF remains small (absolute error of 0.15-0.38 eV). \n",
"\n",
"\n",
"### Including relaxation with wave-function methods\n",
"\n",
"Considering the X-ray absorption and emission spectrum of water, the error of ADC(*n*) calculations with respect to experimental measurements can be considered:\n",
"\n",
"```{note}\n",
"pyscf version - to be updated.\n",
"```\n",
"\n",
"\n",
"```python\n",
"# Create pyscf mol object\n",
"mol = gto.Mole()\n",
"mol.atom = water_xyz\n",
"mol.basis = \"unc-6-31G\"\n",
"mol.build()\n",
"\n",
"# Perform restricted SCF calculation\n",
"scf_gs = scf.RHF(mol)\n",
"scf_gs.kernel()\n",
"\n",
"# Perform CVS-ADC calculations\n",
"xas_adc1 = adcc.cvs_adc1(scf_gs, n_singlets=2, core_orbitals=1)\n",
"xas_adc2 = adcc.cvs_adc2(scf_gs, n_singlets=2, core_orbitals=1)\n",
"xas_adc2x = adcc.cvs_adc2x(scf_gs, n_singlets=2, core_orbitals=1)\n",
"xas_adc3 = adcc.cvs_adc3(scf_gs, n_singlets=2, core_orbitals=1)\n",
"\n",
"# Perform unrestricted SCF calculation\n",
"scf_gs = scf.UHF(mol)\n",
"scf_gs.kernel()\n",
"\n",
"# Copy molecular orbitals\n",
"mo0 = copy.deepcopy(scf_gs.mo_coeff)\n",
"occ0 = copy.deepcopy(scf_gs.mo_occ)\n",
"\n",
"# Create 1s core-hole by setting alpha_0 population to zero\n",
"occ0[0][0] = 0.0\n",
"\n",
"# Perform unrestricted SCF calculation with MOM constraint\n",
"scf_ion = scf.UHF(mol)\n",
"scf.addons.mom_occ(scf_ion, mo0, occ0)\n",
"scf_ion.kernel()\n",
"\n",
"# Perform ADC calculations\n",
"xes_adc1 = adcc.adc1(scf_ion, n_states=2)\n",
"xes_adc2 = adcc.adc2(scf_ion, n_states=2)\n",
"xes_adc2x = adcc.adc2x(scf_ion, n_states=2)\n",
"xes_adc3 = adcc.adc3(scf_ion, n_states=2)\n",
"```\n",
"\n",
"Plotting the errors compared to experiment (given as a dashed line), and the variations in oscillator strength:\n",
"\n",
"```python\n",
"plt.figure(figsize=(10, 6))\n",
"plt.subplot(221)\n",
"plt.title(\"XAS energy\")\n",
"expt, rel = 534.00, 0.37\n",
"plt.plot(0, au2ev * xas_adc1.excitation_energy[0] + rel - expt, \"bo\")\n",
"plt.plot(1, au2ev * xas_adc2.excitation_energy[0] + rel - expt, \"bo\")\n",
"plt.plot(2, au2ev * xas_adc2x.excitation_energy[0] + rel - expt, \"bo\")\n",
"plt.plot(3, au2ev * xas_adc3.excitation_energy[0] + rel - expt, \"bo\")\n",
"plt.plot([-0.5, 3.5], [0, 0], \"k--\", zorder=0)\n",
"plt.xlim((-0.5, 3.5))\n",
"plt.ylabel(\"Error wrt experiment [eV]\")\n",
"plt.xticks([0, 1, 2, 3], (\"ADC(1)\", \"ADC(2)\", \"ADC(2)-x\", \"ADC(3/2)\"))\n",
"\n",
"plt.subplot(222)\n",
"plt.title(\"XAS oscillator strength\")\n",
"plt.plot(0, xas_adc1.oscillator_strength[0], \"bo\")\n",
"plt.plot(1, xas_adc2.oscillator_strength[0], \"bo\")\n",
"plt.plot(2, xas_adc2x.oscillator_strength[0], \"bo\")\n",
"plt.plot(3, xas_adc3.oscillator_strength[0], \"bo\")\n",
"plt.xlim((-0.5, 3.5))\n",
"plt.ylabel(\"Absolute oscillator strength\")\n",
"plt.xticks([0, 1, 2, 3], (\"ADC(1)\", \"ADC(2)\", \"ADC(2)-x\", \"ADC(3/2)\"))\n",
"\n",
"plt.subplot(223)\n",
"plt.title(\"XES energy\")\n",
"expt, rel = 527.10, 0.37\n",
"plt.plot(0, -au2ev * xes_adc1.excitation_energy[0] + rel - expt, \"bo\")\n",
"plt.plot(1, -au2ev * xes_adc2.excitation_energy[0] + rel - expt, \"bo\")\n",
"plt.plot(2, -au2ev * xes_adc2x.excitation_energy[0] + rel - expt, \"bo\")\n",
"plt.plot(3, -au2ev * xes_adc3.excitation_energy[0] + rel - expt, \"bo\")\n",
"plt.plot([-0.5, 3.5], [0, 0], \"k--\", zorder=0)\n",
"plt.xlim((-0.5, 3.5))\n",
"plt.ylabel(\"Error wrt experiment [eV]\")\n",
"plt.xticks([0, 1, 2, 3], (\"ADC(1)\", \"ADC(2)\", \"ADC(2)-x\", \"ADC(3/2)\"))\n",
"\n",
"plt.subplot(224)\n",
"plt.title(\"XES oscillator strength\")\n",
"plt.plot(0, xes_adc1.oscillator_strength[0], \"bo\")\n",
"plt.plot(1, xes_adc2.oscillator_strength[0], \"bo\")\n",
"plt.plot(2, xes_adc2x.oscillator_strength[0], \"bo\")\n",
"plt.plot(3, xes_adc3.oscillator_strength[0], \"bo\")\n",
"plt.xlim((-0.5, 3.5))\n",
"plt.ylabel(\"Absolute oscillator strength\")\n",
"plt.xticks([0, 1, 2, 3], (\"ADC(1)\", \"ADC(2)\", \"ADC(2)-x\", \"ADC(3/2)\"))\n",
"\n",
"plt.tight_layout()\n",
"plt.show()\n",
"```\n",
"\n",
"```{figure} ../../img/xray/xray_adc.svg\n",
"---\n",
"name: xray_adc\n",
"scale: 6%\n",
"---\n",
"```\n",
"\n",
"As can be seen, the smallest discrepancy for XAS is observed for CVS-ADC(2)-x, and for XES it is here for ADC(3/2). In general, it has been observed that CVS-ADC(2)-x performs best for [XAS](https://aip.scitation.org/doi/10.1063/1.4921841), and ADC(2) for [XES](https://pubs.acs.org/doi/10.1021/acs.jctc.8b01046) (in terms of absolute energies - for relative features ADC(2)-x performs a bit [better](https://doi.org/10.26434/chemrxiv-2022-vllwh-v2)). As ADC(3/2) is correct in higher order of perturbation theory and performs well for valence excitations, it could be expected to also perform best for core properties. The fact that this is not the case has [been](https://aip.scitation.org/doi/10.1063/1.4921841) [discussed](https://pubs.acs.org/doi/10.1021/acs.jctc.8b01046) as being due to balancing in relaxation effects, which are highly influential for core properties and vary significantly between XAS and XES.\n",
"\n",
"Furthermore, we see that the oscillator strengths vary greatly for XAS $-$ in particular when comparing CVS-ADC(1) to higher orders of theory $-$ but not very much for XES.\n",
"\n",
"\n",
"### Including relaxation with DFT-based methods\n",
"\n",
"When using a $\\Delta$KS approach, relaxation is explicitly accounted for, as discussed [above](sec:xray_delta_rel). By contrast, TDDFT is a single-excitation theory which cannot properly account for relaxation (which requires doubly-excited configuration in a response theory framework), although it can be partially accounted for through some separate correction scheme. The error introduced by a lack of relaxation will furthermore be counteracted by the self-interaction error, as is discussed [below](sec:xray_topics_sie).\n",
"\n",
"\n",
"\n",
"(sec:xray_topics_rel)=\n",
"## Relativity\n",
"\n",
"Relativistic effects are siginificantly stronger for core than valence MOs, on account of the strong potential experienced by the core orbitals. These effects mainly influences calculated properties in two separate manners:\n",
"\n",
"1. Scalar relativistic effects, which are always present and lead to a contraction (and thus increase in binding energy) of core orbitals. Higher lying orbitals (especially for $l>1$) can become decontracted due to increasing screening from contracted MOs, but this effect is typically much smaller than for, *e.g.*, 1s orbitals\n",
"\n",
"2. Spin-orbit coupling, which breaks degeneracies for MOs of $l>0$ and is thus negligible for, *e.g.*, the $K$- and $L_1$-edges (as the splitting for valence MOs is much smaller), but very important for, *e.g.*, the $L_{2,3}$-edge\n",
"\n",
"\n",
"### Scalar relativistic effects\n",
"\n",
"Including scalar relativistic effects is relatively straightforward, as a number of Hamiltonians are available which can include these effects in a 1-component framework. Furthermore, the effects are typically quite stable over the separate elements, and can thus often be considered as a common absolute off-set in energy. It should be noted that the (scalar) relativistic effects are technically not additive with electron correlation, such that a proper consideration should account for both simultaneously. In practise, this discrepancy is small and relativistc shifts from lower levels of theory are typically sufficient.\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"remove-cell"
]
},
"source": [
"\n",
"(sec:xray_topics_rel_ip)=\n",
"#### For ionization energies\n",
"\n",
"The scalar relativistic effects for ionization energies can be estimated from MO energies or using $\\Delta$SCF:\n",
"\n",
"```python\n",
"# Containers for atomic number, MO energy shift, delta-SCF shift, and delta-SCF IP\n",
"z, mo_diff, delta_diff, delta_ie = [], [], [], []\n",
"\n",
"# Consider even-numbered elements fron helium to argon\n",
"for i in np.arange(2, 20, 2):\n",
" # Atomic calculation\n",
" struct = str(i) + \" 0.0 0.0 0.0\"\n",
" # Create pyscf mol object\n",
" mol = gto.Mole()\n",
" mol.atom = struct\n",
" mol.basis = \"unc-6-31G\"\n",
" mol.build()\n",
" # Perform non-relativistic unrestricted SCF calculation\n",
" scf_nr = scf.UHF(mol)\n",
" scf_nr.kernel()\n",
" # Perform relativistic unrestricted SCF calculation\n",
" scf_rel = scf.UHF(\n",
" mol\n",
" ).x2c1e() # Scalar relativistic effects with the X2C Hamiltonian\n",
" scf_rel.kernel()\n",
" # Append atomic number\n",
" z.append(i)\n",
" # Append MO energy shift\n",
" mo_diff.append(au2ev * (scf_nr.mo_energy[0][0] - scf_rel.mo_energy[0][0]))\n",
" # Copy molecular orbitals and change population of 1s\n",
" mo0 = copy.deepcopy(scf_nr.mo_coeff)\n",
" occ0 = copy.deepcopy(scf_nr.mo_occ)\n",
" occ0[0][0] = 0.0\n",
" # Perform non-relativistic core-hole calculation\n",
" ion_nr = scf.UHF(mol)\n",
" scf.addons.mom_occ(ion_nr, mo0, occ0)\n",
" ion_nr.kernel()\n",
" # Perform relativistic core-hole calculation (can use the same gues wave function)\n",
" ion_rel = scf.UHF(mol).x2c1e()\n",
" scf.addons.mom_occ(ion_rel, mo0, occ0)\n",
" ion_rel.kernel()\n",
" # Append delta energy and relativistic shift\n",
" delta_diff.append(\n",
" au2ev * ((ion_rel.e_tot - scf_rel.e_tot) - (ion_nr.e_tot - scf_nr.e_tot))\n",
" )\n",
" delta_ie.append(au2ev * (ion_rel.e_tot - scf_rel.e_tot))\n",
"```\n",
"\n",
"\n",
"```python\n",
"plt.figure(figsize=(10, 4))\n",
"# Plot absolute scalar relativistic effects\n",
"plt.subplot(121)\n",
"plt.title(\"Absolute relativistic effects\")\n",
"plt.plot(z, mo_diff)\n",
"plt.plot(z, delta_diff)\n",
"plt.xticks(np.arange(2, 20, 2))\n",
"plt.legend((\"From MO energies\", r\"From $\\Delta$SCF\"))\n",
"plt.xlabel(\"Atomic number\")\n",
"plt.ylabel(\"Absolute shift [eV]\")\n",
"\n",
"# Plot percentage change\n",
"plt.subplot(122)\n",
"plt.title(\"Percentage change\")\n",
"plt.plot(z, 100 * np.array(mo_diff) / np.array(delta_ie))\n",
"plt.plot(z, 100 * np.array(delta_diff) / np.array(delta_ie))\n",
"plt.xticks(np.arange(2, 20, 2))\n",
"plt.xlabel(\"Atomic number\")\n",
"plt.ylabel(\"Relative shift [% of total energy]\")\n",
"plt.tight_layout()\n",
"plt.savefig(\"foo.svg\")\n",
"plt.show()\n",
"```\n",
"\n",
"```{figure} ../../img/xray/scalar_rel.svg\n",
"---\n",
"name: xray_scalar_rel\n",
"scale: 6%\n",
"---\n",
"```\n",
"\n",
"Two things are immediately clear:\n",
"\n",
"1. The shift in 1s ionization energy is strickly non-linear with respect to the atomic number, and using a simple model for a hydrogen-like system the main effects are expected to scale as $Z^4$. Fitting to a power function, we obtain $Z^{4.6}$ (with $R^2=0.997$) for the $\\Delta$SCF results, so some deviation from the hydrogenic systems\n",
"\n",
"2. The relativistic shifts for MO energies and $\\Delta$SCF IEs are almost identical, with a discrepancy that increases with $Z$ from 0.00 to 0.31 eV. Performing similar calculations with $\\Delta$CCSD is expected to yield very similar values, such that the scalar relativistic effects for the 1s MO of the second and third row can be estimated from MO energies\n",
"\n",
"\n",
"#### For spectra (XAS)\n",
"\n",
"For spectrum calculations, scalar relativistic effects can be included merely by using an appropriate Hamiltonian for the SCF calculation, as here illustrated for the X-ray absorption spectrum of water:\n",
"\n",
"```python\n",
"# Create pyscf mol object\n",
"mol = gto.Mole()\n",
"mol.atom = water_xyz\n",
"# using a decontracted basis set to better describe contracted relativistic 1s\n",
"mol.basis = \"unc-6-31G\"\n",
"mol.build()\n",
"\n",
"# Perform non-relativistic restricted SCF and CVS-ADC(2) calculations\n",
"scf_nr = scf.RHF(mol)\n",
"scf_nr.kernel()\n",
"adc_nr = adcc.cvs_adc2(scf_nr, n_singlets=5, core_orbitals=1)\n",
"\n",
"# Perform relativistic restricted SCF and CVS-ADC(2) calculations\n",
"scf_rel = scf.RHF(mol).x2c1e()\n",
"scf_rel.kernel()\n",
"adc_rel = adcc.cvs_adc2(scf_rel, n_singlets=5, core_orbitals=1)\n",
"\n",
"plt.figure(figsize=(6, 3))\n",
"# Relativistic and non-relativistc spectra with custom broadening\n",
"x, y = au2ev * adc_nr.excitation_energy, adc_nr.oscillator_strength\n",
"xi, yi = lorentzian(x, y, 534, 542, 0.01, 0.5)\n",
"plt.plot(xi, yi)\n",
"x, y = au2ev * adc_rel.excitation_energy, adc_rel.oscillator_strength\n",
"xi, yi = lorentzian(x, y, 534, 542, 0.01, 0.5)\n",
"plt.plot(xi, yi)\n",
"plt.legend((\"Non-relativistic\", \"Scalar relativistic\"))\n",
"plt.tight_layout()\n",
"plt.show()\n",
"```\n",
"\n",
"We see that the spectrum is shifted upwards by $\\sim$0.38 eV, due to the increase in core-electron binding energy. The same approach can be used for including these effects also for X-ray emission spectrum calculations."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Spin-orbit coupling\n",
"\n",
"Including the spin-orbit effects is more complicated, and can either be\n",
"\n",
"- Included by using 4-component relativistic theory, where they do not need to explicitly addressed\n",
"\n",
"- Included in a non-relativistic approach through a perturbative scheme, or through some *ad hoc* correction \n",
"\n",
"At the present stage, including spin-orbit effects is beyond the scope of this tutorial.\n",
"\n",
"\n",
"(sec:xray_topics_basis)=\n",
"## Basis set considerations\n",
"\n",
"The basis set requirements of X-ray spectrum calculations depend on which spectroscopy is considered\n",
"\n",
"- XPS and XES probe occupied states, and thus need a good description of those\n",
"\n",
"- XAS and RIXS probe excited states, and thus need an improved description of this region\n",
"\n",
"In all cases the relaxation due to the creation of a core-hole needs to be accommodated, which yields requirements of reasonable flexibility of the core and inner valence region. This is not usually the case for standard basis sets, as they are typically constructed with valence properties in mind and thus have a minimal or close to minimal description of the core region. A number of approaches for improving such valence-focused basis sets have been developed by:\n",
"\n",
"1. Augmenting the Dunning basis sets with core-polarizing functions, *e.g.* cc-pV*n*Z $\\rightarrow$ cc-pCV*n*Z\n",
"\n",
"2. Adding extra flexibility by performing full or partial decontraction, such as\n",
" 1. Decontracting 1s of the probed element (*e.g.* u6-311++G\\*\\*, which is based on 6-311++G\\*\\*)\n",
" 2. Fully decontracting the basis set (*e.g.* un6-311++G\\*\\*, which is based on 6-311++G\\*\\*)\n",
"\n",
"3. Using basis functions from the next element, as inspired by the $Z+1$ approximation\n",
"\n",
"As an illustration, we consider the IE of water with different basis sets and levels of theory\n",
"\n",
"*Toggle for u6-311G\\*\\*, un6-311G\\*\\*, and 6-311G\\*\\* (Z+1) basis sets.*\n",
"\n",
"```{toggle}\n",
"```python\n",
"from pyscf import gto\n",
"\n",
"u6311gss = {'O': gto.basis.parse('''\n",
"O S\n",
" 8588.500 1.000000 \n",
"O S\n",
" 1297.230 1.000000 \n",
"O S\n",
" 299.2960 1.000000 \n",
"O S\n",
" 87.37710 1.000000 \n",
"O S\n",
" 25.67890 1.000000 \n",
"O S\n",
" 3.740040 1.000000 \n",
"O SP\n",
" 42.11750 0.113889 0.0365114\n",
" 9.628370 0.920811 0.237153\n",
" 2.853320 -0.00327447 0.819702\n",
"O SP\n",
" 0.905661 1.000000 1.000000\n",
"O SP\n",
" 0.255611 1.000000 1.000000\n",
"O D\n",
" 1.292 1.000000\n",
"'''), 'H': '6-311G**'}\n",
"un6311gss = {'O': gto.basis.parse('''\n",
"O S\n",
" 8588.500 1.000000 \n",
"O S\n",
" 1297.230 1.000000 \n",
"O S\n",
" 299.2960 1.000000 \n",
"O S\n",
" 87.37710 1.000000 \n",
"O S\n",
" 25.67890 1.000000 \n",
"O S\n",
" 3.740040 1.000000 \n",
"O SP\n",
" 42.11750 1.000000 1.000000\n",
"O SP\n",
" 9.628370 1.000000 1.000000\n",
"O SP\n",
" 2.853320 1.000000 1.000000\n",
"O SP\n",
" 0.905661 1.000000 1.000000\n",
"O SP\n",
" 0.255611 1.000000 1.000000\n",
"O D\n",
" 1.292 1.000000\n",
"'''), 'H': '6-311G**'}\n",
"z6311gss = {'O': gto.basis.parse('''\n",
"O S\n",
" 8588.500 0.00189515\n",
" 1297.230 0.0143859\n",
" 299.2960 0.0707320\n",
" 87.37710 0.2400010\n",
" 25.67890 0.5947970\n",
" 3.740040 0.2808020\n",
"O SP\n",
" 42.11750 0.113889 0.0365114\n",
" 9.628370 0.920811 0.237153\n",
" 2.853320 -0.00327447 0.819702\n",
"O SP\n",
" 0.905661 1.000000 1.000000\n",
"O SP\n",
" 0.255611 1.000000 1.000000\n",
"O D\n",
" 1.292 1.000000\n",
"O S\n",
" 11427.10 0.00180093\n",
" 1722.350 0.0137419\n",
" 395.7460 0.0681334\n",
" 115.1390 0.2333250\n",
" 33.60260 0.5890860\n",
" 4.919010 0.2995050\n",
"O SP\n",
" 55.44410 0.114536 0.0354609\n",
" 12.63230 0.920512 0.237451\n",
" 3.717560 -0.00337804 0.820458\n",
"O SP\n",
" 1.165450 1.000000 1.000000\n",
"O SP\n",
" 0.321892 1.000000 1.000000\n",
"O D\n",
" 1.750 1.000000\n",
"'''), 'H': '6-311G**'}"
]
},
{
"cell_type": "markdown",
"metadata": {
"tags": [
"hide-input"
]
},
"source": [
"```{note}\n",
"pyscf version - to be changed\n",
"```\n",
"\n",
"\n",
"```python\n",
"# Investigated basis sets\n",
"basis_sets = [\n",
" \"6-311G**\",\n",
" u6311gss,\n",
" un6311gss,\n",
" z6311gss,\n",
" \"cc-pVDZ\",\n",
" \"cc-pCVDZ\",\n",
" \"cc-pVTZ\",\n",
" \"cc-pCVTZ\",\n",
" \"cc-pVQZ\",\n",
" \"cc-pCVQZ\",\n",
"]\n",
"\n",
"\n",
"# Containers for IE calculated with four different delta calculations: HF, MP2, CCSD, and B3LYP\n",
"ip_hf, ip_mp2, ip_ccsd, ip_b3lyp, n_bas = [], [], [], [], []\n",
"\n",
"for basis in basis_sets:\n",
" # Create pyscf mol object\n",
" mol = gto.Mole()\n",
" mol.atom = water_xyz\n",
" mol.basis = basis\n",
" mol.build()\n",
" # Unrestricted SCF calculation of neutral\n",
" scfres = scf.UHF(mol)\n",
" scfres.kernel()\n",
" # Copy molecular orbitals and change population of 1s\n",
" mo0 = copy.deepcopy(scfres.mo_coeff)\n",
" occ0 = copy.deepcopy(scfres.mo_occ)\n",
" occ0[0][0] = 0.0\n",
" # Unrestricted HF calculation of core-hole and appending delta_HF IP\n",
" scfion = scf.UHF(mol)\n",
" scf.addons.mom_occ(scfion, mo0, occ0)\n",
" scfion.kernel()\n",
" ip_hf.append(au2ev * (scfion.energy_tot() - scfres.energy_tot()))\n",
" # CCSD calculations of neutral and core-hole and appending delta_CCSD IP\n",
" ccsd_res = cc.UCCSD(scfres).run()\n",
" ccsd_ion = cc.UCCSD(scfion).run(max_cycle=500)\n",
" ip_ccsd.append(au2ev * (ccsd_ion.e_tot - ccsd_res.e_tot))\n",
" # Appending delta_MP2 IP\n",
" ip_mp2.append(\n",
" au2ev * (ccsd_ion.e_hf + ccsd_ion.emp2 - (ccsd_res.e_hf + ccsd_res.emp2))\n",
" )\n",
" # Unrestricted KS-DFT calculation of neutral\n",
" scfres = scf.UKS(mol)\n",
" scfres.xc = \"b3lyp\"\n",
" scfres.kernel()\n",
" # Copy molecular orbitals and change population of 1s\n",
" mo0 = copy.deepcopy(scfres.mo_coeff)\n",
" occ0 = copy.deepcopy(scfres.mo_occ)\n",
" occ0[0][0] = 0.0\n",
" # Unrestricted KS-DFT calculation of core-hole and appending delta_DFT IP\n",
" scfion = scf.UKS(mol)\n",
" scf.addons.mom_occ(scfion, mo0, occ0)\n",
" scfion.xc = \"b3lyp\"\n",
" scfion.kernel()\n",
" ip_b3lyp.append(au2ev * (scfion.energy_tot() - scfres.energy_tot()))\n",
" # Append basis set size\n",
" n_bas.append(len(scfres.mo_occ[0]))\n",
"\n",
"plt.figure(figsize=(10, 5))\n",
"# Plot IEs, with error as a function of basis set size\n",
"plt.plot(n_bas, ip_hf, \"r*\")\n",
"plt.plot(n_bas, ip_mp2, \"bv\")\n",
"plt.plot(n_bas, ip_ccsd, \"d\", color=\"orange\")\n",
"plt.plot(n_bas, ip_b3lyp, \"go\")\n",
"plt.legend((\"HF\", \"MP2\", \"CCSD\", \"B3LYP\", \"Expt\", \"Expt - 0.37 eV\"))\n",
"\n",
"# Experimental value as a dotted line, minus relativistic effect as dashed line\n",
"plt.plot([min(n_bas) - 5, max(n_bas) + 5], [539.9, 539.9], \"k:\")\n",
"plt.plot([min(n_bas) - 5, max(n_bas) + 5], [539.9 - 0.37, 539.9 - 0.37], \"k--\")\n",
"\n",
"# Basis set label as x-ticks\n",
"basis_set_labels = [\n",
" \"6-311G**\",\n",
" \"u6-311G**\",\n",
" \"un6-311G**\",\n",
" \"6-311G** (Z+1)\",\n",
" \"cc-pVDZ\",\n",
" \"cc-pCVDZ\",\n",
" \"cc-pVTZ\",\n",
" \"cc-pCVTZ\",\n",
" \"cc-pVQZ\",\n",
" \"cc-pCVQZ\",\n",
"]\n",
"plt.xticks(n_bas, basis_set_labels, rotation=70, fontsize=10)\n",
"\n",
"# Labels and setting x-interval\n",
"plt.ylabel(\"Ionization energy [eV]\")\n",
"plt.xlim((min(n_bas) - 5, max(n_bas) + 5))\n",
"plt.tight_layout()\n",
"plt.show()\n",
"```\n",
"\n",
"\n",
"\n",
"```{figure} ../../img/xray/basis_xps.svg\n",
"---\n",
"name: basis_xps\n",
"---\n",
"```\n",
"\n",
"At the HF and KS-DFT levels of theory we see that some of the smaller basis sets have fortuitous error cancellations, such that both 6-311G\\*\\* and cc-pVTZ yield results close to experimental values. However, augmenting these basis sets with more functions in the core region brings us away from experimental measurements (note that the calculations do not include relativistic effects, so the comparison would be toward experiment minus relativistic shift, *i.e.* the dashed line).\n",
"\n",
"With MP2 and CCSD we note very similar results for the smaller basis sets, which then deviate significantly for cc-pVTZ and above. $\\Delta$CCSD estimates using core-polarized cc-pVTZ and cc-pVQZ are well in line with experimental measurements. Augmenting the 6-311G\\*\\* basis sets with a decontracted core region (u6-311G\\*\\*) yields some improvement at a relatively low increase in computational cost, and further decontracting the full set (un6-311G\\*\\*) or adding basis functions from $Z+1$ (6-311G\\*\\* (Z+1)) provides some additional improvement.\n",
"\n",
"\n",
"(sec:xray_topics_sie)=\n",
"## Self-interaction error in DFT\n",
"\n",
"For a single-electron systems the following equality should hold for the two-electron terms:\n",
"\n",
"$$\n",
"J[\\rho] + E_{xc} [\\rho] = 0\n",
"$$\n",
"\n",
"This cancellation is preserved within HF theory, but for any approximate functionals in KS-DFT the equality will no longer hold, resulting in the potentially significant errors due to spurious self-interactions. Reformulated, the Coulomb self-repulsion of each electron included in the Coulomb operator is exactly canceled by the nonlocal exchange self-attraction in Hartree-Fock, but this is no longer the case when the exchange operator is replaced by an approximate exchange functional.\n",
"\n",
"As such, KS-DFT possess a self-interaction error (SIE), which shows a strong dependence on how compact the respective MOs are. This means that the self-interaction error is significantly more influential for core properties than for valence properties, and the tendency to decontract core orbitals due to the SIE has an opposite effect on transition energies as compared to relaxation effects.\n",
"\n",
"Let us consider the ionization energy of water, as calculated with MO energies and $\\Delta\\textrm{SCF}$ using two different exchange-correlation functionals with variable HF exchange: \n",
"\n",
"```{note}\n",
"pyscf version - to be changed (need more xc-functionals in VeloxChem)\n",
"```\n",
"\n",
"```{note}\n",
"Using MO energies for TDHF and TDDFT yields the same IE as when addressing core-transitions to extremely diffuse basis functions (due to the lack of relaxation), and is thus representative of the error for such XAS calculations.\n",
"```\n",
"\n",
"```python\n",
"# Create pyscf mol object\n",
"mol = gto.Mole()\n",
"mol.atom = water_xyz\n",
"mol.basis = \"6-31G\"\n",
"mol.build()\n",
"\n",
"# Investigated HF exchange and result containers\n",
"hf_exc = [0.0, 0.25, 0.50, 0.75, 1.0]\n",
"hf_blyp, hf_pbe = [[], []], [[], []]\n",
"\n",
"\n",
"def dft_delta(mol, xc):\n",
" # Neutral calculation\n",
" dft_res = scf.UKS(mol)\n",
" dft_res.xc = xc\n",
" dft_res.kernel()\n",
" # Core-hole calculation\n",
" mo0 = copy.deepcopy(dft_res.mo_coeff)\n",
" occ0 = copy.deepcopy(dft_res.mo_occ)\n",
" occ0[0][0] = 0.0\n",
" dft_ion = scf.UKS(mol)\n",
" dft_ion.xc = xc\n",
" scf.addons.mom_occ(dft_ion, mo0, occ0)\n",
" dft_ion.kernel()\n",
" return -au2ev * dft_res.mo_energy[0][0], au2ev * (\n",
" dft_ion.energy_tot() - dft_res.energy_tot()\n",
" )\n",
"\n",
"\n",
"for h_x in hf_exc:\n",
" # BxLYP (Slater exchange 0.08, except when HF i 1.00)\n",
" if h_x == 1.0:\n",
" s_x = 0.00\n",
" else:\n",
" s_x = 0.08\n",
" b_x = 1.00 - h_x - s_x\n",
" v_e, v_l = 0.19, 0.81\n",
" xc = f\"{h_x:} * HF + {s_x:} * Slater + {b_x:} * B88, {v_l:} * LYP + {v_e:} * VWN\"\n",
" E_mo, E_delta = dft_delta(mol, xc)\n",
" hf_blyp[0].append(E_mo)\n",
" hf_blyp[1].append(E_delta)\n",
" # PBEx\n",
" p_x = 1.0 - h_x\n",
" xc = f\"{h_x:} * HF + {p_x:} * PBE, PBE\"\n",
" E_mo, E_delta = dft_delta(mol, xc)\n",
" hf_pbe[0].append(E_mo)\n",
" hf_pbe[1].append(E_delta)\n",
"```\n",
"\n",
"```python\n",
"plt.figure(figsize=(6, 4))\n",
"# DFT ionization energies, including scalar relativistic shift\n",
"plt.plot(hf_exc, np.array(hf_blyp[0]) + 0.37, \"bs--\")\n",
"plt.plot(hf_exc, np.array(hf_blyp[1]) + 0.37, \"ro--\")\n",
"plt.plot(hf_exc, np.array(hf_pbe[0]) + 0.37, \"d--\", color=\"orange\")\n",
"plt.plot(hf_exc, np.array(hf_pbe[1]) + 0.37, \"g^--\")\n",
"\n",
"# Experimental value and legend\n",
"plt.plot([-0.1, 1.1], [539.9, 539.9], \"k-\")\n",
"plt.legend(\n",
" (r\"BxLYP ($\\Delta$SCF)\", \"BxLYP (MO)\", r\"PBEx ($\\Delta$SCF)\", \"PBEx (MO)\", \"Expt\")\n",
")\n",
"\n",
"plt.xlim((-0.025, 1.025))\n",
"plt.ylim((508, 562.0))\n",
"plt.xlabel(\"HF exchange [%]\")\n",
"plt.ylabel(\"Ionization energy [eV]\")\n",
"plt.tight_layout()\n",
"plt.show()\n",
"```\n",
"\n",
"\n",
"```{figure} ../../img/xray/sie.svg\n",
"---\n",
"name: xray_sie\n",
"scale: 6%\n",
"---\n",
"```\n",
"\n",
"We see that the $\\Delta$SCF errors are quite constant between the two base functionals and amounts of HF exchange, while the MO energies strongly scale with HF exchange, but not with the base functional. The total MO error is a function of two effects:\n",
"\n",
"1. Lack of relaxation, which generally leads to an overestimation of IEs and core-excitation energies\n",
"\n",
"2. Self-interaction, which by contrast leads to an underestimation of these energies\n",
"\n",
"For a typical exchange-correlation functional ($\\sim$20-25% HF exchange) the SIE thus dominates, and by tuning the (global) HF exchange to around 60% we achieve IE estimates in good agreement to experiment. We note that this is only acheieved by cancellation of the above two errors, and does not mean that we have an appropriate description of the core-excitation processes. Furthermore, a good exchange-correlation functional should possess good energies *and* densities, and focusing on only energies can push us away from such a goal. Nevertheless, a number of tailored functionals with good absolute energies have been constructed, tuning the HF exchange either in a global or range-separated manner. Note that the parameterization will have to depend on the spectroscopy, probed edge, and row of the periodic table, as the exact balance in relaxation terms and SIE varies.\n",
"\n",
"\n",
"\n",
"(sec:xray_topics_chloc)=\n",
"## Core-hole localization\n",
"\n",
"For $\\Delta$ SCF calculations (either for estimating IEs or as a initial state for XES), a core-hole is created at some atomic site. For many systems, this is unproblematic, but if chemically equivalent atoms are present, the MOs are delocalized over these sites. Performing a straightforward $\\Delta$SCF calculation will then results in a delocalized core-hole, as illustrated below for ethene:\n",
"\n",
"```{note}\n",
"pyscf version - to be changed\n",
"```\n",
"\n",
"```python\n",
"ethene_xyz = \"\"\"\n",
"C0 0.6660048050 0.0000000000 -0.0000000000\n",
"C1 -0.6660048050 -0.0000000000 0.0000000000\n",
"H 1.2278735620 0.9228358077 -0.0000000000\n",
"H 1.2278735620 -0.9228358077 -0.0000000000\n",
"H -1.2278735620 0.9228358077 0.0000000000\n",
"H -1.2278735620 -0.9228358077 0.0000000000\n",
"\"\"\"\n",
"\n",
"\n",
"def Koopman_vs_delta(molecule, basis, ecp=None):\n",
" # Create pyscf mol object\n",
" mol = gto.Mole()\n",
" mol.atom = molecule\n",
" mol.basis = basis\n",
" if ecp:\n",
" mol.ecp = ecp\n",
" mol.build()\n",
" # Perform unrestricted SCF calculation\n",
" scf_gs = scf.UHF(mol)\n",
" scf_gs.kernel()\n",
" # Copy molecular orbitals\n",
" mo0 = copy.deepcopy(scf_gs.mo_coeff)\n",
" occ0 = copy.deepcopy(scf_gs.mo_occ)\n",
" # Create 1s core-hole by setting alpha_0 population to zero\n",
" occ0[0][0] = 0.0\n",
" # Perform unrestricted SCF calculation with MOM constraint\n",
" scf_ion = scf.UHF(mol)\n",
" scf.addons.mom_occ(scf_ion, mo0, occ0)\n",
" scf_ion.kernel()\n",
" # Return MO energy and deltaSCF\n",
" return -au2ev * scf_gs.mo_energy[0][0], au2ev * (\n",
" scf_ion.energy_tot() - scf_gs.energy_tot()\n",
" )\n",
"\n",
"\n",
"# All-electron calculation of IP\n",
"delta_all = Koopman_vs_delta(ethene_xyz, \"unc-6-31G\")\n",
"\n",
"# IE calculation with ECP on one carbon\n",
"delta_ecp = Koopman_vs_delta(ethene_xyz, \"unc-6-31G\", {\"C0\": \"stuttgart\"})\n",
"```\n",
"\n",
"```python\n",
"plt.figure(figsize=(5, 3))\n",
"# IE estimates with all-electron (delocalized core-hole)\n",
"plt.plot([1, 2], [delta_all[0] + 0.11, delta_all[1] + 0.11], \"ro\")\n",
"\n",
"# IE estimates with ECP (localized core-hole)\n",
"plt.plot([1, 2], [delta_ecp[0] + 0.11, delta_ecp[1] + 0.11], \"bs\")\n",
"\n",
"# Experiment and figure settings\n",
"plt.plot([0.5, 2.5], [290.80, 290.80], \"k--\")\n",
"plt.ylabel(\"Ionization energy [eV]\")\n",
"plt.xticks([1, 2], (\"MO\", r\"$\\Delta$SCF\"))\n",
"plt.xlim((0.5, 2.5))\n",
"plt.tight_layout()\n",
"plt.show()\n",
"```\n",
"\n",
"\n",
"```{figure} ../../img/xray/ch_loc.svg\n",
"---\n",
"name: xray_ch_loc\n",
"scale: 6%\n",
"---\n",
"```\n",
"\n",
"As can be seen, the errors from MO energies are very similar for the all-electron calculation (with a delocalized core-hole) and using ECP (delocalized core-hole). For the $\\Delta$SCF calculation, however, only approximately half the relaxation energy is accounted for using a delocalized core-hole. Such a discrepancy for delocalized calculations of IEs has been noted [previously](https://doi.org/10.1063/1.1676850), and it has been [associated](https://doi.org/10.1063/1.433763) with missing correlation effects which are required to properly correct the delocalized calculation. Within a full CI scheme, delocalization will no longer be an issue (and it will capture the small shift in energy between the gerade and ungerade core-hole).\n",
"\n",
"Moving from theory, in experiment the localizatio/delocalization of core-holes is primarily relevant for diatomic systems, as other compounds will possess vibrational effects which break the symmetry. For diatomic systems the localized/delocalized nature of the core-hole has been seen to depend on the details of the measurement, such that both can be probed [ref].\n",
"\n",
"As such, for practical calculation we [recommend](https://doi.org/10.1063/5.0088195) that core-holes are localized for IE and XES calculations. For XAS it does not provide any issues for linear response methods (provided that both gerade and ungerade MOs are kept in the CVS space/equivalent), but for state-specific approaches this will again be an issue (*e.g.* STEX and TP-DFT).\n",
"\n",
"\n",
"\n",
"(sec:xray_topics_tailcvs)=\n",
"## Tailored CVS\n",
"\n",
"\n",
"For large systems, it may become relevant to tailor the CVS space such that each chemically unique atom is probed in turn. While this requires explicit consideration of each set of unique MOs separately, it serves dual purposes which can make it worth it:\n",
"\n",
"1. It lowers the dimensionality of the involved matrices, and thus decrease computational costs\n",
"\n",
"2. The density-of-states decrease, which means that the same energy region can be resolved with a few states per unique MO, rather than having to converge a large number of eigenstates simultaneously\n",
"\n",
"\n",
"While this technically does not fulfill one of foundations on which the CVS approximation is built (the MOs are not energetically well-separated), the spatial separation means that it will typically work. As an example, we consider ethene (which has delocalized MOs) and vinylfluoride (which does not):\n",
"\n",
"```{note}\n",
"pyscf version - to be changed\n",
"```\n",
"\n",
"```python\n",
"vinyl_xyz = \"\"\"\n",
" C 0.000000 -0.246412 -1.271068\n",
" C 0.000000 0.457081 -0.154735\n",
" F 0.000000 -0.119195 1.052878\n",
" H 0.000000 0.272328 -2.210194\n",
" H 0.000000 -1.319906 -1.249847\n",
" H 0.000000 1.530323 -0.095954\n",
"\"\"\"\n",
"\n",
"print(\"Ethene\")\n",
"# Create pyscf mol object\n",
"mol = gto.Mole()\n",
"mol.atom = ethene_xyz\n",
"mol.basis = \"6-31G\"\n",
"mol.build()\n",
"# Perform restricted SCF calculation\n",
"scf_gs = scf.RHF(mol)\n",
"scf_gs.kernel()\n",
"# Perform CVS-ADC calculations\n",
"ethene_adc1_0 = adcc.cvs_adc1(scf_gs, n_singlets=4, core_orbitals=[0])\n",
"ethene_adc1_1 = adcc.cvs_adc1(scf_gs, n_singlets=4, core_orbitals=[1])\n",
"ethene_adc1_01 = adcc.cvs_adc1(scf_gs, n_singlets=4, core_orbitals=[0, 1])\n",
"ethene_adc2_0 = adcc.cvs_adc2(scf_gs, n_singlets=4, core_orbitals=[0])\n",
"ethene_adc2_1 = adcc.cvs_adc2(scf_gs, n_singlets=4, core_orbitals=[1])\n",
"ethene_adc2_01 = adcc.cvs_adc2(scf_gs, n_singlets=4, core_orbitals=[0, 1])\n",
"\n",
"print(\"Vinylfluoride\")\n",
"# Create pyscf mol object\n",
"mol = gto.Mole()\n",
"mol.atom = vinyl_xyz\n",
"mol.basis = \"6-31G\"\n",
"mol.build()\n",
"# Perform restricted SCF calculation\n",
"scf_gs = scf.RHF(mol)\n",
"scf_gs.kernel()\n",
"# Perform CVS-ADC calculations\n",
"vinyl_adc1_1 = adcc.cvs_adc1(scf_gs, n_singlets=4, core_orbitals=[1])\n",
"vinyl_adc1_2 = adcc.cvs_adc1(scf_gs, n_singlets=4, core_orbitals=[2])\n",
"vinyl_adc1_12 = adcc.cvs_adc1(scf_gs, n_singlets=4, core_orbitals=[0, 1, 2])\n",
"vinyl_adc2_1 = adcc.cvs_adc2(scf_gs, n_singlets=4, core_orbitals=[1])\n",
"vinyl_adc2_2 = adcc.cvs_adc2(scf_gs, n_singlets=4, core_orbitals=[2])\n",
"vinyl_adc2_12 = adcc.cvs_adc2(scf_gs, n_singlets=4, core_orbitals=[0, 1, 2])\n",
"```\n",
"\n",
"\n",
"```python\n",
"plt.figure(figsize=(10, 6))\n",
"# Ethene with CVS-ADC(1)\n",
"plt.subplot(221)\n",
"plt.title(\"Ethene\")\n",
"ethene_adc1_0.plot_spectrum(label=\"MO #0\")\n",
"ethene_adc1_1.plot_spectrum(label=\"MO #1\")\n",
"ethene_adc1_01.plot_spectrum(label=\"MOs #0-1\")\n",
"plt.legend()\n",
"plt.ylabel(\"CVS-ADC(1)\")\n",
"\n",
"# Vinylfluoride with CVS-ADC(1)\n",
"plt.subplot(222)\n",
"plt.title(\"Vinylfluoride\")\n",
"vinyl_adc1_1.plot_spectrum(label=\"MO #1\")\n",
"vinyl_adc1_2.plot_spectrum(label=\"MO #2\")\n",
"vinyl_adc1_12.plot_spectrum(label=\"MO #1-2\")\n",
"plt.legend()\n",
"\n",
"# Ethene with CVS-ADC(2)\n",
"plt.subplot(223)\n",
"ethene_adc2_0.plot_spectrum()\n",
"ethene_adc2_1.plot_spectrum()\n",
"ethene_adc2_01.plot_spectrum()\n",
"plt.ylabel(\"CVS-ADC(2)\")\n",
"\n",
"# Vinylfluoride with CVS-ADC(2)\n",
"plt.subplot(224)\n",
"vinyl_adc2_1.plot_spectrum()\n",
"vinyl_adc2_2.plot_spectrum()\n",
"vinyl_adc2_12.plot_spectrum()\n",
"\n",
"plt.tight_layout()\n",
"plt.show()\n",
"```\n",
"\n",
"\n",
"\n",
"```{figure} ../../img/xray/adc_cvs.svg\n",
"---\n",
"name: xray_adc_cvs\n",
"scale: 6%\n",
"---\n",
"```\n",
"\n",
"We see that considering each MO in turn works well for vinylfluoride, where the two carbon 1s are localized, but not for ethene, where they are delocalized. Including or not including fluorine here makes little difference, as that MO is both energetically and spatially well separated from the carbon 1s.\n",
"\n",
"Note that the delocalization [might](https://doi.org/10.1063/5.0088195) not be immediatelly apparent, as just visually looking at the molecule and concluding which atoms are unique might not be sufficient. Care should be taken for larger systems, and explicitly looking at the MOs is advised. Alternatively, using ECPs to localize the core-hole and then restricting the CVS space accordingly will avoid this issue.\n",
" \n",
"\n",
"(sec:xray_topics_nonaufbau)=\n",
"## Post-HF on a non-Aufbau reference\n",
"\n",
"\n",
"```{note}\n",
"- Include a discussion on the potential issues in using post-HF methods for a non-*Aufbau* reference state, as can be seen by noting the MP2 energy denominator (and potential singularities, and discusses [here](https://doi.org/10.1039/D2CP00584K))\n",
"- Include explicit example:\n",
" 1. Go over all permutations (aa, ab, bb), saving the highest MO index which yields denominators below 0.1 au. Just have a list of virtual indices, adding when an index is not already found there. Do this in one shot, instead of saving one virtual at a time.\n",
" 2. Compare this to the scheme in adcc, noting that the present (simpler) scheme is likely just as useful - despite not ordering indices by how problematic they are.\n",
" 3. Use this to calculate the ADC(2)/cc-pCVQZ spectra of ammonia, with and without frozen virtuals.\n",
"```\n",
"\n",
"\n",
"\n",
"(sec:xray_cons_spur_val)=\n",
"## Spurious valence mixing in XAS\n",
"\n",
"Within a damped response theory framework or when performing explicit real-time propagation, a CVS-like scheme which separates core-transitions from the valence-continuum is usually not performed. If such calculations are run with an atom-centered basis set (as most are), this can then lead to spurious mixing with valence-excited states. What happens is that resonances in the valence-continuum occur close to core-excitations, leading to unphysical mixing which affects excitation energies and (more significantly) intensities. These issues occur both for methods featuring only single-excited configurations (*e.g.* TDDFT) and higher-excited configurations (*e.g.* ADC) and need to be avoided.\n",
"\n",
"These issues have been observed for, *e.g.*, $L_{2,3}$-edge of silicon {cite}`xas4crttddft, ledge2016`, and are usually more common at lower energies. They can be identified by studying the amplitudes and the $t^2$ values for post-HF methods, and can be avoided by chosing appropriate basis sets, freezing certain virtual orbitals, or by re-introducing the CVS approximation (or some other flavour) within the CPP/RT framework. The latter approach is recommended at least for CPP-ADC, as converging the response vectors has otherwise been seen to be difficult {cite}`gator`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"celltoolbar": "Tags",
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.15"
}
},
"nbformat": 4,
"nbformat_minor": 4
}