We investigate the performance of time-dependent density functional theory (TDDFT) for reproducing high-level reference x-ray absorption spectra of liquid water and water clusters. For this, we apply the integrated absolute difference (IAD) metric, previously used for x-ray emission spectra of liquid water [T. Fransson and L. G. M. Pettersson, J. Chem. Theory Comput. 19, 7333–7342 (2023)], in order to investigate which exchange–correlation (xc) functionals yield TDDFT spectra in best agreement to reference, as well as to investigate the suitability of IAD for x-ray absorption spectroscopy spectrum calculations. Considering highly asymmetric and symmetric six-molecule clusters, it is seen that long-range corrected xc-functionals are required to yield good agreement with the reference coupled cluster (CC) and algebraic-diagrammatic construction spectra, with 100% asymptotic Hartree–Fock exchange resulting in the lowest IADs. The xc-functionals with best agreement to reference have been adopted for larger water clusters, yielding results in line with recently published CC theory, but which still show some discrepancies in the relative intensity of the features compared to experiment.
INTRODUCTION
Being vital to all known forms of life, water has long been the subject of extensive investigations on its structure, dynamics, and behavior under different conditions.1 Through this, a number of unusual properties have been identified,2 e.g., a decrease in viscosity under pressure (when supercooled), high surface tension, and minimum in the compressibility at 46 °C.2–4 These anomalous behaviors are most pronounced under supercooled conditions, but persist even at ambient temperatures (e.g., a density increase upon melting). Despite comprehensive investigations, a consensus on their origin remains elusive, but it is expected to be found in the intricate details of the structure and dynamics of the hydrogen-bonding network.1,5 Initially believed to fluctuate around a tetrahedral motif, more recent experimental6–9 and simulation10–12 evidence suggests that water can exist in two distinct liquid forms with a liquid–liquid phase transition in the deeply supercooled and pressurized region of the phase diagram. The present debate contrasts a tetrahedral model with a continuum of distorted hydrogen bonds13,14 and the two-state model where ambient liquid water is viewed as a single-phase region beyond the termination of their coexistence line. The two-state perspective suggests a predominance of close-packing of molecules [high-density liquid (HDL)] under ambient conditions, favored by entropy, with local fluctuations into tetrahedral environments [low-density liquid (LDL)], favored by enthalpy. This model aligns with thermodynamics, offering quantitative predictions of water’s thermodynamic properties.15,16
The development of contemporary x-ray radiation sources has enabled a number of different methods of investigation, capable of investigating, for example, local molecular structure, occupied states, and unoccupied states.17–21 Due to the significant differences in core-orbital energies, particularly in the inner-core region (1s), these methods are element-specific, allowing for the exploration of local electronic and atomic structures.20,21 Included here is x-ray absorption spectroscopy (XAS), which probes the absorption cross section of the system, shedding light on unoccupied states and local structure, depending on if the measured region is around or above the ionization energy. Note that, for a more complete understanding, it is important to correlate observables from different measurements, which can provide complementary insights for understanding complicated electronic and geometric structures.22–26
The x-ray absorption spectrum of water has been investigated extensively,19,27–32 including the effects of temperature, variations in pH, isotope substitution, different solutes, as well as comparisons with the bulk and surfaces of different phases of ice and the behavior when adsorbed on surfaces. What is needed in order to reach a fundamental understanding of the structure and dynamics of the H-bonding network is mainly theoretical considerations. The spectrum of liquid water is typically divided into three main parts: a distinct pre-edge peak centered at ∼535 eV, a prominent main-edge peak around 537–538 eV, and a strong and extended post-edge peak spanning 540–542 eV.19,32,33 The interpretation of these spectroscopic features in relation to various hydrogen-bond configurations is well established, although there remains an ongoing debate concerning the extent of distortions involved.13,25,31–47 In the isolated water molecule, the lowest unoccupied molecular orbitals have a large spatial extent and they are polarized toward the hydrogen atoms. Due to this, their shape and orbital energies are sensitive to the formation of donated H-bonds, making XAS a very sensitive tool to probe the donated H-bond local structure. Molecules with asymmetrically donated H-bonds, i.e., one well-defined and the other weak or broken, contribute to the pre-edge, while tetrahedrally H-bonded molecules contribute to a strong post-edge.19,32 Many attempts have been made to simulate the water XAS spectrum using different theoretical approaches, such as transition potential density functional theory (TP-DFT),24,31,41,48–51 the Bethe–Salpeter equation (BSE) formalism,39,47 time-dependent density functional theory (TDDFT),25,50–52 and coupled cluster (CC) theory.43,50 These approaches generate the correct assignment for the general spectrum features, but none of them can fully reproduce the XAS water spectrum in all spectrum regions (pre-, main-, and post-edge).19,49 In addition, it would be very beneficial to combine several different spectroscopies on the same structures, in order to achieve a holistic view of the structure and dynamics of water.24
The modeling of x-ray absorption spectra is challenging on account of a number of factors, such as pronounced relaxation effects stemming from the creation of a core-hole, the necessity to treat valence, Rydberg, and continuum excited states in a balanced way, as well as strong relativistic effects. Still, a plethora of different methods have been developed and utilized,20 encompassing semiempirical, density-based, and wave function-based approaches.53–69 The core transitions in XAS are embedded in a continuum of valence excitations, and several different approaches in which this complication can be circumvented are used. Focusing on methods grounded in linear response theory, the weak coupling between core- and valence-excited states can be leveraged by neglecting such terms, leading to the core–valence separation (CVS) scheme.70,71 Such a strategy has been effectively implemented for various methods, including time-dependent Hartree–Fock and density functional theory (TDHF and TDDFT),54–56,72–74 equation-of-motion coupled cluster theory (EOM-CCSD),61,75 the algebraic-diagrammatic construction (ADC) scheme for the polarization propagator,63,64,70,71 and transition potential density functional theory (TP-DFT).53 Demonstrated to yield negligible errors, the CVS approximation is a reliable choice for directly targeting core-excited states.76,77 Among these methodologies, TDDFT stands out for its computational efficiency; however, it is afflicted by issues of lacking relaxation and self-interaction, which are particularly problematic in calculations of core transitions.78,79 The latter issue arises from spurious self-interactions in the final state, often resulting in x-ray absorption features manifesting at too low energies. Nevertheless, by employing scalar shifts of spectral features or utilizing exchange–correlation (xc) functionals tailored for x-ray properties,55,80 TDDFT remains a popular choice for XAS simulations,20,55–57,72 and efforts to mitigate the self-interaction error are ongoing, including the utilization of electron-affinity TDDFT81 and correction schemes based on many-body perturbation theory.73
Here, we investigate the performance of different methods for calculating x-ray absorption spectra of water, focusing on the ability of time-dependent density functional theory (TDDFT) to reproduce reference spectra obtained by high-quality calculations using the ADC-scheme for the polarization propagator and equation-of-motion coupled cluster singles and doubles (EOM-CCSD) theory. A number of different xc-functionals have been considered, evaluating their capability to replicate ab initio reference spectra for both highly structured (tetrahedral) and disordered (asymmetric) six-molecule water clusters. A measure termed the integrated absolute difference (IAD) is utilized, which probes the correspondence between two spectra.23,82–84 We previously adopted this metric to evaluate the performance of different methodologies for calculating the x-ray emission spectrum (XES) of liquid water, where a correlation between the amount of Hartree–Fock exchange and the ability to simultaneously reproduce reference spectra for highly terahedral and asymmetric clusters was noted.85 The present contribution is a continuation of that work, considering XAS and contrasting this to the results for XES.
The rest of this study is outlined as follows: we first discuss the different structures that are considered, followed by a brief review of the methods of calculating x-ray absorption spectra used here. This is followed by an overview of the IAD and by the computational details. Moving to the Results and Discussion section, we investigate the performance of the reference methods, illustrate the spectra and resulting IADs of a few select xc-functionals, and then perform a more systematic investigation and analysis of the performance of 21 different xc-functionals. Finally, the size of water clusters for XAS calculations is considered, the “best fit” xc-functionals are used to compare the results using multilevel coupled cluster theory43 and experiment,33 and the impact of diffuse basis functions is investigated.
METHODOLOGY
In the present work, we focus on the theoretical modeling of XAS by calibrating approximate methods (that can be applied to larger models) against high-level quantum chemical approaches that are restricted to smaller systems. It is important to note that it is not clear if the inability to, as of yet, obtain good theoretical x-ray absorption spectra of liquid water is due to deficiencies in the molecular dynamics simulations or in the spectrum calculations. Here, water clusters were taken from Ref. 25, where structures were characterized in terms of H-bond symmetry/asymmetry. Specifically, we chose structures representing highly tetrahedral environments (indicative of low-density local environments) and asymmetric configurations (representative of high-density local environments). The former feature relatively short (<2.68 Å) and well-defined H-bonds, while the latter exhibit longer (>2.80 Å) and asymmetric H-bonds and were considered in order to investigate the performance of TDDFT for distinctly dissimilar structural motifs, as has been seen to be important for XES.85 In addition, we conduct calculations on gas-phase water, considering the equilibrium geometry of the isolated water molecule. Finally, structures from Ref. 43 are considered, in order to compare results obtained with the “best fit” exchange–correlation functionals to spectra from multilevel coupled cluster theory.
For wave function-based spectrum calculations, single- and multireference methods can be used, where electron relaxation—resulting from the reduced screening of the nuclei due to the creation of a core-hole, as well as interactions with the excited electron—is accounted for through electron correlation by the use of (at least) doubly excited configurations.20,86,87 Available schemes include single-reference methods, such as coupled cluster (CC) theory, the algebraic diagrammatic construction (ADC) approach, and multireference methods, such as RASSCF, RASPT2, and MR-CC.60–67,87,88 Among correlated ab initio (post-Hartree–Fock) methods for excited states, we employ the algebraic diagrammatic construction (ADC) scheme for the polarization propagator and coupled cluster (CC) theory. Within ADC, the polarization propagator undergoes expansion in a hierarchical series, with the poles and residues of its spectral representation corresponding to excitation energies and transition amplitudes.89–92 This hierarchy of ADC methods is derived by truncating the perturbation expansion at a desired order, with implementations available for ADC(1), ADC(2), ADC(2)-x, and ADC(3), along with recent advancements enabling ADC(4) excitation energy calculations.93 The ADC(3) implementation utilized here is accurate to third order in transition energies but only to second order in property gradients and is thus also referred to as ADC(3/2). Alternatively, a hierarchy of post-HF methods is accessible through the CC approach, which can be explored either in an equation-of-motion (EOM) or in a linear response (LR) formalism. The results presented here are compared against reference values obtained using EOM-CCSD and ADC(2)-x, two methods known for their strong agreement with experimental data while maintaining a robust theoretical foundation.60,61,63,64,87,94–96 These methods offer accurate relative features, albeit with some element-dependent shifts in absolute energy.87 We stress that trends and relative features are often more important than absolute energies, which can simply be compensated for by shifting to experiment.
Moving from the reliable but costly post-HF methods, density functional theory (DFT) provides high computational efficiency, but at the cost of lower reliability—in particular for cases where suitable exchange–correlation (xc) functionals have yet to be identified. Originally formulated to capture both accurate densities and energies, contemporary DFT functionals often prioritize energy accuracy over density precision, leading to smaller errors in energy at the expense of larger errors in density.97 As such, the development and utilization of ab initio wave function methods for computational spectroscopy remain crucial, especially for novel processes and applications. Focusing on the use of linear response TDDFT, low computational cost is achieved, but the issue of self-interaction plagues x-ray spectrum calculations and results in substantial inaccuracies in absolute energies.20,55,78,79 This arises from spurious self-interactions, which decontract the high-density core orbitals, leading to too high orbital energies and thereby underestimating transition energies. These inaccuracies are partially canceled out by the lack of relaxation effects in TDDFT, but most xc-functionals yield too low XAS energies. Functionals tailored for use in x-ray spectroscopies have been developed, but it has been noted that better accuracy in absolute energies does not necessarily lead to better accuracies in relative energies.87 In fact, the improved absolute energies are largely a result of customized amounts of Hartree–Fock exchange, which, in practice, leads to a situation where errors due to self-interaction are canceled out by errors due to a lack of relaxation—for example, a pure generalized gradient approximation (GGA) functional tends to underestimate carbon K-edge transition energies by ∼18 eV, compared to the typical overestimation of about 9 eV by TDHF.87,98
COMPUTATIONAL DETAILS
The equilibrium geometry of the isolated water molecule was calculated at the MP2/cc-pVTZ99 level of theory, using the Q-Chem 5.2100 software package. For the water cluster calculations, six water molecules were included in each cluster, using the ten asymmetric and ten symmetric clusters from Ref. 85 (there referred to as HDL- and LDL-like), originally presented in Ref. 25. Additional spectrum calculations were performed using structures from Ref. 43, considering clusters with 32 water molecules. Unless stated otherwise, property calculations were run using the IGLO-III basis set for the central oxygen (7s6p2d) and IGLO-II for the hydrogens (3s1p).101,102 To localize the core-transition to that from the central water molecule,103 an effective core potential (ECP) of the Stuttgart–Cologne type104 was used for the surrounding oxygens, and the corresponding double-ζ basis set was converted to a triple-ζ level (3s3p1d, uncontracting the third PGTO in the first CGTO and adding a d-function with an exponent of 0.7).50 In order to improve the description of the excited states, the basis set of the central oxygen was augmented with additional 9s9p diffuse functions (even-tempered basis set with the first two exponents amounting to 1.238 and 0.632). A convolution of the obtained energies and intensities using a Lorentzian broadening function was performed to facilitate the analysis, using a half-width at half-maximum (HWHM) of 0.3 eV, except where otherwise stated. The coupled cluster and TDDFT calculations were carried out using Q-Chem 5.2,100 utilizing the frozen core (fc) implementation for CVS-EOM-CCSD61 and the Tamm–Dancoff approximation for TDDFT.105 The CVS-ADC results were obtained using the adcc software package,106 using SCF results obtained from PySCF.107,108
RESULTS AND DISCUSSION
Performance of reference methods
In order to assess the performance of the reference methods, the CVS-ADC(2)-x and fc-CVS-EOM-CCSD spectra of gas phase water and of 20 six-molecule clusters are shown in Fig. 1, compared to experiment.33 For brevity, we will henceforth denote CVS-ADC(2)-x as ADC(2)-x and fc-CVS-EOM-CCSD as CCSD. In addition, TDDFT spectra using the two exchange–correlation functionals, which give best agreement with the reference methods, are included in this figure, featuring CAM-QTP02 (shortened to QTP02) and CAM100%-B3LYP (shortened to CAM100)—see below for details. Regarding the clusters, ten disordered and ten tetrahedral structures were considered individually, representing local asymmetric and symmetric environments, respectively, and the spectra for each selection are presented. We stress that this aspect of the analysis is not aimed at reproducing experimental results, as the chosen structures are limited in number and small in scale. Instead, the clusters were selected specifically for their clearly distinct tetrahedral and asymmetric local structures, allowing for an investigation into the performance of different theoretical models across disparate structural motifs—such an analysis has previously been demonstrated to be important when considering x-ray emission spectra of liquid water.85 Therefore, the comparison with experimental measurements serves merely as an illustrative exercise. The theoretical results have been shifted to overlap with the 4a1-feature of the gas phase (∼534 eV), with an additional positive shift of 1.5 eV for the clusters. These shifts account for effects such as missing relaxation, self-interaction error (for TDDFT), basis set incompleteness, and the lack of relativistic effects (which would here shift the absolute energies by ∼0.3–0.4 eV).60,109 The additional shift for the clusters is a result of missing environment effects, which are poorly accounted for when using clusters of only six molecules. The computed gas phase spectra are overall quite similar, albeit with some difficulties in reproducing the position of the first Rydberg (2b1) state—experimentally this is found 3.0 eV above 4a1,33 but the shift is slightly lower for ADC(2)-x, and even more so for TDDFT. The theoretical spectra for the region above the ionization energy (IE) of gas phase water (539.9 eV110) have some very intense features, which are primarily a result of the discretized continuum,77 and are reduced if a more extended atom-centered basis set is used.49 The spectra of the six-molecule clusters are significantly different from that of liquid water, but are similar between the two reference methods—the asymmetric clusters yield a strong pre-edge and a slightly more pronounced main-edge, and the symmetric clusters yield a stronger post-edge, in particular for the region ∼3 eV above the post-edge.
Left: Computed x-ray absorption spectra of gas phase water, compared to experiment.33 Theoretical results have been shifted in energy by 11.22, 14.87, 0.63, and 1.16 eV for (from top to bottom) QTP02, CAM100, CCSD, and ADC(2)-x, respectively, to match experiment (bottom). Right: x-ray absorption spectra of ten asymmetric (solid line) and ten symmetric (dashed line) six-water clusters, compared to experimental measurements of liquid water.33 Theoretical results have been shifted by an additional 1.5 eV, when moving from the isolated molecule to the clusters. The experimental spectrum of liquid water is included only for illustrative purposes, as the theoretical results are obtained using small and select types of clusters; the dotted vertical lines indicate the approximate positions of experimental features.
Left: Computed x-ray absorption spectra of gas phase water, compared to experiment.33 Theoretical results have been shifted in energy by 11.22, 14.87, 0.63, and 1.16 eV for (from top to bottom) QTP02, CAM100, CCSD, and ADC(2)-x, respectively, to match experiment (bottom). Right: x-ray absorption spectra of ten asymmetric (solid line) and ten symmetric (dashed line) six-water clusters, compared to experimental measurements of liquid water.33 Theoretical results have been shifted by an additional 1.5 eV, when moving from the isolated molecule to the clusters. The experimental spectrum of liquid water is included only for illustrative purposes, as the theoretical results are obtained using small and select types of clusters; the dotted vertical lines indicate the approximate positions of experimental features.
Spectra from selected xc-functionals
In Fig. 2, we illustrate the spectra and IADs of the asymmetric and symmetric clusters, as obtained using CCSD and TDDFT with five different exchange–correlation (xc) functionals, and compared to ADC(2)-x results. The broadened spectra have been shifted to overlap with the pre-edge in the summed ADC(2)-x spectrum of the asymmetric clusters (using the same energy shift for the symmetric structures), and the IADs are calculated over the energy range in this figure. For the TDDFT results, xc-functionals have been selected to represent a range of different functional types: global hybrid (B3LYP111), general range-separated hybrid (ωB97112), and range-separated hybrids tailored for core transitions (CAM-QTP00,113 SRC1-R1,80 and CAM100%-B3LYP114). The Tamm–Dancoff approximation105 was applied for all TDDFT calculations, which only results in minor changes in spectrum features77—this was verified in additional test calculations (results not shown). For brevity, CAM-QTP0X and CAM100%-B3LYP are here referred to as QTP0X and CAM100, respectively (X = 0, 1, 2).
X-ray absorption spectra and integrated absolute differences (IADs) of EOM-CCSD and TDDFT with five different xc-functionals, compared to ADC(2)-x reference spectra. Spectra have been shifted to overlap with the pre-edge peak of the asymmetric clusters (vertical dashed line), and the IADs have been calculated over the illustrated energy region. IADs calculated for the summed spectra (large markers) and each spectrum individually (small markers). The applied energy shifts amount to, from top to bottom: 13.93, 19.74, −0.46, −1.49, 13.42, and −1.86 eV, respectively.
X-ray absorption spectra and integrated absolute differences (IADs) of EOM-CCSD and TDDFT with five different xc-functionals, compared to ADC(2)-x reference spectra. Spectra have been shifted to overlap with the pre-edge peak of the asymmetric clusters (vertical dashed line), and the IADs have been calculated over the illustrated energy region. IADs calculated for the summed spectra (large markers) and each spectrum individually (small markers). The applied energy shifts amount to, from top to bottom: 13.93, 19.74, −0.46, −1.49, 13.42, and −1.86 eV, respectively.
Focusing on the IADs of the summed spectra, the values range from 0.12 to 0.38, with similar values for the asymmetric and symmetric clusters. This can be compared to the x-ray emission spectra discussed in Fig. 2 of Ref. 85, where a similar analysis resulted in IADs of about 0.1–0.6, with potentially large discrepancies between the asymmetric and symmetric structures. Any such large shifts in IAD for the different structures are missing here, and we generally obtain good agreement between the IAD numerical values and a visual comparison between the spectrum features. We note that larger IADs are mainly associated with differences in the intensity and position of features ∼9–10 eV above the pre-edge, which would correspond to features ∼3 eV above the post-edge region when comparing to experiment (see Fig. 1). Considering the IADs of the individual spectra, these are generally larger than for the summed spectra, as a result of fluctuation canceling out in the summed spectra. The regions that primarily contribute to the IAD can be seen in Fig. 2, where the largest differences tend to be attributed to intense features at ∼543 eV (using the energy scale of this figure), in particular for the symmetric clusters, and features around 533–536 eV for the asymmetric structures. A more detailed study of the discrepancies of individual methods can be carried out by looking at the difference spectra, or splitting the total energy region into different intervals and calculating the IAD of each in turn.
IADs of different xc-functionals
A more systematic overview of the IADs of different methods is provided in Fig. 3, where the adopted pre-edge position and the IADs calculated with reference to ADC(2)-x and CCSD using three different broadening parameters (0.2, 0.3, and 0.4 eV HWHM) are reported. Results using ADC(n), CCSD, and TDDFT with a number of different xc-functionals are included, where the xc-functionals have been organized into different sets: long-range corrected functionals (CAM100,50,87,114,115 CAM-B3LYP,116 LRC-ωPBE,117 and ωB97112), the QTP family113,118,119 (LC-QTP, CAM-QTP00, CAM-QTP01, and CAM-QTP02), short-range corrected functionals tailored for XAS (SRC1-R1 and SRC2-R1),80 and global hybrids (PBE0,120 B97,121 and B3LYP with varying amounts of HF exchange).111 The different broadening parameters yield similar trends, with IADs increasing with decreasing HWHM,85 and slightly larger IADs for the symmetric structures. The dependence of the IAD on the broadening was also considered in Ref. 85, which reported relatively preserved trends, but differing absolute values, for HWHMs from 0.01 to 1.00 eV. Large HWHM values smoothen out any differences in spectra, while too small broadenings yield overly sharp features—a broadening of ∼0.2–0.3 eV provides a good compromise, as well as a qualitative agreement with experiment, and 0.3 eV has been used in this study. Furthermore, different broadening schemes can be applied, such as a Voigt profile, which can account for both Gaussian and Lorentzian broadening processes; here, Lorentzian broadening was applied.
Integrated absolute difference (IAD) of summed spectra of ten asymmetric and ten symmetric clusters, as obtained using the ADC(n) hierarchy and for a number of xc-functionals and compared to ADC(2)-x and CCSD reference spectra. The estimated asymmetric pre-edge peak position (peak pos., in eV) is used to align all spectra, and the IADs were calculated using broadenings of 0.2, 0.3, and 0.4 eV (HWHM), as indicated on each horizontal bar in the order 0.4, 0.3, and 0.2 from left to right.
Integrated absolute difference (IAD) of summed spectra of ten asymmetric and ten symmetric clusters, as obtained using the ADC(n) hierarchy and for a number of xc-functionals and compared to ADC(2)-x and CCSD reference spectra. The estimated asymmetric pre-edge peak position (peak pos., in eV) is used to align all spectra, and the IADs were calculated using broadenings of 0.2, 0.3, and 0.4 eV (HWHM), as indicated on each horizontal bar in the order 0.4, 0.3, and 0.2 from left to right.
For the wave function methods, we note that ADC(2)-x and CCSD are in best agreement in terms of integrated absolute difference, but ADC(2) and CCSD are closest in absolute energy. ADC(1) shows an overestimation in absolute energy due to the lack of relaxation, and ADC(3) also yields quite high transition energies. These results are in line with previous observations,87 where CCSD and ADC(2)-x are generally in good agreement, while ADC(2) (and CC2) shows larger discrepancies, and ADC(3) is more substantially off. The poor performance of ADC(3) is due to fortuitous error compensation for ADC(2)-x, which is broken for ADC(3).64 Additional calculations with CIS have also been performed, which give very similar results as ADC(1)—these methods have the same accuracy in terms of energies, but not in intensities.89
For TDDFT, previous studies have demonstrated that long-range corrected functionals are needed to obtain good agreement with ADC(2)-x and CCSD,87 featuring error spreads in energy and intensity of 0.2–0.3 eV and ∼10%, as compared to 0.3–0.6 eV and ∼20% for a pure hybrid. Such trends are seen in the present data, with the long-range corrected functionals featuring lower IADs [in particular when compared to ADC(2)-x]. Our recent study on x-ray emission spectra of water clusters showed that the amount of short-range HF exchange is most important for reproducing reference spectra and 1b1 energy-splits, and there is a dependence on the relative asymmetric/symmetric IADs with respect to the short-range exchange.85 However, such trends are not present here, where we see that the asymmetric and symmetric IADs of the same xc-functional are similar, and the amount of short-range HF exchange is not the main factor in deciding the agreement with reference spectra. While B3LYP with different amounts of HF exchange yields IADs ranging from small to large (with a minimum around 50%, which is also in the region where absolute energies are in best agreement with experiment94), the best performing xc-functionals are CAM100 and members of the QTP family (LC-QTP, QTP01, and QTP02). The parameters of these functionals are provided in Table I, where CAM-B3LYP and QTP00 are also included. These xc-functionals all use HF and B88 for exchange and LYP and VWM for correlation, and the well-performing functionals all have pure HF exchange in the asymptotic limit.
Parameters of long-range corrected xc-functionals, where the exchange is divided into short-range (SR) and long-range (LR) terms, and the slope of the range-correcting error function is given by the μ parameter (expressed in bohr−1).
. | . | Exchange . | Correlation . | ||||
---|---|---|---|---|---|---|---|
Functional . | μ . | HF (SR) . | HF (LR) . | B88 (SR) . | B88 (LR) . | LYP . | VWM . |
CAM-B3LYP | 0.330 | 0.19 | 0.65 | 0.81 | 0.35 | 0.81 | 0.19 |
CAM100 | 0.330 | 0.19 | 1.00 | 0.81 | 0.00 | 0.81 | 0.19 |
QTP00 | 0.290 | 0.54 | 0.91 | 0.46 | 0.09 | 0.80 | 0.20 |
QPT01 | 0.310 | 0.23 | 1.00 | 0.77 | 0.00 | 0.80 | 0.20 |
QTP02 | 0.335 | 0.28 | 1.00 | 0.72 | 0.00 | 1.00 | 0.00 |
LC-QTP | 0.475 | 0.00 | 1.00 | 1.00 | 0.00 | 1.00 | 0.00 |
. | . | Exchange . | Correlation . | ||||
---|---|---|---|---|---|---|---|
Functional . | μ . | HF (SR) . | HF (LR) . | B88 (SR) . | B88 (LR) . | LYP . | VWM . |
CAM-B3LYP | 0.330 | 0.19 | 0.65 | 0.81 | 0.35 | 0.81 | 0.19 |
CAM100 | 0.330 | 0.19 | 1.00 | 0.81 | 0.00 | 0.81 | 0.19 |
QTP00 | 0.290 | 0.54 | 0.91 | 0.46 | 0.09 | 0.80 | 0.20 |
QPT01 | 0.310 | 0.23 | 1.00 | 0.77 | 0.00 | 0.80 | 0.20 |
QTP02 | 0.335 | 0.28 | 1.00 | 0.72 | 0.00 | 1.00 | 0.00 |
LC-QTP | 0.475 | 0.00 | 1.00 | 1.00 | 0.00 | 1.00 | 0.00 |
Additional test calculations have been performed using a BHandHLYP functional with different fractions of HF exchange, as well as with 100% HF in the asymptotic limit, using μ = 0.33 bohr−1—this corresponds to CAM100 with different amounts of short-range HF (SRHF). For this, the IAD generally decreases for SRHF < 50% and increases for SRHF > 50%. Adding 100% asymptotic HF on the PBE0 xc-functional (25% SRHF) yields similar results, with the IAD decreasing to approximately half the original values. This could explain why QTP00 yields poor IADs [in particular when compared to ADC(2)-x], even though it has 0.91 long-range HF—it has 0.54 SRHF, which could lead to larger discrepancies. We further note that the SRC family has relatively large fractions of short-range HF exchange,80 which could potentially help explain why the SRC xc-functionals yield a bit higher IADs than, for example, the QTP family.
The precise value of the integrated absolute difference has previously been noted to depend on the broadening scheme and parameter, adopted energy alignment and energy region, and the number of structures.85 In addition to these factors, different IADs are obtained when considering asymmetric and symmetric structures separately, or when calculating it from the combined spectra—using a broadening of 0.3 eV, we here obtain IADs that are 17% ± 4% smaller when comparing the combined spectra (asymmetric+symmetric) to the average IAD of asymmetric and symmetric separately. Another factor that should be considered is the use of individual area-normalization of each spectrum in turn, as compared to using a common area-normalization from, for example, the symmetric spectra. Such a scheme could potentially introduce some artifacts if there is a systematic difference in relative intensities of the different spectra. For the structures considered here, this yields only a small effect—taking the fraction of the normalization area of the reference and test method for asymmetric minus that of symmetric only yields a small deviation of 2.3% ± 0.8%. As such, we see that considering each type of structure individually, but with a common energy shift, yields reliable results.
Comparing the IADs obtained here to the XES in Ref. 85, the IAD for XAS is seen to adopt a more narrow range of values—0.15–0.65, compared to 0.10–0.75. This is largely a result of the higher density of states in XAS (compare Fig. 2 here to Fig. 2 in Ref. 85), which yields some cancellation of features—more features distributed over an energy region increase the chances that misalignment in any two peaks partially cancels out.
From the results in Fig. 3, we suggest that CAM100 and QTP02 are the “best” performing xc-functionals for calculations of XAS on water, as compared to ADC(2)-x and CCSD reference spectra for these small water clusters. These xc-functionals are not necessarily expected to give similar reliability in relative features for other elements, or even necessarily for other solutions, but it can be noted that CAM100 has previously been identified as a well-performing functional for 1s → π* transitions,87 as well as providing reasonable spectra of liquid water and ice.50,51 Care should always be taken when using TDDFT, and several xc-functionals should be compared before proceeding with any production calculations.
Calculating XAS of liquid water
Moving to larger clusters, Fig. 4 reports the spectra obtained using CAM100 and QTP02, considering asymmetric and symmetric clusters with 6 or 32 molecules (using the same 20 original structures, but different numbers of included water molecules). Significant differences in spectrum features are noted, as the small clusters do not appropriately incorporate long-range interactions or delocalized excited states. In fact, a previous study has shown that over 2000 water molecules (in a QM/MM framework) are required to fully converge spectral features for individual structures, but that substantially smaller structures can be used within averaging schemes.50 It was furthermore seen that the integrated intensity was preserved for cluster sizes ranging from 32 to 96. Considering ice, the spectrum for clusters ranging from 27 to 92 water molecules has been considered using transition-potential DFT (TP-DFT), showing that some fluctuations exhibited by the smaller clusters are quenched when more molecules are added, leading to smoother spectra.49 Finally, TP-DFT spectra from 174 molecules centered in clusters of 39 molecules have been seen to be very similar (albeit with a difference in absolute energy) to calculations considering the same 174 molecules with periodic boundary conditions, with the full unit cell containing 192 molecules.49 As such, clusters with 32 molecules are deemed to be sufficient, but questions on employed basis sets and spectrum broadening schemes need to be considered.49
X-ray absorption spectra of 10 asymmetric and 10 symmetric water clusters, considering structures with 6 or 32 molecules, calculated using two different xc-functionals, CAM100 (top) and QTP02 (bottom).
X-ray absorption spectra of 10 asymmetric and 10 symmetric water clusters, considering structures with 6 or 32 molecules, calculated using two different xc-functionals, CAM100 (top) and QTP02 (bottom).
In a recent study by Folkestad and co-workers,43 the x-ray absorption spectrum of liquid water was calculated using multi-level65 CCSD and CC3, considering 896 structures from a path integral molecular dynamics (PIMD) simulation. The PIMD simulation used 24 beads and was propagated for 35 ps, using the canonical ensemble with a box with 32 molecules and a density of 1 g/cm3, and employed the strongly constrained and appropriately normed (SCAN) xc-functional.122 Structures were extracted from one bead, considering 28 independent snapshots and extracting clusters for each of the 32 molecules. The structures and resulting spectra are available on Zenodo,123 and we have adopted these structures to calculate the absorption spectrum using the two “best” performing xc-functionals identified in the present work. Two commonly adopted xc-functionals, B3LYP and CAM-B3LYP, have been included for comparison. The resulting peak-normalized spectra are shown in Fig. 5, as compared to multi-level CCSD and CC343 and to experiment.33 The spectra in Figs. 5 and 6 have been convoluted using the broadening scheme from Ref. 43, featuring a Voigt profile with a Lorentzian with HWHM of 0.1 eV and a Gaussian with standard deviation of 0.2 eV—this results in a HWHM of ∼0.6 eV. The TDDFT calculations have been performed for smaller clusters of 32 molecules, and using the basis set described in computational details, while the CCSD/CC3 results from Ref. 43 were obtained using larger clusters (∼100 molecules) and a different basis set (aug-cc-pVTZ for the central molecule, aug-cc-pVDZ for the four nearest neighbors, and cc-pVDZ for the remaining waters). A multi-level scheme was applied, where the central molecule was described using CCSD or CC3, the four nearest neighbors were described using CCSD, and the remaining environment adopted a frozen Hartree–Fock density. The differences in cluster size and basis set are unlikely to affect the final TDDFT spectra to any larger extent, and we expect the presented results to be comparable. The theoretical spectra have been shifted such that the pre-edge feature approximately overlaps with that of experiment, with energy shifts amounting to 11.40, 15.15, 1.50, and 0.00 eV for QTP02, CAM100, CCSD, and CC3, respectively. The CC3 calculation thus yields energies in very good agreement with experiment (although relativistic effects are missing). The B3LYP and CAM-B3LYP spectra lack a distinct pre-edge and have instead been shifted by the same energy as CAM100 (15.15 eV).
X-ray absorption spectrum of liquid water, calculated using four xc-functionals for 896 structures from a PIMD trajectory from Ref. 43, compared to multi-level CCSD and CC3 and experiment. The vertical dotted lines indicate main features in the experimental spectrum. From top to bottom: B3LYP, CAM-B3LYP, QTP02, CAM100, CCSD, MLCC3, and experiment.
X-ray absorption spectrum of liquid water, calculated using four xc-functionals for 896 structures from a PIMD trajectory from Ref. 43, compared to multi-level CCSD and CC3 and experiment. The vertical dotted lines indicate main features in the experimental spectrum. From top to bottom: B3LYP, CAM-B3LYP, QTP02, CAM100, CCSD, MLCC3, and experiment.
X-ray absorption spectrum of liquid water, calculated using the CAM100 xc-functional with three different basis sets (described in the main text). Legend provides the size of the different basis sets.
X-ray absorption spectrum of liquid water, calculated using the CAM100 xc-functional with three different basis sets (described in the main text). Legend provides the size of the different basis sets.
Considering first the results obtained using B3LYP and CAM-B3LYP (which has 65% asymptotic HF exchange), the agreement with experiment is poor, in particular for B3LYP. The addition of long-range correction in CAM-B3LYP improves the spectrum features, with distinct main- and post-edges starting to arise, but still feature overall poor agreement with experiment. As such, we conclude that the use of a global xc-functional, such as B3LYP, does not yield good agreement with XAS for liquid water or water clusters and that the addition of long-range correction requires a large fraction (∼100%) of HF exchange in the asymptotic limit to provide good results. Focusing on the remaining theoretical spectra, they all have relatively weak pre- and post-edges and a main edge that is too sharp and intense. In particular, CCSD yields low intensities for the post-edge region and a very sharp main-edge, and QTP02 also suffers from low intensities at higher energies, potentially due to issues in correctly capturing more delocalized excitations. CC3 and CAM100 yield relatively similar features, although the pre-edge of CC3 is a bit more distinct. The relative positions of the pre- and post-edges are relatively well captured in these four calculated spectra (omitting B3LYP and CAM-B3-LYP), but the first feature at higher energies is about 0.5–1.0 eV lower than the experimental post-edge. The exact features of the post-edge are affected by the adopted cluster size, basis set, and broadening scheme, but this is unlikely to substantially impact the total intensity over this region. Furthermore, the pre-edge is composed of local excitations, which should be relatively well captured by the current cluster sizes and basis sets. The origin of the remaining discrepancies to experiment could be due to inaccuracies in the molecular dynamics simulation, spectrum calculations, or both. A thorough investigation of this question, as well as a breakdown on which structural motifs that are associated with which features, is beyond the scope of the present study and will be approached in combination with reproducing other properties in order to achieve a holistic view of the structure of liquid water.24
In order to investigate the influence of the basis set on the spectrum features (in particular the post-edge), Fig. 6 includes CAM100 results using three different basis sets, obtained for the same structures as in Fig. 5. In B1, the same basis set as above is used (IGLO III for the main oxygen, IGLO II for the hydrogens, a triple-ζ ECP basis set for the other oxygens, and 9s9p diffuse functions added to the main oxygen), which results in 986 basis functions. B2 is inspired by Ref. 50, featuring the same basis set as above for the central oxygen, but with 19s19p diffuse functions (using an even-tempered construction starting with exponents 1.238 and 0.884). The hydrogens of the central molecule and the five closest water molecules are also given the same basis as above, but a double- ζ basis set is used for the outer molecules (2s2p for oxygen and 2s for hydrogen). This yields a smaller basis set, but the features are similar to B1. In order to avoid linear dependencies for B2 and B3, the threshold for removing basis functions is increased from 10−6 to 10−5. Finally, inspired by Ref. 49, B3 includes additional 9s9p diffuse basis functions centered at the oxygens of the four nearest water molecules (as well as 19s19p for the central oxygen). This yields a basis set of intermediate size and smoother features—the two post-edge peaks from B1 and B2 are now replaced by an extended feature with a center-of-max ∼6 eV above the pre-edge peak, in good agreement with experiment. The main-edge peak maximum is furthermore slightly less intense and redshifted by ∼0.15 eV. As such, a segmented basis set similar to B3, featuring diffuse functions centered not only on the main oxygen, but also nearest-neighbors, provides improved agreement with experiment. We note, however, that the main-edge is still significantly sharper than experiment.
CONCLUSIONS
The x-ray absorption spectrum of liquid water is still to be fully replicated by theory, and the origin of the discrepancies—inaccuracies in the underlying structures or in the spectrum calculations—remains unresolved. In the present study, we evaluate the efficacy of TDDFT in computing x-ray absorption spectra, using highly reliable ADC(2)-x and CCSD results as reference, considering both highly asymmetric and symmetric six-molecule clusters. Using the integrated absolute difference (IAD)23,82–85 as a metric, we show that xc-functionals with intermediate (∼0.2) amount of short-range Hartree–Fock (HF) exchange, but with pure HF exchange in the asymptotic long-range limit, yield best agreement with reference spectra. The CAM-QTP02 xc-functional is seen to yield best agreement with CCSD, and CAM100%-B3LYP (CAM-B3LYP with 100% asymptotic HF exchange) when compared to ADC(2)-x.
Considering the IAD, this was originally formulated for comparing experimental spectra of different compounds,23,82–84 but we have recently adopted this metric for evaluating the performance of different methods to calculate x-ray emission spectra of liquid water.85 The present work is a continuation of this effort, where we aim to use several different spectroscopies on the same water models to obtain a holistic view of the structure and dynamics of liquid water.24 The IAD is calculated from the difference between shifted and area-normalized spectra and thus measures agreement in both relative energies and intensities, providing a metric that agrees well with the visual comparison of spectra. It is suitable when dealing with ensembles of structures, and investigating the performance of several methods using different reference methods simultaneously, although care must be taken when aligning spectra, and the precise numerical value will depend on broadening scheme and parameter, the feature that spectra are aligned to, and the energy region for the integration. We believe that this measure is highly suitable for comparative studies of ensembles of structures and encourage further adaptation and evaluation. In the present work, the IAD is seen to adopt a more narrow range of values compared to XES85 (∼0.15–0.65, compared to ∼0.10–0.75), largely as a result of the higher density of states in XAS, which results in broader features.
In order to consider more realistic models of liquid water, clusters with 32 molecules were considered, which show noticeable differences compared to the spectra of the six-molecule models used for the reference calculations, although the main qualitative features are relatively unaffected. Furthermore, we compare the performance of two widely used xc-functionals (B3LYP and CAM-B3LYP) and the two “best” performing xc-functionals to recently published results using multi-level CCSD and CC3,43 as well as with experiment.33 B3LYP and CAM-B3LYP yield poor agreement with experiment, and the remaining theoretical spectra generally show a too intense main-edge peak, weak pre-edge features, and too low intensity in the post-edge (in particular for CCSD and CAM-QTP02). These remaining discrepancies could be due to inaccuracies in molecular dynamics, spectrum calculations, or both. Importantly, our CAM100 spectrum is in good agreement with high-level CC3 calculations, indicating that this is, indeed, a reasonable choice for the x-ray absorption spectrum of liquid water. Furthermore, it is seen that augmenting the basis set with additional diffuse functions centered at the nearest-neighbor oxygens yields improved features in the post-edge region, but the intensity remains too low as compared to the main-edge.
In conclusion, the present results suggest that the CAM100 xc-functional provides spectra in good agreement with higher-level reference methods [ADC(2)-x and CCSD] in terms of spectral shape (for six-molecule clusters), good agreement with MLCC3 calculations43 for larger clusters, and reasonable agreement with experiment.33 The integrated absolute difference (IAD) is seen to provide a suitable metric for evaluating the performance of different methods for calculating x-ray absorption spectra, as well as for x-ray emission spectra.85
ACKNOWLEDGMENTS
We acknowledge the support from the European Research Council (ERC) Advanced Grant under Project No. 101021166–GAS-WAT and Vetenskapsrådet under Grant No. 2020-05538. The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at the National Supercomputer Center (NSC), partially funded by Vetenskapsrådet under Grant No. 2022-06725.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Thomas Fransson: Conceptualization (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Lars G. M. Pettersson: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.