Benchmarking the semi-stochastic CC(P;Q) approach for singlet–triplet gaps in biradicals

There has been great progress in ab initio computational quantum chemistry, but an accurate description of multi-reference (MR) situations, such as molecular potential energy surfaces along bond breaking coordinates, electronic spectra of radical and biradical species, and excited states dominated by two-electron transitions, continues to represent a major challenge. While a traditional way of addressing this challenge has been to use MR approaches that rely on multi-dimensional model spaces spanned by multiple reference determinants,^1–8 in this work we focus on methods based on the single-reference (SR) coupled-cluster (CC) theory^9–13 (cf. Refs. 14 and 15 for selected reviews), in which the ground electronic state is expressed using the exponential wave function ansatz $Ψ=eTΦ$⁠, where $T=∑n=1NTn$ is the cluster operator, N is the total number of correlated electrons in the system, T_n is the n-body (n-particle–n-hole or n-tuply excited) component of T, and $Φ$ is the reference determinant that serves as a Fermi vacuum. By truncating T at various excitation ranks, one obtains the well-known hierarchy of SRCC approximations, including the CC method with singles and doubles (CCSD),^16–19 where T is truncated at T₂, the CC method with singles, doubles, and triples (CCSDT),^20–23 where T is truncated at T₃, the CC method with singles, doubles, triples, and quadruples (CCSDTQ),^24–27 where T is truncated at T₄, and so on. As long as the number of strongly correlated electrons is not too large, the CCSD, CCSDT, CCSDTQ, etc., hierarchy and its extensions to excited states and properties other than energy via the equation-of-motion (EOM)^28–35 and linear response (LR)^36–44 CC frameworks rapidly converge to the exact, full configuration interaction (FCI), limit.¹⁵ As a result, the CCSDT, CCSDTQ, and similar methods and their EOM and LR extensions are capable of accurately describing typical MR situations, such as bond rearrangements in chemical reactions, singlet–triplet gaps in biradicals, and excited states having substantial double excitation character, via particle–hole excitations from a single determinant.

The convergence of the CCSD, CCSDT, CCSDTQ, etc., hierarchy toward FCI in the majority of chemical applications is fast, but to achieve a quantitative description, one has to go beyond the basic CCSD level and face high, often prohibitive, computational costs, such as the iterative $no3nu5$ steps of CCSDT or the even more demanding $no4nu6$ steps of CCSDTQ (n_o and n_u are the numbers of correlated occupied and unoccupied orbitals, respectively). Thus, one of the key challenges in the development of SRCC methodologies has been the incorporation of many-electron correlation effects due to higher-than-two-body components of the cluster operator T without running into the enormous computational costs of the high-level parent methods, such as CCSDT or CCSDTQ, while avoiding the failures of perturbative approximations of the CCSD(T)^45,46 type when bond breaking, biradicals, and other MR situations are examined.^{14,15,47–49}

One of the most promising ways of addressing the above challenge within the SRCC framework is offered by the CC(P;Q) formalism, which allows one to correct the CC and EOMCC energies obtained using conventional as well as unconventional truncations in the cluster and EOM excitation operators for the missing many-electron correlation effects of interest using the moment expansions developed in Refs. 50–53. Focusing on the ground electronic states (or, in general, the lowest states of symmetries of interest, as in the case of many lowest-energy singlet, doublet, or triplet states), the CC(P;Q) methods based on correcting the CC energies obtained with the conventional truncations in T at a given excitation rank reduce to the left-eigenstate completely renormalized (CR) CC approaches, such as CR-CC(2,3),^54–58 in which one corrects the CCSD energies for the leading T₃ correlations, and its higher-order extensions.^53,59–62 The CR-CC(2,3) method improves poor performance of CCSD(T) in situations involving covalent bond breaking^{50,54–57,60,63,64} and certain types of noncovalent interactions,^59,65 while retaining the relatively inexpensive iterative $no2nu4$ and noniterative $no3nu4$ steps of CCSD(T), but neither CR-CC(2,3) nor its CR-CCSD(T)^47,48,66,67 and locally renormalized CCSD(T)⁶⁸ predecessors, nor their CCSD(2)_T,^69–72 CCSD(T)_Λ,^73–75 Λ-CCSD(T),^76,77 and CCSD(T–n)^78,79 counterparts, designed to improve CCSD(T) as well, are accurate enough when the coupling of the low-order T_n components, such as T₁ and T₂, and their higher-order counterparts with n > 2, e.g., T₃, becomes larger.^51,52 One can address this issue by solving for T₁ and T₂ clusters in the presence of the dominant higher-than-doubly excited components of T identified with the help of active orbitals, as in the CCSDt, CCSDtq, and similar schemes,^{27,31,32,49,80–85} and account for the remaining correlation effects of interest (e.g., the T₃ correlations not captured by CCSDt) using the CC(P;Q) moment expansions, but the resulting CC(t;3), CC(t,q;3), CC(t,q;3,4), etc. hierarchy,^{50–53,59,65} while accurately reproducing the parent CCSDT and CCSDTQ energetics at the fraction of the computational costs, requires choosing user- and system-dependent active orbitals. The question arises if the CC(P;Q) methodology can be similarly successful without having to define active orbitals tuned to a given calculation.

To address this question, we have recently proposed an automated way of identifying the leading T_n components with n > 2 within the CC(P;Q) framework by merging the deterministic CC(P;Q) formalism with the stochastic configuration interaction (CI) Quantum Monte Carlo (QMC)^86–90 and CC Monte Carlo (CCMC)^91–94 propagations.^95–97 As shown in Ref. 95, the results of the semi-stochastic CC(P;Q) computations, in which the lists of determinants needed to define higher-than-doubly excited cluster components are extracted from the FCIQMC and CCSDT-MC runs, rapidly converge to the parent high-level CC energetics, represented in Ref. 95 by CCSDT, out of the early stages of the underlying QMC propagations, even when electronic quasi-degeneracies and the associated T₃ correlations become substantial. In more recent studies, it has been demonstrated that the semi-stochastic CC(P;Q) calculations exhibit a similarly fast convergence toward the parent CCSDT and CCSDTQ energies when the FCIQMC propagations are replaced by their less expensive truncated CISDT-MC and CISDTQ-MC analogs⁹⁷ and that the CIQMC-driven CC(P;Q) framework can be very effective in converging the EOMCCSDT energies of excited states, including states of the same symmetry as the ground state, states of different symmetries, and challenging excited states dominated by two-electron transitions⁹⁶ (cf., also, Ref. 98). Our group has also explored the deterministic alternative to the semi-stochastic CC(P;Q) methodology by replacing the CIQMC and CCMC propagations, needed to identify the dominant higher-than-doubly excited components of T, by the sequences of Hamiltonian diagonalizations defining the selected CI approach abbreviated as CIPSI,^99–101 obtaining encouraging results.¹⁰² 

In the present work, we investigate the ability of the semi-stochastic CC(P;Q) formalism, driven by either FCIQMC or CISDTQ-MC propagations, to converge the energies of the lowest singlet and triplet states and the corresponding singlet–triplet gaps, ΔE_S-T, resulting from the high-level CCSDT calculations for the selected organic biradicals, including methylene, cyclobutadiene, cyclopentadienyl cation, and trimethylenemethane, and a prototype magnetic system, represented by the (HFH)⁻ ion. Understanding electronic structure and properties of organic biradicals is of major significance in areas such as chemical reaction dynamics,¹⁰³ molecular magnetism,¹⁰⁴ spintronics,^105,106 nonlinear optics,¹⁰⁷ photochemical pathways,^108–110 and photovoltaics.^111–113 The linear (HFH)⁻ species, in which one simultaneously stretches both H–F bonds, is a model inorganic magnetic system, in which two paramagnetic centers carrying unpaired spins associated with the terminal hydrogen atoms are coupled via a polarizable diamagnetic bridge formed by F⁻.¹¹⁴ It is well established that an accurate determination of the ordering of the low-lying singlet and triplet states and ΔE_S-T values in biradicals, which are among their most important physical characteristics, remains a significant challenge for computational quantum chemistry, even when high-level ab initio wave function methods are used in the calculations.^{52,57,115–126} This is because one has to balance substantial nondynamical correlation effects needed for an accurate description of the low-spin singlet state, which can be either of the open-shell type involving the nearly degenerate singly occupied valence molecular orbitals (MOs) or of the multi-configurational type mixing two closed-shell determinants involving the highest occupied and lowest unoccupied MOs, HOMO and LUMO, respectively [the determinant in which HOMO is doubly occupied and the doubly excited determinant of the (HOMO)² → (LUMO)² type], with the dynamical correlations of the high-spin triplet state, which has a predominantly SR character. If we limit ourselves to the black-box SRCC methodologies, the only methods that can provide a reliable and balanced description of the lowest singlet and triplet states and singlet–triplet gaps in biradicals and magnetic systems are the high-level CC approaches, beginning with CCSDT, and their particle-nonconserving double electron-attachment (DEA) and double ionization potential (DIP) EOMCC counterparts,^{123,125–136} especially those that incorporate the high-rank 4-particle–2-hole (4p2h) and 4-hole–2-particle (4h2p) correlations on top of the CC treatment of the underlying closed-shell cores.^{123,125,126,135,136} It is interesting to examine if our semi-stochastic, CIQMC-driven, CC(P;Q) methodology can be as successful in recovering the CCSDT results for the lowest singlet and triplet states and ΔE_S-T gaps of the biradical and magnetic systems listed above as in the previously reported benchmark studies of bond dissociations,^95,97 chemical reaction pathways,^95,97 and excited electronic states.⁹⁶ 

We recall that the CC(P;Q) formalism is a natural generalization of the biorthogonal moment expansions of Refs. 54–56, which in the past resulted in the CR-CC/EOMCC triples corrections to the CCSD/EOMCCSD energies, such as CR-CC(2,3),^54–58 CR-EOMCC(2,3),^56,137 and δ-CR-EOMCC(2,3),¹³⁸ and their extensions to quadruples,^53,59–62 to unconventional truncations in the cluster and EOM excitation operators. In the case of the ground electronic state, or the lowest-energy state of symmetry other than the ground state that can be treated using the SRCC framework, the CC(P;Q) calculation is a two-step procedure. In the first step, abbreviated as CC(P), we solve the CC amplitude equations in the subspace of the many-electron Hilbert space, $H(P)$⁠, referred to as the P space, spanned by the excited determinants that together with the reference determinant |Φ⟩ dominate the target state |Ψ⟩ of interest, for the suitably truncated form of the cluster operator T consistent with the content of the P space, designated as T^(P), and the corresponding energy E^(P). In the second step, we correct the CC(P) energy E^(P) for the many-electron correlation effects captured with the help of the excited determinants that span another subspace of the Hilbert space, called the Q space and denoted as $H(Q)$⁠, using the noniterative correction δ(P;Q) derived from the CC(P;Q) moment expansion introduced in Refs. 50–52. The final CC(P;Q) energy is determined as

where the explicit formulas for the correction δ(P;Q) in terms of the generalized moments^66,67,139 of the CC(P) equations, which correspond to projections of these equations on the complementary Q-space determinants, and the left eigenstate ⟨Φ|(1 + Λ^(P)) of the similarity-transformed Hamiltonian $H̄(P)=e−T(P)HeT(P)$ in the P space, with Λ^(P) representing the relevant hole–particle deexcitation operator, can be found in Refs. 50–53 (see, also, Refs. 59, 65, 95–97, and 102). In the most complete variant of the CC(P;Q) formalism adopted in this work, one uses the Epstein–Nesbet denominators involving $H̄(P)$ in determining the δ(P;Q) corrections. One could employ the Møller–Plesset denominators instead, but the Epstein–Nesbet form has been shown to be generally more accurate (see, e.g., Refs. 52, 53, 95, and 97).

The main advantage of the CC(P;Q) formalism compared to its CR-CC predecessors is the flexibility in defining the underlying P and Q spaces that allows us to relax the lower-order T₁ and T₂ components of the cluster operator T in the presence of their higher-order counterparts, such as the leading T₃ contributions, which the CR-CC(2,3), CCSD(2)_T, Λ-CCSD(T), and other noniterative triples corrections to CCSD cannot do. As explained in Sec. I, one can incorporate the coupling among the lower-rank T₁ and T₂ and higher-rank T_n components with n > 2 within the CC(P;Q) framework by including the dominant higher-than-doubly excited determinants identified with the help of active orbitals or selected CI diagonalizations in the P space, using the δ(P;Q) corrections to capture the remaining correlations of interest, but this is not what we do in this work. Here, we explore the alternative approach that utilizes the semi-stochastic CC(P;Q) methodology of Refs. 95–97, in which the leading higher-than-doubly excited determinants included in the underlying P spaces are identified using CIQMC runs, whereas the corresponding Q spaces, needed to determine corrections δ(P;Q), are populated by the remaining determinants of interest not captured by CIQMC at a given propagation time τ. Given our interest in converging the CCSDT energies of the lowest singlet and triplet states of biradical systems, we adopt in this study the following CC(P;Q) algorithm:

1. Initiate the desired CIQMC (in this study, FCIQMC or CISDTQ-MC) propagation by placing a certain number of walkers on the appropriate reference determinant |Φ⟩. In calculations for the lowest singlet states of symmetries of interest, we used the relevant restricted Hartree–Fock (RHF) determinants as reference functions. In the case of the lowest triplet states, we used the restricted open-shell Hartree–Fock (ROHF) references.

2. After a certain number of CIQMC time steps, referred to as MC iterations, or, equivalently, at some propagation time τ > 0, extract a list of triply excited determinants captured by CIQMC, i.e., those triples that are populated by at least one walker.

3. Solve the CC(P) amplitude equations and the associated left eigenstate problem involving $H̄(P)$ in the P space spanned by all singly and doubly excited determinants and the subset of triply excited determinants identified in step 2, using the same reference |Φ⟩ as that chosen to initiate the CIQMC run, to obtain the cluster operator T^(P), defined as $T1+T2+T3(MC)$⁠, and its $Λ(P)=Λ1+Λ2+Λ3(MC)$ companion, where the list of triples in $T3(MC)$ and $Λ3(MC)$ is extracted from the CIQMC propagation at time τ. Calculate the energy E^(P).

4. Correct the CC(P) energy E^(P) obtained in step 3 for the remaining T₃ correlations using Eq. (1), where the Q space needed to determine the δ(P;Q) correction consists of the triply excited determinants not captured by the CIQMC propagation at time τ, to obtain the CC(P;Q) energy E^(P+Q).

5. Repeat steps 2–4 at some later CIQMC propagation time τ′ > τ to check the convergence of the CC(P;Q) energies E^(P+Q). Stop the calculation if the consecutive E^(P+Q) values do not change within an a priori defined convergence threshold or if the stochastically determined P space captures a large enough fraction of the triply excited determinants sufficient to achieve the desired accuracy relative to the parent CCSDT approach.

In this work, we explored two types of CIQMC propagations to carry out our semi-stochastic CC(P;Q) computations. For smaller systems, including methylene, (HFH)⁻, and cyclobutadiene, we used FCIQMC, which allows the CIQMC algorithm to explore the entire many-electron Hilbert space. For the two largest biradical species considered in our calculations, namely, cyclopentadienyl cation and trimethylenemethane, we used the truncated CISDTQ-MC approach, in which spawning beyond quadruply excited determinants relative to reference |Φ⟩ is disallowed, reducing computation costs compared to FCIQMC, especially in later stages of CIQMC propagations. As shown in Ref. 97, convergence of the CC(P;Q) energies toward their high-level CC parents, such as CCSDT, is not affected by the type of the CIQMC approach used to identify the relevant higher-than-doubly excited determinants (in this work, the triply excited determinants), so choosing CISDTQ-MC as a substitute for FCIQMC has no effect on our conclusions. Following our earlier studies of the semi-stochastic CC(P;Q) methodology^95–97 and related work,⁹⁸ in all of the CIQMC calculations reported in this article, we relied on the initiator CIQMC (i-CIQMC) algorithm introduced in Ref. 87, based on integer walker numbers, in which only the determinants populated by n_a or more walkers are allowed to spawn new walkers onto empty determinants, but the semi-stochastic CC(P;Q) procedure summarized above is flexible and could be merged with other CIQMC techniques developed in recent years, such as those described in Refs. 89 and 90. While the convergence of the CC(P;Q) energies of the lowest singlet and triplet states of biradicals considered in this study and of the corresponding singlet–triplet gaps obtained with the help of i-CIQMC propagations toward the parent CCSDT results is very fast, we will examine the utility of the newer CIQMC algorithms in the future work.

As explained in Refs. 95–97 (cf., also, Ref. 98), the semi-stochastic CC(P;Q) procedure, as summarized above, offers considerable savings in the computational effort compared to its CCSDT parent. Indeed, the computational times and walker populations characterizing the early stages of CIQMC propagations, which are sufficient to produce enough information for the subsequent CC(P;Q) calculations to recover the target CCSDT energetics to within fractions of a millihartree, are very small compared to the converged CIQMC runs. Next, the CC(P) calculations using small fractions of the triply excited determinants captured in the early stages of the CIQMC propagations are much faster than the parent CCSDT computations using all triples. Finally, the computational times required to determine corrections δ(P;Q) due to the remaining T₃ correlations not captured by the preceding CC(P) runs, which scale no worse than $∼2no3nu4$⁠, are much less than the $no3nu5$ timings associated with full CCSDT iterations.

Another interesting feature of the semi-stochastic CC(P;Q) algorithm, as defined by the above steps 1–5, is a systematic behavior of the E^(P+Q) energies as τ varies from 0 to ∞. Indeed, the initial, τ = 0, CC(P;Q) energy, where the P space is spanned by all singly and doubly excited determinants and the Q space consists of all triples, is identical to that obtained using the CR-CC(2,3) correction to CCSD. On the other hand, when τ = ∞, so that the P space is spanned by all singly, doubly, and triply excited determinants and the Q space is empty, the CC(P;Q) energy E^(P+Q) becomes equivalent to its CCSDT parent. Thus, the CIQMC propagation time τ can be regarded as a continuation parameter connecting the approximate treatment of T₃ clusters, represented by CR-CC(2,3), with their complete description offered by full CCSDT. As τ → ∞, the uncorrected CC(P) energies E^(P) converge to their CCSDT counterparts too, but, as shown in Refs. 95–97, and as illustrated via our calculations of the singlet–triplet gaps discussed in Sec. III, the convergence toward CCSDT is in this case substantially slower, since the initial, τ = 0, CC(P) energy is equivalent to that of CCSD, where T₃ = 0. The main role of corrections δ(P;Q) in the context of the semi-stochastic CC(P;Q) algorithm considered in this study is to accelerate convergence of the underlying CC(P) energies toward the desired CCSDT limit.

As explained in Sec. I, in order to assess the performance of our semi-stochastic, CIQMC-driven, CC(P;Q) methodology in converging the full CCSDT data for the singlet–triplet gaps and the corresponding singlet- and triplet-state energies of biradical systems, we applied it to methylene, (HFH)⁻, cyclobutadiene, cyclopentadienyl cation, and trimethylenemethane. Following our earlier studies of the singlet–triplet gaps in the same systems using the deterministic CC(P;Q)⁵² and DEA/DIP-EOMCC^{123,124,135,136} approaches, we used the aug-cc-pVTZ basis set^140,141 for methylene, the 6-31G(d,p) basis^142,143 for the (HFH)⁻ ion, and the cc-pVDZ basis set¹⁴⁰ for cyclobutadiene, cyclopentadienyl cation, and trimethylenemethane. In the case of methylene and trimethylenemethane, we focused on the ability of the semi-stochastic CC(P;Q) algorithm to converge the adiabatic ΔE_S-T data obtained with CCSDT. When executing the semi-stochastic CC(P;Q) calculations for (HFH)⁻, cyclobutadiene, and cyclopentadienyl cation, we focused on recovering the CCSDT values of the vertical singlet–triplet gaps. Throughout this work, we define ΔE_S-T as E_S − E_T, where E_S and E_T are the electronic energies of the corresponding singlet and triplet states, i.e., the positive ΔE_S-T value implies that triplet is lower in energy.

All of the CC calculations reported in this article were performed using our group’s standalone codes, interfaced with the RHF, ROHF, and integral transformation routines in the GAMESS package,^144,145 which were originally developed in Refs. 50–53 and extended to the stochastically generated P spaces for the use in CC(P) and CC(P;Q) computations in Refs. 95–98. The i-FCIQMC [methylene, (HFH)⁻, and cyclobutadiene] and i-CISDTQ-MC (cyclopentadienyl cation and trimethylenemethane) calculations, needed to generate the lists of triples for the semi-stochastic CC(P) and CC(P;Q) runs, were carried out with the HANDE software.^146,147 Each of the i-FCIQMC and i-CISDTQ-MC propagations was initiated by placing 1500 walkers on the relevant reference determinant. The CIQMC time step δτ and the initiator parameter n_a were set at 0.0001 a.u. and 3, respectively. In all post-Hartree–Fock calculations, the core MOs correlating with the 1s orbitals of the C and F atoms were kept frozen. If the true point group of the biradical system of interest was not Abelian, we used its largest Abelian subgroup, since our CC codes interfaced with GAMESS and the CIQMC routines in HANDE cannot handle non-Abelian symmetries.

We begin the discussion of our results by analyzing the performance of the semi-stochastic CC(P;Q) approach in converging the CCSDT energies of the ground (X³B₁) and first-excited (A¹A₁) states of the methylene/aug-cc-pVTZ system and the adiabatic gap between them. The C_2v-symmetric geometries of CH₂ in the two states, optimized using FCI and the [5s3p/3s] triple zeta basis set of Dunning¹⁴⁸ augmented with two sets of polarization functions (TZ2P), were taken from Ref. 149. The electronically nondegenerate triplet ground state has a predominantly SR nature dominated by the $(1a1)2(2a1)2(1b2)2(3a1)1(1b1)1$ configuration, whereas the first-excited singlet state exhibits a significant MR character requiring a linear combination of the $(1a1)2(2a1)2(1b2)2(3a1)2$ and doubly excited $(1a1)2(2a1)2(1b2)2(1b1)2$ closed-shell determinants for a proper zeroth-order description. Because of these fundamentally different characteristics of the X³B₁ and A¹A₁ states, a well-balanced and accurate treatment of dynamical and nondynamical correlation effects is the key to a reliable description of the singlet–triplet gap in methylene. It is, therefore, unsurprising that one usually resorts to methods of the MRCI^149–154 or MRCC^155–158 type, or to the high-level SRCC theories that account for higher-than-doubly excited clusters in an iterative manner, such as full CCSDT used in Refs. 23 and 52, to accomplish this goal (for other examples of high-level ab initio calculations for the X³B₁ and A¹A₁ states of methylene, see Refs. 123, 135, 136, and 159 and references therein). The CCSDT results for the adiabatic singlet–triplet gap in methylene, which are of interest in the present study, are indeed very accurate. As shown, for example, in Ref. 52, the difference between the adiabatic ΔE_S-T value obtained in the CCSDT/TZ2P calculations and the corresponding FCI result of 11.14 kcal/mol¹⁴⁹ is only 0.11 kcal/mol or 38 cm⁻¹. As demonstrated in Ref. 52 as well, the purely electronic A¹A₁ − X³B₁ separation resulting from the CCSDT computations using the aug-cc-pVTZ basis employed in this work is only about 0.15 kcal/mol (∼50 cm⁻¹) higher than the experimentally derived value of 9.37 kcal/mol reported in Ref. 156, obtained by correcting the vibrationless adiabatic singlet–triplet gap determined in Ref. 160 for the relativistic and nonadiabatic (Born–Oppenheimer diagonal correction) effects estimated in Refs. 161 and 162, respectively. It is, therefore, interesting to examine if the semi-stochastic CC(P;Q) approach advocated in this work is capable of reproducing the high-quality CCSDT/aug-cc-pVTZ data for the X³B₁ and A¹A₁ states of methylene and the adiabatic separation between them.

The results of our FCIQMC-driven CC(P;Q) calculations for the methylene/aug-cc-pVTZ system, reported as errors relative to the parent CCSDT data, and their CC(P) counterparts are shown in Table I and Fig. 1. The reference determinants |Φ⟩ used to initiate the i-FCIQMC propagations and to carry out the CC(P), CC(P;Q), CCSD, CR-CC(2,3), and CCSDT calculations were the ROHF determinant in the case of the X³B₁ state and the RHF determinant for the A¹A₁ state. The subsets of triply excited determinants needed to construct the P spaces used in the CC(P) and CC(P;Q) computations at various propagation times τ were the S_z = 1 triples of the B₁ symmetry captured during the i-FCIQMC run for the X³B₁ state and the S_z = 0 triples of the A₁ symmetry captured during the analogous run for the A¹A₁ state. Following the semi-stochastic CC(P;Q) algorithm described in Sec. II, the Q spaces needed to determine corrections δ(P;Q) were defined as the remaining triples not captured by i-FCIQMC.

Let us start our analysis by examining the CC(P) and CC(P;Q) data at τ = 0, where the P spaces do not contain any triply excited determinants. As shown in Table I, the CC(P) energies of the X³B₁ and A¹A₁ states at τ = 0, which are equivalent to those obtained using conventional CCSD, are above their CCSDT [i.e., τ = ∞ CC(P)] counterparts by 4.187 and 5.918 millihartree, respectively. This translates into a 380 cm⁻¹ or ∼11% error in the adiabatic ΔE_S-T value when compared to the 3328 cm⁻¹ singlet–triplet gap obtained with CCSDT. The situation improves when the CC(P;Q) corrections δ(P;Q) due to T₃ correlation effects, calculated by placing all triply excited determinants in the respective Q spaces, are added to the CC(P) energies. The τ = 0 CC(P;Q) or CR-CC(2,3) energy characterizing the X³B₁ state is only 0.177 millihartree above the parent CCSDT value, which is an error reduction relative to CCSDT compared to the underlying CC(P) result by a factor of ∼24. The δ(P;Q) correction improves the τ = 0 CC(P) energy of the more challenging A¹A₁ state as well, although the difference between the resulting CR-CC(2,3) energy and its CCSDT counterpart, of 0.656 millihartree, is almost four times larger than the analogous difference obtained for the X³B₁ state. As a result, the 105 cm⁻¹ error relative to CCSDT characterizing the adiabatic A¹A₁ − X³B₁ separation obtained in the τ = 0 CC(P;Q) or CR-CC(2,3) calculations, while considerably smaller than the 380 cm⁻¹ obtained in the underlying CC(P) (i.e., CCSD) runs, leaves room for further improvements. One can improve the CR-CC(2,3) energies of the X³B₁ and A¹A₁ states and the gap between them by enriching the P spaces defining the CC(P) calculations with the leading triply excited determinants identified using active orbitals and correcting the resulting CCSDt energies for the remaining T₃ correlations that have not been captured by CCSDt,⁵² but our objective here is to examine if one can accomplish the same, or improve the CC(t;3) results reported in Ref. 52 even further, by turning to the more black-box semi-stochastic CC(P;Q) methodology, in which the dominant triply excited determinants are identified with CIQMC.

The results in Table I and Fig. 1 show that when the τ = 0P spaces are augmented with the subsets of triply excited determinants captured in the i-FCIQMC runs at τ > 0 and, following the CC(P;Q) recipe, the resulting CC(P) energies are corrected for the remaining T₃ correlations, the convergence of the total electronic energies of the X³B₁ and A¹A₁ states and the adiabatic separation between them toward their CCSDT parents is rapid. We can see this already in the early stages of the i-FCIQMC propagations. For example, at τ = 0.8 a.u., i.e., after only 8000 δτ = 0.0001 a.u. MC iterations, the errors in the CC(P;Q) energies of the X³B₁ and A¹A₁ states and the corresponding ΔE_S-T value relative to CCSDT are 0.049 millihartree, 0.106 millihartree, and 13 cm⁻¹, respectively, substantially improving the CR-CC(2,3) [i.e., τ = 0 CC(P;Q)] calculations, which give 0.177 millihartree for the X³B₁ state, 0.656 millihartree for the A¹A₁ state, and 105 cm⁻¹ for ΔE_S-T. This confirms our expectation that the main source of errors in the CR-CC(2,3) computations, especially in the case of the more MR A¹A₁ state, which is characterized by larger T₃ effects, is the use of the unrelaxed T₁ and T₂ amplitudes obtained with CCSD in constructing the correction due to triples. The FCIQMC-based CC(P;Q) calculations at τ = 0.8 a.u., which use as little as 16% of all triply excited determinants to define the P space for the X³B₁ state and only 25% of all triples in the P space for the A¹A₁ state, are also more accurate than the CC(t;3) computations reported in Ref. 52, which produced the 0.130 millihartree, 0.409 millihartree, and 61 cm⁻¹ errors relative to CCSDT for the X³B₁ and A¹A₁ energies and ΔE_S-T, respectively. This is all very promising, especially if we realize that the i-FCIQMC propagations used to generate the lists of triples for our semi-stochastic CC(P;Q) runs, which work so well, are very far from convergence when τ = 0.8 a.u. Indeed, as seen in Table S.1 of the supplementary material, the total numbers of walkers at 8000 δτ = 0.0001 a.u. MC iterations, which are 132 689 in the case of the X³B₁ state and 165 564 for the A¹A₁ state, represent tiny fractions, 2.17% and 1.11%, respectively, of the total walker populations at τ = 20.0 a.u., where we stopped our i-FCIQMC propagations.

As demonstrated in Table I and Fig. 1, the convergence of the energies of the X³B₁ and A¹A₁ states and the gap between them resulting from the FCIQMC-driven CC(P;Q) calculations remains fast at the larger propagation times τ as well. For example, if we allow i-FCIQMC to populate the respective P spaces with about 26%–38% of all triples, which happens after 20 000 δτ = 0.0001 a.u. MC iterations, the CC(P;Q) energies of the X³B₁ and A¹A₁ states and the resulting singlet–triplet gap become practically indistinguishable from the parent CCSDT data, with errors in the total electronic energies and ΔE_S-T being only ∼20 microhartree and 2 cm⁻¹, respectively. As shown in Table S.1 of the supplementary material, walker populations characterizing the X³B₁ and A¹A₁ states produced by i-FCIQMC at 20 000 δτ = 0.0001 a.u. MC time steps are still very small compared to the total numbers of walkers at τ = 20.0 a.u., where we terminated our i-FCIQMC propagations (4.11% for the X³B₁ state and 2.17% in the case of the A¹A₁ state). It is also interesting to note that the more MR A¹A₁ state requires a higher fraction of triply excited determinants in the P space than its SR X³B₁ counterpart to achieve similar accuracy levels in the semi-stochastic CC(P;Q) computations for both states. For example, the i-FCIQMC propagation has to capture about 25% of all triples, for the inclusion in the P space, if we are to reduce errors relative to CCSDT in the CC(P;Q) calculations for the A¹A₁ state to ∼0.1 millihartree. In the case of the X³B₁ state, the analogous fraction of triples is about 10% (cf. Table I). This highlights the importance of balancing the SR triplet state with the more MR singlet state in obtaining accurate ΔE_S-T estimates, which is not a problem for the semi-stochastic CC(P;Q) methodology because the underlying i-FCIQMC wave function sampling is very effective in identifying the dominant higher-than-doubly excited determinants, to be included in the relevant P spaces, and the δ(P;Q) corrections to the CC(P) energies take care of the remaining correlation effects of interest.

Before concluding this subsection and discussing other molecular examples, we would like to comment on the effectiveness of the noniterative corrections δ(P;Q), adopted in the CC(P;Q) formalism, in accelerating convergence of the underlying CC(P) calculations toward CCSDT. The CC(P) and CC(P;Q) error curves shown in Fig. 1 illustrate this best. It is clear from this figure that the CC(P;Q) energies of the X³B₁ and A¹A₁ states [Figs. 1(a) and 1(b)] and the corresponding ΔE_S-T values [Fig. 1(c)] converge to the parent CCSDT data much faster than in the case of the uncorrected CC(P) computations. One can see the same by inspecting the numerical data shown in Table I. In this context, it is worth commenting on the CC(P) and CC(P;Q) results obtained after 8000 MC iterations. In that case, the CC(P;Q) calculations reduce the ∼2.4 millihartree errors relative to CCSDT characterizing the CC(P) energies of the X³B₁ and A¹A₁ states to 0.1 millihartree or less, which is a desired behavior, but the CC(P;Q) ΔE_S-T value is less accurate than that obtained with the uncorrected CC(P). One should not read too much into this though. The fact that the CC(P;Q) calculations at 8000 MC iterations increase the very small 3 cm⁻¹ error obtained with CC(P) to 13 cm⁻¹ is a coincidence arising from the accidental cancellation of errors in the CC(P) total electronic energies obtained at this particular propagation time. Indeed, when the later stages of the i-FCIQMC propagations are considered, the differences between the CC(P) and CCSDT values of ΔE_S-T become increasingly negative, reaching −107 cm⁻¹ at 50 000 MC iterations, before eventually converging to the CCSDT limit, whereas the corresponding CC(P;Q) results display a smooth behavior, rapidly approaching CCSDT. In particular, they reduce the relatively large negative error value obtained for ΔE_S-T in the CC(P) calculations at 50 000 MC iterations to a numerical 0 cm⁻¹. This highlights, once again, the ability of the CC(P;Q) corrections δ(P;Q) to offer a well-balanced description of the lowest singlet and triplet states in methylene, in addition to improving the individual state energies.

Our next example is the linear, D_∞h-symmetric, (HFH)⁻ anion, a prototype magnetic system in which unpaired spins of terminal hydrogen atoms couple to singlet and triplet states via a polarizable diamagnetic bridge of F⁻.¹¹⁴ The energies of the lowest two electronic states of the (HFH)⁻ system, including the singlet ground state $X1Σg+$ and the first-excited triplet state $A3Σu+$⁠, and the vertical gap between them, which is proportional to the magnetic exchange coupling constant J and which should approach zero as both H–F bonds are stretched to infinity, were used in the past to test various quantum chemistry approaches.^{52,53,57,58,114,135,136,163,164} Among them were methods developed in our group, including CR-CC(2,3),^57,58 CR-CC(2,4),⁵³ the DIP-EOMCC approaches with full and active-space treatments of 4h2p correlations of top of CCSD,^135,136 and the active-orbital-based CC(t;3), CC(t,q;3), and CC(t,q;3,4) hierarchy.^52,53 Here, we test the alternative to CC(t;3) offered by the semi-stochastic, FCIQMC-driven, CC(P;Q) algorithm aimed at the CCSDT energetics. As in our previous studies,^{52,53,57,58,135,136} we used the 6-31G(d,p) basis set and several stretches of both H–F bonds, including R_H–F = 1.50, 1.75, 2.00, 2.50, and 4.00 Å, where R_H–F is the distance between the hydrogen and fluorine nuclei.

An accurate computation of the singlet–triplet gap in the (HFH)⁻ system is complicated by the fact that, unlike the $A3Σu+$ state, which is weakly correlated and well represented by a single ROHF determinant, its ground-state counterpart $X1Σg+$ displays a substantial MR character that includes a significant contribution from the doubly excited (HOMO)² → (LUMO)² determinant, in addition to the RHF configuration. The MR character of the $X1Σg+$ state, which is already noticeable at shorter H–F separations and which substantially strengthens as R_H–F increases, can be illustrated by the ratio of the FCI expansion coefficients at the (HOMO)² → (LUMO)² and RHF determinants or the equivalent T₂ cluster amplitude extracted from FCI, which increases, in absolute value, from 0.38 at R_H–F = 1.50 Å to 1.17 at R_H–F = 4.00 Å, when the 6-31G(d,p) basis is employed^57,58 (the HOMO and LUMO have different symmetries, σ_g and σ_u, respectively, so that the HOMO → LUMO T₁ amplitude is zero). As a result of all this, it is difficult to balance the lowest two states of the (HFH)⁻ system in a single quantum chemistry calculation, especially when the SRCC framework using the RHF reference determinant for the $X1Σg+$ state and the ROHF reference for the $A3Σu+$ state is employed. Indeed, as shown in Refs. 52, 57, and 58, the differences between the energies obtained in the CCSD/6-31G(d,p) computations and their FCI counterparts at R_H–F = 1.50 Å are 12.674 millihartree for the $X1Σg+$ state and only 2.628 millihartree when the $A3Σu+$ state is considered. The analogous differences at R_H–F = 2.00 Å are 19.398 and 2.068 millihartree, respectively. The observed large discrepancies between the errors in the CCSD energies for the $X1Σg+$ and $A3Σu+$ states translate into a poor description of the singlet–triplet gaps. One can see this by comparing the ΔE_S-T values resulting from the RHF/ROHF-based CCSD/6-31G(d,p) computations at R_H–F = 1.50, 1.75, 2.00, 2.50, and 4.00 Å with the corresponding FCI data. CCSD/6-31G(d,p) gives −7320, −1838, 1656, 3605, and 230 cm⁻¹, respectively, as opposed to −9525, −4911, −2147, −277, and 0 cm⁻¹ obtained with FCI.^52,57,58 If we are to improve the CCSD results within the SRCC framework, we must turn to higher-level theories, such as the CCSDT approach that interests us in this study,^52,53,163 CCSDTQ,⁵³ or the DIP-EOMCC methodology, especially after incorporating 4h2p correlations.^135,136 The CCSDT method is indeed very accurate, reducing the 2205, 3073, 3804, 3882, and 230 cm⁻¹ errors relative to FCI in the ΔE_S-T values obtained with CCSD/6-31G(d,p) at R_H–F = 1.50, 1.75, 2.00, 2.50, and 4.00 Å to 198, 270, 341, 420, and 58 cm⁻¹, respectively.^52,53,163 It also greatly improves the total electronic energies. Indeed, the differences between the CCSDT and FCI energies of the $X1Σg+$ and $A3Σu+$ states in the entire R_H–F = 1.50–4.00 Å region obtained using the 6-31G(d,p) basis set do not exceed 2.276 and 0.389 millihartree, respectively.^52,53 The analogous differences between the CCSD and FCI energies are as large as 20.546 millihartree for the former state and 2.628 millihartree when the latter state is considered. One can reduce the remaining small errors in the CCSDT results even further or practically eliminate them by using CCSDTQ⁵³ or the DIP-EOMCC approaches with 4h2p contributions,^135,136 but the objective of this study is to assess the performance of our semi-stochastic CC(P;Q) methodology in converging the CCSDT data.

The results of our FCIQMC-driven CC(P;Q)/6-31G(d,p) computations for the $X1Σg+$ and $A3Σu+$ states of the linear (HFH)⁻ system at the H–F distances R_H–F = 1.50, 1.75, 2.00, 2.50, and 4.00 Å and the corresponding ΔE_S-T values, along with the underlying CC(P) data, are reported in Tables II–IV and Fig. 2. In all of our CC(P) and CC(P;Q) computations and the underlying i-FCIQMC runs for the D_∞h-symmetric (HFH)⁻ system, we used the D_2h Abelian subgroup of D_∞h. In particular, the i-FCIQMC calculations for the $X1Σg+$ and $A3Σu+$ states were set up to converge the lowest-energy states of the¹A_g(D_2h) and³B_1u(D_2h) symmetries. As a result, the subsets of triply excited determinants used to construct the P spaces for the subsequent CC(P) and CC(P;Q) computations for the $X1Σg+$ state at the various R_H–F and τ values considered in this work were defined as the S_z = 0 triples of the A_g(D_2h) symmetry captured in the underlying i-FCIQMC propagations. Similarly, the subsets of triply excited determinants used to design the P spaces for the CC(P) and CC(P;Q) calculations for the $A3Σu+$ state were the S_z = 1 triples of the B_1u(D_2h) symmetry extracted from i-FCIQMC. In analogy to all other CC(P;Q) computations performed in this work, the Q spaces used to determine the δ(P;Q) corrections to the CC(P) energies were defined as the remaining triples not captured by the respective i-FCIQMC runs.

As shown in Table II, and in line with our earlier CC(P;Q) work⁵² and the above remarks, the CC(P) energies of the $X1Σg+$ state of (HFH)⁻ obtained at τ = 0, which are identical to those resulting from the conventional CCSD calculations reported in Refs. 52, 57, and 58, are characterized by large errors relative to their τ = ∞, i.e., CCSDT, parents. Indeed, the differences between the τ = 0 and τ = ∞ CC(P) energies for the $X1Σg+$ state increase from 11.412 millihartree at R_H–F = 1.50 Å to more than 17 millihartree at R_H–F = 2.00 and 2.50 Å. These differences become smaller at large H–F separations, represented in our calculations by R_H–F = 4.00 Å, where the D_∞h-symmetric (HFH)⁻ system is essentially dissociated into the stretched hydrogen molecule, which has only two electrons, so that CCSD becomes exact, and the closed-shell fluoride ion, which has the electronic structure of the neon atom and which is characterized by small T_n correlations with n > 2, but they remain large when the R_H–F values are smaller. This should be contrasted with the small, ∼1–2 millihartree, differences between the τ = 0 and τ = ∞ CC(P) energies obtained at all values of R_H–F for the predominantly SR $A3Σu+$ state (see Table III). As already alluded to above, and as shown in Table IV, this imbalance in the description of the $X1Σg+$ and $A3Σu+$ states by the CCSD, i.e., τ = 0 CC(P), calculations gives rise to large errors in the resulting ΔE_S-T values relative to their τ = ∞ (CCSDT) counterparts, which range from 2007 to 3462 cm⁻¹ in the R_H–F = 1.50–2.50 Å region. Once again, these errors become small at large H–F separations, such as R_H–F = 4.00 Å used in this work, where (HFH)⁻ is more or less equivalent to the stretched H₂ and F⁻, resulting in the nearly degenerate singlet and triplet states and the 172 cm⁻¹ difference between the τ = 0 and τ = ∞ CC(P) values of ΔE_S-T, but at shorter H–F distances they are large and comparable to or even larger than the singlet–triplet gap values provided by CCSDT or FCI.

The situation dramatically changes, when the τ = 0 CC(P) or CCSD energies are corrected for T₃ correlations with the help of the noniterative correction δ(P;Q), as in the τ = 0 CC(P;Q) calculations, which are equivalent to the purely deterministic CR-CC(2,3) runs reported in Refs. 52, 53, 57, and 58. As shown in Tables II–IV, the τ = 0 CC(P;Q), i.e., CR-CC(2,3), energies of the $X1Σg+$ and $A3Σu+$ states at the various H–F distances considered in this study and the gaps between them are substantially more accurate than their uncorrected CC(P) (i.e., CCSD) counterparts. For example, the CR-CC(2,3) approach reduces the large, more than 17 millihartree, errors in the CCSD energies of the $X1Σg+$ state relative to their CCSDT [τ = ∞ CC(P) or CC(P;Q)] parents at R_H–F = 2.00 and 2.50 Å to ∼1–3 millihartree. We see similarly significant improvements in the CCSD energies of the $X1Σg+$ state by CR-CC(2,3) at other H–F distances, even at the “easiest” R_H–F = 4.00 Å value, where the triples correction δ(P;Q) is capable of reducing the already small, 1.907 millihartree, difference between the CCSD and CCSDT energies to the much smaller (in absolute value) 0.291 millihartree (see Table II). Consistent with our earlier studies,^52,53,57,58 the CR-CC(2,3) method performs even better when the weakly correlated $A3Σu+$ state is examined, reducing the ∼1–2 millihartree errors in the underlying CCSD energetics relative to our CCSDT target to about 0.2 millihartree (see Table III). As a result of all of these accuracy improvements, the singlet–triplet gap values obtained using CR-CC(2,3) are much closer to their CCSDT parents than their CCSD counterparts, reducing the 2007, 2803, 3462, 3462, and 172 cm⁻¹ errors relative to CCSDT obtained with CCSD at R_H–F = 1.50, 1.75, 2.00, 2.50, and 4.00 Å, respectively, by factors ranging from 6 at R_H–F = 2.50 Å to 72 at R_H–F = 1.50 Å, but, as shown in Table IV [see, also, Ref. 52, where one can find a comparison of the CCSD, CR-CC(2,3), and CCSDT ΔE_S-T data for additional H–F distances], the differences on the order of (−600)–(−300) cm⁻¹ between the CR-CC(2,3) and CCSDT singlet–triplet separations in the intermediate R_H–F = 2.00–3.00 Å region remain. The question arises if one can refine the CR-CC(2,3) results by enriching the P spaces used in the CC(P;Q) calculations, which in CR-CC(2,3) consist of only singles and doubles, with the subsets of triply excited determinants identified by i-FCIQMC propagations.

As shown in Tables II–IV and Fig. 2, once the leading triply excited determinants, captured using i-FCIQMC at τ > 0, are included in the respective P spaces and the δ(P;Q) corrections due to the remaining T₃ correlations are added to the energies obtained in the CC(P) calculations, the resulting CC(P;Q) values of the $X1Σg+$ and $A3Σu+$ energies and vertical gaps between them display very fast convergence toward their CCSDT counterparts. This is already observed when the i-FCIQMC propagation times are short, engaging tiny walker populations that are orders of magnitude smaller than those required to converge the i-FCIQMC runs, and the fractions of the triply excited determinants captured by i-FCIQMC are small. For example, after as few as 2000 δτ = 0.0001 a.u. MC iterations, where τ is only 0.2 a.u. and where, as shown in Table S.2 of the supplementary material, the total walker populations characterizing the underlying i-FCIQMC runs are 0.01%–0.11% of the respective numbers of walkers at τ = 20.0 a.u. [the termination time for our i-FCIQMC propagations for (HFH)⁻], the differences between the CC(P;Q) and CCSDT energies obtained for the strongly correlated $X1Σg+$ state are −0.035 millihartree for R_H–F = 1.50 Å, −0.056 millihartree for R_H–F = 1.75 Å, −0.110 millihartree for R_H–F = 2.00 Å, −0.583 millihartree for R_H–F = 2.50 Å, and −0.025 millihartree for R_H–F = 4.00 Å. Despite using only about 10%–30% of all triply excited determinants in the underlying P spaces, the FCIQMC-based CC(P;Q) energies of the $X1Σg+$ state obtained after 2000 MC iterations reduce the errors relative to CCSDT characterizing the CR-CC(2,3) [i.e., τ = 0 CC(P;Q)] computations in the R_H–F = 1.50–4.00 Å region by factors ranging from 5 to 13 (see Table II). In fact, with an exception of R_H–F = 2.00 and 2.50 Å, they are much more accurate than the results produced by the purely deterministic CC(t;3) analog of the semi-stochastic CC(P;Q) methodology, reported in Refs. 52 and 53. One can observe even more dramatic improvements over CR-CC(2,3) offered by the FCIQMC-driven CC(P;Q) approach, when the propagation time τ increases. For example, after 4000 MC iterations, where the i-FCIQMC propagations are still far from being converged (cf. the total walker populations used by our i-FCIQMC runs relative to the termination time τ = 20.0 a.u. in Table S.2 of the supplementary material) and the fractions of triples included in the stochastically determined P spaces, which range from 12% at R_H–F = 4.00 Å to 56% at R_H–F = 1.50 Å, remain relatively small, the differences between the CC(P;Q) and CCSDT energies obtained for the $X1Σg+$ state at R_H–F = 1.50, 1.75, 2.00, 2.50, and 4.00 Å are −28, −9, −17, −50, and −4 microhartree, respectively, reducing the errors relative to CCSDT that characterize the corresponding CR-CC(2,3) calculations by factors ranging from 12 to 86 [2 to 61 when compared to the CC(t;3) results reported in Refs. 52 and 53]. As shown in Table III, the performance of the FCIQMC-driven CC(P;Q) approach becomes even more impressive when the $A3Σu+$ state, which has a SR character, is examined. After 2000 δτ = 0.0001 a.u. MC time steps, the errors in the CC(P;Q) energies relative to their CCSDT parents obtained for the $A3Σu+$ state at R_H–F = 1.50, 1.75, 2.00, 2.50, and 4.00 Å are only −40, −24, −38, −29, and −14 microhartree, respectively. After 4000 MC iterations, they become −10, −10, −12, −9, and −2 microhartree, respectively. Once again, these are considerable improvements compared to CR-CC(2,3) and CC(t;3) that both give errors on the order of −0.2 millihartree,^52,53 especially if we realize that the fractions of triples captured by the i-FCIQMC runs after 2000 and 4000 MC iterations are relatively small (5%–28% and 5%–49%, respectively) and, as shown in Table S.2 of the supplementary material, the corresponding numbers of walkers represent only about 1%–2% of the total numbers of walkers at τ = 20.0 a.u., where we stopped our i-FCIQMC propagations.

As a consequence of the small errors in the CC(P;Q) total energies characterizing the $X1Σg+$ and $A3Σu+$ states in the early stages of the i-FCIQMC propagations, the resulting singlet–triplet gap values are very accurate as well. This is demonstrated in Table IV, where one can see that after 2000 δτ = 0.0001 a.u. MC iterations, which is, as already explained, a very short propagation time engaging tiny walker populations and small fractions of triples, most of the differences between the CC(P;Q) and CCSDT ΔE_S-T values in the R_H–F = 1.50–4.00 Å region are on the order of a few reciprocal centimeter. The only exception is the semi-stochastic CC(P;Q) run at R_H–F = 2.50 Å, where the −122 cm⁻¹ error relative to CCSDT characterizing the singlet–triplet gap obtained after 2000 MC time steps, while representing a five-fold error reduction compared to CR-CC(2,3), is comparable, in magnitude, to the CCSDT value of ΔE_S-T. This happens because the CC(P;Q) energy of the strongly correlated $X1Σg+$ state obtained after 2000 MC iterations at R_H–F = 2.50 Å differs from its CCSDT parent by −0.583 millihartree, whereas the analogous difference between the CC(P;Q) and CCSDT energies for its weakly correlated $A3Σu+$ companion is only −29 microhartree. This is not a problem though, since by running i-FCIQMC a little longer and capturing about 20% of all triply excited determinants in the relevant P spaces, as is the case when 4000 δτ = 0.0001 a.u. MC time steps are considered, one reduces the differences between the CC(P;Q) and CCSDT energies of the $X1Σg+$ and $A3Σu+$ states to −50 and −9 microhartree, respectively (cf. Tables II and III), so that the 122 cm⁻¹ unsigned error in the CC(P;Q) value of ΔE_S-T relative to CCSDT obtained after 2000 MC iterations decreases to less than 10 cm⁻¹. This is yet another illustration of the ability of the semi-stochastic CC(P;Q) methodology pursued in this work to balance the more MR singlet and weakly correlated triplet states of biradical systems in a single computation at the fraction of the cost of the parent high-level CC calculations. As shown in Table IV, at τ = 0.4 a.u., where the i-FCIQMC propagations are still far from being converged, the FCIQMC-driven CC(P;Q) calculations recover the CCSDT values of the singlet–triplet gaps in (HFH)⁻ at all H–F distances considered in this study to within a few reciprocal centimeter, reaching a 1–2 cm⁻¹ or better accuracy after 6000 MC iterations.

Last, but not least, the results reported in Tables II–IV and Fig. 2 also demonstrate the remarkable efficiency of the δ(P;Q) corrections in accelerating the convergence of the CC(P) energies of the $X1Σg+$ and $A3Σu+$ states and the vertical gaps between them toward CCSDT, independent of the H–F distance considered. Let us, for example, compare the uncorrected CC(P) and corrected CC(P;Q) energies of the $X1Σg+$ and $A3Σu+$ states of (HFH)⁻ at the five H–F separations considered in this work obtained after 2000 MC iterations. In the case of the former, more MR, state, the CC(P;Q) corrections δ(P;Q) reduce the positive 2.601, 3.998, 3.511, 6.586, and 0.412 millihartree errors relative to CCSDT resulting from the CC(P) computations at R_H–F = 1.50, 1.75, 2.00, 2.50, and 4.00 Å to the much smaller negative error values of −0.035, −0.056, −0.110, −0.583, and −0.025 millihartree, respectively. When the latter state, which is characterized by much weaker correlations, is considered, the CC(P;Q) approach reduces the 0.995, 0.826, 0.834, 0.502, and 0.239 millihartree errors obtained with CC(P) to −40, −24, −38, −29, and −14 microhartree, respectively. It is interesting to notice that while the errors characterizing the CC(P) calculations for the $A3Σu+$ state are generally much smaller than their $X1Σg+$ counterparts, and the two states have a substantially different character, the error reductions offered by the CC(P;Q) corrections δ(P;Q), by at least one order of magnitude, apply to both states. As already alluded to above, and as shown in Table IV and Figs. 2(e) and 2(f), where we examine the convergence of the CC(P) and CC(P;Q) ΔE_S-T values toward their CCSDT parents, the noniterative corrections δ(P;Q) are also very effective in improving the balance in the description of the $X1Σg+$ and $A3Σu+$ states by the CC(P) approach and smoothing the convergence of the resulting singlet–triplet gaps toward their CCSDT limits. This can be illustrated by comparing the behavior of the error values relative to CCSDT characterizing the CC(P) calculations of ΔE_S-T at R_H–F = 2.00 Å with their CC(P;Q) counterparts, shown in Table IV. In the former case, the 3462 cm⁻¹ error at τ = 0 decreases, in absolute value, to 8 cm⁻¹ at τ = 0.8 a.u. (8000 MC iterations), to increase to 17 cm⁻¹ at τ = 2.0 a.u. (20 000 MC iterations), to decrease again to a numerical 0 cm⁻¹ at τ = 20.0 a.u. (200 000 MC iterations). Once the CC(P) energies of the $X1Σg+$ and $A3Σu+$ states are corrected using the δ(P;Q) corrections, the unsigned errors in the resulting CC(P;Q) values of ΔE_S-T relative to their CCSDT parent monotonically and rapidly decrease, from 282 cm⁻¹ at τ = 0 to a numerical 0 cm⁻¹ at τ ≥ 0.8 a.u. It is clear from Tables II–IV and Fig. 2 that while both the CC(P) and CC(P;Q) energies converge to the parent CCSDT limit, the latter energies and the gaps between them converge to CCSDT a lot faster.

We now proceed to the examination of the performance of the semi-stochastic CC(P;Q) algorithm in calculations involving medium-sized organic biradicals, starting from two prototypical anti-aromatic systems, cyclobutadiene and cyclopentadienyl cation, both described using the cc-pVDZ basis set. As in the rest of the present study, we are mainly interested in how efficient the CIQMC-driven CC(P;Q) methodology is in recovering the CCSDT energies of the lowest singlet and triplet states and gaps between them. In the case of cyclobutadiene and cyclopentadienyl cation discussed in this subsection, we focus on examining vertical singlet–triplet gaps.

We begin with the FCIQMC-driven CC(P;Q) calculations for cyclobutadiene, in which we adopted the D_4h-symmetric geometry that represents the transition state for the automerization of cyclobutadiene proceeding on the lowest singlet potential, optimized with the MR average-quadratic CC (MR-AQCC) approach^165,166 using the cc-pVDZ basis in Ref. 167. We employed this geometry for two reasons. One of them is the fact that we used the same geometry in our earlier CIQMC- and CCMC-based,^95,97 CIPSI-driven,¹⁰² and active-orbital-based⁵¹ CC(P;Q) calculations for cyclobutadiene, when examining its automerization. Because of this, we could verify the correctness of our FCIQMC-driven CC(P;Q) calculations for the lowest-energy singlet state, which is also the ground state of the system. Another is the observation that the D_4h-symmetric transition-state structure characterizing the automerization of cyclobutadiene is practically identical to the D_4h-symmetric minimum on the lowest triplet surface. Indeed, the MR-AQCC/cc-pVDZ C–C and C–H bond lengths defining the transition state on the ground-state singlet potential differ from those characterizing the triplet minimum optimized using unrestricted CCSD (UCCSD) in Ref. 118 by less than 0.009 and 0.001 Å, respectively.

At the D_4h-symmetric geometry used in our calculations, cyclobutadiene is characterized by the delocalization of four π electrons over four π MOs, which gives rise to the close-lying singlet and triplet states that require a highly accurate treatment of electron correlation effects if we are to obtain a well-balanced description of the two states and the small energy separation between them. One can understand this by examining the valence π network of the D_4h-symmetric cyclobutadiene species, which consists of the doubly occupied nondegenerate a_2u orbital, the doubly degenerate e_g level, in which each component MO is occupied by a single electron, and the nondegenerate b_1u orbital, which in the zeroth-order description of the lowest singlet and triplet states remains empty. The two valence electrons in the degenerate e_g shell can couple to a singlet or a triplet, resulting in the open-shell singlet ground state, $X1B1g$⁠, which has a substantial MR character, and the first excited triplet state, $A3A2g$⁠, which is predominantly SR in nature. In order to balance the substantial nondynamical correlation effects, needed for an accurate description of the low-spin $X1B1g$ state, with the dynamical correlations dominating its high-spin triplet $A3A2g$ companion within a conventional, particle-conserving, SRCC framework and produce reliable ΔE_S-T values for cyclobutadiene, which could compete with the high-accuracy ab initio data reported in Refs. 118, 123, 124, and 167–172, one has to consider robust treatments of the connected triply excited clusters, such as that offered by CCSDT.^168,172 Indeed, full CCSDT, which is the target of this investigation, produces high-quality results for the lowest singlet and triplet states of the D_4h-symmetric cyclobutadiene system and the energy separation between them. For example, the ΔE_S-T value obtained in the CCSDT/cc-pVDZ calculations at the transition-state geometry used in the present study, of −4.8 kcal/mol, is practically identical to the results of the state-of-the-art DEA-EOMCC computations including the high-rank 4p2h correlations on top of CCSD, reported in Refs. 123 and 124, which give −5.0 kcal/mol when the cc-pVDZ basis set is employed (for similar recent observations regarding the reliability of full CCSDT in generating virtually exact singlet–triplet gap values for cyclobutadiene, see Ref. 172). It is, therefore, interesting to explore if the semi-stochastic CC(P;Q) methodology investigated in this work is capable of converging the CCSDT results for the $X1B1g$ and $A3A2g$ states of cyclobutadiene and vertical gap between them out of the early stages of CIQMC propagations.

The results of our FCIQMC-driven CC(P) and CC(P;Q) computations for cyclobutadiene are summarized in Table V and Fig. 3. In all of our calculations, starting with the stochastic i-FCIQMC steps and ending with the deterministic CC(P;Q) and CCSDT runs, we used the D_2h Abelian subgroup of the D_4h point group characterizing the cyclobutadiene’s geometry adopted in this work. Consequently, the i-FCIQMC propagations for the $X1B1g$ and $A3A2g$ states were set up to converge the lowest states of the¹A_g(D_2h) and³B_1g(D_2h) symmetries. Consistent with the CC(P) and CC(P;Q) runs that follow the i-FCIQMC steps and the accompanying CCSD, CR-CC(2,3), and CCSDT computations, the reference determinants used to initiate our i-FCIQMC propagations were the closed-shell, D_2h-adapted, RHF function obtained by placing two electrons on one of the e_g valence orbitals for the lowest-energy¹A_g(D_2h) state and the high-spin ROHF determinant, adapted to D_2h as well, for the lowest³B_1g(D_2h) state. As a result, the lists of triply excited determinants extracted from the i-FCIQMC runs at the various propagation times τ > 0, needed to define the P spaces for the CC(P) and CC(P;Q) computations, consisted of the S_z = 0 triples of the A_g(D_2h) symmetry for the $X1B1g$ state and the S_z = 1 triples of the B_1g(D_2h) symmetry in the case of the $A3A2g$ state. Given our interest in converging the CCSDT energetics, the Q spaces used to construct the δ(P;Q) corrections consisted of the remaining triply excited determinants, absent in the i-FCIQMC wave functions of the $X1B1g$ and $A3A2g$ states at a given τ.

The results shown in Table V and Fig. 3 display several similarities with the previously discussed methylene and (HFH)⁻ cases. One cannot, for example, obtain an accurate description of the more MR singlet ground state and the energy separation between the $X1B1g$ and $A3A2g$ states without incorporating the leading triply excited determinants in the P space. Indeed, when the P space consists of only singly and doubly excited determinants, as in the τ = 0 CC(P) (i.e., CCSD) and CC(P;Q) [i.e., CR-CC(2,3)] calculations, one ends up with the enormous errors in the energies of the $X1B1g$ state relative to their CCSDT parent, which are 47.979 millihartree in the former case and 14.636 millihartree when the latter computation is considered. The τ = 0 CC(P;Q) energy of the $A3A2g$ state is a lot more accurate, reducing the large, 23.884 millihartree, error relative to CCSDT obtained in the underlying CC(P) calculation to −60 microhartree, but this does not help too much. The corresponding CR-CC(2,3) triples correction to CCSD, which neglects the coupling of the low-order T₁ and T₂ clusters with their higher-order T₃ counterpart, is incapable of offering a balanced description of the $X1B1g$ and $A3A2g$ states, so that the resulting singlet–triplet gap is very poor. The 9.2 kcal/mol difference between the ΔE_S-T values obtained in the τ = 0 CC(P;Q) or CR-CC(2,3) and CCSDT calculations is so large that the $X1B1g−A3A2g$ separation predicted by CR-CC(2,3) has a wrong sign compared to its −4.8 kcal/mol CCSDT counterpart, while being nearly identical in magnitude. This difference becomes even larger when the uncorrected τ = 0 CC(P), meaning CCSD, calculations are considered (15.1 kcal/mol).

As shown in Table V and Fig. 3, the situation dramatically changes when the P spaces used in the CC(P) and CC(P;Q) calculations are enriched with the subsets of triply excited determinants captured by the i-FCIQMC propagations. The convergence of the CC(P;Q) energies of the $X1B1g$ and $A3A2g$ states, especially the former ones, and the vertical separations between them is particularly impressive. For example, after as few as 6000 δτ = 0.0001 a.u. MC time steps and i-FCIQMC capturing less than 30% of all triples in the P space, where, as demonstrated in Table S.3 of the supplementary material, the walker population characterizing the i-FCIQMC run for the $X1B1g$ state is only 0.02% of the total number of walkers at τ = 8.0 a.u. (the termination time for our i-FCIQMC propagations for cyclobutadiene), the CC(P;Q) approach reduces the 14.636 millihartree difference between the CR-CC(2,3) and CCSDT energies of the strongly correlated singlet ground state to 2.223 millihartree. While the CR-CC(2,3) description of the $A3A2g$ state, which has a largely SR character, is already excellent, the CC(P;Q) calculation performed after 6000 MC iterations, which uses only 26% of triples in the P space and a tiny walker population that amounts to 0.04% of all walkers at τ = 8.0 a.u. in the underlying i-FCIQMC propagation, improves it too, reducing the small, 60 microhartree, unsigned difference between the CR-CC(2,3) and CCSDT energies to an even smaller 51 microhartree. As a consequence of the above improvements, especially for the $X1B1g$ state, the error relative to CCSDT characterizing the ΔE_S-T value obtained in the FCIQMC-driven CC(P;Q) calculations after 6000 δτ = 0.0001 a.u. MC time steps, where the underlying i-FCIQMC propagations are still in their early stages, is only 1.4 kcal/mol, as opposed to 9.2 kcal/mol obtained at τ = 0 with CR-CC(2,3). The resulting $X1B1g−A3A2g$ energy separation, of −3.4 kcal/mol, has not only the correct sign, but is also very close to the −4.8 kcal/mol value obtained with CCSDT. If we wait a little longer, by executing the extra 2000 MC iterations, so that the i-FCIQMC propagations can capture 34%–39% of all triply excited determinants, we can reduce the already small 2.223 millihartree, 51 microhartree, and 1.4 kcal/mol errors in the CC(P;Q) energies of the $X1B1g$ and $A3A2g$ states and separation between them relative to CCSDT, obtained after 6000 MC time steps, to 0.835 millihartree, 31 microhartree, and 0.5 kcal/mol, respectively. It is clear from Table V and Fig. 3 that the convergence of the semi-stochastic CC(P;Q) results for the lowest-energy singlet and triplet states of cyclobutadiene, especially the $X1B1g$ energies and the $X1B1g−A3A2g$ gap values, which the τ = 0 CC(P;Q) or CR-CC(2,3) calculations describe poorly, toward CCSDT is very fast, even when the underlying i-FCIQMC propagations are far from convergence. It is also apparent from our calculations that the noniterative corrections δ(P;Q) play a significant role in accelerating convergence of the corresponding CC(P) energetics toward CCSDT. As shown, for example, in Table V, the relatively large differences between the uncorrected CC(P) energies of the $X1B1g$ and $A3A2g$ states and vertical gap between them obtained at τ = 0.8 a.u., i.e., after 8000 δτ = 0.0001 a.u. MC iterations, and the corresponding CCSDT data, which exceed 11 and 7 millihartree and 3 kcal/mol, respectively, are reduced to 0.835 millihartree, 31 microhartree, and 0.5 kcal/mol, when the CC(P;Q) approach is employed. We can see similar improvements in the CC(P) energies at other τ values.

Most of the observations regarding the performance of the semi-stochastic CC(P;Q) methodology and its CC(P) counterpart remain valid when the larger cyclopentadienyl cation, which is also the largest molecular system considered in our CC(P)/CC(P;Q) work to date, is examined. Following our previous DEA-EOMCC studies of cyclopentadienyl cation,^123,124 where we investigated the effect of high-order 4p2h correlations on its singlet–triplet gap, we used the D_5h-symmetric geometry corresponding to a minimum on the lowest triplet surface obtained in the UCCSD/cc-pVDZ optimization in Ref. 118. At this geometry, cyclopentadienyl cation is characterized by the delocalization of four π electrons over five π MOs, resulting in the doubly occupied nondegenerate $a2″$ orbital, the doubly degenerate $e1″$ shell, in which each component MO is occupied by a single electron, and the doubly degenerate $e2″$ shell, which in the zeroth-order description of the lowest-energy singlet and triplet states remains empty. In analogy to the previously discussed cyclobutadiene system, the two electrons in the degenerate $e1″$ MOs can couple to a singlet or triplet, but compared to cyclobutadiene, where the lowest-energy singlet state is also a ground state, the state ordering in cyclopentadienyl cation is reversed, so that the lowest triplet, designated as $X3A2′$⁠, is the ground state and the lowest-energy singlet, denoted as $A1E2′$⁠, is the first excited state. Similar to all other examples considered in this work, in order to obtain a well-balanced description of the $X3A2′$ state, which has a SR character dominated by dynamical correlations, and its $A1E2′$ companion, which is an open-shell singlet characterized by significant MR correlations, and obtain an accurate value of ΔE_S-T within a conventional SRCC framework, one must turn to higher-level theories that can offer a robust treatment of T_n clusters with n > 2. Otherwise, as shown in Ref. 118, and as confirmed in our calculations, the results can be very poor. For example, the $A1E2′−X3A2′$ separation in cyclopentadienyl cation resulting from the restricted CCSD calculations using the cc-pVDZ basis, which are equivalent to our τ = 0 CC(P) computations, is about 23 kcal/mol. This is in large disagreement with the most accurate ab initio calculations of the singlet–triplet gap in cyclopentadienyl cation performed to date using the DEA-EOMCC formalism including 3p1h as well as 4p2h correlations on top of the CCSD treatment of the underlying closed-shell core, which give about 14 kcal/mol when the cc-pVDZ basis set is employed^123,124 (for the examples of other high-level SRCC and MRCC calculations of the singlet–triplet gap in cyclopentadienyl cation, see Ref. 118; Ref. 124 also provides the well-converged MR perturbation theory data, which agree with the state-of-the-art DEA-EOMCC computations reported in Refs. 123 and 124). The restricted CCSDT approach, which is the target SRCC method in this study, provides a much better description, reducing the ∼9 kcal/mol error relative to the most accurate DEA-EOMCC calculations with up to 4p2h excitations reported in Refs. 123 and 124 obtained with restricted CCSD to less than 3 kcal/mol, when the cc-pVDZ basis set is employed. It would certainly be interesting to examine if the inclusion of higher-than-triply excited clusters, such as T₄, could further improve the CCSDT description of the singlet–triplet gap in cyclopentadienyl cation, but in this work we focus on the ability of the semi-stochastic, CIQMC-based, CC(P;Q) methodology to improve the CR-CC(2,3) ΔE_S-T values and converge the results of CCSDT computations. We hope to return to the topic of the role of T₄ clusters in describing the singlet–triplet gap in cyclopentadienyl cation in one of our future studies. It may be worth pointing out that the $A1E2′−X3A2′$ gap obtained in the restricted CCSDT/cc-pVDZ calculations, which give ΔE_S-T = 16.7 kcal/mol, is in very good agreement with the 16.1 kcal/mol resulting from the DEA-EOMCC/cc-pVDZ computations truncated at 3p1h excitations.^123,124

The results of our CIQMC-driven CC(P) and CC(P;Q) computations for cyclopentadienyl cation are reported in Table VI and Fig. 4. As already alluded to above, to reduce the computational costs of the CIQMC propagations preceding the CC(P) and CC(P;Q) steps, especially in the later stages of the CIQMC runs that are included in Table VI and Fig. 4 for the completeness of our presentation, we replaced the i-FCIQMC algorithm, which we exploited in our calculations for methylene, (HFH)⁻, and cyclobutadiene, by its truncated i-CISDTQ-MC counterpart. It has been established in Ref. 97 that the replacement of i-FCIQMC by i-CISDTQ-MC, when identifying the leading higher-than-doubly excited determinants for the inclusion in the P spaces used in the semi-stochastic CC(P) and CC(P;Q) runs, has virtually no effect on the rate at which these runs converge the parent SRCC energetics. In analogy to cyclobutadiene, all of our i-CISDTQ-MC, semi-stochastic CC(P) and CC(P;Q), and deterministic CCSD, CR-CC(2,3), and CCSDT computations utilized the largest Abelian subgroup of the D_5h point group characterizing the cyclopentadienyl cation’s structure examined in the present study, which is C_2v. This means that in setting up our calculations for the $X3A2′$ state, we treated it as the lowest state of the³B₂(C_2v) symmetry, whereas the doubly degenerate $A1E2′$ state was represented by its¹A₁(C_2v) component. Similar to cyclobutadiene, and to remain consistent with the CC(P), CC(P;Q), and other SRCC runs for cyclopentadienyl cation carried out in this study, the reference determinant used to initiate the i-CISDTQ-MC propagation for the lowest-energy³B₂(C_2v) state was the triplet ROHF determinant. In the case of the¹A₁(C_2v) component of the $A1E2′$ state, we used the RHF determinant obtained by pairing the two valence electrons in one of the $e1″$ MOs to initiate the corresponding i-CISDTQ-MC run. Consistent with the above description, the subsets of triply excited determinants used to construct the P spaces for the semi-stochastic CC(P) and CC(P;Q) computations for the $X3A2′$ state were the S_z = 1 triples of the B₂(C_2v) symmetry captured by i-CISDTQ-MC. In the case of the $A1E2′$ state, represented, as explained above, by its¹A₁(C_2v) component, we used the S_z = 0 triples of the A₁(C_2v) symmetry identified by the i-CISDTQ-MC propagation set up to converge the lowest A₁(C_2v) state. As usual, the corresponding Q spaces were spanned by the remaining triply excited determinants that were not captured by the i-CISDTQ-MC runs when the lists of P-space triples were created.