Environment Effects on X-Ray Absorption Spectra With Quantum Embedded Real-Time Time-Dependent Density Functional Theory Approaches.

In this work we implement the real-time time-dependent block-orthogonalized Manby-Miller embedding (rt-BOMME) approach alongside our previously developed real-time frozen density embedding time-dependent density functional theory (rt-TDDFT-in-DFT FDE) code, and investigate these methods' performance in reproducing X-ray absorption spectra (XAS) obtained with standard rt-TDDFT simulations, for model systems comprised of solvated fluoride and chloride ions ([X@ ( H 2 O ) 8 - , X = F, Cl). We observe that for ground-state quantities such as core orbital energies, the BOMME approach shows significantly better agreement with supermolecular results than FDE for the strongly interacting fluoride system, while for chloride the two embedding approaches show more similar results. For the excited states, we see that while FDE (constrained not to have the environment densities relaxed in the ground state) is in good agreement with the reference calculations for the region around the K and L1 edges, and is capable of reproducing the splitting of the 1s1 (n + 1)p1 final states (n + 1 being the lowest virtual p orbital of the halides), it by and large fails to properly reproduce the 1s1 (n + 2)p1 states and misses the electronic states arising from excitation to orbitals with important contributions from the solvent. The BOMME results, on the other hand, provide a faithful qualitative representation of the spectra in all energy regions considered, though its intrinsic approximation of employing a lower-accuracy exchange-correlation functional for the environment induces non-negligible shifts in peak positions for the excitations from the halide to the environment. Our results thus confirm that QM/QM embedding approaches are viable alternatives to standard real-time simulations of X-ray absorption spectra of species in complex or confined environments.

1 Introduction

X-ray absorption spectroscopy (XAS) is a powerful technique to probe the structural and electronic properties of molecules from an atomistic picture, since the absorbing photons in the X-ray energy range promote excitations of the core electrons to unoccupied or continuum states. The resulting absorption peaks are called edges in XAS and are labelled according to the origin of the core state, for instance K edge for 1s, L1 for 2s, L2 edge for 2p1/2 and L3 edge for the 2p3/2. The spectra features near these edges are called X-ray absorption near-edge structures (XANES). Both the energy range and the spectral shapes directly provide information on the oxidation state, local symmetry, and coordination environment of a selected analyte in the gas, liquid, or solid phase (de Groot, 2001; Bunker, 2009; Bokhoven and Lamberti, 2016; Zimmermann et al., 2020). For instance, K edges correspond to 1s → (n + 1)p dipole transitions, n + 1, being the first virtual p level, implying that, in a simple picture, the edge position is a direct measure of molecular valence states, thus allowing us to monitor the effect of the local environment on a given atom within an analyte. The interpretation of such environmental interplay calls for electronic structure calculations that allow us to access the atomic and molecular energy levels. More specifically, the theoretical modeling of XAS spectra implies the calculations of core-valence excitations.

Within quantum chemical methods, density functional theory (DFT)-based approaches such as time-dependent density functional theory (TDDFT) in its linear response (LR) (frequency domain) formulation, currently offer the best compromise between cost and accuracy for calculating electronic excitations (Norman and Dreuw, 2018; Besley, 2021). While a brute-force application of TDDFT to XAS would be prohibitively expensive as a large number of states (valence excitations, resonance, etc.) need to be determined before arriving at the energy regions pertaining to the core excitations, the introduction of restricted-excitation window TDDFT (in which one can restrict the calculation to access only configurations in which particular core electrons are excited (Stener et al., 2003; Besley and Asmuruf, 2010; Zhang et al., 2012)) or the complex polarization propagator approach (Ekström et al., 2006; Ekström and Norman, 2006; Jiemchooroj et al., 2007; Villaume et al., 2010; Pedersen et al., 2014; Fahleson et al., 2016; Rinkevicius et al., 2016) (from which one can obtain the spectral profiles for a given range of frequencies of the external perturbations from the imaginary part of the dipole polarizability), has allowed TDDFT to be routinely used in simulating XAS.

An alternative to the frequency domain approaches above that is gaining attention in recent years is that of the real-time formulation of TDDFT (rt-TDDFT) (Goings et al., 2017; Li et al., 2020), in which time-dependent properties (such as electronically excited states) are obtained based on integrating the time-dependent Kohn-Sham (TDKS) equations in time. While the theoretical underpinnings, strengths, and limitations in respect to accuracy are similar to traditional linear response (LR) TDDFT methods for obtaining electronic spectra, rt-TDDFT provides fully time-resolved solutions that can potentially incorporate non-linear effects, and also allows for strong external perturbations. With that, rt-TDDFT is used to compute not only spectroscopic properties including XAS (Lopata et al., 2012; Kadek et al., 2015; Li et al., 2020) but also the time and space-resolved electronic response to arbitrary external stimuli (e.g., electron charge dynamics after laser excitation) (Eberly et al., 1991; Cheng et al., 2006; Keldysh, 2017; Mokkath, 2020).

However, as soon as one wishes to treat molecules surrounded by an environment (e.g., species in solution or in otherwise confined spaces), the structural model for the system of interest might become too large to be treated with plain DFT approaches. In such cases, subsystem or embedding approaches (Gomes and Jacob, 2012) appear as a computationally efficient strategy: they allow for combinations of different levels of theory for the subsystems of interest and their surroundings, thereby reducing the overall computational cost. Furthermore, embedding approaches allow for selectively switching on or off interactions between different subsystems, and thus can offer a powerful way to understand the physics of chemistry of a particular process, when the analysis of a full (supermolecular) calculation may prove much more cumbersome to analyze. This is the case when analyzing electronically excited states of confined systems, which can involve transitions within particular subsystems as well as between the different subsystems–and in core states in particular since core states can be embedded into states representing the continuum.

There have been several propositions to couple rt-TDDFT methods to embedding approaches, perhaps the most widely used ones involving the coupling between a quantum subsystem and a classical environment (Lipparini and Mennucci, 2021), described by continuum models (Pipolo and Corni, 2016; Gil et al., 2019) or classical force fields (QM/MM) (Marques et al., 2003; Morzan et al., 2014; Wu et al., 2017; Parise et al., 2018). Although the obvious advantages are in cost reduction, these approaches may not properly describe specific interactions such as hydrogen bonds (for the continuum models) or rely on the availability of an appropriate classical force field (for QM/MM). Classical approaches in any case will be limited in their ability to properly describe phenomena in which a quantum description of the environment is important (such as charge delocalization, coupled excitations, or excitations across many parts of the systems not confined to a small fragment). The alternative in this case is to use quantum embedding theories (QM/QM) (Gomes and Jacob, 2012; Jacob and Neugebauer, 2014; Wesolowski et al., 2015; Sun and Chan, 2016; Goez et al., 2018), and among the fully quantum mechanical approaches to include environmental effects in the molecular response property, we note the family of subsystem DFT approaches (Jacob and Neugebauer, 2014; Wesolowski et al., 2015), to which the frozen density embedding (FDE) scheme is a member. It corresponds to a partitioning of a given system into a set of subsystems that can be, for instance, all represented within the Kohn-Sham framework, which interact through a local embedding potential. A subsystem DFT formulation of the real-time methodology (rt-TDDFT-in-DFT) has been presented in a seminal work by Pavanello et al. together with its formulation within the FDE framework (Krishtal et al., 2015; Krishtal and Pavanello, 2016; Kumar et al., 2017; Genova et al., 2017).

This initial formulation, based upon plane-wave basis representations for the different subsystems, has been shown to properly capture the coupling in the response of the different subsystems, through the dependency of the time-dependent embedding potential on the time-dependent electron densities of all (or a subset of) subsystems, whenever such a coupling is of importance. It should be noted that such couplings between the response of subsystems to external perturbation can also be taken into account in a frequency-domain formulation, but at the expense of determining second (or higher) order derivatives of the interaction energy (Neugebauer, 2007; Neugebauer, 2009a; Neugebauer, 2009b; Höfener et al., 2012; König and Neugebauer, 2013; Pavanello, 2013). That said, applications of linear-response or real-time TDDFT-in-DFT showed that in many cases, the coupling between the response of different subsystems can be ignored and a so-called “uncoupled” TDDFT-in-DFT approach can yield accurate results (Gomes et al., 2008; Gomes et al., 2013; Olejniczak et al., 2017; De Santis et al., 2020a), provided the coupling between subsystems in the ground state is well described by the embedding potential representing the subsystems’ interaction. FDE-based calculation has been shown to perform well for situations in which there are no strong interactions between subsystems, such as covalent bonds. This makes it possible in general to describe interactions such as hydrogen bonds, though in certain cases the approximations intrinsic to FDE, due to the use of approximate kinetic energy density functionals (KEDFs) (Beyhan et al., 2010; Grimmel et al., 2019) in the description of the embedding potential, prevent it from accurately describing stronger non-covalent interactions (Fux et al., 2010; Bouchafra et al., 2018a). While in such cases, a pragmatic solution is to enlarge the active subsystem, that can be potentially problematic in respect to the increase of computational costs, especially if one is interested in replacing DFT by higher-level approaches such as coupled clusters to describe the subsystem of interest.

Another QM/QM family of embedding approaches closely connected to the subsystem DFT approaches mentioned above involves the use of projection operator techniques (Goodpaster et al., 2010; Goodpaster et al., 2011; Manby et al., 2012; Ding et al., 2017; Tölle et al., 2019a; Tölle et al., 2019b; Lee et al., 2019; Graham et al., 2020; Niemeyer et al., 2020; Scholz et al., 2020), and by foregoing the use of the approximate KEDFs, shows better performance in describing strong interactions. These approaches, in some variants also referred to as Manby-Miller embedding (MME), have been shown to be particularly adept at allowing the fragmentation of a particular system through covalent bonds. More recently, block-orthogonalized MME (BOMME) (Ding et al., 2017) has been introduced to alleviate issues that plagued prior MME variants (Manby et al., 2012). BOMME allows one to treat the target system with a high-level Fock matrix and the remaining degrees of freedom with a less expensive Fock matrix by reducing the quality of the basis set and exchange. A combination of BOMME with rt-TDDFT has been recently proposed by Koh et al. (Koh et al., 2017). They demonstrated that rt-BOMME can capture both intermolecular and intramolecular couplings and their induced effects, namely solvent shifts, on spectra of chromophores.

However, in that work only processes involving valence electrons have been considered. Given the interest of core spectra as a means to characterize species in complex environments, it is of great interest to explore the behavior of rt-BOMME for describing XAS. We note that the same is also the case for FDE or TDDFT-in-DFT since these have also, to the best of our knowledge, not yet been explored for core excitations.

Thus, the main goal of this work is to describe the first investigation of the performance of rt-TDDFT-in-DFT and rt-BOMME for core excitations. To do so, we have extended our recently developed Psi4Numpy-based rt-TDDFT-in-DFT to implement the BOMME and rt-BOMME methods. As discussed below, in this manuscript we shall focus on the K and L1 edges of hydrated halides, since halogenated species and their interaction with species in solution or at interfaces are of particular interest in atmospheric sciences (Finlayson-Pitts, 2013; Pillar-Little et al., 2013; Simpson et al., 2015; Kong et al., 2017; Finlayson‐Pitts, 2019; Bartels-Rausch et al., 2021; Yu and Li, 2021). Here, however, in order to simplify our discussion we have considered relatively simple model systems representing the first hydration shell of the halides ( and ) that nevertheless can gauge the ability of the different embedding methods to describe interactions of varying strengths between halides and their environment. Also, due to the scarce experimental data for XAS on such systems, we shall restrict ourselves to comparisons of two limiting cases: the free ions and the rt-TDDFT calculations on the supermolecular system (which then serve as our benchmarks).

2 Materials and Methods

2.1 Theoretical Background

The Frozen Density (FDE) (Wesolowski and Warshel, 1993; Gomes and Jacob, 2012; Jacob and Neugebauer, 2014; Wesolowski et al., 2015) and Block-Orthogonalized Manby-Miller embedding (BOMME) approaches (Ding et al., 2017), and their extension to the rt-TDDFT framework have been described in previous works (Koh et al., 2017; De Santis et al., 2020a). In this section, after brief recapitulation of the rt-TDDFT method, we will outline analogies and differences of rt-FDE and rt-BOMME approaches.

In rt-TDDFT, the one-electron density matrix (t) representing, in the algebraic approximation, the time-dependent electron density evolves in time according to where (t, t 0) is the matrix representation of the time-evolution operator:

The real-time approach is based on the repeated application of Eq. 1 over a discretized time-domain. Time discretization allows us to devise advantageous representations of to be employed in real computer codes. In this work, we employ the exponential midpoint ansatz, which has been successfully employed in the study of valence and core excitations (Lopata and Govind, 2011; Lopata et al., 2012) Extensive discussion on the computational strategies employed to carry out the time-evolution propagation can be found in the seminal work by Castro and co-workers (Castro et al., 2004). In rt-TDDFT, the Fock matrix is defined as where 0 represents the one-electron operator while is the two-electron term c being the fraction of the Hartree-Fock exchange in the exchange correlation potential V . It is worth noting the Fock matrix appearing in Eq. 3 has an implicit time-dependence due to a time-dependent density matrix, and the explicit time-dependence due to the external potential v ext(t).

An embedding mean-field approach is based on the mapping of two different domains within the total system, into two different-quality levels of theory to be employed in each domain. This can be realized by assigning a high-level Fock matrix to the subsystem to be treated accurately, while letting the remaining part be described by a low-level Fock matrix in a reduced basis set. The block-orthogonalized (BO) partitioning scheme proposed by Ding and co-workers (Ding et al., 2017) relies on a projected basis in place of the conventional atomic-orbital (AO) partitioning to define the high- and low-level components of the system. Such a scheme proved to be suitable to remove the artifacts related to the embedding scheme while keeping the expression of the low-level Fock as simple as:

In Eq. 5, quantities expressed in the block-orthogonalized basis are denoted by tildes, and is the transformation matrix from the non-orthogonal AO basis set to the BO basis set:

The sub-blocks appearing in the transformation matrix are the identity matrices AA and BB having dimensions of n and n , mapping subsystem A and B basis sets respectively, and the projection matrix , in which is the AO overlap between the subsystems. Here and hereafter the AA (BB) block denote the subsystem with high- (low-)level theory. The Fock matrix in the BO basis reads as:

In this context different schemes for the calculation of the exchange term (in High) are available. Following Koh et al. (2017), we adopted the simplest scheme for which takes into account only the exact exchange interaction within the AA block:

In the Frozen Density formulation of DFT the entire system is partitioned into N subsystems, and the total density ρ tot( ) is represented as the sum of electron densities of the various subsystems [i.e., ρ ( ) (a = 1, ‥, N)]. In this work we restrict our consideration to a simplified model in which the total density is partitioned in only two contributions as

The total energy of the system can then be written as with the energy of each subsystem (E [ρ ], with i = I, II) given according to the usual definition in DFT as

In the above expression, is the nuclear potential due to the set of atoms which defines the subsystem and is the related nuclear repulsion energy. T [ρ ] is the kinetic energy of the auxiliary non-interacting system, which is, within the Kohn-Sham (KS) approach, commonly evaluated using the KS orbitals. The interaction energy is given by the expression: with and as the nuclear potentials due to the set of atoms associated with subsystems I and II, respectively. The repulsion energy for nuclei belonging to different subsystems is described by the term. The non-additive contributions are defined as: with X = E xc, T . These terms arise because both exchange-correlation and kinetic energy, in contrast to the Coulomb interaction, are not linear functionals of the density.

The electron density of a given fragment (ρ I or ρ II in this case) can be determined by minimizing the total energy functional (Eq. 10) with respect to the density of the fragment while keeping the density of the other subsystem frozen. This procedure is the essence of the FDE scheme and leads to a set of Kohn-Sham-like equations (one for each subsystem) which are coupled by the embedding potential term , that carries all dependence on the other fragment’s density. In the framework of FDE theory, is explicitly given by where the non-additive exchange-correlation and kinetic energy contributions are defined as the difference between the associated exchange-correlation and kinetic potentials defined using ρ tot( ) and ρ I( ). It is worth noting that only the density for the total system is available so that potentials requiring KS orbitals as input are excluded.

For the exchange-correlation potential, one may make use of accurate density functional approximations and its quality is therefore similar to that of ordinary KS. The potential for the non-additive kinetic term ( , in Eq. 15) is more problematic as it relies on less accurate orbital-free kinetic energy density functionals (KEDFs).

In this context, the Thomas-Fermi (TF) kinetic energy functional (Thomas, 1927) or the GGA functional PW91k (Lembarki and Chermette, 1994), are customarily employed. The reader interested in applicability and shortcomings of the functionals associated to term can refer to Fux et al. (2010) and references therein.

In general, the FDE scheme provides a set of coupled equations for the subsystems that have to be solved iteratively. Typically, the “freeze-and-thaw” (FnT) procedure is employed, meaning that the electron density of the active subsystem is determined, keeping the electron density of the other subsystems frozen, and is then frozen when the electron density of the other subsystems is worked out. The subsystems’ densities are converged by repeatedly applying the procedure.

We conclude this section highlighting the main differences between FDE and BOMME approaches. The FDE approach employing an explicit embedding potential allows us to optimize the subsystem of interest limiting the basis set to the sole “active” basis subset. We have already mentioned that the embedding potential relies on the KEDFs, which are in general less accurate than the exchange-correlation counterpart. On the contrary, in BOMME formulation, the self-consistent calculation is carried out in the supermolecular basis, and the embedding is handled implicitly in the calculation. As far as the coupling between subsystems is concerned, the FDE approach is trivial when attempting to estimate the interaction energy of the subsystem and eventually evaluate the net effect of the environment polarization performing an unrelaxed calculation (keeping the environment frozen). It is important to note that since the total density is obtained as the sum of subsystem densities, the partitioning reflects the mean values of observables. In the BOMME approach, the high-level system (AA block) and its environment (BB block) are optimized on the same footing, thus disentangling them could result in a cumbersome procedure. Nevertheless it could be possible to investigate the contribution of the different domains to the overall value of an observable using localization techniques.

2.2 Computational Details

In the ADF (te Velde et al., 2001) calculations, all of which were performed with version 2019.307 (Autschbach et al., 2019), we have employed the AUG/ATZP basis sets for the halogens, and the single-z without polarization (SZ) basis set for the water molecules (Van Lenthe and Baerends, 2003). In supermolecular calculations, we employed the B3LYP functional. The ADF FDE and FnT calculations were performed via the PyADF scripting framework (Jacob et al., 2011). The halogen subsystem has been calculations with the B3LYP functional and the water molecules with BLYP. Since the use of different density functionals for different subsystems in an FnT calculation is currently not possible from within the ADF implementation (Jacob et al., 2008), a PyADF script to carry out such calculations is provided as part of the dataset accompanying this manuscript (De Santis et al., 2021a), and in this case we employed a convergence criteria on the energy of 1 × 10–6. In all cases, the Thomas-Fermi and BLYP functionals have been used to calculate the non-additive kinetic energy and exchange-correlation contributions to the embedding potentials, respectively. We employed supermolecular integration grids of normal (6.0) accuracy in all calculations.

In the (rt-)BOMME and (rt-)FDE calculations in the Psi4Numpy (Smith et al., 2018) framework, we employed version 1.3.2 of the Psi4 code (Smith et al., 2020) as a computational backend. We have employed the equivalent functionals for the ADF calculations for the halogen and water subsystems (B3LYP and BLYP, respectively). As for basis sets, we employed aug-cc-pVTZ (Dunning, 1989; Kendall et al., 1992; Woon and Dunning, 1994) and STO-3G (Hehre et al., 1969) basis sets for the halogen and the water cluster, respectively. We note that for the Psi4 calculations the basis sets employed are those provided by the code’s own basis set library.

For the real-time simulations, the electronic ground state of the halogen-water complex, calculated in absence of an external electric field, was perturbed by an analytic δ-function pulse with a strength of κ = 5.0 × 10–4 a.u. along the three directions, x, y, z. The induced dipole moment has been collected for 56,000 time steps with a length of 0.025 a.u. per time step, corresponding to 33.9 fs of simulation. The choice of such a fine-grained time grid ensures in principle an observable frequency up to 3419.5 eV in the power spectrum distribution. In the case of the fluorine-water complex, the near-edge structure is located in the range of 665–700 eV, thus the time-dependent dipole moment was down-sampled halving the amount of sampling points. The use of Padé approximant-based Fourier transform allowed us to further reduce the length of the signal to be sampled corresponding to an “effective” dipole moment of 24 and 29 fs for the fluorine- and chlorine-water complex respectively. In both cases prior to Fourier transformation an exponential damping e − with λ = 3.0 × 10–4 was applied.

The code implementing the rt-FDE in the Psi4Numpy framework used in this work is part of the PyBertha package (De Santis et al., 2020a; De Santis et al., 2020b; De Santis et al., 2021b) (revision 3c752072). The code implementing the (rt-)BOMME approach is under version control (Git) but does not yet have a public release version (one is envisaged for 2022). The simulations described in the paper have been carried out with revision 3c4c334b.

The structures employed in the calculations were taken from the structures generated by Bouchafra et al. (Bouchafra et al., 2018a; Bouchafra et al., 2018b) for halogens in 50-water droplets–in particular, snapshot 619 for chloride and snapshot 1 for fluoride–and, for reasons of computational cost, we have only kept the nearest 8 water molecules that correspond to the first solvation shell. This setup is exemplified in Figure 1.

FIGURE 1

Structure for snapshot 1 of the fluoride-water droplet system taken from Bouchafra et al. (Bouchafra et al., 2018a; Bouchafra et al., 2018b) (A) and the model used in this work (B), in which the halide and the eight nearest water molecules making up its first solvation shell were extracted from the 50-water droplet.

We provide as part of the Supplemental Information a comparison of the effect of replacing the aug-cc-pVTZ with the aug-cc-pCVTZ basis set (for BOMME and supermolecule calculations), as well as performance metrics for the real-time simulations.

3 Results and Discussion

We proceed now to the presentation and discussion of our results. Before doing so, we recall that since we are interested in the relative performance of the embedding methods with respect to a calculation on the whole system, and not in a comparison to the experiment, we have opted to disregard both relativistic effects and statistical sampling of different solute configurations (e.g., by considering different snapshots from molecular dynamics simulations) as done by some of us in Bouchafra et al. (Bouchafra et al., 2018a), which we aim to consider in subsequent work. Second, we chose to focus only on transitions from the core s orbitals of the halogens, that is, the 1 s (K edges) for F− and Cl−, and the 2 s (L1 edge) of Cl−, since they provide sufficient information for our method comparisons.

3.1 Ground States

Before investigating the outcome of the real-time propagation of the wavefunctions, it is instructive to analyze the differences between the different models: isolated atoms, embedding approaches (FDE and BOMME), and standard (supermolecular) DFT calculations. To this end, we shall focus on the comparison of core orbital energies, on the one hand, since there is a direct connection between them and how environment effects are incorporated (see discussion on theoretical approaches in Section 2.1), and, on the other hand, their values provide an approximation to the ionization potentials–though the very important effects of wavefunction relaxation will still be missing due to the creation of the core hole.

While core orbitals are naturally rather localized, they are nevertheless quite sensitive to changes in the surroundings of the atom due to the presence of the solvent molecules, as we can see from the comparison of values for the isolation anions and the supermolecular systems. With the embedding approaches we expect the orbital energies to be much closer to the supermolecular values, since they introduce the different interactions (electrostatic, kinetic energy, and exchange-correlation) between the halides and the water molecules, albeit in more or less approximated manners. Consequently, the closer an embedding approach yields orbital energies to supermolecular ones, the better suited it can be considered to replace the supermolecular calculation.

Before we can proceed to a comparison between the FDE and BOMME results shown in Table 1, we should note that our Psi4Numpy-based code, in which both embedding approaches are implemented, does not yet implement the “freeze-and-thaw” (FnT) procedure for the FDE case. While that posed no problem in its first application to the rt-TDDFT-in-DFT simulation of a neutral system (De Santis et al., 2020a), prior work by some of us (Gomes et al., 2008; Gomes et al., 2013; Bouchafra et al., 2018a) has shown that for charged systems such as those considered here, the manner in which the environment density has been constructed is important, and that a relaxation of both the active subsystem and the environment densities via FnT can improve the results over a pure FDE calculation, in which the environment’s density and electrostatic potential have been obtained in the absence of the halides.

TABLE 1

Orbital energies (ϵ, in eV) of the core s orbitals the halogen atom in the clusters, (X = F−, Cl−), obtained with different models: (0) isolated halogen atoms; (1) DFT-in-DFT without the relaxation of the solvent (FDE); (2) DFT-in-DFT with the relaxation of the solvent environment (FnT); and (3) block-orthogonalized Manby-Miller embedding (BOMME). In addition to the energies obtained with embedding, we provide energy differences with respect to reference supermolecular DFT calculations (Δϵ = ϵ sup − ϵ model), represented by Δϵ 0, Δϵ 1, Δϵ 2, and Δϵ 3, respectively. Due to the fact that for technical reasons, in our FDE implementation based on the Psi4Numpy framework, we are currently not able to perform an embedding scheme (2), we provide results for (1) and (2) obtained with the ADF code.

Orbital energies (eV)
Framework	X	Orbital	Iso. (0)	FDE (1)	FnT (2)	BOMME (3)	Δϵ 0	Δϵ 1	Δϵ 2	Δϵ 3
ADF	F−	1s	−659.67	−661.16	−661.59		−2.57	−0.67	0.46
	Cl−	1s	−2753.62	−2755.05	−2755.48		−1.55	−0.12	0.32
		2s	−248.46	−249.92	−250.36		−1.59	−0.12	0.31
Psi4Numpy	F−	1s	−659.88	−661.96	−662.13		−2.42	−0.34		−0.17
	Cl−	1s	−2754.83	−2756.28	−2756.25		−1.50	−0.05		−0.08
		2s	−248.78	−250.27	−250.25		−1.55	−0.06		−0.08

To estimate the effect of relaxing the environment on the orbital energies, and indirectly part of its influence on the simulation of core spectra (the other part coming from the effect on the halide virtual orbitals), we also present in Table 1, results obtained with the ADF code, in which both FDE and FnT calculations have been carried out. We see that the FDE calculations tend to overestimate the effect of the environment, via overall more attractive embedding potentials, reflected in lower orbital energies than the supermolecular case, whereas FnT reverses this trend but overcorrects somewhat and yields energies which are slightly higher than the supermolecular ones. For fluoride, FDE and FnT differ by roughly 1.2 eV, whereas for chloride there is a much smaller important difference, of around 0.4 eV. In addition to being smaller in magnitude, the shift for chloride is roughly the same for both 1 and 2 s orbitals, an observation that is consistent with prior work (Gomes et al., 2013) in which we observed that the embedding potential shifted orbital energies in a nearly constant manner across different occupied orbitals.

Comparing the differences between FDE and supermolecule results between ADF and Psi4Numpy, we see a similar trend in that FDE overestimates the effect of the environment. From the comparison of Δϵ 0 for the two codes, we see that discrepancies of around 0.15 eV (for fluoride) and 0.05 eV (for chloride) can be attributed to differences inherent to the two sets of calculations (Slater vs. Gaussian basis sets, etc.), with values calculated with ADF showing larger discrepancies between isolated and supermolecular calculations than Psi4Numpy. If we correct for these differences, we see that for chloride the Δϵ 1 values are consistent between codes, though for fluoride even taking into account such corrections, non-negligible differences between codes, of around 0.15 eV, remain. From this comparison, we believe that we can conclude that, if we were able to carry out such calculations, the Psi4Numpy FnT Δϵ 2 would likely be of around 0.2–0.3 eV for chloride, and 0.4–0.5 eV for fluoride.

The BOMME results show a similar trend to the FDE in overestimating the effect of the environment with respect to supermolecular results, and that such an overestimation is larger for fluoride than for chloride. The magnitude of such an effect, however, is about half of that of FDE for fluoride (-0.17 eV vs -0.34 eV), and roughly equivalent to that of FDE (less than -0.1 eV) for chloride. The differences between BOMME and FDE are consistent with what is known in the literature between the more reliable behavior of projection-based embedding (such as BOMME) in describing cases in which there are strong interactions between the different subsystems with respect to FDE, which suffers from the limited accuracy of the non-additive kinetic energy density functionals used to calculate the non-additive kinetic energy contribution to the embedding potential (Gomes and Jacob, 2012).

From the discussion above, and assuming that the dominant effects in the electronic spectra would come mainly from the energy differences between the core and low-lying virtuals either on the halogen (for both BOMME and FDE) or on the environment (for BOMME), we can expect to see that BOMME excitation energies would be consistently closer to the supermolecular results than FDE, but that such a difference would decrease for chloride. In the following we shall see to what extent this picture holds true. In any case, for the core orbital energies, the FDE results seem to provide a fortuitous error cancellation that places the relatively simple FDE model on par with the much more sophisticated BOMME.

3.2 Core Excited States

Before discussing the behavior of the different approaches for the core states, we note that in the following we shall focus on the edge region for the K edges of both systems, considering energies spanning a somewhat broad window (20–30 eV higher than the first peaks, in order to have a wider region in which to compare the different models) and the L1 edge of chloride. Furthermore, in the discussion below we shall focus on combined contributions from the x, y, z components of the perturbing field. We present a breakdown of these by component of the perturbing field, along with the spectra for the whole regions under consideration (Figures 2, 3, 5) in the Supplemental Information.

FIGURE 2

Simulated K-edge spectra for the fluoride model system, over the roughly 30 eV interval starting at the free ion edge peak. It should be noted that here the peak heights (in arbitrary units) for each family of models: free ion (=iso), FDE, BOMME, and supermolecule (super) have been scaled, with a height of 1 assigned to the most intense transition.

FIGURE 3

Simulated K-edge spectra for the chloride model system, over the roughly 30 eV interval starting at the free ion edge peak. It should be noted that here the peak heights (in arbitrary units) for each family of models: free ion (=iso), FDE, BOMME, and supermolecule (=super) have been scaled, with a height of 1 assigned to the most intense transition.

FIGURE 5

Simulated L1-edge spectra for the fluoride model system in the edge region (contrary to the K edge, no peaks of appreciable intensity have been observed at higher energies). It should be noted that here the peak heights (in arbitrary units) for each family of models: free ion (=iso), FDE, BOMME, and supermolecule (super) have been scaled, with a height of 1 assigned to the most intense transition.

3.2.1 K Edge

The spectra for the K edges of fluoride and chloride are shown in Figure 4. Starting with the simplest systems, for the isolated anions, we note that, as expected, the K-edge spectra corresponds to transitions from the 1 s orbitals to the first virtual halide p orbitals ((n + 1)p). The second peak in the energy range considered, on the other hand, corresponds to transitions from the 1 s to higher-lying halide p-type orbitals ((n + 2)p).

FIGURE 4

Details on the three main energy ranges for the K-edge spectra of F− (A,C,E) and Cl− (B,D,F) model systems selected from the spectra shown in Figures 2, 3.

Second, at the other extreme we have the supermolecular calculations on the microsolvated anions (which here serve as a benchmark to which the embedding approaches will be compared). The first remarkable difference from the free ions is that there is an environment-induced shift in the first region (Figures 4A,B), which at around 1 eV is fairly similar between systems, but about half of what would be expected from the difference in orbital energies between the free ions and the supermolecular systems (Δɛ 0). This is the first indication that core orbital energies are useful for understanding the K-edge absorption spectra from a semi-quantitative viewpoint at best.

The asymmetrical first hydration shell environment also breaks the atomic symmetry, which has as consequence the lifting of the selection rules for the atom, and introduces differential interactions with the different p orbitals that become occupied in the excited states. As a result of that, for both fluoride and chloride, in the supermolecular calculations we observe four transitions within roughly a 1 eV interval, with spacings of around 0.5 eV between the first three peaks (with the fourth being much closer to the third).

For the region corresponding to the second free ion transition (Figures 4E,F), we observe a similar situation to that of the first, with the environment inducing a symmetry breaking of the higher-lying p orbitals. In the case of fluoride, for the supermolecular calculations we observe five peaks, one rather close to that of the free ion (around 691.5 eV), followed by two other peaks around 2 eV higher (at 694 eV), and two more peaks between 697 and 699 eV. In the case of chloride, we also observe five peaks, one nearest to that of the free ion, another peak around 2775.5 eV, and then three additional peaks between 2778 and 2780 eV. Finally, in between the two free ion peaks (Figures 4C,D), in the supermolecular calculation we have a region that contains several peaks.

Considering now the FDE calculations–and recalling that these correspond to a situation in which the density of the environment has not been relaxed in the presence of the anions–we see, for fluoride, a semi-quantitative agreement with the supermolecular calculation; for the lower energy region (Figures 4A,B), the first peak appears in slightly (around 0.2 eV) lower energies than the supermolecular ones, while the second, third, and fourth peaks appear at slightly higher energies. For chloride, the peak positions are overall closer (0.1 eV or less) to the supermolecular one than for fluoride, but now the energies of the first three peaks are slightly overestimated with respect to the supermolecule, while in the fourth we see a slight underestimation. We note that this is in line with the better agreement between FDE and supermolecular 1 s orbital energies than for fluoride. Such a tendency was also observed by Bouchafra et al. (Bouchafra et al., 2018a), though for valence ionizations. It was shown that an FDE model in which only the halide belonged to the active subsystem was not a good representation for the solvated ion, due to the strong water-fluoride interactions in the valence regions, whereas for chloride (and other heavier halides) this simple FDE model containing no explicit halide-water interactions was quite good. In the higher energy region (Figures 4E,F), we see that the peaks from FDE calculations are also close to that of the free ion and of the supermolecule (around 2774 eV), though for FDE we observe another two peaks, just over 2774 eV (that show a very small splitting), and four others between 2777 and 2778 eV. The behavior of FDE for this higher energy range is, therefore, in stark contrast to the lowest energy range considered, since there is not even qualitative agreement with the supermolecular results.

Having in mind that in the supermolecular case the complete system is allowed to respond to the external perturbation, but that by construction the response of the environment is absent in the FDE case, this discrepancy provides the first indication of the importance of the response of the environment for higher energies. This is further underscored by the fact that the FDE calculations show no peaks in the intermediate energy range considered. Consequently, we can safely say that the intermediate energy range is, in effect, dominated by excitations from the halide to virtuals with strong (if not dominant) contributions from the environment. Furthermore, while the behavior of FDE is in line with the difference between the isolated and FDE orbital energies, in particular for the lower energy part of the spectra, the situation is less clear-cut with respect to a comparison to the supermolecule. We consider the discrepancies in this case to partly arise from the lack of relaxation for fluoride virtual orbitals and partly from the lack of coupling between the response of the subsystems, as discussed below.

Now comparing supermolecular and BOMME calculations, we see that for both systems the latter provides an overall improvement over FDE–already in qualitative terms, with BOMME we are able to capture the contributions to the solvent to the different electronic states, and, furthermore, in all energy ranges BOMME systematically approaches (underestimates) the excitation energies. In quantitative terms, for fluoride, BOMME clearly performs better than FDE; this can already be seen from the orbital energy differences, and in the low-energy range, the differences in absolute are not very large but BOMME does show smaller differences. For the larger energy range, where FDE is not even qualitatively correct, BOMME shows discrepancies of around 0.1 eV or less.

For the intermediate energy range, on the other hand, we see more significant differences between the BOMME and supermolecular energies. We attribute this to the use of a GGA functional for the environment, since GGAs tend to underestimate excitation energies with respect to hybrid functionals (Besley, 2021), and it is precisely in this region that contributions from the environment have a prominent role. While the goal of BOMME is to replace a high-level description of the environment for a lower-level one, our results are the first indication that for core excited states in which the environment plays an important role, the quality of the low-level of theory may matter much more than for valence states.

The trends outlined above for fluoride are also seen for chloride; if for the low-energy region BOMME does not bring about as significant an improvement over FDE as for fluoride, for the other two energy ranges BOMME provides a semi-quantitative agreement due to its systematic behavior. We see, however, that for the higher energy range, we have often energy differences between BOMME and the supermolecule in the order of 1.0 eV, and nearly so for many of the states in the intermediate region. This further underscores the importance of the quality of the density functionals employed for core energies, and in particular deeper cores such as the chloride K edge.

We can now compare these results to what one could expect, the simple argument put forward above that the halide orbital energies would provide a dominant contribution to the excited states. From the orbital energies alone, we would expect that both BOMME and FDE would yield excitation energies larger than supermolecular ones; considering only the low-energy range, for both halides this is at odds with the BOMME results, which always show lower energies than the supermolecule, but the simple orbital picture is more consistent with the FDE results, since for three peaks out of four they appear at slightly higher excitation energies (for fluoride, the first FDE peak appears at a lower energy than the supermolecule, whereas for chloride that happens for the fourth peak). We consider this is yet another piece of evidence of the importance of the virtual orbitals from the environment to characterize the excitation energies obtained with BOMME, since such contributions are absent by construction in the case of FDE.

To deepen the discussion on the origin of differences between embeddding and supermolecular results, it is useful to analyse the molecular orbitals involved in the K edges. In the real-time framework, similarly to LR, the absorption cross-section can be interpreted in terms of occupied-virtual molecular orbital (MO) pairs. This approach has been proposed originally by Repisky and coworkers (Repisky et al., 2015) and implemented in the relativistic code ReSpect (Repisky et al., 2020). The scheme has been slightly reworked by Bruner et al. (Bruner et al., 2016) accelerating the methodology applying the Padé approximants to the Fourier transform of the deconvolution of the induced dipole into molecular orbital pairs. The method gives the MO contribution to the dipole strength function at all frequencies. The relative areas under each peak correspond to MO contribution to an excitation at that frequency, which gives a representation of transitions consistent with the linear response.

In order to have a more visual interpretation of transitions, one can return to transition density plots. We recall that the simplest direct approach (Hofmann and Kümmel, 2012) is to evaluate the transition density at a preset frequency according to: where is the Fourier transform of the time-dependent (TD)-induced density δρ( , t) = ρ( , t) − ρ( , 0). Many research groups have recently contributed to this topic (Rossi et al., 2017; Schelter and Kümmel, 2018; Sinha-Roy et al., 2018; Jornet-Somoza and Lebedeva, 2019), focusing on low-frequency excitations. To the best of our knowledge, these methods have not been applied yet in the high-frequency range.

In the Supplemental Information, we can observe from the isosurface plot, for ω = 668.586 eV, that the surrounding waters contribute to the K edge. At the same time, we should point out that the Fourier transform of the TD-induced dipole (using GNU Octave (Eaton et al., 2020) routines) is of exceptionally poor quality in the frequency range of the fluorine K edge. It can be argued reasonably that also is not of the best quality possible. With that, we have been unable to apply the methodology developed by Schelter and co-workers, in which the accurate value of the oscillator strength is extracted from the refined dipole strength function (DSF) of the TD-induced dipole moment. The lack of high quality DSF makes it difficult to provide an accurate estimation of both the transition dipole moment and transition density. This further prevents us from applying more elaborate methods, namely natural transition orbitals. These difficulties should be addressed in dedicated future works.

In the meantime, in the present work we provide in the Supplemental Information an analysis of an LR-TDDFT calculation on the supermolecular system, performed with the NWChem code (Valiev et al., 2010) (using the same basis sets and density functionals as the Psi4 calculations), from which we can determine that lower-energy K-edge transitions of indeed involve virtual orbitals with contributions from orbitals centered on the oxygen atoms, in line with the qualitative picture we have managed to extract from the isosurface plots of .

3.2.2 L1 Edge

The spectra for the chloride L1 edge region are shown in Figure 5. Unlike the case for the K edge, here we do not show a larger energy range since, for energies between 245 and 260 eV, there are no other peaks with appreciable intensity other than those in the picture. A comparison between the free ion at around 253 eV and the supermolecular peaks here show that the energy shift due to the environment is not as marked (around 0.5 eV higher) but nevertheless sufficient to clearly characterize the interaction with the waters through the splitting of the peaks. In this case, we have a near perfect agreement between the supermolecular and BOMME results in terms of peak positions, at around 253.5 eV but also for low intensity transitions around 251.5 eV. We also observe that, unlike for the K edge, here the BOMME results slightly overestimate the supermolecular ones, in agreement with the orbital energy differences (Δϵ 3) in Table 1. We consider that this point, and the absence of other peaks as in the K edge that would indicate more or fewer important contributions from the environment, the L1 edge is almost exclusively dominated by halide to halide transitions. The FDE results, on the other hand, underestimate the effect of the environment and show almost no difference to the free ion results, apart from the fact that a splitting of the peak, much less significant than that seen for BOMME or the supermolecule, is also seen. This is a further indication that the FDE approach has not properly captured the perturbations to the virtual orbitals of chloride induced by the solvent.

4 Conclusion

In this manuscript we have carried out an investigation, to the best of our knowledge for the first time, of the accuracy of fully quantum mechanical (QM/QM), DFT-based embedding approaches–namely, Frozen Density embedding (FDE) and block-orthogonalized Manby-Miller embedding (BOMME) approaches–in the description of core excitation spectra (XAS), obtained by the real-time propagation of the electron density, for model systems representing the hydration of halide ions, comprising the halide ions (fluoride and chloride) as active subsystems, and the eight water molecules in the first solvation shells as the environment.

We note that the BOMME approach and its real-time variant has been implemented within the Psi4Numpy framework in which some of us had previously implemented the real-time TDDFT-in-DFT FDE method, thereby facilitating a rigorous one-to-one comparison between approaches.

From our comparison of the two embedding methods to reference DFT calculations on the whole model systems, we observe first that the BOMME approach can better describe the fluoride core orbital energies in the fluoride-water system than FDE, due to its better handling of the stronger interactions between subsystems, while for chloride both BOMME and FDE perform rather similarly.

In the case of real-time simulations, we have found that the rt-BOMME approach follows the behavior of the reference rt-TDDFT calculations in a very systematic manner across all energy ranges investigated, and as such can potentially become very useful in the investigation of core spectra of species in confined or complex environments.

We observe that rt-BOMME tends to slightly underestimate the supermolecular results around the edge region for both K and L1 edges, where excitations mostly take place between orbitals belonging to the halides. On the other hand, we observe much more important discrepancies in higher energy regions for which, it turns out, the environment plays a more important role.

We attribute these discrepancies to the fact that in rt-BOMME the environment is described with a lower-accuracy GGA functional, a functional class which tends to underestimate core excitation energies due to larger self-interaction errors in comparison to the hybrid functionals which were employed for the active subsystem, and for the whole system in the supermolecular calculations. Our results call for particular attention, in the case of core spectra, in choosing the density functional for the environment, in order to minimize artifacts in the simulations.

The rt-TDDFT-in-DFT simulations carried out under the constraint that the density of the environment has not been relaxed, have nevertheless shown performances similar to the rt-BOMME and reference rt-TDDFT simulations for the (pre-)edge regions. However, since in our implementation the response of the environment is also lacking, large parts of the spectra are either inaccessible or are not correctly described. We intend to address this issue, and introduce the coupling between the response of the different subsystems, in subsequent work.