Double Ionization Potential Equation-of-Motion Coupled-Cluster Approach with Full Inclusion of 4-Hole–2-Particle Excitations and Three-Body Clusters

The single-reference coupled-cluster (CC) theoryCoester (1958); Coester and Kümmel (1960); Čížek (1966, 1969); Paldus, Čížek, and Shavitt (1972) and its equation-of-motion (EOM) extensions to electronically excitedEmrich (1981a, b); Geertsen, Rittby, and Bartlett (1989); Comeau and Bartlett (1993); Stanton and Bartlett (1993) and electron attached and ionized states Nooijen and Bartlett (1995a, b); Hirata, Nooijen, and Bartlett (2000); Musiał and Bartlett (2003); Gour, Piecuch, and Włoch (2005, 2006); Gour and Piecuch (2006); Nooijen and Snijders (1992, 1993); Stanton and Gauss (1994); Bartlett and Stanton (1994); Musiał, Kucharski, and Bartlett (2003); Musiał and Bartlett (2004); Bomble et al. (2005); Kamiya and Hirata (2006); Nooijen and Bartlett (1997a, b); Wladyslawski and Nooijen (2002); Nooijen (2002); Sattelmeyer, Schaefer, and Stanton (2003); Musiał, Perera, and Bartlett (2011); Musiał, Kucharski, and Bartlett (2011); Kuś and Krylov (2011, 2012); Shen and Piecuch (2013, 2014); Ajala, Shen, and Piecuch (2017); Shen and Piecuch (2021); Gulania et al. (2021); Musiał et al. (2012) have become preeminent methods of quantum chemistry. In this Communication, we focus on the double ionization potential (DIP) EOMCC framework, which allows one to directly determine the ground and excited states of doubly ionized molecular species and which is useful in many applications, such as Auger electron spectroscopy, Ghosh, Vaval, and Pal (2017); Skomorowski and Krylov (2021a, b); Jayadev et al. (2023) singlet–triplet gaps in biradicals, Kuś and Krylov (2011, 2012); Shen and Piecuch (2013, 2014) and strong-field-induced chemical reactivity.Stamm et al. (2025)

In the DIP-EOMCC formalism, the ground ($\mu=0$) and excited ($\mu>0$) states of the target ($N-2$)-electron system are expressed as

where the doubly ionizing operator

with $R_{\mu,(n+2)h\mbox{-}np}$ representing its $(n+2)$-hole–$n$-particle [$(n+2)h$–$np$] components, removes two electrons from the CC ground state

of the underlying $N$-electron species, in which

is the cluster operator and $|\Phi\rangle$ is the reference determinant that serves as a Fermi vacuum. The many-body components of the $T$ and $R_{\mu}^{(-2)}$ operators are given by

and

respectively, where, as usual, indices $i,j,\,\ldots$ ($a,b,\,\ldots$) denote the spinorbitals that are occupied (unoccupied) in $|\Phi\rangle$ and $a^{p}$ ($a_{p}$) represents the fermionic creation (annihilation) operator associated with the spinorbital $|p\rangle$. $M_{R}$ and $M_{T}$ control the DIP-EOMCC theory level, with $M_{R}=N-2$ defining exact calculations and $M_{R}<N-2$ and $M_{T}\leq N$ leading to DIP-EOMCC approximations.

The DIP-EOMCC approaches that have been implemented so far include the basic DIP-EOMCCSD(3$h$-1$p$) method, Wladyslawski and Nooijen (2002); Nooijen (2002); Musiał, Perera, and Bartlett (2011); Kuś and Krylov (2011, 2012); Shen and Piecuch (2013, 2014) defined by $M_{R}=1$ and $M_{T}=2$, which describes 2$h$ and 3$h$-1$p$ correlations on top of the CC calculations with singles and doubles (CCSD),Purvis and Bartlett (1982); Cullen and Zerner (1982); Scuseria et al. (1987); Piecuch and Paldus (1989) its higher-level DIP-EOMCCSD(4$h$-2$p$) extensionShen and Piecuch (2013, 2014) corresponding to $M_{R}=2$ and $M_{T}=2$, which also describes 4$h$-2$p$ correlations on top of CCSD, and the DI-EOMCCSDT schemeMusiał, Perera, and Bartlett (2011) corresponding to $M_{R}=1$ and $M_{T}=3$, which accounts for 2$h$ and 3$h$-1$p$ correlations on top of the CC calculations with singles, doubles, and triples (CCSDT).Noga and Bartlett (1987); Scuseria and Schaefer (1988) While all of these methods and the various approximations to them aimed at reducing computational costs are useful tools for determining DIP energies in molecular systems and singlet–triplet gaps in certain biradicals, our recent numerical tests, including the DIPs of a few small molecules examined in this Communication, indicate that the explicit incorporation of 4$h$-2$p$ correlations, needed to achieve a highly accurate description,Shen and Piecuch (2013, 2014) may not be well balanced with the CCSD treatment of the underlying $N$-electron species, resulting in some cases in loss of accuracy compared to the basic DIP-EOMCCSD(3$h$-1$p$) level.

The main goal of this work is to enrich the existing arsenal of DIP-EOMCC methods and, in particular, address the above concerns by examining the high-level DIP-EOMCC approach with a full treatment of both 4$h$-2$p$ and $T_{3}$ correlations, corresponding to $M_{R}=2$ and $M_{T}=3$ and abbreviated as DIP-EOMCCSDT(4$h$-2$p$). The present study describes the efficient formulation and computer implementation of full DIP-EOMCCSDT(4$h$-2$p$) and its less expensive DIP-EOMCCSD(T)(a)(4$h$-2$p$) approximation, including the factorized and programmable expressions that define them. We illustrate the performance of both methods, as coded for nonrelativistic and spin-free scalar-relativistic Hamiltonians in the open-source CCpy packageCCp available on GitHub, by calculating the vertical DIPs of H₂O, CH₄, BN, Cl₂, Br₂, and HBr and comparing the results with those obtained using the DIP-EOMCCSD(3$h$-1$p$) and DIP-EOMCCSD(4$h$-2$p$) approaches and the reliable, high-accuracy, reference data.

The key step of any DIP-EOMCC calculation is a diagonalization of the similarity-transformed Hamiltonian $\overline{H}_{N}=e^{-T}H_{N}e^{T}=(H_{N}e^{T})_{C}$ associated with the $N$-electron CC ground state in the relevant ($N-2$)-electron subspace of the Fock space corresponding to the content of the ionizing operator $R_{\mu}^{(-2)}$. Here, subscript $C$ designates the connected operator product and $H_{N}=H-\langle\Phi|H|\Phi\rangle=F_{N}+V_{N}$ is the electronic Hamiltonian in the normal-ordered form relative to $|\Phi\rangle$, with $F_{N}$ and $V_{N}$ representing its Fock and two-electron interaction components. In the case of the cluster and ionizing operators defined by Eqs. (4) and (2), respectively, the basis states that span the ($N-2$)-electron subspace of the Fock space used in the DIP-EOMCC calculations are $|\Phi_{ij}\rangle=a_{j}a_{i}|\Phi\rangle$ and $|\Phi_{ijk_{1}\ldots k_{n}}^{\phantom{ab}c_{1}\ldots c_{n}}\rangle=a^{c_{1}}\ldots a^{c_{n}}a_{k_{n}}\ldots a_{k_{1}}a_{j}a_{i}|\Phi\rangle$, $n=1,\ldots,M_{R}$. Assuming that $M_{R}\leq M_{T}$ (a condition required for retaining size intensivity of the resulting double ionization energiesShen and Piecuch (2013, 2014, 2021)), the DIP-EOMCC eigenvalue problem is given by

where $\overline{H}_{N,\text{open}}$ refers to the diagrams of $\overline{H}_{N}$ containing external fermion lines and $\omega_{\mu}^{(N-2)}=E_{\mu}^{(N-2)}-E_{0}^{(N)}$ is the vertical DIP energy representing the difference between the total energy of the ground ($\mu=0$) or excited ($\mu>0$) state of the ($N-2$)-electron target system, denoted as $E_{\mu}^{(N-2)}$, and the ground-state CC energy of the underlying $N$-electron species, $E_{0}^{(N)}=\langle\Phi|H|\Phi\rangle+\langle\Phi|\overline{H}_{N}|\Phi\rangle$. In practice, including this work, the solutions of Eq. (7) are obtained using the Hirao–Nakatsuji generalizationHirao and Nakatsuji (1982) of the Davidson diagonalization algorithmDavidson (1975) to non-Hermitian Hamiltonians.

In the DIP-EOMCCSDT(4$h$-2$p$) approach developed in this work, the CCSDT similarity-transformed Hamiltonian, designated as $\overline{H}_{N}^{(\text{CCSDT})}$, is diagonalized in the $(N-2)$-electron subspace of the Fock space spanned by $|\Phi_{ij}\rangle$, $|\Phi_{ijk}^{\phantom{ab}c}\rangle$, and $|\Phi_{ijkl}^{\phantom{ab}cd}\rangle$. The programmable factorized expressions for the required projections of the left-hand-side of Eq. (7), with $T$ truncated at $T_{3}$ and $R_{\mu}^{(-2)}$ truncated at $R_{\mu,4h\mbox{-}2p}$, as implemented in CCpy, are

and

where we use the Einstein summation convention over repeated upper and lower indices and $\mathscr{A}^{pq}=\mathscr{A}_{pq}=1-(pq)$, $\mathscr{A}^{pqr}=\mathscr{A}^{p/qr}\mathscr{A}^{qr}$, and $\mathscr{A}^{pqrs}=\mathscr{A}^{p/qrs}\mathscr{A}^{qrs}$, with $\mathscr{A}^{p/qr}=1-(pq)-(pr)$ and $\mathscr{A}^{p/qrs}=1-(pq)-(pr)-(ps)$, are index antisymmetrizers. The expressions for the one-body ($\bar{h}_{p}^{q}$) and two-body ($\bar{h}_{pq}^{rs}$) components of the similarity-transformed Hamiltonian as well as the additional intermediates entering Eqs. (9) and (10) are provided in Tables 1 and 2. Equations (8)–(10) imply that the diagonalization step of DIP-EOMCCSDT(4$h$-2$p$) has computational costs identical to those characterizing DIP-EOMCCSD(4$h$-2$p$), which scale as $n_{o}^{4}n_{u}^{4}$ or $\mathscr{N}^{8}$, where $n_{o}$ ($n_{u}$) is the number of occupied (unoccupied) orbitals in $|\Phi\rangle$ and $\mathscr{N}$ is a measure of the system size. However, the overall computational effort associated with the DIP-EOMCCSDT(4$h$-2$p$) approach is considerably higher than that of DIP-EOMCCSD(4$h$-2$p$) since in the DIP-EOMCCSDT(4$h$-2$p$) case, one also has to solve the CCSDT equations for the underlying $N$-electron species, which involve $n_{o}^{3}n_{u}^{5}$ steps, as opposed to the much less expensive $n_{o}^{2}n_{u}^{4}$ ($\mathscr{N}^{6}$) steps of CCSD.

Given the high computational costs of the DIP-EOMCCSDT(4$h$-2$p$) method, we also consider the more practical DIP-EOMCCSD(T)(a)(4$h$-2$p$) scheme, in which we adopt the Møller-Plesset (MP) partitioning of the Hamiltonian and, following Ref. Matthews and Stanton, 2016, incorporate the leading $T_{3}$ correlation effects by correcting the $T_{1}$ and $T_{2}$ clusters obtained with CCSD using the formulas

and

where $D_{1}$ and $D_{2}$ are the usual MP denominators for singles and doubles and

is the lowest-order approximation to $T_{3}$, with $D_{3}$ representing the MP denominator for triples. Once $\tilde{T}_{1}$, $\tilde{T}_{2}$, and $\tilde{T}_{3}$ are determined via Eqs. (11)–(13), the resulting CCSD(T)(a) similarity-transformed Hamiltonian, constructed using the recipe described in Ref. Matthews and Stanton, 2016, is diagonalized in the same way as $\overline{H}^{(\text{CCSDT})}_{N}$ in DIP-EOMCCSDT(4$h$-2$p$) with the help of Eqs. (8)–(10). By eliminating the need for performing CCSDT calculations, the most expensive steps of DIP-EOMCCSD(T)(a)(4$h$-2$p$) scale as $n_{o}^{4}n_{u}^{4}$ rather than $n_{o}^{3}n_{u}^{5}$.

To illustrate the performance of the DIP-EOMCCSDT(4$h$-2$p$) and DIP-EOMCCSD(T)(a)(4$h$-2$p$) methods, as implemented for nonrelativistic and spin-free scalar-relativistic Hamiltonians in the open-source CCpy package,CCp we applied them, along with their DIP-EOMCCSD(3$h$-1$p$) and DIP-EOMCCSD(4$h$-2$p$) counterparts, available in CCpy as well, to the vertical DIPs of two sets of small molecules, for which highly accurate reference data can be found in the literature. The first set, which is the primary focus on our discussion below, consisted of the H₂O, CH₄, and BN molecules, as described by the aug-cc-pVTZ basis set,Dunning (1989); Kendall, Dunning, and Harrison (1992) where we compared our DIP-EOMCC DIPs corresponding to the lowest singlet and triplet states of the $\mathrm{(H_{2}O)^{2+}}$, $\mathrm{(CH_{4})^{2+}}$, and $\mathrm{(BN)^{2+}}$ dications with their near-exact counterparts obtained in Ref. Marie et al., 2024 with the configuration interaction (CI) approach using perturbative selection made iteratively, abbreviated as CIPSI,Huron, Malrieu, and Rancurel (1973); Garniron et al. (2017, 2019) extrapolated to the full CI limit. We also calculated the vertical DIPs of the Cl₂, Br₂, and HBr molecules corresponding to the triplet ground states and low-lying singlet states of $\mathrm{(Cl_{2})}^{2+}$, $\mathrm{(Br_{2})}^{2+}$, and $\mathrm{(HBr)}^{2+}$, comparing our DIP-EOMCC results with the reliable experimental data reported in Refs. McConkey et al., 1994; Fleig et al., 2008; Eland, 2003. In this case, to obtain insights into the convergence of our calculated DIP values with the basis set, we used the cc-pVTZ and cc-pVQZ bases.Dunning (1989); Woon and Dunning (1993); Wilson et al. (1999) Following Ref. Marie et al., 2024, the equilibrium geometries of the ground-state H₂O, CH₄, and BN molecules used in our DIP-EOMCC computations were taken from Ref. Marie and Loos, 2024, where they were optimized with the CC3/aug-cc-pVTZ approach. The equilibrium bond lengths of Cl₂, Br₂, and HBr were taken from Ref. Huber and Herzberg, 1979. In setting up and solving the DIP-EOMCC eigenvalue problems for the $\mathrm{(H_{2}O)^{2+}}$, $\mathrm{(CH_{4})^{2+}}$, $\mathrm{(BN)^{2+}}$, $\mathrm{(Cl_{2})}^{2+}$, $\mathrm{(Br_{2})}^{2+}$, and $\mathrm{(HBr)}^{2+}$ target species and executing the preceding CC computations for their neutral parents, we used the restricted Hartree–Fock (RHF) orbitals of the respective H₂O, CH₄, BN, Cl₂, Br₂, and HBr molecules. The relevant RHF reference determinants and transformed one- and two-electron integrals were generated with the PySCF code,Sun et al. (2018, 2020) with which CCpy is interfaced. The scalar-relativistic effects included in our DIP-EOMCC calculations for the DIPs of Cl₂, Br₂, and HBr were handled using the SFX2C-1e spin-free exact two-component approach of Ref. Cheng and Gauss, 2011, as implemented in PySCF, and the lowest-energy orbitals correlating with the chemical cores of non-hydrogen atoms were frozen in post-RHF steps.

The results of our DIP-EOMCCSD(3$h$-1$p$), DIP-EOMCCSD(4$h$-2$p$), DIP-EOMCCSD(T)(a)(4$h$-2$p$), and DIP-EOMCCSDT(4$h$-2$p$) computations for the vertical DIPs of H₂O, CH₄, and BN corresponding to the lowest singlet and triplet states of the $\mathrm{(H_{2}O)^{2+}}$, $\mathrm{(CH_{4})^{2+}}$, and $\mathrm{(BN)^{2+}}$ dications are summarized in Table 3. The vertical DIPs of Cl₂, Br₂, and HBr corresponding to the ground and low-lying excited states of $\mathrm{(Cl_{2})}^{2+}$, $\mathrm{(Br_{2})}^{2+}$, and $\mathrm{(HBr)}^{2+}$ resulting from our DIP-EOMCCSD(3$h$-1$p$), DIP-EOMCCSD(4$h$-2$p$), DIP-EOMCCSD(T)(a)(4$h$-2$p$), and DIP-EOMCCSDT(4$h$-2$p$) calculations are reported in Table 4. In much of the remaining discussion, we focus on the DIPs of H₂O, CH₄, and BN shown in Table 3, where we compare our DIP-EOMCC data with their near-full-CI counterparts determined with the help of CIPSI in Ref. Marie et al., 2024. While brief remarks about our DIP-EOMCC calculations for the DIPs of Cl₂, Br₂, and HBr summarized in Table 4 are given here as well, a more thorough analysis of these computations is provided in the Supplementary Material. As shown in Table 3, the vertical DIPs obtained in the DIP-EOMCCSD(3$h$-1$p$) calculations using the aug-cc-pVTZ basis set are characterized by significant errors relative to their near-exact counterparts resulting from the CIPSI extrapolations performed in Ref. Marie et al., 2024, which are 0.77 and 0.61 eV for the $X\>^{3}\mathrm{B}_{1}$ and $a\>^{1}\mathrm{A}_{1}$ states of $\mathrm{(H_{2}O)^{2+}}$, 0.32 and 0.31 eV for the $X\>^{3}\mathrm{T}_{1}$ and $a\>^{1}\mathrm{E}$ states of $\mathrm{(CH_{4})^{2+}}$, and 0.44 and 0.35 eV for the $X\>^{3}\Sigma^{-}$ and $a\>^{1}\Delta$ states of $\mathrm{(BN)^{2+}}$. For all of the calculated states of $\mathrm{(H_{2}O)^{2+}}$, $\mathrm{(CH_{4})^{2+}}$, and $\mathrm{(BN)^{2+}}$ considered in Table 3 and for the majority of states of $\mathrm{(Cl_{2})}^{2+}$, $\mathrm{(Br_{2})}^{2+}$, and $\mathrm{(HBr)}^{2+}$ examined in Table 4, especially when a larger cc-pVQZ basis set is employed, the DIP-EOMCSCD(3$h$-1$p$) computations overestimate the corresponding DIPs of H₂O, CH₄, BN, Cl₂, Br₂, and HBr relative to the reference data. In both sets of molecular examples shown in Tables 3 and 4, the poor quality of the DIPs obtained with the DIP-EOMCCSD(3$h$-1$p$) approach seems to be a consequence of neglecting the 4$h$-2$p$ component of $R_{\mu}^{(-2)}$.

Indeed, when $R_{\mu,4h\mbox{-}2p}$ is included in $R_{\mu}^{(-2)}$ via the DIP-EOMCCSD(4$h$-2$p$) approach, the errors in the lowest triplet and singlet DIPs of the water molecule relative to the extrapolated CIPSI data obtained with DIP-EOMCCSD(3$h$-1$p$), of 0.77 and 0.61 eV, reduce in the DIP-EOMCCSD(4$h$-2$p$) calculations to 0.25 and 0.25 eV, respectively. The analogous errors characterizing the DIP-EOMCCSD(3$h$-1$p$) computations for methane, which are 0.32 and 0.31 eV, reduce to 0.17 and 0.18 eV, respectively, when the DIP-EOMCCSD(4$h$-2$p$) method is employed. Unfortunately, the DIP-EOMCCSD(4$h$-2$p$) approach does not always improve the DIP-EOMCCSD(3$h$-1$p$) results. For example, the DIP values obtained in the DIP-EOMCCSD(4$h$-2$p$) calculations for BN are less accurate than their DIP-EOMCCSD(3$h$-1$p$) counterparts, increasing the 0.44 and 0.35 eV errors in the DIP-EOMCCSD(3$h$-1$p$) data for the $X\>^{3}\Sigma^{-}$ and $a\>^{1}\Delta$ states of $\mathrm{(BN)^{2+}}$ relative to their extrapolated CIPSI values to 0.61 and 0.67 eV, respectively, when DIP-EOMCCSD(3$h$-1$p$) is replaced by DIP-EOMCCSD(4$h$-2$p$). In general, the DIPs of H₂O, CH₄, and BN resulting from the DIP-EOMCCSD(4$h$-2$p$) calculations reported in Table 3 lie substantially below the corresponding near-full-CI extrapolated CIPSI data. Similar remarks apply to the DIPs of Cl₂, Br₂, and HBr shown in Table 4. This behavior of DIP-EOMCCSD(4$h$-2$p$) can be attributed to the imbalance between the high-level 4$h$-2$p$ treatment of double ionization and the lower-level CCSD description of the neutral species.

The results of our DIP-EOMCCSDT(4$h$-2$p$) calculations, in which the similarity-transformed Hamiltonian of CCSD is replaced by its CCSDT counterpart, allowing us to treat the $N$- and ($N-2$)-electron species in a more accurate and balanced manner, confirm this observation. As shown in Table 3, the DIPs of H₂O, CH₄, and BN become much more accurate when we move from the DIP-EOMCCSD(4$h$-2$p$) approach to its higher-level DIP-EOMCCSDT(4$h$-2$p$) counterpart. The errors in the DIPs associated with the lowest triplet and singlet dicationic states of H₂O, CH₄, and BN relative to the extrapolated CIPSI reference data, which are 0.25 and 0.25 eV for $\mathrm{(H_{2}O)^{2+}}$, 0.17 and 0.18 eV for $\mathrm{(CH_{4})^{2+}}$, and 0.61 and 0.67 eV for $\mathrm{(BN)^{2+}}$ in the DIP-EOMCCSD(4$h$-2$p$) case, reduce to the minuscule 0.02 and 0.03 eV for $\mathrm{(H_{2}O)^{2+}}$, 0.00 and 0.01 eV for $\mathrm{(CH_{4})^{2+}}$, and 0.01 and 0.00 eV for $\mathrm{(BN)^{2+}}$, respectively, when the DIP-EOMCCSDT(4$h$-2$p$) method is employed. The same is generally true when examining the DIPs of Cl₂, Br₂, and HBr shown in Table 4, where the DIP-EOMCCSDT(4$h$-2$p$) approach offers similar improvements. The only exceptions are the $c\>{{}^{1}}\Sigma_{u}^{-}$ state of $\mathrm{(Cl_{2})}^{2+}$ and the three states of $\mathrm{(HBr)}^{2+}$, for which the DIP-EOMCCSD(4$h$-2$p$)/cc-pVQZ DIPs are already very accurate. We conclude by pointing out that the DIP-EOMCCSD(T)(a)(4$h$-2$p$) method, which offers significant savings in the computational effort compared to full DIP-EOMCCSDT(4$h$-2$p$), reproduces the DIPs of H₂O, CH₄, and BN reported in Table 3 as well as the DIPs of Cl₂, Br₂, and HBr considered in Table 4 resulting from the parent DIP-EOMCCSDT(4$h$-2$p$) calculations to within 0.02 eV. While we will continue testing the DIP-EOMCCSD(T)(a)(4$h$-2$p$) approach against the DIP-EOMCCSDT(4$h$-2$p$) and other high-accuracy data, its excellent performance in this study is encouraging.

In summary, we presented the fully factorized and programmable equations defining the DIP-EOMCCSDT(4$h$-2$p$) approach and the perturbative approximation to it abbreviated as DIP-EOMCCSD(T)(a)(4$h$-2$p$). We incorporated the resulting computer codes into the open-source CCpy package available on GitHub. We applied the DIP-EOMCCSDT(4$h$-2$p$) and DIP-EOMCCSD(T)(a)(4$h$-2$p$) methods and their DIP-EOMCCSD(3$h$-1$p$) and DIP-EOMCCSD(4$h$-2$p$) predecessors to the vertical DIPs of H₂O, CH₄, and BN, as described by the aug-cc-pVTZ basis set, and Cl₂, Br₂, and HBr, using the cc-pVTZ and cc-pVQZ bases and the spin-free two-component SFX2C-1e treatment of the scalar-relativistic effects. We demonstrated that with the exception of the higher-lying $c\>{{}^{1}}\Sigma_{u}^{-}$ state of $\mathrm{(Br_{2})}^{2+}$, the DIP values computed with DIP-EOMCCSDT(4$h$-2$p$) are not only in generally good agreement with the available high-accuracy theoretical or experimental reference data, but also more accurate than those obtained with DIP-EOMCCSD(4$h$-2$p$), which uses CCSD instead of CCSDT to construct the underlying similarity-transformed Hamiltonian. The DIP-EOMCCSD(T)(a)(4$h$-2$p$) method, which avoids the most expensive steps of DIP-EOMCCSDT(4$h$-2$p$), turned out to be similarly effective, recovering the vertical DIPs of H₂O, CH₄, BN, Cl₂, Br₂, and HBr obtained with its DIP-EOMCCSDT(4$h$-2$p$) parent to within 0.02 eV. Our future plans include the development of nonperturbative ways of reducing costs of the DIP-EOMCCSDT(4$h$-2$p$) calculations through the active-space treatments of CCSDTOliphant and Adamowicz (1992); Piecuch, Oliphant, and Adamowicz (1993); Piecuch, Kucharski, and Bartlett (1999) and 4$h$-2$p$ amplitudes Shen and Piecuch (2013, 2014) and the use of frozen natural orbitals, combined with Cholesky decomposition and density fitting techniques, which will also be useful in improving our description of relativistic effects following the four-component methodology of Ref. Surjuse et al., 2022.