Real-time equation-of-motion CC cumulant and CC Green’s function simulations of photoemission spectra of water and water dimer

Coupled-cluster (CC) theory^1–7 has evolved into one of the most accurate formulations to capture complex correlation effects defining a broad class of quantum effects that drive chemical transformations, usually identified with bond-forming and bond-breaking processes.^7,8 Over the last few decades, theoretical efforts closely followed by the development of sophisticated computational models that take advantage of the ever-growing computational power of parallel architectures enabled the establishment of a hierarchy of approximations and novel formulations for closed/open-shell molecules, strong correlation, and large-scale systems. Special attention has been paid to single-reference CC formulations that can adequately describe chemical reactions and topological events of the corresponding ground-state potential energy surfaces, including barriers, avoided crossings, and multiple minima. This effort includes an impressive array of formulations based on the inclusion of high-rank collective effects contributing to the expansion of the ground-state wave function. This progress was made possible to achieve thanks to cornerstone implementations of the now ubiquitous models, such as coupled-cluster with singles and doubles (CCSD),⁵ singles, doubles and triples (CCSDT),^9–11 singles, doubles, triples and quadruples (CCSDTQ) and their perturbative counterparts,^12,13 including CCSD[T],^14,15 CCSD(T),¹⁶ CCSD(TQ), and CCSDT(Q) type approaches.^17,18

For its simplicity, the CCSD(T) method has assumed a unique position in high-accuracy chemical simulations of equilibrium properties and chemical reactions. However, a significant effort has been expended to alleviate serious issues associated with the perturbative nature of the triples in CCSD(T) for molecular configurations away from the equilibrium geometry.^19–29 Different design principles drive a large class of these developments in perturbative and non-perturbative formulations.

Structural changes in molecular systems, such as bond-breaking/forming, can be characterized by analyzing their potential energy surface (PES) or scrutinizing processes/states in various energy regimes sensitive to the geometry changes induced by chemical transformations. Excited state extensions of the CC formalism play an essential role in these studies either through equation-of-motion CC (EOMCC)^30–32 methodologies or CC linear response (LR-CC)^33,34 theory, which are usually formulated in frequency space. These methods have been applied widely in studies of excitation energies, excited-state PESs, excited-state non-adiabatic dynamics, frequency-dependent optical properties, and, more recently, multi-component polaritonic systems.³⁵ Green’s function extensions of the CC formalism (CCGF)^36–47 can also capture excited-state correlation effects needed to describe quasiparticles (QPs) and satellite peaks observed in x-ray photoemission spectra (XPS). It is worth mentioning that the development of the CC methods for the description of satellite peaks was paralleled by enabling vertex corrected GWΓ approaches⁴⁸ and adaptive sampling configuration interaction (ASCI) algorithms⁴⁹ capable of attaining high accuracy in describing multi-configurational electronic states for relatively little computational cost.

Time-dependent CC formulations, which provide an alternative way of describing excited-state processes, have attracted significant attention in the last decade.^50–56 In this context, we have recently developed a real-time equation-of-motion CC (RT-EOM-CC) method^57–60 to compute the one-particle GF based on a CCSD cumulant approach, which provides several advantages compared to other excited-state formulations. For example, the approach leads to an explicit exponential cumulant representation of the Green’s function⁶¹ and, at the same time, allows one to formulate a non-perturbative expression for the cumulant in terms of solutions to a set of coupled, first-order nonlinear differential equations for the CC amplitudes. Additionally, the RT-EOM-CC formulation provides a theoretical platform for including nonlinear corrections to the traditional cumulant approximations, which are usually linear in the self-energy.

In this paper, we compare the performance of the CCGF and RT-EOM-CC formulations when applied to core and valence binding energies of the water and water dimer molecules. In particular, we focus our analysis on identifying unique features of spectral functions, i.e., the location of quasiparticle and satellite peaks (for core and valence regions) as functions of geometry changes in bond-breaking processes. In this regard, we focus on the O–H stretch in H₂O and the bridging proton transfer process in (H₂O)₂. We demonstrate that the impact of the geometry changes is especially strong for the satellite region of the core spectral function. We also investigate the effect of the basis set on the accuracies of ionization potentials obtained with the CCGF and RT-EOM-CC models with singles and doubles (CCSDGF/RT-EOM-CCSD) and compare them with large CCSD(T) calculations.

This paper is organized as follows: In Sec. II, we provide a brief overview of the CCGF, IP-EOMCC, and RT-EOM-CC formulations and approximations used in the calculations for the H₂O and (H₂O)₂ systems. In Sec. III, we discuss the ionization potentials and spectral functions in the valence and core regions. Using RT-EOM-CCSD, we investigate how the bond-breaking process affects the spectral functions calculated for core and valence energy regimes. Special focus is on evaluating the performance of the RT-EOM-CCSD formalism in comparison with other CC formulations. In Sec. III, we also discuss the basis set effects on the accuracy of calculated ionization potentials. Finally, Sec. IV briefly summarizes our results.

In this section, we discuss basic aspects of various CC approaches used here to calculate correlated binding energies in various energy windows. In particular, we focus on the coupled-cluster Green’s function formalism and closely related ionization-potential equation-of-motion coupled-cluster (EOM-CC) approach and on the real-time EOM-CC cumulant formulation.

The retarded part of the Green’s function (advanced part can be developed in an analogous way) is defined by the following matrix elements $GpqR(ω)$⁠:

where E₀ is the corresponding ground-state energy for the N-electron system, η is a broadening factor, $ap/ap†$ are creation/annihilation operators for an electron in the p-th spin-orbital, and |Ψ_g⟩ is the ground state wave function of the N-electron system. The bi-variational CC^32,62 formalism utilizes distinct parameterizations for the bra (⟨Ψ_g|) and ket (|Ψ_g⟩) ground-state wave functions, i.e.,

which leads to the following form of the retarded part of the CC Green’s function (CCGF):^{36–40,42,44,63–65}

The similarity transformed operators $Ā$$(A=H,ap,aq†)$ are defined as $Ā=e−TAeT$⁠. $H̄N$ stands for a normal product form of $H̄$⁠. The numerical algorithms for calculating Eq. (4) employ auxiliary operators X_p(ω),

which satisfy the following equations:

In Eq. (5), the i, j, … (a, …) labels refer to spin-orbital labels occupied (unoccupied) in the N-electron reference function |Φ⟩. Using these operators, the matrix elements can be expressed in a simple form as

The tensors x_i(ω)_p, $xija(ω)p$⁠, etc., in Eq. (5) represent sought after cluster amplitudes. In our implementation of CCGF, we approximate the de-excitation Λ operator by the cluster operator T^†. This step is justified by analyzing low-order perturbative contributions to T and Λ operators. Additionally, this approximation does not affect the location of the CC Green’s function poles. The main numerical effort associated with constructing the retarded CC Green’s function is associated with the need to solve a large number of independent linear equations, which, in turn, can contribute to efficient parallel schemes utilizing multiple levels of parallelism. We implemented the CCGF model with single and double excitations (CCSDGF)⁶⁶ using the parallel tensor library TAMM (Tensor Algebra for Many-body Methods).⁶⁷ The current version of the CCGF code utilizes Cholesky-decomposed two-electron integrals that significantly reduce the memory requirements and inter-node communication.

The ionization potential EOMCC (IP-EOMCC) formulation⁶⁸ utilizes the EOMCC-type wave function expansion for the k-th electronic state |Ψ_k⟩ of the N − 1 electron system,

where |Φ⟩ and T are the reference function and cluster operator of the N-electron system, respectively, and the R_k operator is given by the expansion

R_k can be determined by solving a non-Hermitian eigenvalue problem in the N − 1 Hilbert space,

where $ωkIP$ stands for the k-th ionization potential. It should be noted that the eigenvalues of the IP-EOMCC [Eq. (10)] are identical to the singularities (as functions of the ω parameter for η = 0) of the solutions of Eq. (6). For the present studies, we use the TCE implementation of the IP-EOMCC model with singles and doubles (IP-EOMCCSD).⁶⁹ 

The cumulant approximation, in which the core-hole GF is approximated as the product of the free particle GF and the exponential of a cumulant, has been shown to provide accurate spectral functions for extended systems.^70–73 We have recently combined this form of the GF with the CC method to develop^57–60 a real-time equation-of-motion CC (RT-EOM-CC) approach to compute the core and valence one-electron GF based on a CC approximation to the cumulant. In this approximation, the retarded GF can be expressed as

where c is the index of the excited hole spin-orbital (which in our previous studies was restricted to be a core spin-orbital, but here can be any occupied spin-orbital), ϵ_c is its bare or Fock energy, $ENcorr$ is the correlation energy of the N electron ground state, and $CcR(t)$ is the retarded cumulant associated with c. As shown in Eq. (11), the CC approximation naturally results in a GF with an explicit exponential cumulant form,^61,75–76 which, unlike that obtained in self-energy formulations, is non-perturbative. The expression for the time derivative of the cumulant in terms of the time-dependent CC amplitudes takes the simple form as

where $fpq=ϵpδpq−vpcqc$ are the matrix elements of the N − 1 electron Fock operator, ϵ_p is the energy of spin-orbital p, and $vpqrs=pq‖rs$ are the anti-symmetrized two-particle Coulomb integrals over the generic spin-orbitals p, q, r, s (i.e., designating occupied or virtual spin-orbitals). The time-dependent amplitudes $tij…ab…(t)$ of the CC expansion are obtained from a set of coupled, first-order non-linear differential equations analogous to those used in the solution of the static CCSD equations (see, for instance, Ref. 60). These equations have the initial conditions $tij…ab…(0)=0$⁠, which result in the initial condition $CcR(0)=0$ for the cumulant in Eq. (11). As shown in Eq. (12), the RT-EOM-CC cumulant also naturally includes non-linear contributions, in contrast to traditional cumulant approximations that usually linearly depend on the self-energy. RT-EOM-CC uses two main approximations: (1) the separable approximation $0≃ac†N−1$ (here, $N−1$ is the fully correlated N − 1 electron portion of the N electron exact ground state wave function $0$⁠) and (2) a time-dependent (TD) CC ansatz for this N − 1 electron state, i.e., $N−1,t=Ñ(t)eT(t)ϕ$ (here, $Ñ(t)$ is the normalization factor and T(t) is the TD cluster operator). In contrast with traditional CC formulations, T(t) is defined in the N − 1 electron Fock space. Therefore, the reference determinant $ϕ$ has a hole in level c, i.e., $ϕ=acΦ$⁠, where $Φ$ is the traditional ground state N electron Hartree–Fock (HF) Slater determinant. We have shown⁶⁰ that, at the CCSD level, the RT-EOM-CC method gives accurate core quasiparticle (QP) binding energies with a mean absolute error (MAE) from the experiment of about 0.3 eV for the CH₄, NH₃, H₂O, HF, and Ne systems. This method also gives a good treatment of the many-body satellites with errors in the QP-satellite gap of less than 1 eV. Finally, despite the use of the separable approximation, RT-EOM-CCSD provides accurate valence ionization energies, as shown below and elsewhere.⁷⁷ 

The calculations of the water monomer use experimental geometries⁷⁸ where R(OH) = 0.958 Å and θ(HOH) = 104.48°. The molecule has C_2v symmetry and is oriented in the traditional way, with the atoms lying in the yz plane. This results in the molecular orbital configuration $(1a1)2(2a1)2(1b2)2(3a1)2(1b1)2$⁠. For the water dimer structure, we use the best estimate reported by Lane,⁷⁹ which was obtained using the counterpoise corrected CCSD(T)-F12b method at the complete basis set limit, with corrections for higher-order excitations, core correlation, relativistic effects, and diagonal Born–Oppenheimer effects. For reference, this structure has the following parameters: R(OH)_f = 0.956 85 Å, R(OH)_b = 0.964 14 Å, R(OH)_a = 0.958 43 Å, θ(HOH)_d = 104.854°, θ(HOH)_a = 104.945°, R(O⋯O) = 2.909 16 Å, α = 5.686°, and β = 123.458° (for an explanation of these labels, refer to Fig. 3 of Ref. 79). The proton transfer coordinate ΔR(OH)_b used here is defined as the otherwise rigid displacement of the bridging proton, around its optimal position, following the direction parallel to the axis defined by the O⋯O bond. For recent reviews of experimental and theoretical efforts to characterize the water dimer, see Refs. 80 and 81.

The RT-EOM-CCSD simulations were performed using the Tensor Contraction Engine (TCE)⁸² implementation of RT-EOM-CCSD.⁶⁰ All the simulations used a time step of 0.025 a.u. and frozen virtual orbitals above 10 a.u., when present. The total length of the simulations was 90 a.u. for the valence states of the monomer and dimer and 400 a.u. for the core state of the monomer. The valence simulations of the monomer used the aug-cc-pVTZ basis set,⁸³ while the core ones used the Sapporo-TZP⁸⁴ basis set. The RT-EOM-CC simulations of the dimer used the aug-cc-pVTZ basis set.

The CCGF spectral functions were calculated using the highly scalable parallel tensor library TAMM (Tensor Algebra for Many-body Methods)^67,85 implementation of the CCSDGF formalism.^64,65,66 The valence vertical ionization potentials for the water and water dimer were calculated using the IP-EOMCCSD code available in NWChem.^86,87

All CCSD(T) calculations for the water dimer in cc-pVXZ and aug-cc-pVXZ (X = D, T, Q, 5, 6) were obtained with the TAMM implementation of the CCSD(T) formalism,^67,85,88 with a linear dependence threshold for the basis sets of 10⁻⁵, a SCF convergence cutoff of 10⁻⁸ for the energy, and a CCSD convergence cutoff of 10⁻⁸. In the largest CCSD(T) calculations corresponding to the aug-cc-pV6Z, more than 1100 orbitals were correlated.

Table I shows the IP-EOMCCSD, relativistic IP-EOMCCSD, and RT-EOM-CCSD results for the ionization potentials corresponding to the 1²B₁, 1²B₂, and 1²A₁ states of H₂O⁺ (in the equilibrium structure of H₂O) calculated with various basis sets of the cc-pVXZ⁹¹ and aug-cc-pVXZ (X = D,T,Q,5)⁸³ type. The observed basis set effects can vary with the targeted states. For example, the difference between cc-pVDZ and cc-pV5Z corresponding to the lowest ionization potential (1²B₁) amounts to 0.917 eV. The analogous difference for the 1²B₂ state is equal to 0.609 eV. The utilization of larger aug-cc-pVXZ basis sets significantly reduces the dependence of the calculated IP-EOMCCSD ionization energies on the basis set employed, which is well illustrated by the 1²B₁ ionization potential, where the difference between IP-EOMCCSD results for aug-cc-pVDZ and aug-cc-pV5Z basis sets is reduced to 0.356 eV. It should also be noted that increasing the basis set size in the IP-EOMCCSD simulations of the 1²B₁ state results in overshooting its experimental value (12.621 ± 0.002). The aug-cc-pVTZ basis set provides the best IP-EOMCCSD estimate for this state. In analogy to the IP-EOMCCSD method, the RT-EOM-CCSD simulations using the same aug-cc-pVTZ basis set provide an accurate (12.587 eV) estimate of the experimental value. The RT-EOM-CCSD results also agree with the IP-EOMCCSD results and with the experiment for the aug-cc-pVTZ basis set for the remaining states. Table I also shows results obtained with the fully four-component relativistic IP-EOMCCSD.^94–94 The effects of relativity have been included by using a Dirac–Coulomb–Gaunt Hamiltonian, taken a Gaussian charge density distribution nuclear model to mimic the finite nucleus. The calculations use the dyall. ae4z basis for both H and O, where large and small component basis functions are uncontracted. All electrons are correlated, and none of the virtual spinors are frozen. The computed IPs are found to be very precise in comparison to the experiment. However, the effect of relativity has a negligible role.

Figure 1 (top) shows a comparison of the ionization potential of the 1²A″ state of (H₂O)₂ as a function of proton transfer coordinate calculated with the IP-EOMCCSD and RT-EOM-CCSD approaches and with the ΔCCSD, ΔCCSD(T), and ΔCCSDT methods as energy differences between the N and N − 1 electron systems. Figure 1 (bottom) shows the same ionization potentials relative to those obtained with the ΔCCSDT method. This behavior may be related to the fact that in the IP-EOMCCSD approach, as a consequence of the composite wave function ansatz given by Eq. (8), the energies of the N − 1 electron system are aware of the ground-state correlation effects for the N electron system, which may result in a better balance between correlation effects characterizing ground states of N and N − 1 electron systems. For all geometries considered, the ΔCCSD(T) results are in excellent agreement with the ΔCCSDT ones. At the same time, ΔCCSD(T) results provide systematic improvements of the ΔCCSD ionization potentials. It is also worth noting that the IP-EOMCCSD results are consistently better for the larger separations (using ΔCCSDT results as a reference point) than the ΔCCSD ones.

The effect of the basis set quality on the low-lying ionization energies calculated with the IP-EOMCCSD formalism for the equilibrium structure of the water dimer is shown in Table II. As in the case of the water monomer (see Table I), one can observe a significant impact of the basis set size on the calculated values of the ionization energies. In particular, the same trend of increasing the IP values with the growing size of the basis set can be noted. For example, for the aug-cc-pVXZ (X = D,T,Q), the aug-cc-pVQZ basis is needed to get a satisfactory agreement with the experiment for the IP corresponding to the 1²A″ state. In the case of the 1²A′ state, the IP-EOMCCSD results align with the experimental values. The RT-EOM-CCSD results in the aug-cc-pVDZ basis for 1²A″ are better than the IP-EOMCCSD ionization potential obtained with the same basis set.

In order to assess the basis set size effects on the ionization potentials, we have performed a series of calculations for the lowest ionization potential of the water dimer changes. These results are compiled in Table III. As described above, the ionization potential is calculated as the difference in energy between the CCSD(T) ground state energies of the water dimer and its cation at the optimal geometry. The basis set sizes vary from the smallest set (cc-pVDZ), with 50 basis functions for the water dimer, to the largest (aug-cc-PV6Z) with 1120 basis functions. Interestingly, the augmentation functions seem to significantly improve the convergence of the ionization potential to its basis set limit, as demonstrated by the smaller (0.25 eV) range of IP values for the augmented sets than for the unaugmented ones (1 eV). From Table III, it is evident that even a modest augmented basis set such as aug-cc-pVDZ does a good job compared to a much more saturated unaugmented basis set.

Figure 2 (bottom) shows the ionization energies and QP strengths for the HOMO 1²B₁ and HOMO-1 1²A₁ states of H₂O, computed with the RT-EOM-CCSD method and the aug-cc-pVTZ basis set. For these states, the RT-EOM-CCSD method gives good results, with errors compared to the experiment of only 0.02 and 0.03 eV, respectively. Both 1²B₁ and 1²A₁ states are sensitive to the bond stretching, with a decrease in the ionization energy of 1 eV along the bond dissociation. The ionization energies have very similar behavior over the whole range, yet, for these states, the QP strength shows very different behavior. While the QP strength of the 1²B₁ state decreases monotonically ∼6%, the one for the 1²A₁ varies only ∼2% but goes through a minimum at 1.1 Å. This can be clearly seen in the middle panel of Fig. 3, which shows that the spectral function of the 1²A₁ state is nearly featureless in the satellite region, in accordance with the little weight transfer from the QP. For the 1²B₁ state (Fig. 3, bottom), the ∼6% weight transfer is manifested in a small increase in the satellite weight in the −50 to −30 eV region.

Figure 4 (bottom) shows the results for the ionization of the HOMO 1²A″ and HOMO-1 1²A′ states of (H₂O)₂ as a function of proton transfer coordinate. For reference, the top panel of Fig. 4 shows the potential energy surfaces for the ground and first ionized (1²A″) states computed at the standard CCSD level. Unlike the monomer case, which shows little variation as a function of coordinate, the two lowest ionized states of the dimer vary significantly. The predicted decrease in IP associated with the proton transfer is expected, given the formation of an $(OH)−(H3O)+$ cluster. The IP of this proton-transferred system should be dominated by the ionization of the (OH)⁻ fragment, which has a low IP despite being“solvated” by the (H₃O)⁺ fragment.⁹⁸ Moreover, the first two ionized states of (OH)⁻(H₃O)⁺ are expected to become nearly degenerate, given that the HOMO state of isolated OH⁻ is doubly degenerate (1π). The CCSD PES also predicts an adiabatic ionization energy from (H₂O)₂ to the disproportionated ion (OH)(H₃O)⁺ of 10.83 eV, in reasonable agreement with other estimated values of 10.64 eV.⁹⁹ The QP strength (Fig. 4) shows a small (∼4%) variation over the proton transfer range for the 1²A″ state and goes through a minimum in the region 0.4–0.6 Å. The 1²A′ state shows very similar behavior. The QP strength minima coincides with the transition state region in the ionized state potential energy surface and results from an increase in many-body electron correlation in the bond breaking/reforming region, which manifests itself in the spectral function as a transfer of weight to the satellites (Fig. 3). Besides the QP shift as a function of proton transfer coordinate, which actually results in an overall shift, the spectral function shows very little change, as shown in Fig. 5.

We are interested in the change in satellite features as a function of proton transfer coordinate, ΔR(OH)_b. From Fig. 5, a shoulder satellite peak can be observed on the lhs of the main QP peak at all the studied ΔR(OH)_b’s, and the satellite peak position red-shifts as ΔR(OH)_b gets larger. This observation is also confirmed by independent static CCSDGF calculations using the same basis set for the same (H₂O)₂ systems at the same ΔR(OH)_b’s and over the same energy regime, and the corresponding results are shown in Fig. 6. The only difference between Figs. 5 and 6 is that, in the CCSDGF results, the satellites are more separated from the QP and more refined satellite structure can be observed. At least two shoulder satellites are observed on the lhs of the QP peak, while in the RT-EOM-CCSD results, the refined satellite structure is embedded in the background of the QP peak. These peaks would likely be better resolved if with longer RT-EOM-CCSD calculations. Further wave function analysis of the satellites features shows a π → σ* transition (e.g., from HOMO-1 to LUMO) along with the main ionization process from the HOMO. However, it is worth noting that the satellites predicted by CCSDGF are usually only qualitative, and systematic improvement toward quantitative evaluation would require including higher-order excitation in the CCGF framework as already demonstrated in the previous CCGF studies (see, for example, Ref. 100). On the other hand, previous studies for small- and medium-sized molecules (e.g., H₂O, NH₃, and CH₄)^57–60 suggested that RT-EOM-CCSD is capable of providing better predictions for the satellite positions attributed to its explicit consideration of the time evolution of the correlation in the N − 1 electron state. Given the fact that only singles and doubles are included in both CCSDGF and RT-EOM-CCSD calculations and the satellite peaks are close to each other, we conjecture that the satellite peak on the lhs of the main QP peak also features a π → σ* transition along with the main ionization. It is worth mentioning that the satellite assignment in the RT-EOM-CCSD method is challenging, and new methods in this regard are now under intensive development by the same authors.

Figure 2 (top) shows the ionization energy and QP strength for the core state of the water monomer, computed with the RT-EOM-CCSD method and the Sapporo-TZP basis set. Unlike with the valence states, the core ionization energy is nearly insensitive to the deformation of the bond with a change of only ∼0.7 eV. At the experimental bond length, the agreement with the experimental core ionization energy error of about 0.5 eV is larger than for the valence states yet reasonable. On the other hand, the core ionization QP strength is much more sensitive to the bond length than for the valence states, decreasing by about 10% over the same range. This implies a significant transfer of weight from the QP to the satellite region. This transfer can be seen in Fig. 3 (top), which focuses on the satellite region of the spectral function as a function of bond length. Unlike the QP, which barely changes position, the satellite region 60 eV below it shows remarkable changes in weight distribution. As the bond dissociates, the main feature, 26 eV below the QP, becomes less pronounced and is replaced by a pair of features about 9 and 12 eV below the QP.

In the present studies, we have compared the performance of several high-accuracy coupled-cluster formulations in calculating the ionization potentials and spectral functions of the water monomer and dimer for various geometries corresponding to bond-breaking and proton transfer processes. We established that the ionization potentials calculated with the CCSDGF, IP-EOMCCSD, and RT-EOM-CCSD methods reproduce those calculated with the CCSD(T) and CCSDT approaches for those systems. When compared to the experiment, these methods provide accurate results with average errors on the order of 0.1 eV. The direct comparison and crosscheck between the CCSDGF and RT-EOM-CCSD results also provide insights into the possible many-body excitations responsible for the satellite region above the quasiparticle peak. The peak assignment can be further improved by including higher-order excitation in the CCGF and RT-EOM-CC frameworks, and new methods are currently being explored. We find that the satellite region of the water monomer is particularly sensitive to bond breaking. We have also demonstrated that the CCSD(T) formalism can almost perfectly reproduce the CCSDT ionization potentials for cases where CCSDT calculations were possible for neutral and charged systems. This prompted us to study the effect of basis set quality on the accuracy of the CCSD(T) ionization potentials. We find that for the water dimer, the aug-cc-pVTZ and aug-cc-pV5Z basis sets provide results converged to 0.1 and 0.01 eV, respectively. The largest CCSD(T)/aug-cc-pV6Z calculations for the (H₂O)₂ and $(H2O)2+$ systems involved 1120 orbitals. While both CCSD(T) IPs obtained with the cc-pV6Z and aug-cc-pV6Z basis are almost identical, the discrepancy between calculated IPs obtained with the smallest and the largest basis sets of a given type is significantly larger for the cc-pVXZ series and amounts to 0.89 eV, showing that the augmented basis sets provide faster convergence. Finally, we also reported results with the fully four-component relativistic IP-EOMCCSD method for the lowest three states of the water monomer and found that its predictions are precise in comparison to the available experimental values. However, we found the effects of relativity to be almost negligible for the valence ionization potential values.

This work was supported by the Computational Chemical Sciences Program of the U.S. Department of Energy, Office of Science, BES, Chemical Sciences, Geosciences and Biosciences Division in the Center for Scalable and Predictive methods for Excitations and Correlated phenomena (SPEC) at PNNL, with computational support from NERSC, a DOE Office of Science User Facility, under Contract No. DE-AC02-05CH11231. B.P. also acknowledges the support from the Laboratory Directed Research and Development (LDRD) Program at PNNL.

The authors have no conflicts of interest to declare.

The data that support the findings of this study are available from the corresponding author upon reasonable request.