# Polarizable embedding complex polarization propagator in four- and two-
component frameworks

Joel Creutzberg  [  Erik. D. Hedegård  [<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Explicit embedding methods combined with the complex polarization propagator
(CPP) enable modeling of spectroscopy for increasingly complex systems with a
high density of states. We present the first derivation and implementation of
the CPP in four- and exact two-component (X2C) polarizable embedding (PE)
frameworks. We denote the developed methods PE-4c-CPP and PE-X2C-CPP,
respectively. We illustrate the methods by estimating the solvent effect on
UV-vis and X-ray atomic absorption (XAS) spectra of [Rh(H2O)6]^3+ and
[Ir(H2O)6]^3+ immersed in aqueous solution. We moreover estimate solvent
effects on UV-vis spectra of a platinum complex that can be photo-chemically
activated (in water) to kill cancer cells. Our results clearly show that
inclusion of the environment is required: UV-vis and (to a lesser degree) XAS
spectra can become qualitatively different from vacuum calculations.
Comparison of PE-4c-CPP and PE-X2C-CPP methods shows that X2C essentially
reproduces the solvent effect obtained with the 4c methods.

###### keywords:

Polarizable embedding, complex polarization propagator, relativistic effects

Lund University] Division of Theoretical Chemistry, Lund University, SE- 223
62 Lund, Sweden University of Southern Denmark] Department of Physics,
Chemistry and Pharmacy, University of Southern Denmark, Campusvej 55, 5230
Odense, Denmark [Lund University] Division of Theoretical Chemistry, Lund
University, SE- 223 62 Lund, Sweden IR,NMR,UV

##  1 Introduction

Linear response functions1, 2 or polarization propagators constitute the
foundation for a large part of quantum chemical calculations on electronic
excitation spectra done today. The theory has been developed to encompass many
different wave function parameterizations, but its’ formulation within density
functional theory (DFT) is probably the most popular.3, 4, 5 While still
widely employed, the original formulations suffer from several shortcomings.
One shortcoming is the use of inherently non-relativistic frameworks, which is
insufficient for heavier elements6, particularly for spectroscopies involving
excitation of core-electrons7. Another shortcoming is the divergence for near-
resonant perturbations, leading to numerical instabilities in spectral regions
close to the electronic excitations.

The first shortcoming has been addressed by the development of relativistic
TD-DFT8, 9, 10 ranging from four-component (4c) Dirac–Kohn–Sham11, 12, 13 to
various two-component frameworks14, 15, 13, 16, 17. In parallel developments,
the second shortcoming has also been addressed: it has been shown that the
divergence can be removed by introducing a (phenomenological) finite lifetime
of the excited states.18, 19, 20 The term including the life-time is complex,
leading to complex and real parts of the resulting response functions; hence
the name complex polarization propagator (CPP). Depending on the perturbation
operator, the real and complex parts of the response functions are associated
with different physical (scattering and absorption) phenomena. For calculation
of electronic absorption, the CPP model exploits the imaginary part of the
polarizability tensor, which is proportional to the absorption cross-section.
Calculating the absorption cross-section at a given frequency allows targeting
any window of frequencies, which is another advantage of the CPP model:
conventional response theory employs a bottom-up approach, calculating lower
excitation energies first. This approach restricts conventional response
solvers to low-lying electronic excitations and regions with low density of
states. Accordingly, the CPP model allows modeling of all parts of the
spectrum on the same footings and is therefore often used for core-excitation
spectra.7

Seeing that relativistic effects are expected to be large for core-electrons,
it is not surprising that the CPP model has been extended to a relativistic
framework21, 22 where applications for both K-edge and L-edge spectra 23, 24
have been reported. Particular for the L-edge involving the 2p-orbitals, a
relativistic description is pertinent since spin-orbit coupling splits the
2p-orbitals into p${}_{\frac{1}{2}}$ and p${}_{\frac{3}{2}}$, giving rise to
L${}_{\text{II}}$ and L${}_{\text{III}}$ edges. Generally, the inclusion of
spin-orbit coupling will also entail the inclusion of a (non-degenerate)
triplet manifold, resulting in a high density of states21, 25, giving the CPP
model another advantage in relativistic calculations.

With the development of (non-)relativistic TD-DFT in resonant as well as for
non-resonant spectral regions, we seem to command a broad range of
computational methods for modeling spectroscopy. Still, spectroscopic methods
are used on ever more complex systems, often in the condensed phase, where
interactions with a (potentially heterogeneous) environment may become
important. Since large cluster models quickly become computationally
intractable, efficient models to describe interactions with solvent or other
environment molecules are required. One strategy is to include the environment
as a structure-less continuum26, 27. Another strategy has been to explicitly
include the environment by combining a QM method with fast molecular mechanics
force fields in a QM/MM hybrid.28, 29 We also mention another explicit
embedding model employing frozen QM densities for the environment30, 31, 32,
but our paper focus on classical treatment of the environment.

The explicit treatment of the environment has some advantages, particularly
for strong interactions (e.g. hydrogen bonds) or heterogeneous environments.
Most QM/MM models are non-relativistic and employ a rather simplistic point-
charge description for the environment, allowing only QM parts to be polarized
by the MM part, but not vice versa. These treatments may be sufficient in some
cases, e.g., for optimizing structures. However, molecular properties and
excitation energies (as we focus on here) are known to be more sensitive,
requiring an more advanced description of the environment electrostatics33.

The description of the electrostatic interactions between QM system and
environment can be improved by employing a polarizable environment.28
Development of such polarizable embedding (PE) models dates some years back34,
35 and is still ongoing.36, 37, 38, 39, 40, 41, 42, 33 A number of works have
appeared with PE models combined with linear response theory 43, 44, 45, 46,
47, 42 including TD-DFT.36, 48, 49, 50, 33 More recently, these models have
also included finite life- times51, 52, but only within a non-relativistic
framework. In fact, PE models were only fairly recently extended to
relativistic Hamiltonians.53, 54

In this paper, we continue our development of relativistic PE models for
spectroscopy: building on recent work with a relativistic linear response PE
model53 we here report the derivation and implementation of a new relativistic
PE-CPP model based on four- or exact-two-component (X2C) DFT. The environment
is described by multipole moments and site-polarizabilities49, 42 and we
denote our method PE-4c-CPP or PE-X2C-CPP. The methods are implemented in the
relativistic electronic structure program DIRAC9, 55 and we here use them to
calculate UV-vis and X-ray atomic absorption (XAS) spectra of [Rh(H2O)6]^3+
and [Ir(H2O)6]^3+ in vacuum and immersed in aqueous solution. The two systems
are mainly chosen as toy-models of heavy metals in aqueous environments, where
we see many future applications of the method. We also apply our newly
developed method on the platinum complex trans-trans-
trans-[Pt(N3)2(OH)2(NH3)2], which is a prototype complex for photo-activated
cancer therapy.56, 57. Therefore, its’ UV-vis spectrum has been studied
extensively, although only recently by relativistic methods25, 58 and without
explicit solvent effects.

##  2 Theory

In this section, we introduce the PE model and its’ extension to a CPP model
in a relativistic framework. We exclusively work within the Born–Oppenheimer
approximation in a second-quantization59 Dirac–Coulomb Kohn–Sham (DC–KS)
framework, using Hartree units. The method is also implemented with a
Hartree–Fock (HF) wave function. The derived expressions are generic for non-
relativistic, scalar relativistic (spin-free), X2C- and 4c-KS frameworks, so
differences between non-relativistic and relativistic frameworks will be
highlighted. We use $p,q,r,s$ indices for general orbitals, $i,j,k,l$ for
occupied orbitals and $a,b,c,d$ for unoccupied orbitals.

Relativistic PE models are rare in the literature and we accordingly have kept
some detail concerning this model. A derivation of the 4c-PE model (including
conventional linear response theory) was given in Ref. 53. We will often refer
to this or the to literature on vacuum CPP methods60, 21 to keep focus on the
novel parts of the relativistic PE-4c-CPP and PE-X2C-CPP models.

### The relativistic polarizable embedding method

The PE model works with an energy expression similar to QM/MM models

$\displaystyle
E=E_{\mathrm{QM}}+E_{\mathrm{es}}+E_{\mathrm{ind}}+E_{\mathrm{env}}.$ (1)

The terms $E_{\mathrm{QM}}$ and $E_{\mathrm{env}}$ describe interactions
exclusively within the QM and environment parts, respectively. Since
$E_{\mathrm{env}}$ does not depend on any wave function parameters, it will
not figure in the response equations in the next section and we will not
consider this term further here. The terms $E_{\mathrm{es}}$ and
$E_{\mathrm{ind}}$ describe the interaction between QM and environment: the
former is the electrostatic interactions and the latter the mutual
polarization of the environment and QM parts (this term is missing from most
QM/MM models). We discuss $E_{\mathrm{QM}}$, $E_{\mathrm{es}}$ and
$E_{\mathrm{ind}}$ from Eq. (1) in more detail below.

The QM energy in a 4c-KS framework based on a Dirac–Coulomb Hamiltonian is
given61, 62

$\displaystyle
E_{\mathrm{QM}}=\sum_{pq}h_{pq}D_{pq}+\sum_{pq}j_{pq}D_{pq}+E_{\mathrm{xc}}[\rho]+E_{\mathrm{nn}},$
(2)

where $D_{pq}=\langle 0|\hat{a}^{\dagger}_{p}\hat{a}_{q}|0\rangle$ is an
element of the one-electron reduced density matrix with creation
($\hat{a}^{\dagger}_{p}$) and annihilation ($\hat{a}_{q}$) operators working
on orbitals $p$ and $q$, respectively. The KS reference determinant is denoted
as $|0\rangle$. Further, $\rho$ is the electron density and $h_{pq}$ are
integrals over the one-electron operator, which for a 4c framework reads

$\displaystyle\hat{h}=\hat{h}_{D}+\hat{V}_{\mathrm{ext}}=\left(\begin{array}[]{cc}\bm{0}_{2}&c(\bm{\sigma}\cdot\hat{\mathbf{p}})\\\
c(\bm{\sigma}\cdot\hat{\mathbf{p}})&-2c^{2}\mathbf{I}_{2}\end{array}\right)+\left(\begin{array}[]{cc}V_{\mathrm{ext}}\mathbf{I}_{2}&\bm{0}_{2}\\\
\bm{0}_{2}&V_{\mathrm{ext}}\mathbf{I}_{2}\end{array}\right).$ (7)

For the remaining terms in Eq. (2) and the definition of constants and
operators in Eq. (7), we refer to the definitions in Refs. 62, 61. We define
$E_{\mathrm{xc}}[\rho]$ including exact exchange as in Ref. 61. An analytic
solution can be obtained for the one-electron part of Eq. (2) and we denote
this wave function

$\displaystyle|\psi^{4c}\rangle=\left(\begin{array}[]{c}|\psi^{L}\rangle\\\
|\psi^{S}\rangle\end{array}\right),$ (10)

where $|\psi_{L}\rangle$ and $|\psi_{S}\rangle$ are the large- and small-
component wave functions, respectively. These are themselves two-component
wave functions, yielding the 4c form of $|\psi^{4c}\rangle$ as required by the
structure of Eq. (7). We will use this nomenclature when we discuss X2C
methods below, but for now we focus on the 4c-KS method.

The optimal energy in Eq. (2) is obtained through diagonalization of the KS-
matrix, $\mathbf{F}$, which is equivalent to requiring the off-diagonal
elements to vanish

$\displaystyle F_{ia}=\langle i|\hat{f}_{0}|a\rangle=0.$ (11)

We have here defined the KS-operator in second-quantization

$\displaystyle\hat{f}_{0}=\sum_{pq}f_{pq}\hat{a}^{\dagger}_{p}a_{q}.$ (12)

Expressions for the integrals in $f_{pq}$ can be found in the literature (see,
e.g., Ref. 61); they contain the one electron part, i.e., Eq. (7), and the
Coulomb and exchange–correlation terms. Differences to the non-relativistic
theory manifest in the $4\times 4$ matrix form of the operator in Eq. (7); the
KS reference state, $|0\rangle$, is accordingly on a similar form as Eq. (10).
In a finite-orbital basis, the individual orbitals can be written

$\displaystyle|p\rangle=\sum^{N_{L}}_{\kappa}\left(\begin{array}[]{c}|\chi^{L}_{\kappa}\rangle\\\
0\end{array}\right)C^{L}_{\kappa}+\sum^{N_{s}}_{\lambda}\left(\begin{array}[]{c}0\\\
|\chi^{S}_{\lambda}\rangle\end{array}\right)C^{S}_{\lambda}$ ;
$\displaystyle|\chi^{X}_{i}\rangle=\left(\begin{array}[]{c}|\chi^{X}_{i\alpha}\rangle\\\
|\chi^{X}_{i\beta}\rangle\end{array}\right),$ (19)

where the primitive basis functions $|\chi^{L}_{\kappa}\rangle$ and
$|\chi^{S}_{\lambda}\rangle$, similar to Eq. (10), are two-component spinors
$|\chi^{X}_{i}\rangle$ ($X=L,S$) with $\alpha$\- and $\beta$ parts. The
expansion coefficients $C^{L}_{\kappa}$ and $C^{S}_{\lambda}$ for the large
and small component spinors are in the DIRAC program partitioned according to
a quaternion symmetry scheme, reducing to complex or real algebra, depending
on the spatial symmetry group. However, since the implementation of the PE
model currently does not exploit spatial symmetry, we will always use full
quaternion symmetry.

From Eq. (19) it follows that the matrix representations of the one-electron
Hamiltonian and the KS matrix in Eq. (11) have the following structures

$\displaystyle\mathbf{h}_{D}=\left(\begin{array}[]{cc}\mathbf{h}^{LL}_{D}&\mathbf{h}^{LS}_{D}\\\
\mathbf{h}^{SL}_{D}&\mathbf{h}^{SS}_{D}\end{array}\right)$
$\displaystyle\mathbf{F}=\left(\begin{array}[]{cc}\mathbf{F}^{LL}&\mathbf{F}^{LS}\\\
\mathbf{F}^{SL}&\mathbf{F}^{SS}\end{array}\right),$ (24)

where each element is a $2\times 2$ block with the generic matrix elements
$h^{XY}_{D,\kappa\lambda}=\langle\chi^{X}_{\kappa}|\hat{h}_{D}|\chi^{Y}_{\lambda}\rangle$
and
$F^{XY}_{\kappa\lambda}=\langle\chi^{X}_{\kappa}|\hat{f}^{\mathrm{KS}}_{0}|\chi^{Y}_{\lambda}\rangle$.

We now turn to the PE operators giving rise to the two terms $E_{\mathrm{es}}$
and $E_{\mathrm{ind}}$ in Eq. (1). These operators have been described
earlier63, 45, 51, 42 with few exceptions64, 53 within non-relativistic
frameworks. The electrostatic energy contribution is obtained as the
expectation value of $\hat{V}^{\mathrm{es}}$

$\displaystyle E_{\mathrm{es}}=\langle 0|\hat{V}^{\mathrm{es}}|0\rangle,$ (25)

where the operator $\hat{V}$ is defined

$\displaystyle\hat{V}^{\mathrm{es}}=\sum_{s}\sum_{|k|=0}\frac{(-1)^{|k|}}{k!}\Bigr{(}\mathbf{M}_{s}^{(k)}\sum_{pq}\mathbf{t}^{(k)}_{pq,s}\hat{a}^{\dagger}_{p}\hat{a}_{q}+\mathbf{M}_{s}^{(k)}\sum_{M}Z_{M}\mathbf{T}_{Ms}^{(k)}\Bigr{)}.$
(26)

We have in accordance with previous literature63, 65, 45, 51, 53, 42 used a
three-dimensional multi-index $k=(k_{x},k_{y},k_{z})$ for multipoles
$\mathbf{M}^{(k)}_{s}$ located at sites $s$; thus $k=0$ corresponds to point
charges, $\mathbf{M}^{(0)}_{s}=q_{s}$, $k=1$ corresponds to components of the
dipole vector,
$\mathbf{M}^{(1)}_{s}=\bm{\mu}_{s}=(\mu_{s,x},\mu_{s,y},\mu_{s,z})$, and so
on. This notation is also used in the definition of the integrals
$\mathbf{t}^{(k)}_{pq,s}=-\langle p|\mathbf{T}^{(k)}_{s}|q\rangle$ over
interaction operators, $\mathbf{T}^{(k)}_{j}$, defined as in the monograph by
Stones.66 We construct the operator in Eq. (26) employing the library PELib67
which previously has been interfaced with DIRAC.53 The multipole moments are
obtained as atom-localized multipole moments from the LoProp method68. Since
this requires a separate QM calculation, the computational cost is kept low by
dividing the environment into fragments. Multipole moments are in this way
calculated for each fragment.

The term accounting for the mutual polarization of QM and environment regions
in Eq. (1) is given

$E_{\mathrm{ind}}=-\frac{1}{2}(\bm{\mu}^{\mathrm{ind}})^{T}\langle
0|\hat{\bm{\mathcal{E}}}|0\rangle,$ (27)

where the field vector operator, $\hat{\bm{\mathcal{E}}}$, is arranged as a
vector where each entry contains the field operator for a polarizable site,
$\hat{\bm{\mathcal{E}}}=(\hat{\bm{\mathcal{E}}}_{1},\ldots,\hat{\bm{\mathcal{E}}}_{s},\ldots,\hat{\bm{\mathcal{E}}}_{S})$.
The induced dipole moment is a vector arranged similarly. The field on site
$s$ has contributions from QM electrons and nuclei as well as the multipoles,
leaving the total field operator for site $s$ as

$\displaystyle\hat{\bm{\mathcal{E}}}_{s}=\hat{\bm{\mathcal{E}}^{\mathrm{e}}}_{s}+\bm{\mathcal{E}}^{\mathrm{nuc}}_{s}+\bm{\mathcal{E}}^{\mathrm{es}}_{s}.$
(28)

The electronic and nuclear fields in Eq. (28) are given as

$\displaystyle\hat{\bm{\mathcal{E}}}^{e}_{s}$
$\displaystyle=\sum_{pq}\mathbf{t}^{(1)}_{pq,s}\hat{a}^{\dagger}_{p}\hat{a}_{q}=-\sum_{pq}\bigr{\langle}p|\frac{\mathbf{r}-\mathbf{r}_{s}}{|\mathbf{r}-\mathbf{r}_{s}|^{3}}|q\bigl{\rangle}\hat{a}^{\dagger}_{p}\hat{a}_{q}$
(29) $\displaystyle\bm{\mathcal{E}}_{s}^{n}$
$\displaystyle=-\sum_{M}Z_{M}\mathbf{T}^{(1)}_{Ms}=\sum_{M}\frac{Z_{M}(\mathbf{R}_{M}-\mathbf{r}_{s})}{|\mathbf{R}_{M}-\mathbf{r}_{s}|^{3}}.$
(30)

The individual induced dipole moments $\bm{\mu}^{\mathrm{ind}}$ in Eq. (27)
are obtained by self-consistently solving the matrix equation69
$\bm{\mu}^{\mathrm{ind}}=\mathbf{R}^{Relay}\bm{\mathcal{E}}$ in each self-
consistent-field iteration. The matrix $\mathbf{R}^{Relay}$ is given
explicitly in previous publications.63, 65, 45, 51, 53, 42 Here it is
sufficient to note that it does not contain any wave function parameters and
the diagonal contains polarizabilities (isotropic or an-isotropic) for each
polarizable site. The polarizabilities are obtained with the LoProp method68,
employing the same fragments as for the multipole moments.

From the expressions $E_{\mathrm{es}}$ and $E_{\mathrm{ind}}$ in Eqs. (25) and
(27), we have previously shown53 that the KS-operator due to a polarizable
environment, is modified with the PE-operator, $\hat{v}^{\mathrm{PE}}$:

$\displaystyle\hat{f}^{\mathrm{tot}}$
$\displaystyle=\hat{f}_{0}+\hat{v}^{\mathrm{PE}}$ (31)
$\displaystyle\hat{v}^{\mathrm{PE}}$
$\displaystyle=\hat{V}^{\mathrm{es}}-\bm{\mu}^{\mathrm{ind}}\hat{\bm{\mathcal{E}}}^{\mathrm{e}}=\hat{V}^{\mathrm{es}}-\langle
0|\hat{\bm{\mathcal{E}}}|0\rangle^{T}\,\mathbf{R}^{Relay}\,\hat{\bm{\mathcal{E}}}^{\mathrm{e}}.$
(32)

Contrary to a non-relativistic framework, the $\hat{v}^{\mathrm{PE}}$ operator
will attain a $4\times 4$ matrix structure. The only non-zero matrix elements
of $\hat{v}^{\mathrm{PE}}$ over four-component spinors (Eq. 19) are within the
large-large or small-small blocks.53 Thus, only $\mathbf{F}^{LL}$ and
$\mathbf{F}^{SS}$ blocks are modified in matrix representation of the total KS
operator

$\displaystyle\mathbf{F}^{\mathrm{KS}}=\left(\begin{array}[]{cc}\mathbf{F}^{LL}+\mathbf{v}^{\mathrm{PE},LL}&\mathbf{F}^{LS}\\\
\mathbf{F}^{SL}&\mathbf{F}^{SS}+\mathbf{v}^{\mathrm{PE},SS}\end{array}\right).$
(35)

where $\mathbf{v}^{\mathrm{PE},XY}$ contains the elements,
$\hat{v}^{\text{PE},XY}_{\kappa\lambda}=\langle\chi^{X}_{\kappa}|\hat{v}^{\text{PE}}|\chi^{Y}_{\lambda}\rangle$.

In addition to PE-4c-KS (and HF), our implementation also covers an X2C
Hamiltonian. In X2C methods, the large- and small-component equations
resulting from the form of $|\psi^{4c}\rangle$ in Eq. (10) are decoupled.
Formally, this is done by the decoupling matrix $\mathbf{U}$

$\displaystyle\mathbf{U}^{\dagger}\mathbf{h}_{D}\mathbf{U}=\left(\begin{array}[]{cc}\mathbf{h}_{++}&0\\\
0&\mathbf{h}_{--}\end{array}\right)$
$\displaystyle\mathbf{U}=\left(\begin{array}[]{cc}\mathbf{\Omega}_{+}&-\mathbf{R}^{\dagger}\mathbf{\Omega}_{-}\\\
\mathbf{R}\mathbf{\Omega}_{+}&\mathbf{\Omega}_{-}\end{array}\right),$ (40)

where
$\mathbf{\Omega}_{+}=\left(1+\mathbf{R}^{\dagger}\mathbf{R}\right)^{-\frac{1}{2}}$,
$\mathbf{\Omega}_{-}=\left(1+\mathbf{R}\mathbf{R}^{\dagger}\right)^{-\frac{1}{2}}$,
and $\mathbf{R}$ describes the coupling of the large and small component wave
functions70. The decoupling in Eq. (40) allows us to focus only on the
positive energy solutions

$\displaystyle\mathbf{h}_{++}|\psi^{2c}_{+}\rangle=E|\psi^{2c}_{+}\rangle,$
(41)

where
$|\psi^{2c}_{+}\rangle=\left(1+\mathbf{R}^{\dagger}\mathbf{R}\right)^{\frac{1}{2}}|\psi^{L}\rangle$
can be expanded in the large-component basis functions, thereby reducing the
computational cost significantly. For HF or DFT methods, a transformation of
the two-electron integrals in the the KS matrix (Eq. 11) is in principle
required, but this transformation would make the calculation more expensive
than the corresponding 4c calculation. In the PE-X2C implementation, we follow
the general procedure in DIRAC where $\mathbf{U}$ is constructed from the one-
electron Hamiltonian and two-electron integrals from the large-component basis
are used unmodified70. Leaving the two-electron part untransformed corresponds
to neglecting two-electron spin-spin and spin-orbit contributions, but these
can be very accurately approximated by atomic mean-field corrections71 (see
the Computational Details Section).

The extraction of $\mathbf{U}$ and construction of the two-component (X2C)
wave function in DIRAC occur prior to the (self-consistent) calculations of
the PE contributions in Eq. (32). Hence, the PE contributions in PE-X2C are
always constructed from the one-electron reduced density matrix corresponding
to the X2C reference state. We verify that this works well by comparison of
PE-X2C-CPP and PE-4c-CPP in the Results and Discussion Section. An alternative
route is to use a so-called mean-field X2C strategy72, i.e., solve the PE-4c
KS (or HF) equations and perform the transformation on the final KS (or Fock)
matrix. We leave an investigation of this route for future studies.

We will in the next section provide the theory behind the PE-CPP method, but
without distinguishing between 4c and X2C frameworks, since their equation
structures are equivalent. However, we note that the transformation in Eq.
(40) should be carried over to the property operators to avoid so-called
picture change errors 73, 74, 75. Following the general implementation in the
DIRAC program, all property operators are always transformed.

### The relativistic polarizable embedding complex polarization propagator

Introducing $\hat{v}^{\text{PE}}$ in Eq. (31) will also lead to modified
response equations due to the polarizable environment 36, 48, 49, 50, 33, 51,
52. We now examine these changes within a relativistic PE model that includes
a finite life time.

In response theory we add a time-dependent interaction operator, $\hat{V}(t)$,
as perturbation to the Hamiltonian; the operator defines a periodic, time-
dependent field

$\hat{V}(t)=\sum_{k}e^{-i\omega_{k}t}\sum_{X}\mathcal{E}_{X}(\omega_{k})\hat{X},$
(42)

where $\hat{X}$ is a property operator, $\mathcal{E}_{X}(\omega_{k})$ are
field-strength parameters (defined as Fourier components) with frequencies
$\omega_{k}$. We follow the nomenclature from Ref. 53, labeling fields from
perturbation operators with frequencies, $\omega_{k}$, while components of the
fields from QM region and multipoles at the polarizable sites in the
environment are labeled as in the previous section (cf. Eq. 28).

We will here work with a time-dependent KS reference determinant,
parameterized as

$|\tilde{0}(t)\rangle=e^{-\hat{\kappa}(t)}|0\rangle,$ (43)

where $\kappa(t)$ is the time-dependent orbital-rotation operator. Since we
are here only interested in linear response theory, we keep only linear
parameters in the perturbation expansion $\kappa(t)=\kappa^{(1)}(t)+\cdots$
which we define

$\hat{\kappa}^{(1)}(t)=\sum_{k}\sum_{ai}\Bigr{(}\kappa_{ai}(\omega_{k})\hat{q}^{\dagger}_{ai}-\kappa^{*}_{ai}(-\omega_{k})\hat{q}_{ai}\Bigl{)}e^{-i\omega_{k}t}.$
(44)

The parameters
$\kappa_{ai}(\omega_{k})=\sum_{X}\mathcal{E}_{X}(\omega_{k})\kappa^{\omega_{k}}_{ai}$
are linear in the field-strength and we have defined
$\hat{q}_{ai}=\hat{a}^{\dagger}_{i}\hat{a}_{a}$. Using the quasi-energy
formulation of response theory76, a time-dependent variational principle is
established through the time-averaged quasi-energy

$\\{Q(t)\\}_{T}=\frac{1}{T}\int_{-T/2}^{T/2}\langle\tilde{0}(t)|\big{(}\hat{H}-i\frac{\partial}{\partial
t}\big{)}|\tilde{0}(t)\rangle dt,$ (45)

from which response functions and response equations can be derived76, 2. The
Hamiltonian is given $\hat{H}=\hat{H}^{\mathrm{tot}}+\hat{V}(t)$: the first
part, $\hat{H}^{\mathrm{tot}}$, contains the usual components of the KS-DFT
framework77 corresponding to $E_{\mathrm{QM}}$ in Eq. (2) and additionally
operators corresponding to $E_{\mathrm{es}}$ and $E_{\mathrm{ind}}$ (Eqs.
25–27)53. The second part of $\hat{H}^{\mathrm{tot}}$ is the external
pertubating field in Eq. (42). By inserting $\hat{H}$ in Eq. (45) we can write
the quasi energy77

$\displaystyle\\{Q_{\mathrm{tot}}(t)\\}_{T}=\\{\bar{Q}(t)\\}_{T}+\sum_{k}\sum_{X}\mathcal{E}_{X}(\omega_{k})\\{\bar{Q}_{X}(t)\\}_{T},$
(46)

where $\\{\bar{Q}(t)\\}_{T}$ has origin in the field-independent part of
$\hat{H}$ (and also includes the time-derivative of Eq. 45) and
$\\{\bar{Q}_{X}(t)\\}$ contains the time-dependent part with the property
operator, $\hat{X}$. Following Ref. 53, we can divide both terms according to

$\displaystyle\\{\bar{Q}(t)\\}_{T}=Q_{0}+\\{Q^{\mathrm{PE}}(t)\\}_{T}$
$\displaystyle\\{\bar{Q}_{X}(t)\\}_{T}=\\{Q_{X}(t)\\}_{T}+\\{Q^{\mathrm{PE}}_{X}(t)\\}_{T},$
(47)

where $\\{Q^{\mathrm{PE}}(t)\\}_{T}$ and $\\{Q^{\mathrm{PE}}_{X}(t)\\}_{T}$
are the additional terms, arising due to the (polarizable) environment (see
Ref. 53 for explicit expressions). Note also that the field-dependent term in
Eq. (47) has been augmented with a contribution due to the polarizable
environment; this contribution does not have origin in
$\hat{H}^{\mathrm{tot}}$ and was not included in initial work with the (non-
relativistic) PE-model.39, 49, 44, 47 The contribution arise as the perturbing
field $\hat{V}(t)$ also induces a field, $\tilde{\bm{\mathcal{E}}}(t)$, within
the environment, modifying the total induced field in Eq. (28) to
$\tilde{\bm{\mathcal{E}}}^{\mathrm{tot}}_{s}=\tilde{\bm{\mathcal{E}}}^{\mathrm{e}}_{s}+\bm{\mathcal{E}}^{\mathrm{nuc}}_{s}+\bm{\mathcal{E}}^{\mathrm{es}}_{s}+\tilde{\bm{\mathcal{E}}}_{s}(t)$
within the linear response framework. This effect is denoted the effective
external field (EEF) contribution.78, 79

Response equations can now be obtained through derivatives of the quasi-
energy. In this last part of the derivation, our paper differs from the
previous relativistic PE models.53, 64 It has been shown for exact-state
response theory80 that employing a wave function without reference to the
excited state life-time (as the one in Eq. 43), corresponds to wave function
whose norm is constant in time. This is a somewhat absurd situation since the
corresponding excited states thus have infinite lifetime.80, 2 A physically
more pleasing expression including an excited state life-time can be obtained
by employing a phenomenological, common life-time for all excited states
($\gamma$) in the definition of the wave function. As shown (for an exact-
state) in Refs. 80, 2 this leads identical response equations, but with state-
energy differences $\omega_{k}$ replaced by the complex entity
$z=\omega_{k}+i\gamma$.18, 60, 80, 21, 20 By analogy, we arrive at the PE-CPP
equations by inserting the phase-isolated KS determinant in Eq. (43) with the
orbital rotation operator in Eq. (44) into the quasi energy, Eqs. (46)–(47),
while replacing $\omega\rightarrow z$. This leads to the following response
equations

$\displaystyle\bm{\kappa}^{z}=-\Bigl{(}\bar{\mathbf{E}}^{[2]}-z\mathbf{S}^{[2]}\Bigr{)}^{-1}\bar{\mathbf{E}}^{[1]}_{X}.$
(48)

The bars over the electronic Hessian, $\mathbf{\bar{E}}^{[2]}$, as well as the
property gradient vector, $\mathbf{\bar{E}}^{[1]}_{X}$, signify that they both
contain additional contributions due to the polarizable environment. For the
Hessian, we obtain

$\displaystyle\bar{\mathbf{E}}^{[2]}=\left(\begin{array}[]{cc}\mathbf{A}&\mathbf{B}\\\
\mathbf{B}^{*}&\mathbf{A}^{*}\end{array}\right)$ (51)

with

$\displaystyle A_{ai;bj}$ $\displaystyle=\langle
0|[\hat{q}^{\dagger}_{ai},[\hat{q}_{bj},\hat{f}^{\mathrm{KS}}_{0}]]|0\rangle+\langle
0|[\hat{q}^{\dagger}_{ai},\hat{v}^{\mathrm{Hxc}}]|0\rangle$
$\displaystyle+\langle
0|[\hat{q}^{\dagger}_{ai},[\hat{q}_{bj},\hat{v}^{\mathrm{PE}}]]|0\rangle+\langle
0|[\hat{q}^{\dagger}_{ai},\hat{v}^{\mathrm{ind}}]|0\rangle$ (52)
$\displaystyle B_{ai;bj}$ $\displaystyle=\langle
0|[\hat{q}_{ai},[\hat{q}_{bj},\hat{f}^{\mathrm{KS}}_{0}]]|0\rangle+\langle
0|[\hat{q}_{ai},\hat{v}^{\mathrm{Hxc}}]|0\rangle$ $\displaystyle+\langle
0|[\hat{q}_{ai},[\hat{q}_{bj},\hat{v}^{\mathrm{PE}}]]|0\rangle+\langle
0|[\hat{q}_{ai},\hat{v}^{\mathrm{ind}}]|0\rangle.$ (53)

The two first terms in Eqs. (52)–(53) are the vacuum contributions, whereas
the following two terms are the additional terms originating from the PE part.
The first term contains the KS operator in Eq. (12) while the second term with
the $\hat{v}^{\text{Hxc}}$ operator is the Hartree-exchange-correlation kernel
(note that this term was left out in Eq. 58 of Ref. 53 by mistake). In the
last two terms, we identify the PE operator in Eq. (32) and we have
additionally defined

$\displaystyle\hat{v}^{\mathrm{ind}}$ $\displaystyle=-\langle
0|[\hat{q}_{bj},\hat{\bm{\mathcal{E}}}^{\text{e}}]|0\rangle^{T}\,\mathbf{R}^{Relay}\,\hat{\bm{\mathcal{E}}}^{\mathrm{e}}.$
(54)

Meanwhile, the elements of the property gradient becomes

$\displaystyle\bar{E}^{[1]}_{ai,X}=\langle
0|[\hat{q}_{ai},\hat{X}]|0\rangle+\frac{\text{d}\bm{\mathbf{\mu}}^{\mathrm{ind}}_{\text{ext},X}(\omega_{k})}{\text{d}\bm{\mathcal{E}}_{X}(\omega_{k})}\langle
0|[\hat{q}_{ai},\hat{\bm{\mathcal{E}}}^{\mathrm{e}}]|0\rangle,$ (55)

where the last term, containing
$\bm{\mu}^{\mathrm{ind}}_{\text{ext},X}(\omega_{k})=\mathbf{R}^{Relay}\bm{\mathcal{E}}_{X}(\omega_{k})$,
gives rise to the EEF effect.

Explicit construction of the Hessian in Eq. (48) is computationally expensive,
but is avoided by expanding the solution vector
$\bm{\kappa}^{z}=\sum_{i=1}^{n}\mathbf{b}_{i}a_{i}(z)$ in a set of trial
vectors, $\\{\mathbf{b}_{i}\\}$, meaning that only

$\displaystyle\bm{\sigma}_{i}=\bar{\mathbf{E}}^{[2]}\mathbf{b}_{i},$ (56)

is required. Several algorithms for constructing these trial vectors have been
presented (see, e.g. Ref. 81 and reference therein). We rely on the scheme
developed in DIRAC, where the trial vectors are divided according to their
Hermitian and time-reversal symmetries.82, 13, 21 This scheme has also been
extended to response equations with complex frequencies21 but the addition of
the PE terms will not lead to any changes. Thus, the only modifications of the
vacuum equations are due to the modifications in $\bar{\mathbf{E}}^{[2]}$ and
$\bar{\mathbf{E}}^{[1]}_{X}$, leaving $\bm{\sigma}_{i}$ to be

$\displaystyle\sigma_{ai}$ $\displaystyle=\langle
0|[-\hat{q}_{ai},\tilde{f}_{0}]|0\rangle+\langle
0|[-\hat{q}_{ai},\tilde{v}^{\text{Hxc}}]|0\rangle$ $\displaystyle+\langle
0|[-\hat{q}_{ai},\tilde{v}^{\text{PE}}]|0\rangle+\langle
0|\tilde{\bm{\mathcal{E}}}^{e}|0\rangle^{T}\,\mathbf{R}^{Relay}\,\langle
0|[-q_{ai},\hat{\bm{\mathcal{E}}}^{e}]|0\rangle$
$\displaystyle=\tilde{f}^{\text{KS}}_{0,ai}+\tilde{v}^{\text{Hxc}}_{ai}+\tilde{v}^{\text{PE}}_{ai}+(\tilde{\bm{\mu}}^{\text{e}})^{T}\cdot\bm{\mathcal{E}}^{\text{e}}_{ai}.$
(57)

Note that we have used the nomenclature $\tilde{O}$ for one-index transformed
operators.1

##  3 Computational Details

Structure optimizations of [Rh(H2O)6]^3+ and [Ir(H2O)6]^3+: For the solvated
structures, we employed the Maestro builder83 to construct the metal ions,
immersed in a cubic solvent box (25 $\times$ 25 $\times$ 25 Å) of explicit
water molecules, modeled with TIP3P84. The systems were optimized with QM/MM
using QSite85, 86, 83. The QM regions were comprised of the [Rh(H2O)6]^3+ and
[Ir(H2O)6]^3+ ions, for which we employed the B3LYP functional87, 88, 89
combined with the LACVP∗ basis set.90, 91, 92, 93 The MM region was relaxed
with the truncated Newton method. As an example, we show the optimized
structures of [Ir(H2O)6]^3+ in Figure 1. Vacuum structures where optimized
with same same functional and basis set and resulted in symmetric structures
($C_{i}$ symmetry, although note that this was not used in the calculations).

Figure 1: The solvated sytems used in the PE calcuations for [Rh(H2O)6]^3+ and
[Ir(H2O)6]^3+ (shown with and without solvent, using [Ir(H2O)6]^3+ as
example). The solvated trans-trans-trans-Pt(N3)2(OH)2(NH3) (also shown without
solvent below [Ir(H2O)6]^3+) is treated similar.

Structure of trans-trans-trans-[Pt(N3)(OH)2(NH3)2] in water: Optimization of
the platinum complex (denoted trans-Pt from here) was more challenging than
the simple aqua ions. We employed a similar strategy as for the aqua ions,
although the box-size of explicit water molecules was larger (35 $\times$ 35
$\times$ 35 Å).

First, the structure of the trans-Pt complex was optimized with frozen
surroundings using the QM/MM method implemented in QSite85, 86, 83 (and same
functional/basis set as for the aqua ions). From this system, a sphere of 6 Å
was optimized with the BP8687, 94 functional and the def2-SV(P) basis set95,
96, 97 in ORCA 98 (water molecules within 4 Å of the trans-Pt complex was
allowed to structurally relax; the rest was kept frozen). For computational
efficiency, a pre-optimization with PBEh-3c99 method was done. After this, the
6 Å sphere was re-inserted into the original box and a 10 Å sphere was cut
out. All PE calculations were performed in this 10 Å system.

Polarizable embedding potentials: Atom-centered multipoles (charges, dipoles,
quadrupole moments) and static anisotropic polarizabilities were calculated
with the LoProp68 method in Molcas100 on individual fragments (each water
molecule constitutes a fragment). The setup of the calculations was done using
the PEAS101 script. The calculations employed the B3LYP functional and a
A-6-31PGP basis set (note that this basis set is specially constructed to fit
the LoProp method).

Complex polarization propagator calculations: Calculation of frequency-
dependent electric dipole polarizabilities were carried out on [Rh(H2O)6]^3+,
[Ir(H2O)6]^3+ and the trans-Pt complex with CPP for three cases, namely 1) the
vacuum optimized structures (denoted vacuum), 2) the structures immersed in a
solvent described by a potential (denoted PE), and 3) as case 2), but with the
addition of an effective external field (denoted PE-EEF). We additionally
carried out a set of calculations without a PE potential, using the QM/MM
optimized structures (i.e. they are optimized within solvent). This will allow
us to separate structural and electronic solvent effects. The calculation were
done in the frequency ranges corresponding to the UV-vis and various XAS
regions, using CPP and PE-CPP in a development version of DIRAC.9, 55 For the
UV-vis spectra the CPP calculations were carried out in the ranges of 1.7–8.8
eV for [Rh(H2O)6]^3+, 1.8–8.8 eV for [Ir(H2O)6]^3+ and 2.0-6.8 eV for trans-Pt
with steps of 0.068 eV (0.0025 Hartree). In order to obtain more data points
for the UV-vis spectra of [Rh(H2O)6]^3+ and [Ir(H2O)6]^3+, a second scan was
performed with start and end 0.3 eV before the original scans. The core
spectra for the LII/LIII edges were calculated in the ranges of 3104.0– 3124.2
eV and 2962.0–2975.5 eV for [Rh(H2O)6]^3+, and 12769.9–12784.0 eV and
11136.9–11147.3 eV for [Ir(H2O)6]^3+. Here, the frequency windows used for
each calculation were 0.61 eV (0.0225 Hartree), with steps of 0.068 eV (0.0025
Hartree). For the XAS calculations, channel restrictions102, 103, 104, 23 were
imposed, where only transitions from the 2p orbitals (p${}_{\frac{1}{2}}$ for
the LII-edge and p${}_{\frac{3}{2}}$ for the LIII-edge) were allowed within
the frequency windows. From the resulting polarizabilities,
$\alpha_{\alpha\beta}$ ($\alpha,\beta=x,yz,z$-coordinates), the cross-section
$\sigma(\omega)$ was calculated as

$\sigma(\omega)=\frac{4\pi\omega}{3c}\text{Im}(\alpha_{xx}+\alpha_{yy}+\alpha_{zz}).$
(58)

Here $\omega$ is the frequency of the corresponding polarizabilities, and $c$
is the vacuum speed of light. The QM region in PE-CPP calculations was always
chosen as [M(H2O)6]^3+ (M = Rh, Ir) and [Pt(N3)(OH)2(NH3)2]. We employed the
CAM-B3LYP105 functional with a def2-SV(P) basis set for the ligands, and a
dyall.v2z106 basis set for the metals. All calculations were carried out using
uncontracted basis sets. We used a non-colinear form of the exchange-
correlation kernel as implemented in DIRAC.107 A damping factor ($\gamma$) of
4.557$\times 10^{-3}$ Hartree (1000 cm-1) was used in all calculations.
Relativistic effects were included through the use of both 4c and X2C
Hamiltonians. The 4c calculations were carried out with (SS$|$SS) integrals
replaced by interatomic SS energy corrections.108 For the X2C calculations,
atomic-mean-field calculations were used to obtain a correction for the two-
electron spin-orbit interaction71, 109. To probe the quality of the atomic-
mean-field corrections, we compared X2C with these corrections to full 4c-KS
calculations for the aqua ions in vacuum and the resulting spectra are
essentially identical (data not shown). For trans-Pt, only X2C calculations
were carried out. To ensure that the calculated spectra did not contain
spurious/artificial transitions our assignment of the transitions involved a
careful examination of the response vectors. This included visual inspection
of the molecular orbitals involved in transitions, as well as inspection of
the atomic orbital contributions to the molecular orbitals in question.

##  4 Results and Discussion

###  4.1 Electronic structures of the model systems

Before discussing the UV-vis and X-ray spectra in the two subsections below,
we first make a few comments on the choice of test systems. We chose the aqua
ions since they are among the simplest transition metal ions showing
significant relativistic effects, allowing us to carry out both (PE-)4c and
(PE-)X2C calculations without exhausting our computational resources. Further,
they are inert $\text{d}^{6}$ aqua ions and therefore well-defined in aqueous
solution without fast exchange of water ligands; they have been experimentally
characterized by UV-vis110, 111 (in water) and in case of [Ir(H2O)6]^3+, also
by XAS.112 Further, their (approximately) octahedral symmetry dictates a
d-orbital splitting with a full $t^{6}_{2g}$ set and empty $e_{g}$ orbitals.
Thus, they can be expected to display types of excitations where the solvent
effect is rather different, namely d-d transitions ($t_{2g}\rightarrow e_{g}$)
and ligand-to-metal charge transfer (LMCT) transitions, involving the empty
$e_{g}$ set of d-orbitals. Analyzing the response vectors of the two aqua
complexes for the transitions in the UV-vis region (under 8 eV) showed indeed
that all intense transitions involve empty metal d-orbitals (the occupied
orbitals involved may either be ligand or other d-orbitals). For the trans-Pt
complex, only LMCT excitations are investigated and the nature of these
transitions was found to be mainly from $\pi$ orbitals on the azide ligands
(N3-) to that of the d-orbitals on Pt.

On one hand, the d-d transitions are not expected to be associated with large
re-organization of the electron density during the excitation, and likely show
little (electronic) effect of the solvent. On the other hand, LMCT transitions
can be expected to show large re-organization of the electron density during
excitation and hence larger (electronic) effect of the solvent is expected. We
will investigate this in the following sections and we will in general
separate solvent effects due to structural and electronic perturbation,
although a few words of caution are needed: the dipole-forbidden d-d
transitions can only gain intensity through symmetry lowering. This is partly
possible in the water solvent, but for the vacuum-optimized structures, the
only mechanism to lower the symmetry is through vibronic coupling. Vibrational
effects are not included in our calculations and a fair comparison to the
vacuum-optimized structures is not possible in this case. We therefore only
comment on electronic solvent contributions in this case.

Solvent effects including vibrational effects could be included through
molecular dynamics (MDs) and this would indeed be required for a proper
comparison with experimental results. However, we stress that our goal here is
to quantify different solvent effects (and to investigate if the solvent
effect is reproduced by X2C) rather than to make a detailed comparison with
experiment. In this regard, it should also be noted that only the d-d
transitions are experimentally observed for the aqua ions, since the LMCT
transitions are too high in energy. We include the LMCT in our analyses,
although in some cases they may overlap with transitions for the water solvent
(we discuss the consequences of this for the PE model below). The trans-Pt
complex will allow us to analyze a LMCT transition in the optically
transparent window of the water solvent.

We have also investigated the nature of electronic structure of the
investigated complexes regarding relativistic effects: the d-orbitals in the
two aqua complexes are heavily mixed in $\alpha$ or $\beta$-components, which
is only correctly obtained if spin-orbit coupling is accounted for; see Tables
S1 and S2 in the supporting information (SI). Our analyses further showed the
transitions at the L${}_{\text{II}}$ and L${}_{\text{III}}$-edges to be from
occupied p${}_{\frac{1}{2}}$ and p${}_{\frac{3}{2}}$-orbitals to unoccupied
d-orbitals (as expected) for [Rh(H2O)6]^3+ and [Ir(H2O)6]^3+. Accordingly, the
core-transitions are also metal-based and it could be expected that the
solvent effects are similar (although differences are seen, as will be shown
below). Still, in all cases the transitions in the aqua complexes cannot be
described without the inclusion of spin-orbit coupling. We have previously
investigated the platinum complex in the same manner and came to the same
conclusion.25 We will in the subsections below therefore not discuss non-
relativistic results, but focus on quantifying the effect of either solvation
or the effect of switching from 4c to X2C Hamiltonian.

###  4.2 UV-vis spectra of [Rh(H2O)6]^3+ and [Ir(H2O)6]^3+

The impact of going from vacuum to an aqueous environment on the UV-vis
spectra of the two aqua complexes is shown in Figures 2 and 3. We show the
effect of PE in panel (a) and that of additionally including EEF effects in
panels (b). Selected peak maxima and associated absorption cross-sections,
$\sigma(\omega)$, are given in Tables 1 and 2. Note that we in Figures 2 and 3
only show comparison to vacuum-optimized structures while the Tables also
includes calculations where we employ the structures optimized within the
solvent, but without a PE potential. The corresponding Figures comparing PE
and vacuum calculations with solvent optimized structures are given in the SI
(Figures S1 and S2).

As expected, the two complexes display intense ligand-to-metal charge transfer
excitations at high energies and two lower-lying d-d transitions. The d-d
transitions are not very intense and can not be discerned in Figures 2 and 3
(a magnified figure is provided in the SI, Figure S1 and S2). Although our
results for the d-d transitions correspond well with experiment111, 110 at 3.1
and 3.9 eV for [Rh(H2O)6]^3+ and 3.9 and 4.6–4.7 eV for [Ir(H2O)6]^3+ (see
Tables 1 and 2) we will not discuss them in detail, but rather focus on the
LMCT transitions. These transitions clearly demonstrate impact of solvation
leads to qualitatively different spectra, where both excitation energies at
peak maxima change and the absorption cross sections increase drastically.

(a)

(b)

Figure 2: UV-vis absorption spectra for [Rh(H2O)6]^3+. (a) Comparison of
spectra calculated in vacuum and with environment described through PE, where
(b) compares spectra calculated in environment with PE and EEF effects.

(a)

(b)

Figure 3: UV-vis absorption spectra for [Ir(H2O)6]^3+. (a) Comparison of spectra calculated in vacuum (including vacuum optimization) and with environment described through PE, where (b) compares spectra calculated in environment with PE and EEF effects. Table 1: Selected peak maxima positions $\Delta E$ (in eV) for [Rh(H2O)6]3+ (the corresponding spectra are given in Figure 2). The absorption cross sections, $\sigma(\omega)$, are given in parantheses.. For $\Delta E^{vac}_{vac-opt}$, the lack of vibronic coupling means that the d-d transitions have no intensity and they are therefore not reported. $\Delta E^{vac}_{vac-opt}$ | $\Delta E^{vac}_{solv-opt}$ | $\Delta E^{PE}_{solv-opt}$ | $\Delta E^{EEF}_{solv-opt}$ | Assignment
---|---|---|---|---
- | 7.63 (0.880) | 7.46 (1.348) | 7.46 (1.176) | LMCT
6.03 (0.272) | 7.08 (1.040) | 7.15 (1.205) | 7.15 (1.070) | LMCT
5.76 (0.107) | 6.81 (1.190) | 6.91 (1.221) | 6.91 (1.053) | LMCT
- | 6.54 (0.952, 6.51 (0.948) | 6.51 (0.893) | 6.51 (0.801) | LMCT
- | 3.89 (0.015), 3.85 (0.015) | 3.82 (0.013) | 3.82 (0.012) | d-d
- | 3.14 (0.010) | 3.11 (0.009) | 3.11 (0.008) | d-d
- | 2.97 (0.010) | 2.93 (0.009) | 2.93 (0.008) | d-d

For the case with structures optimized in solvent (i.e. the electronic solvent
contributions), the peak maxima positions and their associated absorption
cross-section are given in the $\Delta E^{vac}_{solv-opt}$ column in Tables 1
and 2. For the LMCT transitions, the pure solvent contributions to the peak
maxima vary between 0.0–0.17 eV for [Rh(H2O)6]^3+ and 0.07–0.17 eV for
[Ir(H2O)6]^3+, showing that the PE contributions are non-negligible.

The absorption cross-sections are also significantly affected: for
[Rh(H2O)6]^3+ at the peak maximum at 7.46 eV (with PE), we see an increase in
$\sigma(\omega)$ from 0.880 a.u. without PE to 1.348 a.u. with PE,
corresponding to an 53% increase. For the other LMCT peaks, the corresponding
increases are 16% and 2% for the peaks at 7.15 and 6.91 eV (with PE) or a
decrease of 6% for the peaks at 6.54 eV (with PE).

For [Ir(H2O)6]^3+, $\sigma(\omega)$ generally increases 20–27% due to the PE
environment: for the transition at 7.95 (with PE), $\sigma(\omega)$ increases
from 0.619 a.u. to 0.787 a.u. (27%), while $\sigma(\omega)$ for the transition
at 7.75 eV (with PE) increases from 0.610 a.u. to 0.733 a.u. (20 %).

The difference between d-d and LMCT transitions are clearly seen in their
response to the (pure) solvent effects (although we again emphasize that
vibronic coupling is not accounted for). The d-d transitions show little
effect of the solvent after the structural solvent perturbation is accounted
for, whereas the LMCT excitations are associated with a much more pronounced
electronic solvent effect (see $\Delta E^{vac}_{solv-opt}$ and $\Delta
E^{PE}_{solv-opt}$ in Tables 1 and 2 as well as Figures S1 and S2).

We also investigated the effect of an EEF, which is shown on Figures 2b and 3b
for [Rh(H2O)6]^3+ and [Ir(H2O)6]^3+, respectively. As expected, the excitation
energies at peak maxima positions are identical, but there is occasionally an
effect on the absorption cross sections: to quantify the effect, we compare
$\sigma(\omega)$ for the peak maxima positions in the $\Delta E^{PE}_{solv-
opt}$ column with the corresponding numbers in the $\Delta E^{EEF}_{solv-opt}$
column in Tables 1 and 2. This comparison shows that $\sigma(\omega)$ for LMCT
transitions are between 13-16 % ([Rh(H2O)6]^3+) and 10–12 % ([Ir(H2O)6]^3+)
larger without EEF. For both complexes, the effect from the EEF is much
smaller for the d-d transitions.

Table 2: Selected peak maxima positions $\Delta E$ (in eV) for [Ir(H2O)6]3+ (the corresponding spectra are given in Figure 3). The absorption cross sections, $\sigma(\omega)$, are given in parantheses. For $\Delta E^{vac}_{vac-opt}$, the lack of vibronic coupling means that the d-d transitions have no intensity and they are therefore not reported.  $\Delta E^{vac}_{vac-opt}$ | $\Delta E^{vac}_{solv-opt}$ | $\Delta E^{PE}_{solv-opt}$ | $\Delta E^{EEF}_{solv-opt}$ | Assignment
---|---|---|---|---
7.58 (0.351) | 7.79 (0.619) | 7.95 (0.787) | 7.95 (0.717) | LMCT
6.90 (0.070) | 7.82 (0.610) | 7.75 (0.733) | 7.75 (0.656) | LMCT
| 7.58 (0.593) | 7.72 (0.729) | 7.72 (0.652) | LMCT
- | 4.25 (0.008), 4.21 (0.008) | 4.18 (0.008) | 4.18 (0.007) | d-d
- | 3.77 (0.008), 3.74 (0.008) | 3.70 (0.007) | 3.70 (0.007) | d-d

With respect to the LMCT excitations, it should be noted that they in some
cases for [Ir(H2O)6]^3+ and for the highest-lying LMCT excitation for
[Rh(H2O)6]^3+ are close to resonant regions of the water solvent: the first
excitation of water is around 7.5 eV113. The use of static polarizabilities
will most likely lead to some error in these regions, but we leave
investigation of this for further studies. A remedy may be to include parts of
the solvent molecules fully in the QM region. We have so far not investigated
QM size effects, since we present both four- and two-component calculations
and the former becomes computationally heavy with extended QM regions. As we
will show in a section below, PE-4c and PE-X2C calculations provide close to
identical description of the solvent effects, which is a good foundation for
using X2C calculations in future investigations of the QM size. Alternatively,
it is possible to employ frequency-dependent polarizabilities, as has been
done recently by several different groups114, 115, 116 (within non-
relativistic frameworks). Still, frequency-dependent polarizabilities are not
without challenges for response-based methods, although it has been pointed
out115 that CPP methods are beneficial in this regard, since the frequency is
known at the input stage. Again, we will not pursue this further in this
paper, but rather turn to a LMCT excitation far from the resonant region of
water, namely the one in the trans-Pt complex investigated in next subsection.

###  4.3 UV-vis spectrum of trans-trans-trans-[Pt(N3)(OH)2(NH3)2]

(a)

(b)

(c)

(d)

Figure 4: UV-vis absorption spectra for trans-trans-
trans-[Pt(N3)2(OH)2(NH3)2].(a) Comparison of spectra calculated in vacuum
(including vacuum optimization) and with its environment described through PE.
(b) Compares spectra calculated in vacuum (with solvent optimized structure)
to that of spectra with PE. (c) Compares the effect of PE with an EEF to that
of only using PE. (d) Compares the effect of PE to that of only including
electrostatic embedding with either multipoles up to quadrupoles (ES) or
charges (ES,q).

The UV-vis spectra for trans-Pt is shown in Figure 4. We have previously
calculated this spectrum in vacuum with conventional response theory.25 We
found that around 80 roots in the response equations were necessary to cover
the UV-vis region; while this is possible for present days response solvers,
it is certainly still high. Meanwhile, the PE-X2C-CPP is not affected by the
high density of states: in Figure 4 we show UV-vis spectra for (a) the
structure optimization in vacuum and (b) with the structure optimized in
solvent. From the figure, we see that including the solvent has a significant
impact on the resulting spectra: the vacuum-optimized structure only shows one
peak at 4.48 eV (with a shoulder at 4.89 eV), whereas the solvent-optimized
structures has two intense peaks in the energy range 4.0–4.5 eV
(experimentally56, the complex has one large transition at 4.35 eV in water).
Interestingly, the single intense peak is restored (now at 4.38 eV), while an
additional (small) peak appears at 5.27 eV when a PE potential is included.
The absorption cross section also changes significantly, as seen from Table 3.
Depending on whether we compare $\Delta E^{vac}_{vac-opt}$ or $\Delta
E^{vac}_{solv-opt}$ to $\Delta E^{PE}_{solv-opt}$ (see Table 3), the
absorption cross section changes from 3.464 a.u. or 1.319 a.u. to 2.779 a.u.,
corresponding to a decrease of 62 % or an increase of 52 % (the latter purely
due to the electronic effect of the solvent).

As in previous section, we also quantify the effect on an EEF, which is shown
in panel (c) of Figure 4: while the effect is smaller than the effect of
moving from vacuum to PE, this effect is not negligible; $\sigma(\omega)$ is
around 12 % larger without EEF.

Table 3: Selected peak maxima positions $\Delta E$ (in eV) for trans-trans-trans-[Pt(N3)2(OH)2(NH3)2] (the corresponding spectra are given in Figure 4. The absorption cross sections, $\sigma(\omega)$, are given in parantheses (in a.u.).  $\Delta E^{vac}_{vac-opt}$ | $\Delta E^{vac}_{solv-opt}$ | $\Delta E^{PE}_{solv-opt}$ | $\Delta E^{EEF}_{solv-opt}$ | Assignment
---|---|---|---|---
4.89 (0.934) | 4.52 (1.740) | 5.27 (0.883) | 5.27 (0.795) | LMCT
4.48 (3.464) | 3.97 (1.319) | 4.38 (2.779) | 4.38 (2.495) | LMCT
Table 4: Selected peak maxima positions $\Delta E$ (in eV) for trans-trans-trans-[Pt(N3)2(OH)2(NH3)2] with increasing sophistication of the PE potential (the corresponding spectra are given in Figure 4(d). The absorption cross sections, $\sigma(\omega)$, are given in parantheses (in a.u.). $\Delta E^{es,q}_{solv-opt}$ | $\Delta E^{es}_{solv-opt}$ | $\Delta E^{PE}_{solv-opt}$ | Assignment
---|---|---|---
5.23 (0.725) | 5.30 (0.781) | 5.27 (0.883) | LMCT
4.28 (2.495) | 4.42 (2.524) | 4.38 (2.779) | LMCT

In addition to the dissection of structural and electron solvent effects, we
have for the trans-Pt complex also dissected the electronic effect from the PE
potential into the effect from multipoles and the additional effect of
polarization. Morevoer, we also test the case where only charges are included.
The spectra are compared in Figure 4 (d) and peak maxima positions together
with corresponding absorption cross sections are given in Table 4. We focus
again on the most intense LMCT transiton around 4.3–4.4 eV. As seen from the
Figure, the additional effect of polarization is rather small for the
excitation energy; in fact the large qualitative charge from the vacuum
spectrum in Figures 4 (a)–(b) are introduced already by including the charges,
although the higher order multiples have some effect on the excitation energy:
adding multipoles increases the excitation energy with 0.14 eV (the effect on
$\sigma(\omega)$ is only 1%). The effect on the excitation energy with respect
to adding polarization is smaller; it amounts to a 0.04 eV decrease. Yet, the
effect on $\sigma(\omega)$ is considerably larger as the absorption cross
section increases from 2.524 a.u. to 2.779 a.u. corresponding to 10 %.

###  4.4 X-ray atomic spectroscopy for [Rh(H2O)6]^3+ and [Ir(H2O)6]^3+

Before discussing the XAS spectra, it should be emphasized that the
calculation of L-edge spectra is met with a number of issues not seen in
corresponding UV-vis spectra calculations. Although the CPP method
(relativistic or non-relativistic) has a number of advantages, it also comes
with disadvantages and PE-CPP inherits both. One disadvantage (inherent to all
response methods) is the lack of core-polarization, while another disadvantage
(specific to the choice of DFT) is the self-interaction error; both of these
issues can lead to large shifts on the transitions compared with experiment7
and we also observe this here. Several schemes to improve the absolute
excitation energies have been proposed (see e.g. the work by South et al. 24),
but since we here only focus on the effect of the environment, these schemes
were not employed. Moreover, it is not uncommon in the calculation of L-edge
spectra to see spurious transitions between valence- and high-lying virtual
orbitals.24, 23 As described more thoroughly elsewhere117, 118, 24, this
artifact is related to the employed local Gaussian basis sets, which are not
appropriate for continuum or near-continuum states. Indeed, our analysis of
the response vectors revealed a small number of such spurious transitions for
a few of the higher-lying L-edge peaks. The issue is not guaranteed to be
resolved with larger (and more diffuse) basis sets24 and this has therefore
not been pursued. Previous investigations have sometimes chosen to focus
exclusively on excitations below the ionization threshold118, 24 and we will
also do that here (we estimate the ionization threshold from the orbital
energy of either p${}_{\frac{1}{2}}$ or p${}_{\frac{3}{2}}$).24 A pragmatic
approach is to invoke channel-restrictions as a remedy to spurious
transitions102, 103, 104, 23 and we have also employed this strategy.

(a)

(b)

(c)

(d)

Figure 5: L-edge spectra of [Rh(H2O)6]^3+. The panels (a) and (b) compares of vacuum and PE calculations for (a) L${}_{\text{II}}$-edge and (b) the L${}_{\text{III}}$-edge, while panels (c) and (d) compares PE calculations environment with and without EEF effects for (c) L${}_{\text{II}}$-edge and (d) the L${}_{\text{III}}$-edge.  Table 5: Selected peak maxima positions $\Delta E$ (in eV) for the L${}_{\text{II}}$\- and L${}_{\text{III}}$-edges in [Rh(H2O)6]3+ (the corresponding spectra are shown in Figure 5). The absorption cross sections, $\sigma(\omega)$, are given in parantheses.  $\Delta E^{vac}_{vac-opt}$ | $\Delta E^{vac}_{solv-opt}$ | $\Delta E^{PE}_{solv-opt}$ | $\Delta E^{EEF}_{solv-opt}$ | Assignment
---|---|---|---|---
3108.28 (0.172) | 3108.40 (0.171) | 3108.40 (0.163) | 3108.40 (0.155) | p${}_{\frac{1}{2}}$$\rightarrow$d (L${}_{\text{II}}$)
2963.75 (0.432) | 2963.81 (0.408) | 2963.81 (0.404) | 2963.81 (0.381) | p${}_{\frac{3}{2}}$$\rightarrow$d (L${}_{\text{III}}$)

(a)

(b)

(c)

(d)

Figure 6: L-edge spectra of [Ir(H2O)6]^3+. The panels (a) and (b) compares of vacuum and PE calculations for (a) L${}_{\text{II}}$-edge and (b) the L${}_{\text{III}}$-edge, while panels (c) and (d) compares PE calculations environment with and without EEF effects for (c) L${}_{\text{II}}$-edge and (d) the L${}_{\text{III}}$-edge. Table 6: Selected peak maxima positions $\Delta E$ (in eV) for the L${}_{\text{II}}$\- and L${}_{\text{III}}$-edges in [Ir(H2O)6]3+ (the corresponding spectra are shown in Figure 6). The absorption cross sections, $\sigma(\omega)$, are given in parentheses.  $\Delta E^{vac}_{vac-opt}$ | $\Delta E^{vac}_{solv-opt}$ | $\Delta E^{PE}_{solv-opt}$ | $\Delta E^{EEF}_{solv-opt}$ | Assignment
---|---|---|---|---
12771.6 (0.024) | 12771.8 (0.022) | 12771.9 (0.017) | 12771.9 (0.016) | p${}_{\frac{1}{2}}$$\rightarrow$d (L${}_{\text{II}}$)
11139.2 ( 0.090) | 11139.4 (0.085) | 11139.4 (0.066) | 11139.4 (0.062) | p${}_{\frac{3}{2}}$$\rightarrow$d (L${}_{\text{III}}$)

The L${}_{\text{II}}$ and L${}_{\text{III}}$ spectra are shown in Figures 5
and 6 while peak-maxima positions are provided in Tables 5 and 6 for
[Rh(H2O)6]^3+ and [Ir(H2O)6]^3+, respectively. As for the UV-vis spectra, the
table also show results, $\Delta E^{vac}_{solv-opt}$, for spectra calculated
in vacuum with underlying solvent-optimized structures. The corresponding
Figures are given in the SI (Figures S3 and S4). The calculated XAS spectra
contain for both edges (and complexes) one main peak, corresponding to either
a p${}_{\frac{1}{2}}$$\rightarrow$d or a p${}_{\frac{3}{2}}$$\rightarrow$d
transition located on the metal. This transition corresponds to the main peak
in the experimental spectrum112 of [Ir(H)2O)6]^3+. As expected, the CPP value
is shifted from the experiment, both due to the lack of orbital-polarization
effects and self-interaction error in the DFT functional.24, 23, 7 We note
that two other peaks appear in the experimental spectrum at 16–17 eV and 44 eV
above the main peak.112 None of these peaks are within the frequency window
for our L${}_{\text{III}}$-edge calculations, wheres we for the
L${}_{\text{II}}$-edge included up to around 12 eV above the main peak. While
we do see small appearances of peaks around 5-10 eV above the main peak, we do
not include these in our discussion as they are just around the ionization
threshold (12782 eV for L${}_{\text{II}}$ in the PE calculations).

(a)

(b)

(c)

(d)

Figure 7: PE-4c- and PE-X2C-CPP spectra in UV-vis and XAS regions for
[Ir(H2O)6]^3+. (a) Full UV-vis spectrum and (c) focuses on the d-d
transitions. The XAS spectra are shown in (b) for the L${}_{\text{III}}$ (X2C
spectra shifted 12.9 eV) and (d) L${}_{\text{II}}$ (X2C spectra shifted 30.61
eV).

Focusing now on the main transition, this excitation is located entirely on
orbitals on the metal. In relation to any potential solvent effect, we
therefore anticipated a case similar to d-d transitions discussed above, i.e.,
only benign effects from the solvent. For the [Rh(H2O)6]^3+ complex, this is
indeed the case. However, for the [Ir(H2O)6]^3+ complex, the change in the
absorption cross-section is significant: we obtain a decrease in
$\sigma(\omega)$ from 0.024 a.u. to 0.017 with PE and to 0.016 a.u. with EEF
for the L${}_{\text{II}}$-edge. The corresponding numbers for the
L${}_{\text{III}}$-edge are decreases of 0.090 a.u. to 0.066 a.u. (PE) or
0.062 a.u. (EEF), cf. Table 6. Thus, $\sigma(\omega)$ decreases 29-33 % for
the L${}_{\text{II}}$-edge and 27-31% for the L${}_{\text{III}}$-edge. By
comparing calculations without PE in the solvated structure ($\Delta
E^{vac}_{solv-opt}$ in Table 6) to the pure vacuum calculation ($\Delta
E^{vac}_{vac-opt}$), we see that the structural effect is rather small for the
absorption cross section. Thus, the change is rather introduced by the
electronic (PE) part, $\Delta E^{PE}_{solv-opt}$ in Table 6, showing that the
environment electrostatics can be important, also for L-edge XAS spectroscopy.
Note also that EEF effects only account for 11-12 % of the total change.

###  4.5 Comparing four- and exact two-component calculations

In this last section, we investigate the impact of using a two-component
Hamiltonian (X2C) on the solvated UV-vis and XAS spectra. This is done by
comparing the solvated spectra calculated with PE-X2C-CPP to their PE-4c-CPP
counterparts. The results are shown in Figures 7 for [Ir(H2O)6]^3+ (results
for [Rh(H2O)6]^3+ are similar and provided in the SI, Figure S5). The UV-vis
spectra are practically identical, whereas the change of Hamiltonian leads to
shifts in absolute excitation energies in the core spectra. However, applying
this shift as done in Figure 7, the spectra become identical, showing the
PE-X2C-CPP method works extremely well.

##  5 Conclusion

We have presented the derivation and implementation of the CPP method in a
relativistic PE framework. The implementation is illustrated by calculations
of both UV-vis and XAS spectra of solvated aqua-complexes [Rh(H2O)6]^3+ and
[Ir(H2O)6]^3+. We next move to the larger complex trans-trans-
trans-[Pt(N3)2(OH)2(NH3)2], where we also show that inclusion of the
environment can qualitatively change the obtained spectra. In this case, we
show that the largest electronic effect is from the charges, but
polarizabilities in the environment have non-negligible effect on the
absorption cross section.

In all cases, the complexes required a relativistic method to correctly
capture their electronic structures, while treatment of the environment is
indispensable for a correct physical description of the UV-vis and XAS
spectra. Finally, our calculations show that a PE-X2C-CPP model essentially
reproduces corresponding PE-4c-CPP spectra (this was done for the aqua ions);
this highlights that the solvent effects are described well on the two-
component level.

While our current target systems are experimentally known, we mostly consider
them toy models. We are currently investigating the dynamic effects of the
solvent, which is not included in our present calculations. This is mainly due
to the large computational demands for these types of calculations and the
challenge with lack of accurate force field parameters for transition metal
systems. Another line of work we will pursue is to extend the method to
electronic circular dichroism (ECD) spectroscopy.

The authors thank The Villum Foundation, Young Investigator Program (grant no.
29412), the Swedish Research Council (grant no. 2019-04205), and Independent
Research Fund Denmark (grant no. 0252-00002B) for support. The computations
were performed on computer resources provided by the Swedish National
Infrastructure for Computing (SNIC) at Lunarc (Lund University) and HPC2N
(Umeå University).

The supporting information contains Tables with Mulliken ($\alpha$ and
$\beta$) populations of the metal d-orbitals and energies of metal 2p
orbitals. Moreover, the SI contains spectra for the solvent optimized hexaaqua
ions without PE and a zoom in of the region with d-d transitions.

## References

  * Olsen and Jørgensen 1985 Olsen, J.; Jørgensen, P. Linear and non-linear response functions for an exact state and for an MCSCF state. _J. Chem. Phys._ 1985, _82_ , 3235–3264
  * Helgaker et al. 2012 Helgaker, T.; Coriani, S.; Jørgensen, P.; Kristiansen, K.; Olsen, J.; Ruud, K. Recent advances in wave function based methods of molecular-property calculations. _Chem. Rev._ 2012, _112_ , 543–631
  * Runge and Gross 1984 Runge, E.; Gross, E. K. U. Density-functional theory for time-dependent systems. _Phys. Rev. Lett._ 1984, _52_ , 997–1000
  * 4 Casida, M. E. _Recent advances in density functional methods, Part I_ ; World Scientific, Singapore
  * Bauernschmitt and Ahlrichs 1996 Bauernschmitt, R.; Ahlrichs, R. Treatment of electronic excitations within the adiabatic approximation of time dependent density functional theory. _Chem. Phys. Lett._ 1996, _256_ , 454–464
  * Pyykkö 2012 Pyykkö, P. Relativistic effects in chemistry: More common than you thought. _Ann. Rev. of Phys. Chem._ 2012, _63_ , 45–64
  * Norman and Dreuw 2018 Norman, P.; Dreuw, A. Simulating X-ray spectroscopies and calculating core-excited states of molecules. _Chem. Rev._ 2018, _118_ , 7208–7248
  * Liu and Xiao 2018 Liu, W.; Xiao, Y. Relativistic time-dependent density functional theories. _Chem. Soc. Rev_ 2018, _47_ , 4481–4509
  * Saue et al. 2020 Saue, T. et al. The DIRAC code for relativistic molecular calculations. _J. Chem. Phys._ 2020, _152_ , 204104
  * Repiský et al. 2020 Repiský, M.; Komorovský, S.; Kadek, M.; Konecny, L.; Ekström, U.; Malkin, E.; Kaupp, M.; Ruud, K.; Malkina, O. L.; Malkin, V. G. ReSpect: Relativistic spectroscopy DFT program package. _J. Chem. Phys._ 2020, _152_ , 184101
  * Gao et al. 2004 Gao, J.; Liu, W.; Song, B.; Liu, C. Time-dependent four-component relativistic density functional theory for excitation energies. _J. Chem. Phys._ 2004, _121_ , 6658–6666
  * Gao et al. 2005 Gao, J.; Zou, W.; Liu, W.; Xiao, Y.; Peng, D.; Song, B.; Liu, C. Time-dependent four-component relativistic density-functional theory for excitation energies. II. The exchange-correlation kernel. _J. Chem. Phys._ 2005, _123_ , 054102
  * Bast et al. 2009 Bast, R.; Jensen, H. J. Aa.; Saue, T. Relativistic adiabatic time-dependent density functional theory using hybrid functionals and noncollinear Spin magnetization. _Int. J. Quantum Chem._ 2009, _109_ , 2091–2112
  * Wang et al. 2005 Wang, F.; Ziegler, T.; Van Lenthe, E.; Van Gisbergen, S.; Baerends, E. J. The calculation of excitation energies based on the relativistic two-component zeroth-order regular approximation and time-dependent density-functional with full use of symmetry. _J. Chem. Phys._ 2005, _122_ , 204103
  * Peng et al. 2005 Peng, D.; Zou, W.; Liu, W. Time-dependent quasirelativistic density-functional theory based on the zeroth-order regular approximation. _J. Chem. Phys._ 2005, _123_ , 144101
  * Kühn and Weigend 2013 Kühn, M.; Weigend, F. Implementation of two-component time-dependent density functional theory in TURBOMOLE. _J. Chem. Theory Comput._ 2013, _19_ , 5341–5348
  * Egidi et al. 2016 Egidi, F.; Goings, J. J.; Frisch, M. J.; Li, X. Direct atomic-orbital-based relativistic two-component linear response method for calculating excited-state fine structures. _J. Chem. Theory Comput._ 2016, _12_ , 3711–3718
  * Norman et al. 2001 Norman, P.; Bishop, D. M.; Jensen, H. J. A.; Oddershede, J. Near-resonant absorption in the time-dependent self-consistent field and multiconfigurational self-consistent field approximations. _J. Chem. Phys._ 2001, _115_ , 10323–10334
  * Norman et al. 2004 Norman, P.; Ruud, K.; Helgaker, T. Density-functional theory calculations of optical rotatory dispersion in the nonresonant and resonant frequency regions. _J. Chem. Phys._ 2004, _120_ , 5027–5035
  * Kauczor and Norman 2014 Kauczor, J.; Norman, P. Efficient calculations of molecular linear response properties for spectral regions. _J. Chem. Theory Comput._ 2014, _10_ , 2449–2455
  * Villaume et al. 2010 Villaume, S.; Saue, T.; Norman, P. Linear complex polarization propagator in a four-component Kohn–Sham framework. _J. Chem. Phys._ 2010, _133_ , 064105
  * Konecny et al. 2019 Konecny, L.; Repisky, M.; Ruud, K.; Komorovsky, S. Relativistic four-component linear damped response TDDFT for electronic absorption and circular dichroism calculations. _J. Chem. Phys._ 2019, _151_ , 194112
  * Fransson et al. 2016 Fransson, T.; Burdakova, D.; Norman, P. K- and L-edge X-ray absorption spectrum calculations of closed-shell carbon, silicon, germanium, and sulfur compounds using damped four-component density functional response theory. _Phys. Chem. Chem. Phys._ 2016, _18_ , 13591–13603
  * South et al. 2016 South, C.; Shee, A.; Mukherjee, D.; Wilson, A. K.; Saue, T. 4-Component relativistic calculations of L3 ionization and excitations for the species UO${}^{2+}_{2}$, OUN+ and UN2. _Phys. Chem. Chem. Phys._ 2016, _18_ , 21010–21023
  * Creutzberg and Hedegård 2020 Creutzberg, J.; Hedegård, E. D. Investigating the influence of relativistic effects on absorption spectra for platinum complexes with light-activated activity against cancer cells. _Phys. Chem. Chem. Phys._ 2020, _22_ , 27013–27023
  * Tomasi et al. 2005 Tomasi, J.; Mennucci, B.; Cammi, R. Quantum mechanical continuum solvation models. _Chem. Rev._ 2005, _105_ , 2999–3093
  * Di Remigio et al. 119 Di Remigio, R.; Bast, R.; Frediani, L.; Saue, T. Four-component relativistic calculations in solution with the polarizable continuum model of solvation: Theory, implementation, and application to the Group 16 dihydrides H2X (X = O, S, Se, Te, Po). _J. Phys. Chem. A_ 119, _119_ , 5061–5077
  * Warshel and Levitt 1976 Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. _J. Mol. Biol._ 1976, _103_ , 227–249
  * Senn and Thiel 2009 Senn, H. M.; Thiel, W. QM/MM methods for biomolecular systems. _Angew. Chem. Int. Ed._ 2009, _48_ , 1198–1229
  * Wesolowski and Warshel 1993 Wesolowski, T. A.; Warshel, W. Frozen density functional approach for ab initio calculations of solvated molecules. _J. Phys. Chem._ 1993, _97_ , 8050–8053
  * Gomes et al. 2008 Gomes, A. S. P.; Jacob, C. R.; Visscher, L. Calculation of local excitations in large systems by embedding wave-function theory in density-functional theory. _Phys. Chem. Chem. Phys_ 2008, _10_ , 5353–5362
  * Höfener et al. 2012 Höfener, S.; Gomes, A. S. P.; Visscher, L. Molecular properties via a subsystem density functional theory formulation: A common framework for electronic embedding. _J. Chem. Phys._ 2012, _136_ , 044104
  * Bondanza et al. 2020 Bondanza, M.; Nottoli, M.; Cupellini, L.; Lipparini, F.; Mennucci, B. Polarizable embedding QM/MM: The future gold standard for complex (bio)systems? _Phys. Chem. Chem. Phys._ 2020, _22_ , 14433–14448
  * Thompson and Schenter 1995 Thompson, M. A.; Schenter, G. K. Excited states of the bacteriochlorophyll b dimer of Rhodopseudomonas viridis: A QM/MM study of the photosynthetic reaction center that includes MM polarization. _J. Phys. Chem._ 1995, _99_ , 6374–6386
  * Thompson 1996 Thompson, M. A. QM/MMpol: A consistent model for solute/solvent polarization. Application to the aqueous solvation and spectroscopy of formaldehyde, acetaldehyde, and acetone. _J. Phys. Chem._ 1996, _100_ , 14492–14507
  * Jensen et al. 2003 Jensen, L.; van Duijnen, P. T.; Snijders, J. G. A discrete solvent reaction field model for calculating molecular linear response properties in solution. _J. Chem. Phys._ 2003, _119_ , 3800–3809
  * Jensen et al. 2003 Jensen, L.; van Duijnen, P. T.; Snijders, J. G. A discrete solvent reaction field model within density functional theory. _J. Chem. Phys._ 2003, _118_ , 514–521
  * Gordon et al. 2007 Gordon, M. S.; Slipchenko, L.; Li, H.; Jensen, J. H. Chapter 10 The effective fragment potential: A general method for predicting intermolecular interactions. _Annu. Rep. Comput. Chem._ 2007, _3_ , 177–193
  * Söderhjelm et al. 2009 Söderhjelm, P.; Husberg, C.; Strambi, A.; Olivucci, M.; Ryde, U. Protein influence on electronic spectra modeled by multipoles and polarizabilities. _J. Chem. Theory Comput._ 2009, _5_ , 649–658
  * Dziedzic et al. 2016 Dziedzic, J.; Mao, Y.; Shao, Y.; Ponder, J.; Head-Gordon, T.; Head-Gordon, M.; Skylaris, C. K. TINKTEP: A fully self-consistent, mutually polarizable QM/MM approach based on the AMOEBA force field. _J. Chem. Phys._ 2016, _145_ , 124106
  * Hedegård and Reiher 2016 Hedegård, E. D.; Reiher, M. Polarizable embedding density matrix renormalization group. _J. Chem. Theory Comput._ 2016, _12_ , 4242–4253
  * Steinmann et al. 2019 Steinmann, C.; Reinholdt, P.; Nørby, M. S.; Kongsted, J.; Olsen, J. M. H. Response properties of embedded molecules through the polarizable embedding model. _Int. J. Quantum Chem._ 2019, _119_ , e25717
  * Kongsted et al. 2003 Kongsted, J.; Osted, A.; Mikkelsen, K. V.; Christiansen, O. Linear response functions for coupled cluster/molecular mechanics including polarization interactions. _J. Chem. Phys._ 2003, _118_ , 1620–1633
  * Sneskov et al. 2011 Sneskov, K.; Schwabe, T.; Kongsted, J.; Christiansen, O. The polarizable coupled cluster method. _J. Chem. Phys._ 2011, _134_ , 104108
  * Hedegård et al. 2013 Hedegård, E. D.; List, N. H.; Jensen, H. J. A.; Kongsted, J. The multi-configuration self-consistent field method within a polarizable embedded framework. _J. Chem. Phys._ 2013, _139_ , 044101
  * List et al. 2013 List, N. H.; Jensen, H. J. A.; Kongsted, J.; Hedegård, E. D. A unified framework for the polarizable embedding and continuum methods within multiconfigurational self-consistent field theory. _Adv. Quant. Chem._ 2013, _66_ , 195–238
  * Hedegård et al. 2015 Hedegård, E. D.; Olsen, J. M. H.; Knecht, S.; Kongsted, J.; Jensen, H. J. A. Polarizable embedding with a multiconfiguration short-range density functional theory linear response method. _J. Chem. Phys._ 2015, _114_ , 1102–1107
  * Yoo et al. 2008 Yoo, S.; Zahariev, F.; Sok, S.; Gordon, M. S. Solvent effects on optical properties of molecules: A combined time-dependent density functional theory/effective fragment potential approach. _J. Chem. Phys._ 2008, _129_ , 144112
  * 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
  * Loco et al. 2016 Loco, D.; Polack, É.; Caprasecca, S.; Lagardère, L.; Lipparini, F.; Piquemal, J. P.; Mennucci, B. A QM/MM approach using the AMOEBA polarizable embedding: From ground state energies to electronic excitations. _J. Chem. Theory Comput._ 2016, _12_ , 3654–3661
  * Pedersen et al. 2014 Pedersen, M. N.; Hedegård, E. D.; Olsen, J. M. H.; Kauczor, J.; Norman, P.; Kongsted, J. Damped response theory in combination with polarizable environments: The polarizable embedding complex polarization propagator method. _J. Chem. Theory Comput._ 2014, _10_ , 1164–1171
  * Nørby et al. 2017 Nørby, M. S.; Olsen, J. M. H.; Steinmann, C.; Kongsted, J. Modeling electronic circular dichroism within the polarizable embedding approach. _J. Chem. Theory Comput._ 2017, _13_ , 4442–4451
  * Hedegård et al. 2017 Hedegård, E. D.; Bast, R.; Kongsted, J.; Olsen, H. J. A., J. M. H. Jensen Relativistic polarizable embedding. _J. Chem. Theory Comput._ 2017, _13_ , 2870–2880
  * Krause et al. 2016 Krause, K.; Bauer, M.; Klopper, W. Approaching phosphorescence lifetimes in solution: The two-component polarizable-embedding approximate coupled-cluster method. _J. Chem. Theory Comput._ 2016, _12_ , 2853–2860
  * 55 DIRAC, a relativistic ab initio electronic structure program, Release DIRAC21 (2021), written by R. Bast, A. S. P. Gomes, T. Saue, L. Visscher, and H. J. Aa. Jensen, with contributions from I. A. Aucar, V. Bakken, K. G. Dyall, S. Dubillard, U. Ekström, E. Eliav, T. Enevoldsen, E. Faßhauer, T. Fleig, O. Fossgaard, L. Halbert, E. D. Hedegård, T. Helgaker, B. Helmich–Paris, J. Henriksson, M. Iliaš, Ch. R. Jacob, S. Knecht, S. Komorovský, O. Kullie, J. K. Lærdahl, C. V. Larsen, Y. S. Lee, N. H. List, H. S. Nataraj, M. K. Nayak, P. Norman, G. Olejniczak, J. Olsen, J. M. H. Olsen, A. Papadopoulos, Y. C. Park, J. K. Pedersen, M. Pernpointner, J. V. Pototschnig, R. di Remigio, M. Repisky, K. Ruud, P. Sałek, B. Schimmelpfennig, B. Senjean, A. Shee, J. Sikkema, A. Sunaga, A. J. Thorvaldsen, J. Thyssen, J. van Stralen, M. L. Vidal, S. Villaume, O. Visser, T. Winther, and S. Yamamoto
  * Mackay et al. 2006 Mackay, F. S.; Woods, J. A.; Moseley, H.; Ferguson, J.; Dawson, A.; Parsons, S.; Sadler, P. J. A photoactivated trans-diammine platinum complex as cytotoxic as cisplatin. _Chem. Eur. J._ 2006, _12_ , 3155–3161
  * Bednarski et al. 2006 Bednarski, P. J.; Grünert, R.; Zielzki, M.; Wellner, A.; Mackay, F. S.; Sadler, P. J. Light-activated destruction of cancer cell nuclei by platinum diazide complexes. _Chem. Biol._ 2006, _13_ , 61–67
  * Freitag and González 2021 Freitag, L.; González, L. The role of triplet states in the photodissociation of a platinum azide complex by a density matrix renormalization group method. _J. Phys. Chem. Lett._ 2021, _12_ , 4876–4881
  * Helgaker et al. 2004 Helgaker, T.; Jørgensen, P.; Olsen, J. _Molecular Electronic-Structure Theory_ ; Wiley, 2004
  * Norman et al. 2005 Norman, P.; Bishop, D. M.; Jensen, H. J. A.; Oddershede, J. Nonlinear response theory with relaxation: The first-order hyperpolarizability. _J. Chem. Phys._ 2005, _123_ , 194103
  * Saue and Helgaker 2002 Saue, T.; Helgaker, T. Four-component relativistic Kohn–Sham theory. _J. Comput. Chem._ 2002, _23_ , 814–823
  * Salek et al. 2002 Salek, P.; Vahtras, O.; Helgaker, T.; Ågren, H. Density-Functional Theory of linear and nonlinear time-dependent molecular properties. _J. Chem. Phys._ 2002, _117_ , 9630–9645
  * Olsen et al. 2010 Olsen, J. M.; Aidas, K.; Mikkelsen, K. V.; Kongsted, J. Solvatochromic shifts in uracil: A combined MD-QM/MM study. _J. Chem. Theory and Comput._ 2010, _6_ , 249–256
  * Krause and Klopper 2016 Krause, K.; Klopper, W. Communication: A simplified coupled-cluster Lagrangian for polarizable embedding. _J. Chem. Phys._ 2016, _144_ , 041101
  * Olsen and Kongsted 2011 Olsen, J. M.; Kongsted, J. Molecular properties through polarizable embedding. _Adv. Quant. Chem._ 2011, _61_ , 107–143
  * Stone 2002 Stone, A. J. _The theory of intermolecular forces_ ; Oxford University Press, 2002
  * Olsen 2016 Olsen, J. M. H. PElib, The polarizable embedding library, version 1.2.4. https://gitlab.com/pe-software/pelib-public, 2016
  * 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
  * Applequist et al. 1972 Applequist, J.; Carl, J. R.; Fung, K.-K. An atom dipole interaction model for molecular polarizability. Application to polyatomic molecules and determination of atom polarizabilities. _J. Am. Chem. Soc._ 1972, _94_ , 2952–2960
  * Saue 2011 Saue, T. Relativistic hamiltonians for chemistry: A primer. _Chem. Phys. Chem._ 2011, _12_ , 3077–3094
  * Heß et al. 1996 Heß, B. A.; Marian, C. M.; Wahlgren, U.; Gropen, O. A mean-field spin-orbit method applicable to correlated wavefunctions. _Chem. Phys. Lett._ 1996, _251_ , 365–371
  * Sikkema et al. 2009 Sikkema, J.; Visscher, L.; Saue, T.; Iliaš, M. The molecular mean-field approach for correlated relativistic calculations. _J. Chem. Phys._ 2009, _131_ , 124116
  * Baerends et al. 1990 Baerends, E. J.; Schwarz, W. H. E.; Schwerdtfeger, P.; Snijders, J. G. Relativistic atomic orbital contractions and expansions: magnitudes and explanations. _J. Phys. B: At. Mol. Opt._ 1990, _23_ , 3225–3240
  * Kellö and Sadlej 1998 Kellö, V.; Sadlej, A. J. Picture change and calculations of expectation values in approximate relativistic theories. _Int. J. Quantum Chem._ 1998, _68_ , 159–174
  * Dyall 2000 Dyall, K. G. Relativistic electric and magnetic property operators for two-component transformed Hamiltonians. _Int. J. Quantum Chem._ 2000, _78_ , 412–421
  * Christiansen et al. 1998 Christiansen, O.; Jørgensen, P.; Hättig, C. Response functions from Fourier component variational perturbation theory applied to a time-averaged quasienergy. _Int. J. Quantum Chem._ 1998, _68_ , 1–52
  * Salek et al. 2005 Salek, P.; Helgaker, T.; Saue, T. Linear response at the 4-component relativistic density-functional level: application to the frequency-dependent dipole polarizability of Hg, AuH and PtH2. _Chem. Phys._ 2005, _311_ , 187–201
  * Jensen et al. 2005 Jensen, L.; Swart, M.; van Duijnen, P. T. Microscopic and macroscopic polarization within a combined quantum mechanics and molecular mechanics model. _J. Chem. Phys._ 2005, _122_ , 034103
  * List et al. 2016 List, N. H.; Jensen, H. J. A.; Kongsted, J. Local electric fields and molecular properties in heterogeneous environments through polarizable embedding. _Phys. Chem. Chem. Phys._ 2016, _18_ , 10070–10080
  * Kristiansen et al. 2009 Kristiansen, K.; Kauczor, J.; Kjærgaard, T.; Jørgensen, P. Quasienergy formulation of damped response theory. _J. Chem. Phys._ 2009, _131_ , 044112
  * Kauczor et al. 2011 Kauczor, J.; Jørgensen, P.; Norman, P. On the efficiency of algorithms for solving Hartree–Fock and Kohn–Sham response equations. _J. Chem. Theory Comput._ 2011, _7_ , 1610–1630
  * Saue and Jensen 2003 Saue, T.; Jensen, H. J. A. Linear response at the 4-component relativistic level: Application to the frequency-dependent dipole polarizabilities of the coinage metal dimers. _J. Chem. Phys._ 2003, _118_ , 522–536
  * qsi 2021 Schrödinger Release 2021-1: QSite, Schrödinger, LLC, New York, NY). 2021
  * Jorgensen et al. 1983 Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. _J. Chem. Phys._ 1983, _79_ , 926–935
  * Philipp and Friesner 1999 Philipp, D. M.; Friesner, R. A. Mixed ab initio QM/MM modeling using frozen orbitals and tests with alanine dipeptide and tetrapeptide. _J. Comput. Chem._ 1999, _20_ , 1442–1457
  * Murphy et al. 2000 Murphy, R. B.; Philipp, D. M.; Friesner, R. A. A mixed quantum mechanics/molecular mechanics (QM/MM) method for large-scale modeling of chemistry in protein environments. _J. Comput. Chem._ 2000, _21_ , 1442–1457
  * Becke 1988 Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. _Phys. Rev. A_ 1988, _38_ , 3098–3100
  * Lee et al. 1988 Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. _Phys. Rev. B_ 1988, _37_ , 785–789
  * Becke 1993 Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. _J. Chem. Phys._ 1993, _98_ , 5648–5652
  * Hehre et al. 1972 Hehre, W. J.; Ditchfield, R.; Pople, J. A. Self-consistent molecular orbital methods. XII. Further extensions of Gaussian-type basis sets for use in molecular orbital studies of organic molecules. _J. Chem. Phys._ 1972, _56_ , 2257–2261
  * Ditchfield et al. 1971 Ditchfield, R.; Hehre, W. J.; Pople, J. A. Self-consistent molecular-orbital methods. IX. An extended Gaussian-type basis for molecular-orbital studies of organic molecules. _J. Chem. Phys._ 1971, _54_ , 724–728
  * Hariharan and Pople 1973 Hariharan, P. C.; Pople, J. A. The influence of polarization functions on molecular orbital hydrogenation energies. _Theor. Chim. Acta_ 1973, _28_ , 213–222
  * Hay and Wadt 1985 Hay, P. J.; Wadt, W. R. Ab initio effective core potentials for molecular calculations. Potentials for the transition metal atoms Sc to Hg. _J. Chem. Phys._ 1985, _82_ , 270–283
  * Perdew 1986 Perdew, J. P. Density-functional approximation for the correlation energy of the inhomogeneous electron gas. _Phys. Rev. B_ 1986, _33_ , 8822
  * Andrae et al. 1990 Andrae, D.; Häußermann, U.; Dolg, M.; Stoll, H.; Preuß, H. Energy-adjusted ab initio pseudopotentials for the second and third row transition elements. _Theor. Chim. Acta_ 1990, _77_ , 123–141
  * Schäfer et al. 1992 Schäfer, A.; Horn, H.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets for atoms Li to Kr. _J. Chem. Phys._ 1992, _97_ , 2571–2577
  * Weigend and Ahlrichs 2005 Weigend, F.; Ahlrichs, R. Energy-adjusted ab initio pseudopotentials for the second and third row transition elements. _Phys. Chem. Chem. Phys._ 2005, _7_ , 3297–3305
  * Neese 2018 Neese, F. Software update: the ORCA program system, version 4.0. _WIREs Comput. Mol. Sci._ 2018, _8_ , 1327–133
  * Grimme et al. 2015 Grimme, S.; Brandenburg, J. G.; Bannwarth, C.; Hansen, A. Consistent structures and interactions by density functional theory with small atomic orbital basis sets. _J. Chem. Phys._ 2015, _143_ , 054107
  * Fdez. Galván et al. 2019 Fdez. Galván, I. et al. OpenMolcas: From source code to insight. _J. Chem. Theory Comput._ 2019, _15_ , 5925–5964
  * Olsen 2014 Olsen, J. H. Development of quantum chemical methods towards rationalization and optimal design of photoactive proteins. Ph.D. thesis, Faculty of Science, University of Southern Denmark, Campusvej 55 Odense M, 2014
  * Stener et al. 2003 Stener, M.; Fronzoni, G.; de Simone, M. Time dependent density functional theory of core electrons excitations. _Chem. Phys. Lett._ 2003, _373_ , 115–123
  * DeBeer George et al. 2008 DeBeer George, S.; Petrenko, T.; Neese, F. Time-dependent density functional calculations of ligand K-edge X-ray absorption spectra. _Inorg. Chim. Acta_ 2008, _361_ , 965–972
  * Besley and Noble 2007 Besley, N. A.; Noble, A. Time-dependent density functional theory study of the X-ray absorption spectroscopy of acetylene, ethylene, and benzene on Si(100). _J. Phys. Chem. C_ 2007, _111_ , 3333–3340
  * 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
  * Dyall and Gomes 2009 Dyall, K. G.; Gomes, A. S. P. Revised relativistic basis sets for the 5d elements Hf–Hg. _Theor. Chem. Acc._ 2009, _125_ , 97
  * Bast et al. 2009 Bast, R.; Saue, T.; Henriksson, J.; Norman, P. Role of noncollinear magnetization for the first-order electric-dipole hyperpolarizability at the four-component Kohn-Sham density functional theory level. _J. Chem. Phys._ 2009, _130_ , 024109
  * Visscher 1997 Visscher, L. Approximate molecular relativistic Dirac-Coulomb calculations using a simple Coulombic correction. _Theor. Chem. Acc._ 1997, _98_ , 68–70
  * 109 AMFI, an atomic mean-field spin-orbit integral program (1999), Stockholm University, Sweden. Written by B. Schimmelpfennig
  * Brorson et al. 1996 Brorson, M.; Dyxenburg, M. R.; Galsbøl, F.; Simonsen, K.; Liptay, G.; Trautwein-Fritz, K.; Stackhouse, J.; Ciglic, M. I.; Haugg, M.; Trabesinger-Rüf, N.; Weinhold, E. G. Ligand field analysis of the electronic spectra of [M(H2O)6]3+, [M(NH3)6]3+, [M(en)3]3+, [M(chxn)3]3+ and [M(tacn)2]3+ complexes (M = Co${}^{\text{III}}$, Rh${}^{\text{III}}$, Ir${}^{\text{III}}$). Spectrochemical strength and nephelauxetic effect of aqua, ammine and amine ligands. _Acta Chem. Scand._ 1996, _50_ , 289–293
  * Gajhede et al. 1993 Gajhede, M. S.; Skov, K.; Købberøe, L. The preparation and crystal structure of fac-[Ir(tacn)Cl3] (tacn = 1,4,7-triazacyclononane) and the preparation of Cs[Ir(H2O)6](SO4)2 $\cdot$ 6H2O. _Acta Chem. Scand._ 1993, _47_ , 271–274
  * Carrera et al. 2007 Carrera, F.; Torrico, F.; Richens, D. T.; Muñoz-Páez, A.; Martínez, J. M.; Pappalardo, R. R.; Marcos, E. S. Combined experimental and theoretical approach to the study of structure and dynamics of the most inert aqua ion [Ir(H2O)6]3+ in aqueous solution. _J. Phys. Chem. B_ 2007, _111_ , 8223–8233
  * Rubio et al. 2008 Rubio, M.; Serrano-Andrés, L.; Merchán, M. Excited states of the water molecule: Analysis of the valence and Rydberg character. _J. Chem. Phys._ 2008, _128_ , 104305
  * Harczuk et al. 2015 Harczuk, I.; Vahtras, O.; Ågren, H. Frequency-dependent force fields for QMMM calculations. _Phys. Chem. Chem. Phys._ 2015, _17_ , 7800–7812
  * Nørby et al. 2017 Nørby, M. S.; Vahtras, O.; Norman, P.; Kongsted, J. Assessing frequency-dependent site polarisabilities in linear response polarisable embedding. _Mol. Phys._ 2017, _115_ , 39–47
  * Wildman et al. 2019 Wildman, A.; Donati, G.; Lipparini, F.; Mennucci, B.; Li, X. Nonequilibrium environment dynamics in a frequency-dependent polarizable embedding model. _J. Chem. Theory Comput._ 2019, _15_ , 43–51
  * Cacelli et al. 1991 Cacelli, I.; Carravetta, V.; Rizzo, A.; Moccia, R. The calculation of photoionisation cross sections of simple polyatomic molecules by L2 methods. _Phys. Rep._ 1991, _205_ , 283–351
  * Ekström and Norman 2006 Ekström, U.; Norman, P. X-ray absorption spectra from the resonant-convergent first-order polarization propagator approach. _Phys. Rev. A_ 2006, _74_ , 042722