# X-ray Absorption Spectra for Aqueous Ammonia and Ammonium: Quantum
Mechanical versus Molecular Mechanical Embedding Schemes

Sarai Dery Folkestad  Alexander C. Paul  Regina Paul (née Matveeva)
Department of Chemistry, Norwegian University of Science and Technology, NTNU,
7491 Trondheim, Norway  Peter Reinholdt  Department of Physics, Chemistry and
Pharmacy, University of Southern Denmark, SDU, Campusvej 55, 5230 Odense,
Denmark  Sonia Coriani  Department of Chemistry, Technical University of
Denmark, DTU, Kemitorvet Bldg 207, 2800 Kongens Lyngby, Denmark  Michael
Odelius  Department of Physics, Stockholm University, 10691 Stockholm, Sweden
Henrik Koch  Department of Chemistry, Norwegian University of Science and
Technology, NTNU, 7491 Trondheim, Norway<EMAIL_ADDRESS>
###### Abstract

The X-ray absorption (XA) spectra of aqueous ammonia and ammonium are computed
using a combination of coupled cluster singles and doubles (CCSD) with
different quantum mechanical and molecular mechanical embedding schemes.
Specifically, we compare frozen Hartree–Fock (HF) density embedding,
polarizable embedding (PE), and polarizable density embedding (PDE).
Integrating CCSD with frozen HF density embedding is possible within the CC-
in-HF framework, which circumvents the conventional system-size limitations of
standard coupled cluster methods. We reveal similarities between PDE and
frozen HF density descriptions, while PE spectra differ significantly. By
including approximate triple excitations, we also investigate the effect of
improving the electronic structure theory. The spectra computed using this
approach show an improved intensity ratio compared to CCSD-in-HF. Charge
transfer analysis of the excitations shows the local character of the pre-edge
and main-edge, while the post-edge is formed by excitations delocalized over
the first solvation shell and beyond.

###### keywords:

American Chemical Society, LaTeX

These authors contributed equally to this work. These authors contributed
equally to this work. These authors contributed equally to this work.
IR,NMR,UV

##  1 Introduction

Over the last fifty years, ammonia (NH3) and ammonium (NH${}_{4}^{+}$) have
raised the interest of the scientific community for various reasons: One
reason is the controversy surrounding their solvation structure.1 Another
point of interest is the rapid rotation of ammonia in liquid water, which has
been under debate for several years.2 Ammonia and ammonium are also known in
the context of wastewater treatment as the primary forms of inorganic
nitrogen. 3 Moreover, ammonia solution has emerged as a potentially effective
CO2 absorbent, 4, 5 an important feature in times of increasing environmental
pollution.

Ammonia and ammonium have been the subject of numerous experimental and
theoretical studies. 6, 7, 8, 9, 10, 1, 11, 12, 13, 14 Among these are
investigations using X-ray absorption (XA) spectroscopy.1, 12 When simulating
XA spectra of liquids and solutions, it is crucial to account for solvent
effects, as the electronic spectra of a molecule in solution are directly
influenced by a solvent. When calculating such spectra at a fully quantum
mechanical level, including a large number of solvent molecules poses a
challenge. Therefore, hybrid methods, which involve treating solute and
solvent at different levels of theory, have been actively used to describe
aqueous solutions. These methods combine high-level quantum mechanical (QM)
calculations with more approximate models such as polarizable continuum (PCM),
molecular mechanics (MM), and Hartree–Fock (HF) or Kohn-Sham density
embedding.15, 16

Recently, Reinholdt et al.12 presented XA spectra of NH3 and NH${}_{4}^{+}$ in
water computed using the combination of coupled cluster (CC) methods with
polarizable embedding (PE). Coupled cluster singles and approximate doubles
(CC2) and coupled cluster singles and doubles (CCSD) were the methods used.
The spectra obtained using CCSD showed satisfactory results and a notable
improvement compared to conventionally used transition-potential density
functional theory (TP-DFT).17 However, as the authors themselves pointed out,
it would be beneficial to employ multilevel coupled cluster methods and/or
extend the number of water molecules in the QM region to validate their
conclusions.

In this work, we model the XA spectra of aqueous NH3 and NH${}_{4}^{+}$ using
polarizable density embedding (PDE)-CCSD and coupled cluster methods embedded
in a frozen Hartree–Fock (HF) density (CC-in-HF),18, 19, 20 and compare our
results to the PE-CCSD spectra from Ref. 12. The CC-in-HF calculations were
performed with a development version of the $e^{\mathcal{T}}$ program,19 while
the PDE-CCSD calculations were carried out with the Dalton program21 (version
2022).

##  2 Theory

The coupled cluster wave function is obtained through the exponential ansatz,

$\displaystyle|\mathrm{CC}\rangle=e^{T}|\mathrm{R}\rangle~{},$ (1)

where the exponential of the cluster operator, $T$, acts on a reference
determinant (typically the Hartree–Fock state). The cluster operator generates
excitations of the reference, and can be written in terms of excitation
operators, $\tau_{\mu}$, and corresponding amplitudes, $t_{\mu}$:

$\displaystyle
T=\sum_{\mu}t_{\mu}\tau_{\mu}=T_{1}+T_{2}+T_{3}+\ldots+T_{N_{e}}~{}.$ (2)

For instance, single, double, and triple excitations of the reference are
generated by $T_{1}$, $T_{2}$, and $T_{3}$, respectively.

The coupled cluster energy is obtained by projecting the Schrödinger equation
onto the reference determinant,

$\displaystyle\langle\mathrm{R}\hskip 1.00006pt|\hskip
1.00006pte^{-T}He^{T}\hskip 1.00006pt|\hskip
1.00006pt\mathrm{R}\rangle=E_{\textrm{CC}}~{},$ (3)

and the cluster amplitudes are determined by projection onto the space of
excited determinants $\langle\mu|=\langle\textrm{R}|\tau_{\mu}^{\dagger}$ ,

$\displaystyle\Omega_{\mu}=\langle\mu\hskip 1.00006pt|\hskip
1.00006pte^{-T}He^{T}\hskip 1.00006pt|\hskip 1.00006pt\mathrm{R}\rangle=0~{}.$
(4)

In practice, the cluster operator must be truncated at a certain excitation
level, which defines the hierarchy of standard coupled cluster models. Among
the most commonly used models is the coupled cluster singles and doubles
(CCSD) model, where $T=T_{1}+T_{2}$. In the spin-adapted closed-shell
formulation of the theory, we have

$\displaystyle T_{1}=\sum_{ai}t^{a}_{i}E_{ai},\quad
T_{2}=\frac{1}{2}\sum_{aibj}t^{ab}_{ij}E_{ai}E_{bj},$ (5)

where $E_{ai}$ is a singlet excitation operator, $t^{a}_{i}$ and $t^{ab}_{ij}$
are the cluster amplitudes. We use indices $\\{i,j,\ldots\\}$ for occupied
orbitals, $\\{a,b,\ldots\\}$ for virtual orbitals, and $\\{p,q,\ldots\\}$ to
denote general orbitals.

Excitation energies and transition moments can be calculated with linear
response theory or the equation-of-motion22, 23, 24 (EOM) framework. In EOM
coupled cluster theory, the ground and excited states are expressed as the
eigenvectors of the similarity-transformed Hamiltonian,

$\displaystyle\boldsymbol{\bar{H}}=\begin{pmatrix}0&\boldsymbol{\eta}^{T}\\\
0&\boldsymbol{A}\\\ \end{pmatrix}+E_{\mathrm{CC}}\boldsymbol{I},$ (6)

where $A_{\mu\nu}=\langle\mu\hskip 1.00006pt|\hskip
1.00006pt[\bar{H},\tau_{\nu}]\hskip 1.00006pt|\hskip
1.00006pt\mathrm{HF}\rangle$ is the Jacobian matrix and
$\eta_{\nu}=\langle\mathrm{HF}\hskip 1.00006pt|\hskip
1.00006pt[\bar{H},\tau_{\nu}]\hskip 1.00006pt|\hskip
1.00006pt\mathrm{HF}\rangle$. The ground state amplitudes have been determined
from eq (4). The $k$’th state is determined by

$\displaystyle\boldsymbol{\bar{H}}\boldsymbol{R}^{k}$
$\displaystyle=(E_{\mathrm{CC}}+\omega_{k})\boldsymbol{R}^{k}$ (7)
$\displaystyle\boldsymbol{\bar{H}}^{T}\boldsymbol{L}^{k}$
$\displaystyle=(E_{\mathrm{CC}}+\omega_{k})\boldsymbol{L}^{k},$ (8)

where $\omega_{k}$ is the corresponding excitation energy (for $k>0$). Since
$\boldsymbol{\bar{H}}$ is non-Hermitian, the left, $\boldsymbol{L}^{k}$, and
right, $\boldsymbol{R}^{k}$, eigenvectors differ, and so do the state vectors,

$\displaystyle|k\rangle=\sum_{\mu\geq
0}e^{T}R^{k}_{\mu}|\mu\rangle,\quad\langle k|=\sum_{\mu\geq
0}L^{k}_{\mu}\langle\mu|e^{-T}.$ (9)

We need both left and right state vectors to calculate, e.g., oscillator
strengths:25, 26

$\displaystyle f_{k}=\frac{2}{3}\omega_{k}\langle 0\hskip 1.00006pt|\hskip
1.00006pt\boldsymbol{\mu}\hskip 1.00006pt|\hskip 1.00006ptk\rangle\cdot\langle
k\hskip 1.00006pt|\hskip 1.00006pt\boldsymbol{\mu}\hskip 1.00006pt|\hskip
1.00006pt0\rangle,$ (10)

where $0$ denotes the coupled cluster ground states.

Core excitations often have double excitation character, and accurate modeling
relies on the inclusion of triple excitations in $T$. While CCSD scales as
$\mathcal{O}(N^{6})$ with system size, the exact inclusion of triple
excitations, as in CCSDT, will yield $\mathcal{O}(N^{8})$-scaling equations, a
cost that is prohibitive for most molecular systems. Triple excitations are,
therefore, often included perturbatively at an $\mathcal{O}(N^{7})$-cost, as
in the CC327, 28 model.

In this work, we calculate the XA spectra of aqueous ammonia and ammonium.
Accurately simulating XA spectra of solutions or liquids requires many
representative geometries and a high-level electronic structure description of
both solute and neighboring solvent molecules. Hence, a large number of
calculations on a sizable molecular system must be carried out, rendering an
$\mathcal{O}(N^{7})$-scaling model such as CC3 unfeasible.

###  2.1 Multilevel coupled cluster theory

Multilevel coupled cluster (MLCC) theory29 can be used to reduce the cost of
coupled cluster calculations while retaining high accuracy in intensive
properties. In MLCC, the highest-order excitations in the cluster operator are
restricted to an active space. For instance, in MLCCSDT, $T_{3}$ is confined
to an active orbital space, whereas $T_{1}$ and $T_{2}$ include single and
double excitations between all occupied and virtual orbitals in the molecular
system. The ground and excited states are solved analogously to the standard
coupled cluster models. However, the projection space is constrained to
correspond to the excitations included in $T$. In this paper, we use MLCC3,30,
31 a multilevel analog of CC3 with triple excitations restricted to an active
space.

A similar active space strategy can also define a Hartree–Fock density
embedding scheme: by restricting all cluster amplitudes to the active orbital
space, we get the CC-in-HF15, 32, 18 model. For example, in the CCSD-in-HF
model, both $T_{1}$ and $T_{2}$ are restricted. Since high-scaling equations
are limited to the active orbital space, we can use CC-in-HF on large
molecular systems. Unlike QM/MM approaches, the CC-in-HF models offer a single
wave function for the whole molecular system under consideration. Thereby, we
include Pauli-repulsion between the target (active region) and environment
(rest of the system) regions. In the CC-in-HF scheme, an effective Fock matrix
is used to include the interactions with the environment in the coupled
cluster calculation:

$\displaystyle
f^{\text{eff}}_{\mu\nu}=h_{\mu\nu}+\sum_{\gamma\delta}\Big{(}(\mu\nu|\gamma\delta)-\frac{1}{2}(\mu\delta|\gamma\nu)\Big{)}D^{\mathrm{target}}_{\gamma\delta}+\sum_{\gamma\delta}\Big{(}(\mu\nu|\gamma\delta)-\frac{1}{2}(\mu\delta|\gamma\nu)\Big{)}D^{\mathrm{env}}_{\gamma\delta}.$
(11)

Here, $h_{\mu\nu}$ are the one-electron integrals in the atomic orbital (AO)
basis, $(\mu\nu|\gamma\delta)$ are the electron repulsion integrals in
Mulliken ordering, and $D^{\mathrm{target}}_{\gamma\delta}$ and
$D^{\mathrm{env}}_{\gamma\delta}$ are the active and environment densities,
respectively. The total density,

$\displaystyle{D}_{\mu\nu}={D}^{\mathrm{target}}_{\mu\nu}+{D}^{\mathrm{env}}_{\mu\nu},$
(12)

is optimized at the Hartree–Fock level of theory, and the partitioning of the
orbital space occurs between the Hartree–Fock and coupled cluster
calculations. Note that all AOs—both of the target region and the
environment—can contribute to the active molecular orbitals, i.e., the AO
indices $\mu$ and $\nu$ are not restricted to the target region in eq (11).
Conversely, AOs centered on the active region can contribute to the
environment density. A reduced dimension of the effective Fock matrix is only
obtained through transformation to the active MO basis.

Finally, we can combine the MLCC and the CC-in-HF approaches to obtain the so-
called MLCC-in-HF model. This can be done by introducing two concentric active
spaces, where all $t$-amplitudes are restricted to the outer space, and the
high-order $t$-amplitudes are limited to the inner space.18

Successful application of the CC-in-HF and MLCC approaches depends on the
partitioning of the orbital space. We need to choose a suitable type of
Hartree–Fock orbitals and determine the size of the active space. For the XA
spectra of ammonia and ammonium, we select the outer active orbital layer
using localized orbitals. This ensures that the solute and a fixed number of
solvent molecules are included in the correlated part of the calculation. Any
kind of localized Hartree–Fock orbitals can be used. In this work, we use
Cholesky occupied orbitals15 and orthonormalized projected atomic orbitals33,
34 (PAOs) for the virtual space. The Cholesky orbitals result from a partial
Cholesky factorization of the Hartree–Fock density in the atomic orbital
basis:

$\displaystyle{D}_{\mu\nu}=2\sum_{i}C^{i}_{\mu}C^{i}_{\nu}+{D}^{\mathrm{env}}_{\mu\nu}.$
(13)

The elements of the Cholesky vector, $\boldsymbol{C}^{i}$, are the orbital
coefficients of the active occupied orbital $i$. The partial factorization is
obtained by restricting the pivots of the decomposition to the atomic orbitals
of the target region. Furthermore, to ensure a tight active space, we exclude
pivoting elements corresponding to diffuse orbitals with exponents smaller or
equal to 0.6.

For the virtual space, active orbitals are constructed by projecting the
occupied orbitals out of the atomic orbitals centered on the target region,

$\displaystyle\boldsymbol{C}^{\text{PAO}}=\boldsymbol{I}-\frac{1}{2}\boldsymbol{DS}^{\mathrm{target}}.$
(14)

The overlap matrix $\boldsymbol{S}^{\mathrm{target}}$ is rectangular and
contains the columns of the AO overlap matrix that correspond to AOs centered
on the active atoms. Before the coupled cluster calculation, we remove linear
dependencies and orthonormalize the remaining virtual orbitals.

For the inner active space, used in the MLCC3 calculations, we use correlated
natural transition orbitals35 (CNTOs). The CNTOs are constructed from CCSD (or
CCSD-in-HF) excitation vectors, $\boldsymbol{R}^{k}$, by diagonalizing the
unit-trace matrices

$\displaystyle M^{k}_{ij}=$
$\displaystyle\sum_{a}R^{k}_{ai}R^{k}_{aj}+\frac{1}{2}\sum_{abk}(1+\delta_{ai,bk}\delta_{ij})R^{k}_{aibk}R^{k}_{ajbk}$
(15)

and

$\displaystyle N^{k}_{ab}=$
$\displaystyle\sum_{i}R^{k}_{ai}R^{k}_{bi}+\frac{1}{2}\sum_{ijc}(1+\delta_{ai,cj}\delta_{ab})R^{k}_{aicj}R^{k}_{bicj},$
(16)

and transforming the occupied and virtual blocks of the orbital coefficient
matrix with the resulting eigenvector matrices. Using CNTOs, a compact
representation of the excitation captured by $\boldsymbol{R}^{k}$ can be
achieved. To obtain an orbital space that targets a collection of excited
states, the matrices $\boldsymbol{M}=\sum_{k}\boldsymbol{M}^{k}$ and
$\boldsymbol{N}=\sum_{k}\boldsymbol{N}^{k}$ are diagonalized. The resulting
eigenvector matrices are used as transformation matrices for the occupied and
virtual orbital coefficients, respectively.

###  2.2 Polarizable density embedding

The polarizable density embedding36 scheme shares similarities with the
Hartree–Fock density embedding outlined above, and is an extension of
polarizable embedding (PE) approach by Olsen et al.37 In PDE, contributions
from the environment are added to the Fock matrix of the target system
according to

$\boldsymbol{f}^{\text{eff}}=\boldsymbol{f}+\boldsymbol{v}^{\text{env},\text{PDE}}.$
(17)

The embedding operator contains three terms, accounting for the electrostatic
(es), non-electrostatic repulsion (rep), and induction (ind) interactions

$\boldsymbol{v}^{\text{env},\text{PDE}}=\boldsymbol{v}^{\text{es}}+\boldsymbol{v}^{\text{rep}}+\boldsymbol{v}^{\text{ind}}.$
(18)

The electrostatic term includes the effects of the nuclear-electron attraction
and the electron-electron repulsion between electrons in the target system and
the nuclei and electrons of the environment; the latter of which is contained
through frozen fragment densities, $\boldsymbol{D}^{f}$, of the environment
molecules:

$v^{\mathrm{es}}_{\mu\nu}=-\sum_{f}\sum_{m\in
f}Z_{m}v_{\mu\nu}(\mathbf{R}_{m})+\sum_{f}\sum_{\gamma\delta\in
f}D^{f}_{\gamma\delta}(\mu\nu|\gamma\delta).$ (19)

Here, $Z_{m}$ is the charge of nucleus $m$ in fragment $f$, $v_{\mu\nu}$ is a
nuclear attraction integral, $D_{\gamma\delta}^{f}$ is the density matrix of
environment fragment $f$, and $(\mu\nu|\gamma\delta)$ are (intermolecular)
electron repulsion integrals.

Similarly to the CC-in-HF approach, the PDE model retains the “exact”
electron-electron repulsion from the electronic densities of the environment
rather than approximating this interaction with a distributed multipole model,
as is done in the simpler PE models. However, unlike in CC-in-HF, the
environment density in PDE is made up of frozen fragment densities that are
calculated in isolation36; in CC-in-HF the total density
($\boldsymbol{D}=\boldsymbol{D}^{\mathrm{target}}+\boldsymbol{D}^{\mathrm{env}}$)
is optimized at the Hartree–Fock level of theory. As a result, the total
Hartree–Fock density ($\boldsymbol{\rho}=\frac{1}{2}\boldsymbol{D}$) is
idempotent in CC-in-HF, which is not the case for PDE.

Compared to PE, the PDE potential also includes Pauli-repulsion through a
Huzinaga-Cantu-style38 repulsion operator

$v^{\mathrm{rep}}_{\mu\nu}=-\sum_{f}\sum_{\gamma\delta\in
f}\sum_{i}^{N_{\mathrm{occ}},f}\varepsilon_{i}^{f}C^{f}_{\gamma
i}C^{f}_{\delta i}S_{\mu\gamma}S_{\nu\delta},$ (20)

which penalizes the wave function overlap between the target region and the
environment fragment wave functions. The repulsion operator depends on the
environment fragment orbital energies $\varepsilon_{i}^{f}$, MO coefficients
$C_{\gamma i}^{f}$, and the intermolecular overlap integrals $S_{\mu\gamma}$.
This can be compared to the exchange contribution to the target system Fock
matrix in CC-in-HF (see eq (11)).

Finally, the PDE models include mutual polarization effects between the target
region and the environment which is described by an induced dipole model, as
in PE37. The induction operator appears as

$v^{\mathrm{ind}}_{\mu\nu}=-\sum_{m=1}^{N}\sum_{\alpha=x,y,z}{\mu}^{\mathrm{ind}}_{m,\alpha}{t}_{\mu\nu,\alpha}^{(1)}({\mathbf{R}}_{m}).$
(21)

where ${\mu}^{\mathrm{ind}}_{m,\alpha}$ is a component of an induced dipole
moment and ${t}_{\mu\nu,\alpha}^{(1)}$ are matrix elements of the electric
field operator from the QM region. The induced dipoles are obtained as
solutions to the linear problem

$\boldsymbol{\mathbf{\mu}}_{m}^{\mathrm{ind}}\left(\mathbf{F}_{\mathrm{tot}}\right)=\boldsymbol{\alpha}_{m}\mathbf{F}_{\mathrm{tot}}(\mathbf{R}_{m})=\boldsymbol{\alpha}_{m}\left(\mathbf{F}(\mathbf{R}_{m})+\sum_{m^{\prime}\neq
s}\mathbf{T}_{mm^{\prime}}^{(2)}\boldsymbol{\mathbf{\mu}}_{m^{\prime}}^{\mathrm{ind}}\right),$
(22)

where $\mathbf{F}$ is the electric field, and $\mathbf{T}^{(2)}$ is the
dipole-dipole interaction tensor. Thus, the polarizability
$\boldsymbol{\alpha}_{m}$ creates an induced dipole in response to the
electric fields from the QM density, the static fields from all other
environment fragments, and the fields from every other induced dipole moment
in the environment. As a result, the induced dipoles depend on the QM density,
which in turn depends on the values of the induced dipoles, ensuring mutual
polarization between the environment and the QM density.

###  2.3 Charge transfer analysis of excitations

When a molecular system absorbs a photon, the electrons rearrange themselves,
causing a change in density,

$\displaystyle\Delta\boldsymbol{D}_{k}=\boldsymbol{D}_{k}-\boldsymbol{D}_{0},$
(23)

where $\boldsymbol{D}_{k}$ is the density of an excited state and
$\boldsymbol{D}_{0}$ the ground state density. The trace of the one-electron
density equals the number of electrons. Therefore, the trace of the density
difference, $\Delta\boldsymbol{D}_{k}$, must be zero. Using, for example,
Löwdin population analysis, we can transform the density difference from the
delocalized molecular orbital (MO) basis to the local atomic orbital (AO)
basis,

$\displaystyle\Delta\boldsymbol{\rho}^{\mathrm{AO}}_{k}=\boldsymbol{S}^{\frac{1}{2}}\,\boldsymbol{C}\,\Delta\boldsymbol{D}_{k}\,\boldsymbol{C}^{T}\,\boldsymbol{S}^{\frac{1}{2}},$
(24)

where $\boldsymbol{S}$ is the AO overlap matrix and $\boldsymbol{C}$ contains
the MO-coefficients. Since
$\boldsymbol{C}^{T}\boldsymbol{S}\boldsymbol{C}=\boldsymbol{I}$, the trace of
the density difference remains zero in the AO basis. We can now partition the
trace of $\Delta\boldsymbol{\rho}^{\mathrm{AO}}_{k}$ into contributions from
the AOs on different fragments of the molecular system. For instance, for the
calculation of XA spectra of a solvent or liquid, we can partition the trace
with respect to (A) AOs on the the solute, (B) AOs on the water molecules
closest to the solute, and (C) the rest of the AOs:

$\displaystyle\begin{split}\mathrm{Tr}(\Delta\boldsymbol{\rho}^{\mathrm{AO}}_{k})&=\sum_{\alpha}[\Delta\boldsymbol{\rho}^{\mathrm{AO}}_{k}]_{\alpha\alpha}\\\
&=\sum_{\alpha\in
A}|\boldsymbol{\rho}\boldsymbol{D}^{\mathrm{AO}}_{k}]_{\alpha\alpha}+\sum_{\alpha\in
B}|\boldsymbol{\rho}\boldsymbol{D}^{\mathrm{AO}}_{k}]_{\alpha\alpha}+\sum_{\alpha\in
C}|\boldsymbol{\rho}\boldsymbol{D}^{\mathrm{AO}}_{k}]_{\alpha\alpha}\\\
&=0.\end{split}$ (25)

The trace of a single subsystem is then interpreted as the number of electrons
detached from or attached to this subsystem, depending on the sign. Such
analysis has been reported before for TDDFT using Mulliken population analysis
in Ref. 39 and more generally in Ref. 40.

##  3 Results and discussion

Using PDE-CCSD and CC-in-HF, we describe the solute and four closest water
molecules with coupled cluster theory, whereas the remaining solvent molecules
are described at a lower level of theory. These methods serve to reduce the
high cost associated with applying coupled cluster methods to large systems
without compromising the accuracy in intensive molecular properties. The
calculations are performed on the 195 spherical cluster geometries studied in
Ref. 12. For the PDE36, electron densities at the CAMB3LYP41/aug-cc-pVTZ level
are computed for each molecule except the solute and four closest water
molecules. The embedding potential is constructed from these fragment
densities and atomic polarizabilities, obtained with the LoProp42 method. The
spectra are calculated from a CCSD wave function optimized in the presence of
the embedding potential. For CC-in-HF, the orbitals localized on the solute
and the four closest water molecules are included in the coupled cluster
region. The remaining orbitals are described using a frozen HF density.15, 19,
18 We use two versions of CC-in-HF: CCSD-in-HF and MLCC3-in-HF (multilevel
coupled cluster singles and doubles with perturbative triples),30, 31 which
are described in more detail in Ref. 20. In all calculations the core-valence
separation approximation 43, 44, 45 is used to target core excitations.

Figure 1:  XA spectra of ammonia (panels A, C) and ammonium (panels B, D) in
water clusters for different embedding methods. Solute and four closest water
molecules are described with 6-311++G∗∗ basis, while remaining solvent
molecules are treated at the PE, PDE, and HF/6-31G∗∗ levels of theory. The PE-
CCSD, PDE-CCSD, CCSD-in-HF, and MLCC3-in-HF results have been shifted by
$-1.9\text{\,}\mathrm{eV}$, $-2.8\text{\,}\mathrm{eV}$,
$-2.5\text{\,}\mathrm{eV}$, $-1.5\text{\,}\mathrm{eV}$, respectively. The PE-
CCSD data was re-plotted based on calculations from Ref. 12, and the
experimental data was digitized from Ref. 1 using Web Plot Digitizer46. Figure
2:  Comparison of XA spectra computed using Pople (Panels A, B) and Dunning
(Panels C, D) basis sets. XA spectra of ammonia (Panels A, C) and ammonium
(Panels C, D) in water clusters at the CCSD-in-HF and MLCC3-in-HF levels of
theory. In the upper row 6-311++G∗∗ was used in the CCSD space and 6-31G∗∗ for
the remaining solvent, while in the lower row aug-cc-pVTZ was used for the
solute, aug-cc-pVDZ for the four closest water molecules and cc-pVDZ for the
remaining solvent. Constant energy shifts of $-2.5\text{\,}\mathrm{eV}$ and
$-2.0\text{\,}\mathrm{eV}$ were applied for CCSD-in-HF and
$-1.5\text{\,}\mathrm{eV}$, $-1.0\text{\,}\mathrm{eV}$ for MLCC3-in-HF for the
Pople and Dunning bases, respectively.

To compare with the study of Reinholdt et al., 12 we performed calculations
using the same basis set (6-311++G**) in the coupled cluster region. The
remaining water molecules were described using either PDE or HF embedding with
the 6-31G** basis set. The resulting spectra for both ammonia and ammonium are
presented in Figure 1. We use a constant broadening with Voigt profiles of
$0.2\text{\,}\mathrm{eV}$ Lorentzian full width at half maximum (FWHM) and
Gaussian standard deviation (yielding a total FWHM of $\sim
0.59\text{\,}\mathrm{eV}$) to obtain smooth theoretical spectra. The average
spectra are normalized such that the area underneath the curve equals that of
the experiment in the ranges $400\text{\,}\mathrm{eV}$ –
$408\text{\,}\mathrm{eV}$ and $402.5\text{\,}\mathrm{eV}$ –
$411\text{\,}\mathrm{eV}$ for NH3 and NH${}_{4}^{+}$, respectively. The
theoretical spectra for ammonia were shifted to align approximately with the
pre-edge of the experiment (see Figure 1). For ammonium the same shifts were
applied.

Comparing the theoretical spectra of ammonia and ammonium, shown in Figure 1,
we observe that the spectral shapes of the PE-CCSD spectra agree better with
experiment compared to both PDE-CCSD and CCSD-in-HF. Since both the PDE and HF
embedding schemes offer a higher level of theoretical description of the
solvent molecules, the favourable agreement of PE-CCSD with experiment is
likely due to cancellation of errors. For instance, a density-based embedding
scheme includes the effects of Pauli-repulsion, which is not the case for PE.
We note that PDE-CCSD and CCSD-in-HF exhibit very similar spectral features.
Both methods overestimate the ratio between the main-edge (at
$402.8\text{\,}\mathrm{eV}$) and the post-edge ($403\text{\,}\mathrm{eV}$ to
$409\text{\,}\mathrm{eV}$) of the ammonia spectrum (see Panels A and C of
Figure 1), compared to experiment and PE-CCSD.

The main-edges in PE-CCSD and PDE-CCSD exhibit a shoulder at approximately
$404.5\text{\,}\mathrm{eV}$, a feature not observed for CC-in-HF. It has been
shown for a water dimer that PDE is overestimating the repulsion energy
compared to HF.36 Therefore, in CC-in-HF, the transitions in this energy
region are probably shifted to lower energies so that the shoulder merges with
the main-edge. With the selected normalization scheme, the intensity of the
post-edge region is underestimated by PDE and HF-embedding, except for a very
broad feature at $\sim 408\text{\,}\mathrm{eV}$. This broad band is shifted to
lower energies for CCSD-in-HF compared to PDE-CCSD. In contrast to PE-CCSD, we
do not observe a spectral feature at the ionization limit at approximately
$407\text{\,}\mathrm{eV}$. Similar to the study by Reinholdt et al.,12 we do
not treat the continuum differently than the other regions of the spectrum.
Therefore, we conclude that this spectral feature is more likely an artifact
of the PE description, rather than the result of the discretized description
of the continuum.

Inclusion of approximate triple excitations in MLCC3-in-HF compared to CCSD-
in-HF yields an improved intensity ratio between the main- and post-edge in
the ammonia spectrum, resulting in a smoother transition between the two
regions. The absolute peak positions in MLCC3-in-HF are more accurate than in
CCSD-in-HF.

For aqueous ammonium (see panels B and D of Figure 1), the important spectral
characteristics – the shoulder between $403\text{\,}\mathrm{eV}$ and
$404\text{\,}\mathrm{eV}$, the main-edge at $405.7\text{\,}\mathrm{eV}$, and
the post-edge between $407\text{\,}\mathrm{eV}$ and $412\text{\,}\mathrm{eV}$
– are reproduced with PDE-CCSD and CCSD-in-HF. However, compared to PE-CCSD,
the transition from main- to post-edge is characterized by a significant drop
in intensity. This intensity decrease is more pronounced for PDE-CCSD, because
the broad post-edge feature is shifted to higher energies compared to CCSD-in-
HF.

Similar to what we observed for ammonia, the relative intensities of the main-
and post-edge are improved going from CCSD-in-HF to MLCC3-in-HF for ammonium:
in the MLCC3-in-HF spectrum, the main-edge intensity is reduced and the post-
edge intensity is increased. Additionally, we observe a red-shift of the post-
edge feature from CCSD-in-HF to MLCC3-in-HF. While the PE-CCSD spectrum agrees
best with the experiment, it contains a rather sharp feature in the range
$\qtyrange{408}{410}{}$ which is much broader in the other embedding schemes.
This additional peak resembles the one in ammonia, and further supports that
it might be an artifact of the polarizable embedding.

To examine the basis set effect on the spectra, we test various basis set
combinations for the solute, the closest four waters and the remaining solvent
molecules on an arbitrarily chosen snapshot. The results are shown in Figures
S1 and S2 in the Supporting Information. A combination of aug-cc-pVTZ on the
solute, aug-cc-pVDZ on the four closest water molecules, and cc-pVDZ on the
remaining solvent molecules was chosen as a compromise between accuracy of the
basis and computational efficiency. Using these basis sets, the XA spectra of
both ammonia and ammonium were recomputed and compared with previous spectra
(see Figure 2). For both ammonia and ammonium, the change of basis set leaves
the pre- and main-edges of the spectra unchanged. However, for ammonia, the
post-edge feature is broadened and red-shifted when the Dunning basis sets are
used. This change in spectral shape shows that the post-edge feature of
ammonia observed using the Pople basis is not a real resonance, but rather an
artifact of a limited basis set. For ammonium, on the other hand, the post-
edge feature around $408.5\text{\,}\mathrm{eV}$ does not shift, and
compression of the spectrum with increased basis set size is only observed
above $409\text{\,}\mathrm{eV}$.

Figure 3:  XA spectrum plotted together with the charge transfer number (CT)
calculated for the solute molecules (ammonia, left column) and (ammonium right
panel). In the lower panel, each excitation has been assigned a color based on
its charge transfer character. In the upper panel, we show the broadened
spectrum where the color gradient indicates the charge transfer character
averaged for all transitions within a bin of the size equal to the FWHM.

We have performed a charge transfer analysis for aqueous ammonia and ammonium
solutions with MLCC3-in-HF. This analysis is based on the calculation of the
so-called charge transfer number for each excited state of the solute
molecule.39, 40, 47, 48 The charge transfer number shows the significance of
the interaction with the solvent. A negative value indicates the removal of
electrons from the solute while a positive value implies their addition.
Further details on the charge transfer analysis can be found in Supporting
Information and Ref. 20. The results are presented in Figure 3. As described
above, and evident from the similarity in the respective values for NH3 and
NH${}_{4}^{+}$, the charge transfer number acts as an indicator of
environmental influence rather than actual electron transfer. From Figure 3 we
see that the excitations in the pre- and main-edge in the XA spectra of both
molecular systems are localized on the respective solute molecules. In the
post-edge regions, we observe an increasing charge transfer indicating larger
influence from the solvent molecules. This reflects the more diffuse nature of
excitations in these regions.

In Figure S3 in the Supporting Information we show the charge transfer from
the solute and four surrounding water molecules to the remaining solvent. We
see that the excitations responsible for generating the pre-edge, main-edge,
and the first half of the post-edge are well described by the four water
molecules in the coupled cluster region. An increasing degree of charge
transfer in the second half of the post-edge is observed. However, it is
important to note that the charge transfer analysis is only performed on the
orbitals within the coupled cluster region. Consequently, effects beyond the
four nearest water molecules can only be included through the diffuse basis
functions and at the Hartree–Fock level.

To ensure the comparability with the results by Reinholdt et al.,12 we have
included only four water molecules in the coupled cluster region. However, it
is natural to increase the number of solvent molecules in the coupled cluster
region to examine the potential improvement in the spectra. This is
particularly interesting in the case of ammonium, where the number of water
molecules in the first solvation shell is known to be 5.2.49 We have
calculated the spectra of aqueous ammonium and ammonia with six water
molecules in the coupled cluster region for a subset of 30 snapshots. The
results are presented in Figures S4 and S5 in the Supporting Information.
Although indicating that increasing the coupled cluster region could improve
the imbalance between the intensities of the main and post edges, no
significant changes in the spectra are observed. These results are in
agreement with the charge transfer analysis (see Figure S3), indicating that
the CC region is sufficiently large. However, a further improvement of the
electronic structure, achievable through a higher level of coupled cluster
theory, or improved quality of the structures can further reduce the
discrepancy between theory and experiment. The quality of the structures can
be, for example, improved by accounting for nuclear quantum effects in the
molecular dynamics simulations, which play an important role in capturing the
delicate nature of hydrogen bonds.50 From a recent study on the XA of liquid
water,20 we have observed both a redistribution of intensity from the main- to
the post-edge and a broadening effect from the inclusions of NQEs.

##  4 Conclusions

We have compared the performance of the different embedding schemes —PE-CCSD,
PDE-CCSD, CCSD-in-HF and MLCC3-in-HF—in simulating the XA spectra. Our
comparison reveals similarities between PDE-CCSD and CCSD-in-HF, while PE-CCSD
results differ noticeably. Given that PDE-CCSD represents a theoretical
refinement of PE-CCSD36, 51 and CC-in-HF methods are fully quantum mechanical,
we conclude that the close agreement of PE-CCSD with experiment is due to
error cancellation. Approaches that treat the environment through a density-
type embedding will avoid issues relating to the lack of Pauli-repulsion, as
is a concern with PE.

A charge transfer analysis demonstrates the local character of the excitations
responsible for the pre- and main-edge, in contrast to the more diffuse nature
of post-edge excitations in both solutions. Exploring the impact of an
increased number of solvent molecules in the coupled cluster region did not
provide any significant improvement. Therefore, we conclude that our methods
offer a robust description of the electronic structure of the solute, with
potential for improvement by using state-of-the-art path integral molecular
dynamics trajectories or through higher-level electronic structure models.

This work has received funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 Research and Innovation Program (grant
agreement No. 101020016). S.D.F., A.C.P, and H.K. acknowledge funding from the
Research Council of Norway through FRINATEK (project No. 275506). M.O.
acknowledges funding from the Swedish Research Council (grant agreement No.
2021-04521), and from the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie (grant agreement No.
860553). S.C. acknowledges support from the Independent Research Fund
Denmark–Natural Sciences (grant No. 7014-00258B). We acknowledge computing
resources through UNINETT Sigma2 –the National Infrastructure for High
Performance Computing and Data Storage in Norway (project No. NN2962k).

The Supporting Information contains an examination of the basis set
requirements on one arbitrarily selected snapshot, charge transfer analysis on
the solute and the four closest solvent molecules, XA spectra obtained using
various broadening schemes (Gaussian and Lorentzian), comparison of XA spectra
obtained when four vs. six water molecules are included in the coupled cluster
region. Geometries are available from https://doi.org/10.5281/zenodo.10390676.

##  S1 Basis Set Study

Figure S1:  Comparison of different basis sets for the XA spectrum of ammonia
in water at the MLCC3-12/72 level of theory using a single snapshot of the
dynamics simulation (step 2200). Individual excitations have been broadened
using Voigt profiles with $0.2\text{\,}\mathrm{eV}$ Lorentzian fwhm and
$0.2\text{\,}\mathrm{eV}$ Gaussian standard deviation (total fwhm $\sim
0.59\text{\,}\mathrm{eV}$).  Figure S2:  Comparison of different basis sets
for the XA spectrum of ammonium in water at the MLCC3-12/72 level of theory
using a single snapshot of the dynamics simulation (step 2200). Individual
excitations have been broadened using Voigt profiles with
$0.2\text{\,}\mathrm{eV}$ Lorentzian fwhm and $0.2\text{\,}\mathrm{eV}$
Gaussian standard deviation (total fwhm $\sim 0.59\text{\,}\mathrm{eV}$).

##  S2 Charge Transfer Analysis

Figure S3:  Spectrum plotted together with the charge transfer number for
charge transfer from the central molecule and the 4 closest water molecules
(ammonia (left) and ammonium (right)). The color map is identical to the one
in Figure 3 of the main paper.

##  S3 Increasing the number of solvent molecules in CC region

Figure S4: The effect of including 6 closest water molecules in CC-in-HF
calculations for NH3. Average over 30 snapshots. Voigt profiles with
Lorentzian fwhm of 0.2 eV and 0.2 Gaussian standard deviation. XA spectra of
ammonia in water clusters at the CCSD-in-HF and MLCC3-in-HF level of theory.
The spectra have been shifted by $-2\text{\,}\mathrm{eV}$ and
$-1\text{\,}\mathrm{eV}$, respectively.  Figure S5:  The effect of including 6
closest water molecules in CC-in-HF calculations for NH${}_{4}^{+}$. Average
over 30 snapshots. Voigt profiles with Lorentzian fwhm of 0.2 eV and 0.2
Gaussian standard deviation. XA spectra of ammonium in water clusters at the
CCSD-in-HF and MLCC3-in-HF level of theory. The spectra have been shifted by
$-2\text{\,}\mathrm{eV}$ and $-1\text{\,}\mathrm{eV}$, respectively.

##  S4 Convergence from CCSD to CC3

Figure S6:  Convergence from CCSD to CC3 by increasing the number of
occupied/virtual CNTOs included in the MLCC3 orbital space. The spectra are
obtained for a single snapshot of the dynamics simulation of ammonia (step
2200). Individual excitations have been broadened using Voigt profiles with
$0.2\text{\,}\mathrm{eV}$ Lorentzian fwhm and $0.2\text{\,}\mathrm{eV}$
Gaussian standard deviation (total fwhm $\sim 0.59\text{\,}\mathrm{eV}$).
Figure S7:  Convergence from CCSD to CC3 by increasing the number of
occupied/virtual CNTOs included in the MLCC3 orbital space. The spectra are
obtained for a single snapshot of the dynamics simulation of ammonium (step
2200). Individual excitations have been broadened using Voigt profiles with
$0.2\text{\,}\mathrm{eV}$ Lorentzian fwhm and $0.2\text{\,}\mathrm{eV}$
Gaussian standard deviation (total fwhm $\sim 0.59\text{\,}\mathrm{eV}$).

##  S5 XA spectra using various broadening schemes

Figure S8:  XA spectra of NH3 (top) and NH${}^{+}_{4}$ (bottom) in water
clusters at the CCSD-in-HF and MLCC3-in-HF levels of theory. aug-cc-pVTZ was
used for the solute, aug-cc-pVDZ for the four closest water molecules and cc-
pVDZ for the remaining solvent. Broadening: Lorentzian profiles with
$0.2\text{\,}\mathrm{eV}$ fwhm. Constant energy shifts of
$-2.0\text{\,}\mathrm{eV}$ and $-1.1\text{\,}\mathrm{eV}$ were applied for
CCSD-in-HF and MLCC3-in-HF respectively.  Figure S9:  XA spectra of NH3 (top)
and NH${}^{+}_{4}$ (bottom) in water clusters at the CCSD-in-HF and MLCC3-in-
HF levels of theory. aug-cc-pVTZ was used for the solute, aug-cc-pVDZ for the
four closest water molecules and cc-pVDZ for the remaining solvent.
Broadening: Gaussian profiles with $0.4\text{\,}\mathrm{eV}$ fwhm. Constant
energy shifts of $-2.0\text{\,}\mathrm{eV}$ and $-1.1\text{\,}\mathrm{eV}$
were applied for CCSD-in-HF and MLCC3-in-HF respectively.

## References

  * Ekimova et al. 2017 Ekimova, M.; Quevedo, W.; Szyc, Ł.; Iannuzzi, M.; Wernet, P.; Odelius, M.; Nibbering, E. T. Aqueous Solvation of Ammonia and Ammonium: Probing Hydrogen Bond Motifs with FT-IR and Soft X-Ray Spectroscopy. _J. Am. Chem. Soc._ 2017, _139_ , 12773–12783
  * Perrin and Gipe 1986 Perrin, C. L.; Gipe, R. K. Rotation, Solvation, and Hydrogen Bonding of Aqueous Ammonium Ion. _J. Am. Chem. Soc._ 1986, _108_ , 1088–1089
  * Ronan et al. 2021 Ronan, E.; Aqeel, H.; Wolfaardt, G. M.; Liss, S. N. Recent Advancements in the Biological Treatment of High Strength Ammonia Wastewater. _World J. Microbiol. Biotechnol._ 2021, _37_ , 158
  * Chen et al. 2012 Chen, H.; Dou, B.; Song, Y.; Xu, Y.; Wang, X.; Zhang, Y.; Du, X.; Wang, C.; Zhang, X.; Tan, C. Studies on Absorption and Regeneration for CO2 Capture by Aqueous Ammonia. _Int. J. Greenhouse Gas Control_ 2012, _6_ , 171–178
  * Feng et al. 2021 Feng, D.; Guo, D.; Zhang, Y.; Sun, S.; Zhao, Y.; Chang, G.; Guo, Q.; Qin, Y. Adsorption-Enrichment Characterization of CO2 and Dynamic Retention of Free NH3 in Functionalized Biochar with H2O/NH3· H2O Activation for Promotion of New Ammonia-Based Carbon Capture. _J. Chem. Eng._ 2021, _409_ , 128193
  * Pullman and Armbruster 1974 Pullman, A.; Armbruster, A. An Ab Initio Study of The Hydration and Ammoniation of Ammonium Ions. _Int. J. Quantum Chem._ 1974, _8_ , 169–176
  * Rode and Tanabe 1988 Rode, B. M.; Tanabe, Y. Simulation of Preferential Cation Solvation in Aqueous Ammonia. _J. Chem. Soc. Faraday Trans._ 1988, _84_ , 1779–1788
  * Brugé et al. 1999 Brugé, F.; Bernasconi, M.; Parrinello, M. Ab initio simulation of rotational dynamics of solvated ammonium ion in water. _J. Am. Chem. Soc._ 1999, _121_ , 10883–10888
  * Bacelo 2002 Bacelo, D. E. Theoretical study of microscopic solvation of ammonia in water clusters: NH3 (H2O)n, $n$ = 3, 4. _J. Phys. Chem. A_ 2002, _106_ , 11190–11196
  * Weinhardt et al. 2011 Weinhardt, L.; Weigand, M.; Fuchs, O.; Bär, M.; Blum, M.; Denlinger, J.; Yang, W.; Umbach, E.; Heske, C. Nuclear dynamics in the core-excited state of aqueous ammonia probed by resonant inelastic soft x-ray scattering. _Phys. Rev. B_ 2011, _84_ , 104202
  * Guo et al. 2020 Guo, J.; Zhou, L.; Zen, A.; Michaelides, A.; Wu, X.; Wang, E.; Xu, L.; Chen, J. Hydration of NH${}_{4}^{+}$ in Water: Bifurcated Hydrogen Bonding Structures and Fast Rotational Dynamics. _Phys. Rev. Lett._ 2020, _125_ , 106001
  * Reinholdt et al. 2021 Reinholdt, P.; Vidal, M. L.; Kongsted, J.; Iannuzzi, M.; Coriani, S.; Odelius, M. Nitrogen K-Edge X-ray Absorption Spectra of Ammonium and Ammonia in Water Solution: Assessing the Performance of Polarizable Embedding Coupled Cluster Methods. _J. Phys. Chem. Lett._ 2021, _12_ , 8865–8871
  * Carter-Fenk and Head-Gordon 2022 Carter-Fenk, K.; Head-Gordon, M. On the Choice of Reference Orbitals for Linear-Response Calculations of Solution-Phase K-Edge X-Ray Absorption spectra. _Phys. Chem. Chem. Phys._ 2022, _24_ , 26170–26179
  * Patkar et al. 2022 Patkar, D.; Bharati Ahirwar, M.; Deshmukh, M. M. A Tug of War between the Self-and Cross-associating Hydrogen Bonds in Neutral Ammonia-Water Clusters: Energetic Insights by Molecular Tailoring Approach. _ChemPhysChem_ 2022, _23_ , e202200476
  * Sanchez de Meras et al. 2010 Sanchez de Meras, A. M.; Koch, H.; Cuesta, I. G.; Boman, L. Cholesky decomposition-based definition of atomic subsystems in electronic structure calculations. _J. Chem. Phys._ 2010, _132_ , 204105
  * Manby et al. 2012 Manby, F. R.; Stella, M.; Goodpaster, J. D.; Miller III, T. F. A Simple, Exact Density-Functional-Theory Embedding Scheme. _J. Chem. Theory Comput._ 2012, _8_ , 2564–2568
  * Triguero et al. 1998 Triguero, L.; Pettersson, L.; Ågren, H. Calculations of near-edge x-ray-absorption spectra of gas-phase and chemisorbed molecules by means of density-functional and transition-potential theory. _Phys. Rev. B_ 1998, _58_ , 8097
  * Folkestad et al. 2021 Folkestad, S. D.; Kjønstad, E. F.; Goletto, L.; Koch, H. Multilevel CC2 and CCSD in Reduced Orbital Spaces: Electronic Excitations in Large Molecular Systems. _J. Chem. Theory Comput._ 2021, _17_ , 714–726
  * Folkestad et al. 2020 Folkestad, S. D.; Kjønstad, E. F.; Myhre, R. H.; Andersen, J. H.; Balbi, A.; Coriani, S.; Giovannini, T.; Goletto, L.; Haugland, T. S.; Hutcheson, A.; Høyvik, I.-M.; Moitra, T.; Paul, A. C.; Scavino, M.; Skeidsvoll, A. S.; Tveten, Å. H.; Koch, H. eT 1.0: An open source electronic structure program with emphasis on coupled cluster and multilevel methods. _J. Chem. Phys._ 2020, _152_ , 184103
  * Folkestad et al. 2023 Folkestad, S. D.; Paul, A. C.; Paul, R.; Coriani, S.; Odelius, M.; Iannuzzi, M.; Koch, H. Enhanced Understanding of X-Ray Absorption in Liquid Water Using Novel Coupled Cluster Methodologies. _arXiv preprint arXiv:2308.09495_ 2023,
  * Aidas et al. 2014 Aidas, K.; Angeli, C.; Bak, K. L.; Bakken, V.; Bast, R.; Boman, L.; Christiansen, O.; Cimiraglia, R.; Coriani, S.; Dahle, P.; Dalskov, E. K.; Ekström, U.; Enevoldsen, T.; Eriksen, J. J.; Ettenhuber, P.; Fernández, B.; Ferrighi, L.; Fliegl, H.; Frediani, L.; Hald, K.; Halkier, A.; Hättig, C.; Heiberg, H.; Helgaker, T.; Hennum, A. C.; Hettema, H.; Hjertenæs, E.; Høst, S.; Høyvik, I.-M.; Iozzi, M. F.; Jansík, B.; Jensen, H. J. Aa.; Jonsson, D.; Jørgensen, P.; Kauczor, J.; Kirpekar, S.; Kjærgaard, T.; Klopper, W.; Knecht, S.; Kobayashi, R.; Koch, H.; Kongsted, J.; Krapp, A.; Kristensen, K.; Ligabue, A.; Lutnæs, O. B.; Melo, J. I.; Mikkelsen, K. V.; Myhre, R. H.; Neiss, C.; Nielsen, C. B.; Norman, P.; Olsen, J.; Olsen, J. M. H.; Osted, A.; Packer, M. J.; Pawlowski, F.; Pedersen, T. B.; Provasi, P. F.; Reine, S.; Rinkevicius, Z.; Ruden, T. A.; Ruud, K.; Rybkin, V. V.; Sałek, P.; Samson, C. C. M.; de Merás, A. S.; Saue, T.; Sauer, S. P. A.; Schimmelpfennig, B.; Sneskov, K.; Steindal, A. H.; Sylvester-Hvid, K. O.; Taylor, P. R.; Teale, A. M.; Tellgren, E. I.; Tew, D. P.; Thorvaldsen, A. J.; Thøgersen, L.; Vahtras, O.; Watson, M. A.; Wilson, D. J. D.; Ziolkowski, M.; Ågren, H. The Dalton quantum chemistry program system. _WIREs Comput. Mol. Sci._ 2014, _4_ , 269–284
  * Sekino and Bartlett 1984 Sekino, H.; Bartlett, R. J. A linear response, coupled-cluster theory for excitation energy. _Int. J. Quantum Chem._ 1984, _26_ , 255–265
  * Geertsen et al. 1989 Geertsen, J.; Rittby, M.; Bartlett, R. J. The equation-of-motion coupled-cluster method: Excitation energies of Be and CO. _Chem. Phys. Lett._ 1989, _164_ , 57–62
  * Bartlett 2012 Bartlett, R. J. Coupled-cluster theory and its equation-of-motion extensions. _Wiley Interdiscip. Rev. Comput. Mol. Sci._ 2012, _2_ , 126–138
  * Stanton and Bartlett 1993 Stanton, J. F.; Bartlett, R. J. The equation of motion coupled-cluster method. A systematic biorthogonal approach to molecular excitation energies, transition probabilities, and excited state properties. _J. Chem. Phys._ 1993, _98_ , 7029–7039
  * Koch et al. 1994 Koch, H.; Kobayashi, R.; Sanchez de Merás, A.; Jørgensen, P. Calculation of size-intensive transition moments from the coupled cluster singles and doubles linear response function. _J. Chem. Phys._ 1994, _100_ , 4393–4400
  * Koch et al. 1997 Koch, H.; Christiansen, O.; Jørgensen, P.; Sanchez de Merás, A. M.; Helgaker, T. The CC3 model: An iterative coupled cluster approach including connected triples. _J. Chem. Phys._ 1997, _106_ , 1808–1818
  * Paul et al. 2020 Paul, A. C.; Myhre, R. H.; Koch, H. New and Efficient Implementation of CC3. _J. Chem. Theory Comput._ 2020, _17_ , 117–126
  * Myhre et al. 2014 Myhre, R. H.; Sanchez de Meras, A. M.; Koch, H. Multi-level coupled cluster theory. _J. Chem. Phys._ 2014, _141_ , 224105
  * Myhre and Koch 2016 Myhre, R. H.; Koch, H. The multilevel CC3 coupled cluster model. _J. Chem. Phys._ 2016, _145_ , 044111
  * Paul et al. 2022 Paul, A. C.; Folkestad, S. D.; Myhre, R. H.; Koch, H. Oscillator Strengths in the Framework of Equation of Motion Multilevel CC3. _J. Chem. Theory Comput._ 2022, _18_ , 5246–5258
  * Folkestad and Koch 2020 Folkestad, S. D.; Koch, H. Equation-of-Motion MLCCSD and CCSD-in-HF Oscillator Strengths and Their Application to Core Excitations. _J. Chem. Theory Comput._ 2020, _16_ , 6869–6879
  * Pulay 1983 Pulay, P. Localizability of dynamic electron correlation. _Chem. Phys. Lett._ 1983, _100_ , 151–154
  * Sæbø and Pulay 1993 Sæbø, S.; Pulay, P. Local treatment of electron correlation. _Annu. Rev. Phys. Chem._ 1993, _44_ , 213–236
  * Høyvik et al. 2017 Høyvik, I.-M.; Myhre, R. H.; Koch, H. Correlated natural transition orbitals for core excitation energies in multilevel coupled cluster models. _J. Chem. Phys._ 2017, _146_
  * Olsen et al. 2015 Olsen, J. M. H.; Steinmann, C.; Ruud, K.; Kongsted, J. Polarizable Density Embedding: A New QM/QM/MM-Based Computational Strategy. _J. Phys. Chem. A_ 2015, _119_ , 5344–5355
  * Olsen et al. 2010 Olsen, J. M.; Aidas, K.; Kongsted, J. Excited States in Solution through Polarizable Embedding. _J. Chem. Theory Comput._ 2010, _6_ , 3721–3734
  * Huzinaga and Cantu 1971 Huzinaga, S.; Cantu, A. A. Theory of separability of Many-electron systems. _J. Chem. Phys._ 1971, _55_ , 5543–5549
  * Richard and Herbert 2011 Richard, R. M.; Herbert, J. M. Time-Dependent Density-Functional Description of the 1La State in Polycyclic Aromatic Hydrocarbons: Charge-Transfer Character in Disguise? _J. Chem. Theory Comput._ 2011, _7_ , 1296–1306
  * Luzanov and Zhikol 2010 Luzanov, A. V.; Zhikol, O. A. Electron invariants and excited state structural analysis for electronic transitions within CIS, RPA, and TDDFT models. _Int. J. Quantum Chem._ 2010, _110_ , 902–924
  * Yanai et al. 2004 Yanai, T.; Tew, D. P.; Handy, N. C. A new hybrid exchange–correlation functional using the Coulomb-attenuating method (CAM-B3LYP). _Chem. Phys. Lett._ 2004, _393_ , 51–57
  * Gagliardi et al. 2004 Gagliardi, L.; Lindh, R.; Karlström, G. Local properties of quantum chemical systems: The LoProp approach. _J. Chem. Phys._ 2004, _121_ , 4494–4500
  * Cederbaum 1987 Cederbaum, L. S. Many-body theory of multiple core holes. _Phys. Rev. A_ 1987, _35_ , 622–631
  * Coriani and Koch 2015 Coriani, S.; Koch, H. Communication: X-ray absorption spectra and core-ionization potentials within a core-valence separated coupled cluster framework. _J. Chem. Phys._ 2015, _143_ , 181103
  * Coriani and Koch 2016 Coriani, S.; Koch, H. Erratum: “Communication: X-ray absorption spectra and core-ionization potentials within a core-valence separated coupled cluster framework” [J. Chem. Phys. 143, 181103 (2015)]. _J. Chem. Phys._ 2016, _145_ , 149901
  * Rohatgi 2022 Rohatgi, A. Webplotdigitizer: Version 4.6. 2022; https://automeris.io/WebPlotDigitizer, accessed November 15, 2023
  * Plasser and Lischka 2012 Plasser, F.; Lischka, H. Analysis of Excitonic and Charge Transfer Interactions from Quantum Chemical Calculations. _J. Chem. Theory Comput._ 2012, _8_ , 2777–2789
  * Plasser et al. 2014 Plasser, F.; Wormit, M.; Dreuw, A. New tools for the systematic analysis and visualization of electronic excitations. I. Formalism. _J. Chem. Phys._ 2014, _141_ , 024106
  * Aydin et al. 2020 Aydin, F.; Zhan, C.; Ritt, C.; Epsztein, R.; Elimelech, M.; Schwegler, E.; Pham, T. A. Similarities and differences between potassium and ammonium ions in liquid water: a first-principles study. _Phys. Chem. Chem. Phys._ 2020, _22_ , 2540–2548
  * Ceriotti et al. 2016 Ceriotti, M.; Fang, W.; Kusalik, P. G.; McKenzie, R. H.; Michaelides, A.; Morales, M. A.; Markland, T. E. Nuclear Quantum Effects in Water and Aqueous Systems: Experiment, Theory, and Current Challenges. _Chem. Rev._ 2016, _116_ , 7529–7550
  * Hrsǎk et al. 2018 Hrsǎk, D.; Olsen, J. M. H.; Kongsted, J. Polarizable Density Embedding Coupled Cluster Method. _J. Chem. Theory Comput._ 2018, _14_ , 1351–1360