# Non-standard calculations#

In this chapter we will consider some non-standard calculations. Some can be useful but need care for practical considerations (i.e. transient X-ray spectroscopy), and some are more illustrations of alternative (but not recommended) ways of modelling spectra. Be careful before directly adopting these approaches for some practical problem.

import copy
import matplotlib.pyplot as plt
import numpy as np
from pyscf import gto, mp, scf
import veloxchem as vlx

def lorentzian(x, y, xmin, xmax, xstep, 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 / 2.0) / (
(xi[i] - x[k]) ** 2 + (gamma / 2.0) ** 2
)
return xi, yi

au2ev = 27.211386


## Transient X-ray absorption spectroscopy#

An increasingly large field of study is to use a pump-probe protocol, where some chemical reaction is initiated by some pump (typically a laser), and the dynamics are then probed using some X-ray spectroscopy. The modeling of such process is quite new, and we will here present some protocols of how such processes can be considered.

### Valence-excited reference state#

By using an excited state as the references state, the X-ray absorption spectra can be considered in a manner similar to XES:

1. Perform a ground state calculation

2. Use the result as initial guess for a second calculation, where the occupations are changed to emulate that of the desired (valence) excited state

3. Perform a wave function optimization with above configuration, using a constraint such as MOM to avoid a collapse to the ground state

4. Calculate the X-ray absorption spectra of this reference state

When considering (singlet) valence-excitations of a closed-shell molecule, two different reference state spin-muliticiplicity configurations are possible:

• Low-spin open-shell reference (LSOR), which moves an electron within the same multiplicity

• High-spin open-shell reference (HSOR), which flips the spin of the moved electron, and thus forms a triplet reference state (when considering singlet excitations of a closed-shell system)

Note that the LSOR wave function will be heavily spin contaminated, as seen below. This has been claimed to lead to minimal errors, but care should be taken. Further note that converging to a reference state which represents the actual valence excitation is non-trivial, and may indeed be all but impossible. As such, this two-step approach is useful, but has some issues. Finally, if using correlated methods, there are some concerns that the observed relative energies are off from experimental values, which could relate to issues with using a non-Aufbau reference state.

As an example we consider the first valence-excited state of water:

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

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

print('GS:')
scf_gs = scf.UHF(mol)
scf_gs.kernel()

print('LSOR:')
mo0 = copy.deepcopy(scf_gs.mo_coeff)
occ0 = copy.deepcopy(scf_gs.mo_occ)
occ0 = 0.0
occ0 = 1.0
scf_lsor = scf.UHF(mol)
scf_lsor.kernel()

print('HSOR:')
mo0 = copy.deepcopy(scf_gs.mo_coeff)
occ0 = copy.deepcopy(scf_gs.mo_occ)
occ0 = 0.0
occ0 = 1.0
scf_hsor = scf.UHF(mol)
scf_hsor.kernel()

GS:

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

LSOR:

converged SCF energy = -75.7142931989804  <S^2> = 1.0127856  2S+1 = 2.2474747

HSOR:
converged SCF energy = -75.7278282670684  <S^2> = 2.0038371  2S+1 = 3.002557

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

-75.72782826706842


As can be seen, the $$\langle S^2 \rangle$$ of the LSOR reference is close to 1, and thus represents an almost perfect mix of a singlet and triplet.

Calculating the valence-excitation energies and comparing it to energy differences obtained from the different reference states:

adc_gs = adcc.adc2(scf_gs, n_states=5)

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

  1    10       0.22058  543ms  [0.37446036 0.40081598 0.44565934 0.46679031 0.48399895]
2    20     0.0030002  143ms  [0.27713767 0.3056213  0.35895717 0.3715342  0.39063453]

  3    30     9.786e-05  180ms  [0.27640472 0.30463815 0.35766488 0.37047079 0.38929738]
=== Converged ===
Number of matrix applies:    30
Total solver time:           875.266ms

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

  1     8        1.6222  479ms  [20.26123431 20.31992537 20.32302691 20.3464984  20.77599022]
2    16      0.059729  107ms  [19.77470975 19.80230856 19.84312303 19.87818817 20.5373859 ]

  3    24     0.0025223  181ms  [19.75718509 19.78334186 19.84056229 19.86585919 20.53732181]
4    32    0.00064877  150ms  [19.75656768 19.78236576 19.84024749 19.86584883 20.53732181]

  5    38    2.0368e-05  145ms  [19.756435   19.78211629 19.84024749 19.86584882 20.53732181]
=== Converged ===
Number of matrix applies:    38
Total solver time:             1s  73ms

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

  1     8        1.4197  442ms  [19.65025265 20.45690115 20.50695083 20.60577589 20.87460994]
2    16       0.10186  109ms  [19.48097485 19.99916885 20.03363736 20.13668576 20.59740388]

  3    24     0.0025102  167ms  [19.48081058 19.97028865 20.02430938 20.10063975 20.59579242]
4    32    0.00010515  163ms  [19.48081058 19.96947945 20.02423709 20.09930154 20.59579241]
=== Converged ===
Number of matrix applies:    32
Total solver time:           891.704ms

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

  1     9        1.3402   1.2s  [19.65310247 20.44001891 20.46846981 20.50328038 20.87239002]
2    18      0.081792  159ms  [19.48153505 19.99518847 20.03915758 20.06588885 20.59952219]

  3    27     0.0012024  165ms  [19.48137709 19.98562519 20.0121162  20.05726204 20.59826671]
4    36    3.3226e-05  173ms  [19.48137709 19.98556281 20.01155346 20.05707419 20.59826671]
=== Converged ===
Number of matrix applies:    36
Total solver time:             1s 682ms

+--------------------------------------------------------------+
| adc2                                        any ,  converged |
+--------------------------------------------------------------+
|  #        excitation energy     osc str    |v1|^2    |v2|^2  |
|          (au)           (eV)                                 |
|  0     0.2764047      7.521355   0.0000    0.9634   0.03655  |
|  1     0.3046381      8.289626   0.0139    0.9591   0.04088  |
|  2     0.3576649      9.732557   0.0000    0.9686   0.03135  |
|  3     0.3704708      10.08102   0.0000     0.967   0.03295  |
|  4     0.3892974      10.59332   0.0000    0.9636   0.03642  |
+--------------------------------------------------------------+

-8.230663923526496

-7.85749980669949


The first allowed transition energy is relatively close to the $$\Delta$$MP2 results (as we use the MP2 energies), in particular for the LSOR. Resulting spectra can now be plotted as:

plt.figure(figsize=(6,3))
plt.legend()
plt.show() We see the occurance of a lower-energy feature, representing the 1s to SOMO transition, as well as an upward shift of remaining features.

• To be added

## MOM approaches#

Note

• Decide if to include

• IMOM versus MOM (currently only IMOM)

• Thanks to Anna Kristina and Josefine

As an illustration, the MOM approach from pyscf has here been adapted to IMOM, i.e. considering overlap to the initial rather than most recent wave function. This has been done for both restricted (som simultaneously moving 2 electrons at a time) and unrestricted wave functions.

def imom(mf, mocoeffs, occnum):
"""Use initial maximum overlap method to determine occupation number for each orbital in every
iteration. It can be applied to unrestricted as well as restricted open-shell systems."""

# Copy initial MO coefficients for occupied orbitals
# This can be adapted for MOM / IMOM - need to look into
if isinstance(mf, dft.rks.RKS) or isinstance(mf, scf.hf.RHF):
coef_occ_a = mocoeffs[:, occnum > 0]
elif isinstance(mf, dft.uks.UKS) or isinstance(mf, scf.uhf.UHF):
coef_occ_a = mocoeffs[:, occnum > 0]
coef_occ_b = mocoeffs[:, occnum > 0]

# Write new get_occ routine for DFT kernel
def get_occ(mo_energy=None, mo_coeff=None):
# Get new MO energy and coefficients
if mo_energy is None:
mo_energy = mf.mo_energy
if mo_coeff is None:
mo_coeff = mf.mo_coeff

# Restricted
if isinstance(mf, dft.rks.RKS) or isinstance(mf, scf.hf.RHF):
mo_occ = np.zeros_like(occnum)
nocc_a = int(np.sum(occnum) / 2)  # number of occupied orbitals

# Construct s_a = MO_coeff(old)*S*MO_coeff(new)
SC = np.matmul(mf.get_ovlp().T, mo_coeff)
s_a = np.matmul(coef_occ_a.T, SC)

# Get indices for orbitals with largest p = sum_i s_a(ij)s_a(ij)
idx_a = np.argsort(np.einsum("ij,ij->j", s_a, s_a))[::-1]

# Set occupation for selected MOs
mo_occ[idx_a[:nocc_a]] = 2.0

elif isinstance(mf, dft.uks.UKS) or isinstance(mf, scf.uhf.UHF):
mo_occ = np.zeros_like(occnum)
nocc_a = int(np.sum(occnum))  # number of occupied alpha orbitals
nocc_b = int(np.sum(occnum))  # number of occupied beta orbitals

# Construct s_a = MO_coeff(old)*S*MO_coeff(new)
SC_a = np.matmul(mf.get_ovlp().T, mo_coeff)
s_a = np.matmul(coef_occ_a.T, SC_a)
# Likewise for s_b
SC_b = np.matmul(mf.get_ovlp().T, mo_coeff)
s_b = np.matmul(coef_occ_b.T, SC_b)

# Get indices for orbitals with largest p_a = sum_i s_a(ij)s_a(ij)
idx_a = np.argsort(np.einsum("ij,ij->j", s_a, s_a))[::-1]
# Likewise for p_b
idx_b = np.argsort(np.einsum("ij,ij->j", s_b, s_b))[::-1]

# Set occupation for selected MOs
mo_occ[idx_a[:nocc_a]] = 1.0
mo_occ[idx_b[:nocc_b]] = 1.0

return mo_occ

# Redefine get_occ for DFT kernel as new function
mf.get_occ = get_occ

return mf


As examples, this can be used to calculate the single- and double-ionization energies of water:

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

# Index of core MO
core = 0

# Initial restricted/unrestricted SCF
dft_re = dft.RKS(mol)
dft_re.xc = "b3lyp"
dft_re.kernel()
dft_ur = dft.UKS(mol)
dft_ur.xc = "b3lyp"
dft_ur.kernel()

# Single ionization energy
mo0 = copy.deepcopy(dft_ur.mo_coeff)
occ0 = copy.deepcopy(dft_ur.mo_occ)
occ0[core] = 0.0
dft_ion = dft.UKS(mol)
dft_ion.xc = "b3lyp"
imom(dft_ion, mo0, occ0)
dft_ion.kernel()

# Double ionization energy
mo0 = copy.deepcopy(dft_re.mo_coeff)
occ0 = copy.deepcopy(dft_re.mo_occ)
occ0[core] = 0.0
dft_ion = dft.RKS(mol)
dft_ion.xc = "b3lyp"
imom(dft_ion, mo0, occ0)
dft_ion.kernel()

print(
"Single ionization energy:",
np.around(au2ev * (dft_ion.energy_tot() - dft_ur.energy_tot()), 2),
"eV",
)
print(
"Double ionization energy:",
np.around(au2ev * (dft_ion.energy_tot() - dft_ur.energy_tot()), 2),
"eV",
)

res = """
converged SCF energy = -76.3480474320894
converged SCF energy = -76.3480474320839  <S^2> = 2.5504043e-11  2S+1 = 1
converged SCF energy = -56.4125421194362  <S^2> = 0.75441021  2S+1 = 2.0044054
converged SCF energy = -32.9492316840399
Single ionization energy: 542.47 eV
Double ionization energy: 1180.94 eV
"""
print(res)

converged SCF energy = -76.3480474320894
converged SCF energy = -76.3480474320839  <S^2> = 2.5504043e-11  2S+1 = 1
converged SCF energy = -56.4125421194362  <S^2> = 0.75441021  2S+1 = 2.0044054
converged SCF energy = -32.9492316840399
Single ionization energy: 542.47 eV
Double ionization energy: 1180.94 eV


## Impact of channel restriction for TDDFT#

The influence of channel-restriction/CVS schemes for TDDFT can be considered for smaller systems by diagonalizing the full-space and CVS-space matrices, thus constructing the total global spectrum. This will not be possible for larger systems/basis sets, or for correlated method, where other CVS relaxation schemes have been used instead.

For a discussion of the underlying theory, see related sections, and here we primarily note that we will construct and resolve the RPA, obtaining transition energies as poles by solving the eigenvalue problem: $$$E^{}X = {\lambda} S^{}X$$$

where

$\begin{split} {\lambda} =\begin{pmatrix} \lambda & 0\\ 0 & -\lambda \end{pmatrix} \end{split}$

By rewriting the linear response equation, we obtain the oscillator strength as:

$\begin{split} (T^{\mu})^{2} = \sum_{\beta=x,y,z}{\mu_{\beta}^{}}^\dagger X_n X_n^{\dagger}\mu_{\beta}^{}\\ f_{osc}^{if} = \frac{2}{3}\cdot \omega_{if}\cdot(T^{\mu})^2 \end{split}$

This will now be done for water, where the ground state is calculated as:

basis = "6-31G"
vlx_bas = vlx.MolecularBasis.read(vlx_mol, basis)

scf_drv = vlx.ScfRestrictedDriver()
scf_settings = {"conv_thresh": 1.0e-6}
method_settings = {"xcfun": "b3lyp"}
scf_drv.update_settings(scf_settings, method_settings)
scf_drv.compute(vlx_mol, vlx_bas)


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

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

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

* Info * Molecular grid with 35950 points generated in 0.00 sec.

* Info * Overlap matrix computed in 0.00 sec.

* Info * Kinetic energy matrix computed in 0.00 sec.

* Info * Nuclear potential matrix computed in 0.00 sec.


* Info * Orthogonalization matrix computed in 0.00 sec.


* Info * SAD initial guess computed in 0.00 sec.


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


* Info * Overlap matrix computed in 0.00 sec.

* Info * Kinetic energy matrix computed in 0.00 sec.

* Info * Nuclear potential matrix computed in 0.00 sec.


* Info * Orthogonalization matrix computed in 0.00 sec.



Iter. |    Kohn-Sham Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change
--------------------------------------------------------------------------------------------
1       -76.384877806788    0.0000000000      0.05814398      0.01079242      0.00000000

                  2       -76.384691793132    0.0001860137      0.07420609      0.01375856      0.03127560

                  3       -76.385185991153   -0.0004941980      0.00038521      0.00008834      0.01651042

                  4       -76.385186023175   -0.0000000320      0.00007657      0.00002027      0.00029093

                  5       -76.385186023934   -0.0000000008      0.00000194      0.00000051      0.00002930

                  6       -76.385186023935   -0.0000000000      0.00000014      0.00000003      0.00000214


*** SCF converged in 6 iterations. Time: 0.45 sec.


               Spin-Restricted Kohn-Sham:
--------------------------
Total Energy                       :      -76.3851860239 a.u.
Electronic Energy                  :      -85.5413307434 a.u.
Nuclear Repulsion Energy           :        9.1561447194 a.u.
------------------------------------
Gradient Norm                      :        0.0000001389 a.u.

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

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

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

Molecular Orbital No.   2:
--------------------------
Occupation: 2.000 Energy:   -1.01751 a.u.
(   1 O   1s  :     0.21) (   1 O   2s  :    -0.46) (   1 O   3s  :    -0.45)
(   1 O   1p0 :     0.15) (   2 H   1s  :    -0.15) (   3 H   1s  :    -0.15)

Molecular Orbital No.   3:
--------------------------
Occupation: 2.000 Energy:   -0.52917 a.u.
(   1 O   1p+1:     0.52) (   1 O   2p+1:     0.23) (   2 H   1s  :    -0.27)
(   2 H   2s  :    -0.15) (   3 H   1s  :     0.27) (   3 H   2s  :     0.15)

Molecular Orbital No.   4:
--------------------------
Occupation: 2.000 Energy:   -0.35661 a.u.
(   1 O   2s  :    -0.20) (   1 O   3s  :    -0.37) (   1 O   1p0 :    -0.54)
(   1 O   2p0 :    -0.38)

Molecular Orbital No.   5:
--------------------------
Occupation: 2.000 Energy:   -0.29214 a.u.
(   1 O   1p-1:     0.65) (   1 O   2p-1:     0.50)

Molecular Orbital No.   6:
--------------------------
Occupation: 0.000 Energy:    0.05639 a.u.
(   1 O   2s  :    -0.15) (   1 O   3s  :    -1.09) (   1 O   1p0 :     0.29)
(   1 O   2p0 :     0.44) (   2 H   2s  :     0.94) (   3 H   2s  :     0.94)

Molecular Orbital No.   7:
--------------------------
Occupation: 0.000 Energy:    0.14625 a.u.
(   1 O   1p+1:    -0.42) (   1 O   2p+1:    -0.75) (   2 H   2s  :    -1.28)
(   3 H   2s  :     1.28)

Molecular Orbital No.   8:
--------------------------
Occupation: 0.000 Energy:    0.80982 a.u.
(   1 O   1p+1:    -0.20) (   1 O   2p+1:    -0.53) (   2 H   1s  :    -0.97)
(   2 H   2s  :     0.67) (   3 H   1s  :     0.97) (   3 H   2s  :    -0.67)

Molecular Orbital No.   9:
--------------------------
Occupation: 0.000 Energy:    0.87512 a.u.
(   1 O   2s  :     0.18) (   1 O   3s  :    -0.28) (   1 O   1p0 :    -0.78)
(   1 O   2p0 :     0.50) (   2 H   1s  :    -0.71) (   2 H   2s  :     0.59)
(   3 H   1s  :    -0.71) (   3 H   2s  :     0.59)

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



### Full-space results#

For the full-space results, we construct the Hessian and property gradients:

lrs = vlx.LinearResponseEigenSolver()
lrs.update_settings(scf_settings, method_settings)
rsp_drv = vlx.LinearResponseSolver()
rsp_drv.update_settings(scf_settings, method_settings)

# Electronic Hessian
E2 = lrs.get_e2(vlx_mol, vlx_bas, scf_drv.scf_tensors)

# Property gradients for dipole operator
"electric dipole", "x", vlx_mol, vlx_bas, scf_drv.scf_tensors
)
"electric dipole", "y", vlx_mol, vlx_bas, scf_drv.scf_tensors
)
"electric dipole", "z", vlx_mol, vlx_bas, scf_drv.scf_tensors
)

# Dimension
c = int(len(E2) / 2)

# Overlap matrix
S2 = np.identity(2 * c)
S2[c : 2 * c, c : 2 * c] *= -1

* Info * Molecular grid with 35950 points generated in 0.01 sec.

* Info * Processing Fock builds... (batch size: 160)
* Info *   batch 1/1


Calculating the excitation energies and oscillator strengths:

# Set up and solve eigenvalue problem
Sinv = np.linalg.inv(S2)  # for clarity - is identical
M = np.matmul(Sinv, E2)
eigs, X = np.linalg.eig(M)

# Reorder results
idx = np.argsort(eigs)
eigs = np.array(eigs)[idx]
X = np.array(X)[:, idx]

# Compute oscillator strengths
fosc = []
for i in range(int(len(eigs) / 2)):
j = i + int(len(eigs) / 2)  # focus on excitations
Xf = X[:, j]
Xf = Xf / np.sqrt(np.matmul(Xf.T, np.matmul(S2, Xf)))
tm = np.dot(Xf, V1_x) ** 2 + np.dot(Xf, V1_y) ** 2 + np.dot(Xf, V1_z) ** 2
fosc.append(tm * 2.0 / 3.0 * eigs[j])


### CVS-space results#

The dimensions of the CVS-space:

# Number of virtuals
nocc = vlx_mol.number_of_alpha_electrons()
nvirt = vlx.MolecularBasis.get_dimensions_of_basis(vlx_bas, vlx_mol) - nocc
n = nocc * nvirt

# CVS space
res_mo = 1
res_indx = 0
c = res_mo * nvirt


Hessian and property gradients:

## Construct CVS objects

# Define starting index for deexcitation
c_int_deex = n + res_indx

# CVS Hessian
E2_cvs = np.zeros((2 * c, 2 * c))
E2_cvs[0:c, 0:c] = E2[res_indx:c, res_indx:c]
E2_cvs[0:c, c : 2 * c] = E2[res_indx:c, c_int_deex : c_int_deex + c]
E2_cvs[c : 2 * c, 0:c] = E2[c_int_deex : c_int_deex + c, res_indx:c]
E2_cvs[c : 2 * c, c : 2 * c] = E2[
c_int_deex : c_int_deex + c, c_int_deex : c_int_deex + c
]

# CVS overlap matrix
S2_cvs = np.identity(2 * c)
S2_cvs[c : 2 * c, c : 2 * c] *= -1

# CVS property gradients
V1_cvs_x = np.zeros(2 * c)
V1_cvs_x[0:c] = V1_x[res_indx:c]
V1_cvs_x[c : 2 * c] = V1_x[c_int_deex : c_int_deex + c]
V1_cvs_y = np.zeros(2 * c)
V1_cvs_y[0:c] = V1_y[res_indx:c]
V1_cvs_y[c : 2 * c] = V1_y[c_int_deex : c_int_deex + c]
V1_cvs_z = np.zeros(2 * c)
V1_cvs_z[0:c] = V1_z[res_indx:c]
V1_cvs_z[c : 2 * c] = V1_z[c_int_deex : c_int_deex + c]


Calculating energies and intensities:

# Set up and solve eigenvalue problem
Sinv = np.linalg.inv(S2_cvs)  # for clarity - is identical
M = np.matmul(Sinv, E2_cvs)
eigs_cvs, X_cvs = np.linalg.eig(M)

# Reorder results
idx = np.argsort(eigs_cvs)
eigs_cvs = np.array(eigs_cvs)[idx]
X_cvs = np.array(X_cvs)[:, idx]

# Compute oscillator strengths
fosc_cvs = []
for i in range(int(len(eigs_cvs) / 2)):
j = i + int(len(eigs_cvs) / 2)  # focus on excitations
Xf_cvs = X_cvs[:, j]
Xf_cvs = Xf_cvs / np.sqrt(np.matmul(Xf_cvs.T, np.matmul(S2_cvs, Xf_cvs)))
tm_cvs = (
np.dot(Xf_cvs, V1_cvs_x) ** 2
+ np.dot(Xf_cvs, V1_cvs_y) ** 2
+ np.dot(Xf_cvs, V1_cvs_z) ** 2
)
fosc_cvs.append(tm_cvs * 2.0 / 3.0 * eigs_cvs[j])


### Comparison#

Comparing the obtained spectra:

plt.figure(figsize=(6,3))
x, y = eigs[int(len(eigs) / 2) :], fosc
x1, y1 = lorentzian(x, y, 518 / au2ev, 525 / au2ev, 0.001, 0.3 / au2ev)
plt.plot(x1 * au2ev, y1)
x, y = eigs_cvs[int(len(eigs_cvs) / 2) :], fosc_cvs
x2, y2 = lorentzian(x, y, 518 / au2ev, 525 / au2ev, 0.001, 0.3 / au2ev)
plt.plot(x2 * au2ev, y2)
plt.plot(x2 * au2ev, 10*(y1 - y2))
plt.legend(("Full-space", "CVS-space", "Difference x10"))
plt.tight_layout()
plt.show() We note a small difference between the two approaches. This difference has some dependence on the basis set, but will continue to be small for the K-edge.

Note

To be added: comparison to CVS functionalities in VeloxChem (currently unavailable).

## Ionization energies from XAS#

Can model the ionization energy by targetting excitations to extremely diffuse molecular orbital. Here: consider both spectrum and IE of water, at a CVS-ADC(2)-x/6-31G[+diffuse] level of theory.

• To be added

## Ionization energies from XES#

Can model the ionization energy by targetting decays from an extremely diffuse molecular orbital, with a core-electron already moved to it. The need of preparing this state makes the approach rather useless, but it can still provide insight into relaxation effects of XES. Here: consider the IE of water, at a CVS-ADC(2)/6-31G[+diffuse] level of theory.

• To be added

## XES considered in reverse#

Instead of modeling XES by looking at the (negative) eigenstates of a core-hole reference state, it is possible to use a valence-hole reference state, and then consider core-excitations into this hole. This thus models the XES process more from the point of view of XAS, i.e. in reverse. However, this approach presupposes:

1. That a single-particle picture works quite well, and that specific valence-hole will yield reasonable representations of the valence-hole left following decay into a core-hole

2. That the valence-hole can be converged well, and that the different (non-orthogonal) reference states from each calculation are sufficiently consistent with each other

Furthermore, the calculation is made more tedious, as

1. It requires one explicit calculation per state

As such, this approach is by not suggested for practical calculations, but it does illustrate the ease of which such exotic calculations can be run, as well as teach us something about the relaxation of core-transitions.

Considering the X-ray emission spectra of water, we create valence-holes in the four frontier valence MOs and construct a spectrum from this:

# Containers for resulting energies and intensities
xes_rev_E = []
xes_rev_f = []

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

# Copy molecular orbital coefficients
mo0 = copy.deepcopy(scf_res.mo_coeff)

# Reference calculation from core-hole
occ0 = copy.deepcopy(scf_res.mo_occ)
occ0 = 0.0
scf_ion = scf.UHF(mol)
scf_ion.kernel()
xes_E = -au2ev * adc_xes.excitation_energy

for i in range(4):
occ0 = copy.deepcopy(scf_res.mo_occ)
occ0[4 - i] = 0.0
# Perform unrestricted SCF calculation with MOM constraint
scf_ion = scf.UHF(mol)
scf_ion.kernel()
# Perform ADC calculation of first four states

plt.figure(figsize=(6, 3))