Benchmarking the Semi-Stochastic CC(P;Q) Approach for Singlet-Triplet Gaps in Biradicals

We recently proposed a semi-stochastic approach to converging high-level coupled-cluster (CC) energetics, such as those obtained in the CC calculations with singles, doubles, and triples (CCSDT), in which the deterministic CC($P$;$Q$) framework is merged with the stochastic configuration interaction Quantum Monte Carlo propagations [J. E. Deustua, J. Shen, and P. Piecuch, Phys. Rev. Lett. 119, 223003 (2017)]. In this work, we investigate the ability of the semi-stochastic CC($P$;$Q$) methodology to recover the CCSDT energies of the lowest singlet and triplet states and the corresponding singlet-triplet gaps of biradical systems using methylene, ${\rm (HFH)}^{-}$, cyclobutadiene, cyclopentadienyl cation, and trimethylenemethane as examples.


I. INTRODUCTION
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 |Ψ = e T |Φ , where T = N n=1 T n 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][17][18][19] where T is truncated at T 2 , the CC method with singles, doubles, and triples (CCSDT), [20][21][22][23] where T is truncated at T 3 , the CC method with singles, doubles, triples, and quadruples (CCSDTQ), [24][25][26][27] where T is truncated at T 4 , 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 equationof-motion (EOM) [28][29][30][31][32][33][34][35] and linear response (LR) [36][37][38][39][40][41][42][43][44] CC frameworks rapidly converge to the exact, full configuration interaction (FCI), limit. 15 As a result, the CCSDT, CCSDTQ, and similar methods and their EOM and LR a) Corresponding author; e-mail: piecuch@chemistry.msu.edu. 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 n 3 o n 5 u steps of CCSDT or the even more demanding n 4 o n 6 u 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-thantwo-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][48][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][55][56][57][58] in which one corrects the CCSD energies for the leading T 3 correlations, and its higher-order extensions. 53,[59][60][61][62] The CR-CC(2,3) method improves poor performance of CCSD(T) in situations involving covalent bond breaking 50,[54][55][56][57]60,63,64 and certain types of noncovalent interactions, 59,65 while retaining the relatively inexpensive iterative n 2 o n 4 u and noniterative n 3 o n 4 u steps of CCSD(T), but neither CR-CC(2,3) nor its CR-CCSD(T) 47,48,66,67 and locally renormalized CCSD(T) 68 predecessors, nor their CCSD(2) T , 69-72 CCSD(T) Λ , [73][74][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 1 and T 2 , and their higher-order counterparts with n > 2, e.g., T 3 , becomes larger. 51,52 One can address this issue by solving for T 1 and T 2 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][81][82][83][84][85] and account for the remaining correlation effects of interest (e.g., the T 3 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][51][52][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][87][88][89][90] and CC Monte Carlo (CCMC) [91][92][93][94] propagations. [95][96][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-thandoubly 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 3 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 97 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 96 (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][100][101] obtaining encouraging results. 102 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, 103 molecular magnetism, 104 spintronics, 105,106 nonlinear optics, 107 photochemical pathways, [108][109][110] and photovoltaics. [111][112][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 − . 114 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 highlevel ab initio wave function methods are used in the calculations. 52,57,[115][116][117][118][119][120][121][122][123][124][125][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) 2 → (LUMO) 2 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 singlettriplet 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][126][127][128][129][130][131][132][133][134][135][136] especially those that incorporate the high-rank 4-particle-2-hole (4p2h) and 4-hole-2particle (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, CIQMCdriven, 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. 96

II. THEORY
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][55][56][57][58] CR-EOMCC(2,3), 56,137 and δ-CR-EOMCC(2,3), 138 and their extensions to quadruples, 53,[59][60][61][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 manyelectron 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 similaritytransformed Hamiltonian H (P ) = e −T (P ) He T (P ) in the P space, with Λ (P ) representing the relevant hole-particle deexcitation operator, can be found in Refs. [50][51][52][53]also,Refs. 59,65,[95][96][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 1 and T 2 components of the cluster operator T in the presence of their higher-order counterparts, such as the leading T 3 contributions, which the CR-CC(2,3), CCSD(2) T , Λ-CCSD(T), and other noniterative triples corrections to CCSD cannot do. As explained in the Introduction, one can incorporate the coupling among the lower-rank T 1 and T 2 and higher-rank T n components with n > 2 within the CC(P ;Q) frame-work 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: 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 , and its Λ (P ) = Λ 1 + Λ 2 + Λ 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 3 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 manyelectron 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 higherthan-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][96][97] and related work, 98 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][96][97]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 3 correlations not captured by the preceding CC(P ) runs, which scale no worse than ∼2n 3 o n 4 u , are much less than the n 3 o n 5 u 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 3 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 Section 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 3 = 0. The main role of corrections δ(P ; Q) in the context of the semistochastic CC(P ;Q) algorithm considered in this study is to accelerate convergence of the underlying CC(P ) energies toward the desired CCSDT limit.

III. RESULTS
As explained in the Introduction, in order to assess the performance of our semi-stochastic, CIQMCdriven, 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 singlettriplet gaps in the same systems using the deterministic CC(P ;Q) 52 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 140 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 transfor-mation 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.

A. Methylene
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 3 B 1 ) and first-excited (A 1 A 1 ) states of the methylene/aug-cc-pVTZ system and the adiabatic gap between them. The C 2v -symmetric geometries of CH 2 in the two states, optimized using FCI and the [5s3p/3s] triple zeta basis set of Dunning 148 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 (1a 1 ) 2 (2a 1 ) 2 (1b 2 ) 2 (3a 1 ) 1 (1b 1 ) 1 configuration, whereas the first-excited singlet state exhibits a significant MR character requiring a linear combination of the (1a 1 ) 2 (2a 1 ) 2 (1b 2 ) 2 (3a 1 ) 2 and doubly excited (1a 1 ) 2 (2a 1 ) 2 (1b 2 ) 2 (1b 1 ) 2 closed-shell determinants for a proper zeroth-order description. Because of these fundamentally different characteristics of the X 3 B 1 and A 1 A 1 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][150][151][152][153][154] or MRCC [155][156][157][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 3 B 1 and A 1 A 1 states of methylene, see Refs. 123, 135, 136, and 159 and references therein). The CCSDT results for the adiabatic singlettriplet 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 149 is only 0.11 kcal/mol or 38 cm −1 . As demonstrated in Ref. 52 as well, the purely electronic A 1 A 1 − X 3 B 1 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 −1 ) 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 semistochastic CC(P ;Q) approach advocated in this work is capable of reproducing the high-quality CCSDT/aug-cc-pVTZ data for the X 3 B 1 and A 1 A 1 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 3 B 1 state and the RHF determinant for the A 1 A 1 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 1 symmetry captured during the i-FCIQMC run for the X 3 B 1 state and the S z = 0 triples of the A 1 symmetry captured during the analogous run for the A 1 A 1 state. Following the semi-stochastic CC(P ;Q) algorithm described in Section 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 3 B 1 and A 1 A 1 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 −1 or ∼11% error in the adiabatic ∆E S-T value when compared to the 3328 cm −1 singlet-triplet gap obtained with CCSDT. The situation improves when the CC(P ;Q) corrections δ(P ; Q) due to T 3 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 3 B 1 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 1 A 1 state as well, although the difference between the resulting CR-CC(2,3) energy and its CCSDT counterpart, of 0.656 millihartree, is almost 4 times larger than the analogous difference obtained for the X 3 B 1 state. As a result, the 105 cm −1 error relative to CCSDT characterizing the adiabatic A 1 A 1 − X 3 B 1 separation obtained in the τ = 0 CC(P ;Q) or CR-CC(2,3) calculations, while considerably smaller than the 380 cm −1 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 3 B 1 and A 1 A 1 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 3 correlations that have not been captured by CCSDt, 52 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 τ = 0 P 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 3 correlations, the convergence of the total electronic energies of the X 3 B 1 and A 1 A 1 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 3 B 1 and A 1 A 1 states and the corresponding ∆E S-T value relative to CCSDT are 0.049 millihartree, 0.106 millihartree, and 13 cm −1 , respectively, substantially improving the CR-CC(2,3) [i.e., τ = 0 CC(P ;Q)] calculations, which give 0.177 millihartree for the X 3 B 1 state, 0.656 millihartree for the A 1 A 1 state, and 105 cm −1 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 1 A 1 state, which is characterized by larger T 3 effects, is the use of the unrelaxed T 1 and T 2 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 3 B 1 state and only 25% of all triples in the P space for the A 1 A 1 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 −1 errors relative to CCSDT for the X 3 B 1 and A 1 A 1 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 132689 in the case of the X 3 B 1 state and 165564 for the A 1 A 1 state, represent tiny fractions, 2.17% and 1.11%, respec-tively, 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 3 B 1 and A 1 A 1 states and the gap between them resulting from the FCIQMCdriven 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 20000 δτ = 0.0001 a.u. MC iterations, the CC(P ;Q) energies of the X 3 B 1 and A 1 A 1 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 −1 , respectively. As shown in Table  S.1 of the supplementary material, walker populations characterizing the X 3 B 1 and A 1 A 1 states produced by i-FCIQMC at 20000 δτ = 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 3 B 1 state and 2.17% in the case of the A 1 A 1 state). It is also interesting to note that the more MR A 1 A 1 state requires a higher fraction of triply excited determinants in the P space than its SR X 3 B 1 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 1 A 1 state to ∼0.1 millihartree. In the case of the X 3 B 1 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 semistochastic 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 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 ) en-ergies of the X 3 B 1 and A 1 A 1 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 −1 error obtained with CC(P ) to 13 cm −1 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 −1 at 50000 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 50000 MC iterations to a numerical 0 cm −1 . 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.

B. (HFH) −
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 − . 114 The energies of the lowest two electronic states of the (HFH) − system, including the singlet ground state X 1 Σ + g and the first-excited triplet state A 3 Σ + 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), 53 the DIP-EOMCC approaches with full and active-space treatments of 4h2p correlations of top of CCSD, 135,136 and the active-orbitalbased 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 A 3 Σ + u state, which is weakly correlated and well represented by a single ROHF determinant, its groundstate counterpart X 1 Σ + g displays a substantial MR character that includes a significant contribution from the doubly excited (HOMO) 2 → (LUMO) 2 determinant, in addition to the RHF configuration. The MR character of the X 1 Σ + 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) 2 → (LUMO) 2 and RHF determinants or the equivalent T 2 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 1 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 X 1 Σ + g state and the ROHF reference for the A 3 Σ + 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 X 1 Σ + g state and only 2.628 millihartree when the A 3 Σ + u state is considered. The analogous differences at R H-F = 2.00 A are 19.398 and 2.068 millihartree, respectively. The observed large discrepancies between the errors in the CCSD energies for the X 1 Σ + g and A 3 Σ + 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 A with the corresponding FCI data. CCSD/6-31G(d,p) gives −7320, −1838, 1656, 3605, and 230 cm −1 , respectively, as opposed to −9525, −4911, −2147, −277, and 0 cm −1 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, 53 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 −1 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 −1 , respectively. 52,53,163 It also greatly improves the total electronic energies. Indeed, the differences between the CCSDT and FCI energies of the X 1 Σ + g and A 3 Σ + 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 53 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 X 1 Σ + g and A 3 Σ + 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 X 1 Σ + g and A 3 Σ + u states were set up to converge the lowest-energy states of the 1 A g (D 2h ) and 3 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 X 1 Σ + 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 A 3 Σ + 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 52 and the above remarks, the CC(P ) energies of the X 1 Σ + 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 X 1 Σ + 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 A 3 Σ + u state (see Table III). As already alluded to above, and as shown in Table IV, this imbalance in the description of the X 1 Σ + g and A 3 Σ + 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 cm −1 to 3462 cm −1 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 2 and F − , resulting in the nearly degenerate singlet and triplet states and the 172 cm −1 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 3 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 X 1 Σ + g and A 3 Σ + 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 X 1 Σ + 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 X 1 Σ + 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 A 3 Σ + 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 −1 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 −1 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 3 correlations are added to the energies obtained in the CC(P ) calculations, the resulting CC(P ;Q) values of the X 1 Σ + g and A 3 Σ + 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 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 A, remain relatively small, the differences between the CC(P ;Q) and CCSDT energies obtained for the X 1 Σ + 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 A 3 Σ + 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 A 3 Σ + 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 X 1 Σ  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 −1 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 X 1 Σ + 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 A 3 Σ + 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 X 1 Σ + g and A 3 Σ + u states to −50 and −9 microhartree, respectively (cf. Tables II and III), so that the 122 cm −1 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 −1 . 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 −1 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 X 1 Σ + g and A 3 Σ + 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 X 1 Σ + g and A 3 Σ state are generally much smaller than their X 1 Σ + 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  Once the CC(P ) energies of the X 1 Σ + g and A 3 Σ + 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 −1 at τ = 0 to a numerical 0 cm −1 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.

C. Cyclobutadiene and cyclopentadienyl cation
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 4hsymmetric 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, 102 and active-orbitalbased 51 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, X 1 B 1g , which has a substantial MR character, and the first excited triplet state, A 3 A 2g , which is predominantly SR in nature. In order to balance the substantial nondynamical correlation effects, needed for an accurate description of the lowspin X 1 B 1g state, with the dynamical correlations dominating its high-spin triplet A 3 A 2g 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,[167][168][169][170][171][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 singlettriplet 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 X 1 B 1g and A 3 A 2g 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 X 1 B 1g and A 3 A 2g states were set up to converge the lowest states of the 1 A g (D 2h ) and 3 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 1 A g (D 2h ) state and the high-spin ROHF determinant, adapted to D 2h as well, for the lowest 3 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 X 1 B 1g state and the S z = 1 triples of the B 1g (D 2h ) symmetry in the case of the A 3 A 2g 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 X 1 B 1g and A 3 A 2g 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 X 1 B 1g and A 3 A 2g 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 X 1 B 1g 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 A 3 A 2g 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 1 and T 2 clusters with their higher-order T 3 counterpart, is incapable of offering a balanced description of the X 1 B 1g and A 3 A 2g 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 X 1 B 1g − A 3 A 2g 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). 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 X 1 B 1g and A 3 A 2g 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 X 1 B 1g 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 A 3 A 2g 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 shown in
As a consequence of the above improvements, especially for the X 1 B 1g state, the error relative to CCSDT characterizing the ∆E S-T value obtained in the FCIQMCdriven 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 X 1 B 1g − A 3 A 2g 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 X 1 B 1g and A 3 A 2g 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 semistochastic CC(P ;Q) results for the lowest-energy singlet and triplet states of cyclobutadiene, especially the X 1 B 1g energies and the X 1 B 1g − A 3 A 2g 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 X 1 B 1g and A 3 A 2g 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 a ′′ 2 orbital, the doubly degenerate e ′′ 1 shell, in which each component MO is occupied by a single electron, and the doubly degenerate e ′′ 2 shell, which in the zeroth-order description of the lowest-energy singlet and triplet states remains empty. In analogy to the previously discussed cy-clobutadiene system, the two electrons in the degenerate e ′′ 1 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 X 3 A ′ 2 , is the ground state and the lowest-energy singlet, denoted as A 1 E ′ 2 , is the first excited state. Similar to all other examples considered in this work, in order to obtain a well-balanced description of the X 3 A ′ 2 state, which has a SR character dominated by dynamical correlations, and its A 1 E ′ 2 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 A 1 E ′ 2 − X 3 A ′ 2 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 The restricted CCSDT approach, which is the target SRCC method in this study, provides a much better description, reducing the approximately 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 higherthan-triply excited clusters, such as T 4 , 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 4 clusters in describing the singlet-triplet gap in cyclopentadienyl cation in one of our future studies. It may be worth pointing out that the A 1 E ′ 2 − X 3 A ′ 2 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 X 3 A ′ 2 state, we treated it as the lowest state of the 3 B 2 (C 2v ) symmetry, whereas the doubly degenerate A 1 E ′ 2 state was represented by its 1 A 1 (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 3 B 2 (C 2v ) state was the triplet ROHF determinant. In the case of the 1 A 1 (C 2v ) component of the A 1 E ′ 2 state, we used the RHF determinant obtained by pairing the two valence electrons in one of the e ′′ 1 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 X 3 A ′ 2 state were the S z = 1 triples of the B 2 (C 2v ) symmetry captured by i-CISDTQ-MC. In the case of the A 1 E ′ 2 state, represented, as explained above, by its 1 A 1 (C 2v ) component, we used the S z = 0 triples of the A 1 (C 2v ) symmetry identified by the i-CISDTQ-MC propagation set up to converge the lowest A 1 (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.
Our calculations for cyclopentadienyl cation, summarized in Table VI and Fig. 4, demonstrate that the CC(P ;Q) energies of the X 3 A ′ 2 and A 1 E ′ 2 states and vertical gaps between them display fast convergence toward the respective CCSDT values with the propagation time τ . This is particularly apparent in the case of the CC(P ;Q) energies of the more MR A 1 E ′ 2 state and the A 1 E ′ 2 − X 3 A ′ 2 separation, which cannot be accurately described if the underlying P spaces contain only singly and doubly excited determinants. Indeed, the CR-CC(2,3) energy of the A 1 E ′ 2 state, which is equivalent to the τ = 0 CC(P ;Q) value, is much more accurate than the result of the associated CC(P ) or CCSD calculation, which produces the enormous error relative to CCSDT exceeding 38 millihartree, but the substantial, > 6 millihartree, difference with the CCSDT energy remains. The situation for the SR X 3 A ′ 2 state, where the CR-CC(2,3) approach reduces the nearly 29 millihartree error relative to CCSDT obtained in the CCSD calculations to ∼0.2 millihartree, is a lot better, but this does not help the resulting ∆E S-T value, which differs from its CCSDT counterpart by almost 4 kcal/mol (almost a quarter of the CCSDT value of ∆E S-T ). The discrepancy between the errors in the CR-CC(2,3) energies of the X 3 A ′ 2 and A 1 E ′ 2 states is simply too large. Clearly, one needs to incorporate some triples in the corresponding P spaces, especially in the case of the more challenging A 1 E ′ 2 state. Once the τ = 0 P spaces are augmented with the leading triply excited determinants identified by the i-CISDTQ-MC propagations and the noniterative corrections δ(P ; Q) are added to the CC(P ) energies to estimate the effects of the remaining T 3 correlations, we observe smooth convergence of the resulting CC(P ;Q) energetics toward their respective CCSDT limits. This includes significant improvements in the poor description of the A 1 E ′ 2 state and the A 1 E ′ 2 − X 3 A ′ 2 separation by CR-CC(2,3). As shown in Table VI Table S.4 of the supplementary material), the 6.245 millihartree and 3.8 kcal/mol errors in the energy of the A 1 E ′ 2 state and the ∆E S-T value relative to CCSDT obtained with CR-CC(2,3) reduce in the CC(P ;Q) calculations to 2.248 millihartree and 1.3 kcal/mol, respectively. By running i-CISDTQ-MC a little longer and capturing about 50-60% of all triples in the relevant P spaces, as is the case after 20000 MC iterations, where the walker populations compared to τ = 8.0 a.u. are still tiny, the errors in the CC(P ;Q) values of the A 1 E ′ 2 energy and ∆E S-T relative to their CCSDT parents drop down by an order of magnitude compared to 10000 MC iterations, to 0.217 millihartree and 0.1 kcal/mol, respectively. Although the excellent description of the predominantly SR X 3 A ′ 2 state by the CR-CC(2,3) approach hardly needs any improvement, the i-CISDTQ-MC-driven CC(P ;Q) calculations are helping here too, reducing the 0.245 millihartree difference between the CR-CC(2,3) and CCSDT energies to 0.108 millihartree after 10000 MC iterations (26 microhartree when the number of MC iterations is increased to 20000). As anticipated, the uncorrected CC(P ) energies of the X 3 A ′ 2 and A 1 E ′ 2 states converge to the respective CCSDT limits too, but they do it at a much slower pace than their CC(P ;Q) counterparts. A comparison of the results of the CC(P ) and CC(P ;Q) calculations for the A 1 E ′ 2 − X 3 A ′ 2 gap shown in Table VI and Fig. 4 (c) may create an impression as if the noniterative corrections δ(P ; Q) offer very little, but this would be misleading. The relatively fast convergence of the CC(P ) values of ∆E S-T toward their CCSDT parent in the early stages of the underlying i-CISDTQ-MC propagations, which compares well with that observed in the corresponding CC(P ;Q) computations, is a result of the fortuitous cancellation of large errors characterizing the CC(P ) energies of the X 3 A ′ 2 and A 1 E ′ 2 states. Since no other system examined in this study displays similar error cancellations, and since costs of computing corrections δ(P ; Q), which offer major error reductions in the individual CC(P ) energies, while accelerating their convergence toward the SRCC target, are low, we recommend using the δ(P ; Q)-corrected CC(P ;Q) energetics.

D. Trimethylenemethane
Our final example is trimethylenemethane, a fascinating non-Kekulé hydrocarbon examined as early as in 1948 173 and 1950, 174 in which four valence π electrons are delocalized over four closely spaced π-type orbitals. Assuming the D 3h symmetry, which is the symmetry of the minimum-energy structure on the ground-state triplet surface of trimethylenemethane, the four MOs of this system's valence π network consist of the nondegenerate 1a ′′ 2 orbital, the doubly degenerate 1e ′′ shell, and the nondegenerate 2a ′′ 2 orbital. If one adopts the C 2v symmetry, relevant to the low-lying singlet states, which is also the largest Abelian subgroup of D 3h exploited in our CCSD, CR-CC(2,3), CCSDT, and CIQMC-driven CC(P ) and CC(P ;Q) computations, the nondegenerate 1a ′′ 2 and 2a ′′ 2 orbitals in a D 3h description become the 1b 1 and 3b 1 MOs, respectively, whereas the degenerate 1e ′′ shell splits into the 1a 2 and 2b 1 components. The first experimental identification of trimethylenemethane dates back to 1966, 175 a definitive experimental verification, using electron paramagnetic resonance, of its triplet ground state was accomplished already in 1976, 176 and the electronic structure of trimethylenemethane has been well understood for decades (cf., e.g., Ref. 177 and references therein), but an accurate characterization of its triplet ground state and lowlying singlet states and energy separations between them continues to present a significant challenge to quantum chemistry approaches. 52,115,116,123,135,136, The D 3h -symmetric triplet ground state, designated as X 3 A ′ 2 (in a C 2v description adopted in this study, X 3 B 2 ), which is dominated by the |{core}(1a ′′ 2 ) 2 (1e ′′ 1 ) 1 (1e ′′ 2 ) 1 | configuration (in C 2v , |{core}(1b 1 ) 2 (1a 2 ) 1 (2b 1 ) 1 |), is relatively easy to describe, but the next two states, which are the nearly degenerate singlets stabilized by the Jahn-Teller distortion that lifts their exact degeneracy in a D 3h description, are not. The lower of the two singlets, which is characterized by a C s -symmetric minimum that can be approximated by a twisted C 2v structure and which is, therefore, usually designated as the A 1 B 1 state, is an open-shell singlet that emerges from the |{core}(1b 1 ) 2 (1a 2 ) 1 (2b 1 ) 1 | configuration. The second singlet, labeled as the B 1 A 1 state, is a C 2v -symmetric multi-configurational state dominated by the |{core}(1b 1 ) 2 (1a 2 ) 2 | and |{core}(1b 1 ) 2 (2b 1 ) 2 | closed-shell determinants. The A 1 B 1 state, although lower in energy compared to its B 1 A 1 counterpart, has not been observed experimentally due to unfavorable Frank-Condon factors, 193,202 so we do not consider it in this work. However, the second singlet, B 1 A 1 , has been detected in photoelectron spectroscopy experiments reported in Refs. 202 and 203, which located it at 16.1 ± 0.1 kcal/mol above the X 3 A ′ 2 ground state. Thus, following our previous deterministic, active-orbital-based, CC(P ;Q) work 52 and the state-of-the-art DEA-and DIP-EOMCC computations with up to 4p2h and 4h2p excitations reported in Refs. 123, 135, and 136, in carrying out the CIQMC-driven CC(P ) and CC(P ;Q) calculations discussed in this subsection and executing the accompanying CCSD, CR-CC(2,3), and CCSDT runs, we focused on the D 3h -symmetric triplet ground state, X 3 A ′ 2 , the C 2v -symmetric B 1 A 1 singlet, and the adiabatic gap between them, adopting the geometries of the two states optimized using the spin-flip density functional theory (SF-DFT) and the 6-31G(d) basis in Ref. 115. In analogy to other organic biradicals discussed in this article, we employed the cc-pVDZ basis set, so that the parent CCSDT computations, needed to judge the performance of our semi-stochastic CC(P ) and CC(P ;Q) methods, and the more expensive CC(P ) and CC(P ;Q) calculations employing large, near-100%, fractions of triples in the relevant P spaces (captured in the later stages of the underlying CIQMC propagations) were not too difficult to execute on the computers available to us. As shown in our earlier deterministic CC(P ;Q) work, 52 in which we tested the active-orbital-based CC(t;3) method, which recovers the CCSDT energetics to within small fractions of kilocalorie per mole, and as confirmed by the authors of Ref. 199, who managed to perform the CCSDT/cc-pVTZ calculations, the use of a larger cc-pVTZ basis changes the adiabatic B 1 A 1 − X 3 A ′ 2 gap by about 0.5-1 kcal/mol, i.e., the use the cc-pVDZ basis is sufficient to draw meaningful conclusions regarding the performance of the semi-stochastic CC(P ) and CC(P ;Q) approaches.
While the main goal of this study is to examine the efficiency of the CIQMC-driven CC(P ;Q) approaches in converging the CCSDT energetics, it is worth pointing out that the parent CCSDT calculations using the ROHF reference determinant for the X 3 A ′ 2 state and the RHF reference for the more strongly correlated B 1 A 1 state, in spite of their SR character, are capable of producing a reasonably accurate description of the adiabatic B 1 A 1 − X 3 A ′ 2 separation in trimethylenemethane. Indeed, the purely electronic B 1 A 1 − X 3 A ′ 2 gap, designated, in analogy to other singlet-triplet gaps considered in this work, as ∆E S-T , resulting from the ROHF/RHFbased CCSDT/cc-pVDZ computations using the SF-DFT/6-31G(d) geometries of the X 3 A  Table VII as the τ = 0 CC(P ) and CC(P ;Q) data, respectively]. While much of the 3.6 kcal/mol difference between the electronic B 1 A 1 −X 3 A  199, the CCSDT value of the adiabatic B 1 A 1 − X 3 A ′ 2 gap increases with the basis set too, to 23.1 kcal/mol when the cc-pVTZ basis is employed, which reinforces our view that without accounting for T 4 correlations one cannot bring the results of conventional SRCC computations to a close agreement with the experimentally derived data. While the examination of the role of T 4 clusters, basis set, geometries of the X 3 A ′ 2 and B 1 A 1 states employed in the calculations, ∆ZPVE corrections, etc. would certainly be interesting, it would also be outside the scope of the present study. Thus, in the remainder of this subsection, we return to the analysis of the performance of the CIQMC-driven CC(P ;Q) approach and its CC(P ) counterpart, especially their ability to converge the parent CCSDT energetics when the cc-pVDZ basis is employed.
The results of our semi-stochastic CC(P )/cc-pVDZ and CC(P ;Q)/cc-pVDZ computations for the X 3 A ′ 2 and B 1 A 1 states of trimethylenemethane and the adiabatic gap between them, along with the associated CCSD, CR-CC(2,3), and CCSDT data, are summarized in Table VII and Fig. 5. As in the case of cyclopentadienyl cation, to reduce the computational costs of the underlying CIQMC propagations, especially in their later stages, we resorted to the truncated i-CISDTQ-MC approach. In analogy to cyclobutadiene and cyclopentadienyl cation, we terminated our i-CISDTQ-MC propagations after 80000 δτ = 0.0001 a.u. MC time steps, where the differences between the CC(P ;Q) and CCSDT energies of the X 3 A ′ 2 and B 1 A 1 states fall below 1 microhartree. Consistent with the CC(P ), CC(P ;Q), and other SRCC calculations for trimethylenemethane reported in Table VII and Fig. 5, we used the ROHF determinant to initiate the i-CISDTQ-MC propagation for the D 3h -symmetric X 3 A ′ 2 (in C 2v , X 3 B 2 ) state and the RHF determinant to initiate the i-CISDTQ-MC run for the C 2v -symmetric B 1 A 1 state. The lists of triply excited determinants captured by the i-CISDTQ-MC runs at the various times τ > 0, needed to construct the P spaces for the CC(P ) and CC(P ;Q) computations, were the S z = 1 triples of the B 2 (C 2v ) symmetry in the case of the X 3 A ′ 2 state and the S z = 0 triples of the A 1 (C 2v ) symmetry when considering the B 1 A 1 state. The remaining triples not captured by i-CISDTQ-MC defined the corresponding Q spaces.
It is clear from the results presented in Table VII and Fig. 5 that the semi-stochastic CC(P ;Q) approach is very effective in converging the parent CCSDT energetics characterizing the X 3 A ′ 2 and B 1 A 1 states of trimethylenemethane and the adiabatic gap between them. It offers substantial improvements in the results of the CR-CC(2,3) calculations in the early stages of the underlying i-CISDTQ-MC propagations, especially when the multiconfigurational B 1 A 1 state and the adiabatic B 1 A 1 − X 3 A ′ 2 separation ∆E S-T , which are poorly described by CR-CC(2,3), are examined, while greatly accelerating the convergence of the CC(P ) energies toward CCSDT. Indeed, after 6000 δτ = 0.0001 a.u. MC iterations, which is a very short propagation time engaging only ∼0.1% of the total walker populations at τ = 8.0 a.u., where we terminated our i-CISDTQ-MC runs (cf. Table S.5 of the supplementary material), and i-CISDTQ-MC capturing as little as 14-17% of all triply excited determinants, the semi-stochastic CC(P ;Q) methodology reduces the 13.370 millihartree difference between the CR-CC(2,3) and CCSDT energies of the B 1 A 1 state and the 8.1 kcal/mol error in the CR-CC(2,3) value of the B 1 A 1 − X 3 A ′ 2 gap relative to CCSDT to 1.260 millihartree and 0.6 kcal/mol, respectively, which is a chemical accuracy regime. Interestingly, the i-CISDTQ-MCbased CC(P ;Q) value of ∆E S-T obtained after 6000 MC iterations matches the quality of the B 1 A 1 − X 3 A ′ 2 gap resulting from the fully deterministic CC(P ;Q) calculations using the CC(t;3) approach, which give a 0.5 kcal/mol error relative to CCSDT when the cc-pVDZ basis set is employed. 52 After the additional 4000 MC time steps, where the i-CISDTQ-MC propagations for the X 3 A ′ 2 and B 1 A 1 states are still very far from convergence and where the fractions of triples captured by i-CISDTQ-MC increase to about 30%, the small errors relative to CCSDT characterizing the i-CISDTQ-MCbased CC(P ;Q) values of the energy of the B 1 A 1 state and ∆E S-T at τ = 0.6 a.u. drop down by factors of 4-6, to 0.314 millihartree and 0.1 kcal/mol, respectively, illustrating how rapid the convergence of the CIQMCdriven CC(P ;Q) calculations toward the parent SRCC data can be. While the CR-CC(2,3) description of the X 3 A ′ 2 state, which has a SR character, is much better than in the case of its strongly correlated B 1 A 1 counterpart, the semi-stochastic CC(P ;Q) computations offer great improvements in this case too. They are, for example, capable of reducing the ∼0.4 millihartree difference between the CR-CC(2,3) and CCSDT energies to a 0.1 millihartree level after 10000 MC iterations and i-CISDTQ-MC capturing less than 30% of all triples. In analogy to all other molecular examples considered in this article, the uncorrected CC(P ) values of the energies of the X 3 A ′ 2 and B 1 A 1 states and separation between them converge to their CCSDT limits too, but it is clear from Table VII and Fig. 5 that they do it at a much slower rate than their CC(P ;Q) counterparts. This can be illustrated by comparing the errors relative to CCSDT characterizing the CC(P ) and CC(P ;Q) energies of the X 3 A ′ 2 and B 1 A 1 states and separation between them obtained after 6000 MC iterations. They are more than 11 millihartree, about 21 millihartree, and almost 6 kcal/mol, respectively, in the former case and only 0.253 millihartree, 1.260 millihartree, and 0.6 kcal/mol, when the CC(P ) energies are corrected for the remaining T 3 correlations using the CC(P ;Q) approach. As explained in Section II, the CC(P ) energies converge to CCSDT more slowly than their δ(P ; Q)-corrected CC(P ;Q) counterparts, since the initial, τ = 0, CC(P ) calculation for a given electronic state is equivalent to CCSD, where T 3 = 0. The CIQMC-driven CC(P ;Q) calculations start from CR-CC(2, 3), which provides information about T 3 clusters via noniterative corrections to CCSD. This once again emphasizes the benefits of using corrections δ(P ; Q) in the context of the semi-stochastic CC(P ;Q) work.

IV. CONCLUSIONS
One of the most promising ideas in the area of converging high-level SRCC energetics without having to resort to the very expensive methods such as CCSDT or CCSDTQ has been the CC(P ;Q) formalism, originally introduced in Refs. 50-52, in which one solves the CC amplitude equations in a suitably defined subspace of the many-electron Hilbert space, called the P space, and then improves the resulting CC(P ) energies using the a posteriori corrections δ(P ; Q) determined with the help of another subspace of the Hilbert space, referred to as the Q space. In addition to conventional choices of the P and Q spaces, which result in methods such as CR-CC(2,3), 54-58 one can consider various unconventional ways of setting up these spaces in the CC(P ;Q) calculations. [50][51][52][53]59,65,[95][96][97]102 In this work, we have focused on the semi-stochastic formulation of the CC(P ;Q) formalism, introduced in Ref. 95 and further developed in Refs. 96 and 97, in which the leading higher-than-doubly excited determinants incorporated in the P space, needed to solve the CC(P ) equations, are identified, in an automated fashion, by the stochastic wave function propagations employing the CIQMC methodology of Refs. 86-90. The specific objective of this study has been the investigation of the effectiveness of the semi-stochastic CC(P ;Q) approaches driven by i-FCIQMC and i-CISDTQ-MC propagations, which were used to capture the leading triply excited determinants for the inclusion in the respective P spaces, in converging the singlettriplet gaps and the underlying total electronic energies resulting from the high-level CCSDT computations. Molecular systems that have been used to examine the ability of the semi-stochastic CC(P ;Q) methodology to recover the CCSDT energetics of the low-lying singlet and triplet states included the methylene, cyclobutadiene, cyclopentadienyl cation, and trimethylenemethane biradicals and a prototype magnetic system represented by the (HFH) − ion.
Cyclopentadienyl cation and trimethylenemethane are the largest polyatomic species that have been used in the semi-stochastic CC(P ;Q) calculations reported to date.
An accurate determination of the singlet-triplet gaps in the above five systems requires a well-balanced treatment of substantial nondynamical correlation effects, needed for a reliable description of the low-spin singlet states that have a manifestly MR character, and dynamical correlations of the high-spin triplet states, which are SR in nature. Within the conventional, particleconserving, black-box SRCC framework, the only methods that can do this in a robust manner are the highlevel CC approaches with a full treatment of T n correlations with n > 2, beginning with CCSDT. The basic CCSD approximation, which ignores the T n cluster components with n > 2 altogether, and the noniterative T 3 corrections to CCSD, including even CR-CC(2,3), which can handle selected MR situations, such as single bond breaking, fail. Putting aside the quasi-perturbative character of the majority of the existing triples corrections to CCSD, a significant part of the problem is their inability to capture the coupling between the low-rank T 1 and T 2 clusters with their higher-rank T 3 counterpart, which cannot be ignored when the low-lying singlet states of biradicals are considered. We have shown in this study that by relaxing T 1 and T 2 amplitudes in the presence of the T 3 component defined using the list of the leading triply excited determinants identified by the CIQMC propagations and using the δ(P ; Q) corrections to account for the remaining T 3 correlations that are not described by the underlying CC(P ) calculations, the semi-stochastic CC(P ;Q) methodology provides an efficient mechanism for incorporating the coupling among the T 1 , T 2 , and T 3 clusters relevant to a reliable description of the singlettriplet gaps in biradicals, while capturing the vast ma-jority of the many-electron correlation effects included in the parent CCSDT computations. In particular, we have demonstrated that the semi-stochastic CC(P ;Q) calculations are capable of reaching millihartree or submillihartree accuracy levels relative to the parent CCSDT results, including total electronic energies of the lowest singlet and triplet states and adiabatic as well as vertical gaps between them, with small (typically about 20-30%; sometimes even less) fractions of triply excited determinants captured in the early stages of the underlying CIQMC propagations and with tiny walker populations that are orders of magnitude smaller than the total numbers of walkers required to converge them. We have also demonstrated the vital role of the noniterative corrections δ(P ; Q) in accelerating and, in the case of the singlet-triplet gap values, smoothing convergence of the corresponding CC(P ) energetics toward CCSDT.
The numerical results and analyses reported in this article encourage us to pursue the semi-stochastic CC(P ;Q) methodology even further. Putting aside the need for improving the computational efficiency of our CIQMCdriven CC(P ;Q) codes, so that we could consider larger molecular problems than those considered in this study, it would be interesting to examine if our semi-stochastic CC(P ;Q) computations, including those discussed in this work, could take advantage of the recent advances in CIQMC, such as the adaptive CIQMC algorithm of Refs. 89 and 90 that might replace the i-CIQMC approaches that we have used so far. As implied by the remarks made in Section III, some of the biradicals considered in this work might benefit from an accurate description of T 3 as well as T 4 correlations. It would, therefore, be interesting to examine if one could improve the results presented in this work by using the semi-stochastic CC(P ;Q) approaches that aim at converging the CCSDTQ energetics, such as those discussed in Ref. 97. Last, but not least, we have started developing extensions of the semi-stochastic CC(P ) and CC(P ;Q) approaches, pursued in Refs. 95 and 97 and this study, and their EOMCC analogs 96,98 to the particle-nonconserving EOMCC models, especially EA/IP-EOMCC and DEA/DIP-EOMCC. As pointed out in this article, and as shown in Refs. [123][124][125]135, and 136, the deterministic DEA/DIP-EOMCC approaches with the explicit inclusion of 4p2h and 4h2p components of the corresponding electron attaching and ionizing operators on top of CCSD can provide a nearly exact description of the singlet-triplet in biradical species, so it will be interesting if their semistochastic counterparts can do the same.

DATA AVAILABILITY
The data that support the findings of this study are available within the article and its supplementary material. I. Convergence of the CC(P ) and CC(P ;Q) energies of the X 3 B1 and A 1 A1 states of methylene, as described by the aug-cc-pVTZ basis set, and of the corresponding adiabatic singlet-triplet gaps toward their parent CCSDT values. The geometries of the X 3 B1 and A 1 A1 states, optimized in the FCI calculations using the TZ2P basis set, were taken from Ref. 149. The P spaces used in the CC(P ) and CC(P ;Q) calculations were defined as all singly and doubly excited determinants and subsets of triply excited determinants extracted from the i-FCIQMC propagations with δτ = 0.0001 a.u. The Q spaces used to determine the CC(P ;Q) corrections consisted of the triply excited determinants not captured by the corresponding i-FCIQMC runs. The i-FCIQMC calculations preceding the CC(P ) and CC(P ;Q) steps were initiated by placing 1500 walkers on the ROHF (X 3 B1 state) and RHF (A 1 A1 state) reference determinants and the na parameter of the initiator algorithm was set at 3. In all post-Hartree-Fock calculations, the lowest core orbital was kept frozen and the spherical components of d and f orbitals were employed throughout. TABLE II. Convergence of the CC(P ) and CC(P ;Q) energies of the X 1 Σ + g state of (HFH) − , as described by the 6-31G(d,p) basis set, at selected H-F distances RH-F toward their parent CCSDT values. The P spaces used in the CC(P ) and CC(P ;Q) calculations were defined as all singly and doubly excited determinants and subsets of triply excited determinants extracted from the i-FCIQMC propagations with δτ = 0.0001 a.u. The Q spaces used to determine the CC(P ;Q) corrections consisted of the triply excited determinants not captured by the corresponding i-FCIQMC runs. The i-FCIQMC calculations preceding the CC(P ) and CC(P ;Q) steps were initiated by placing 1500 walkers on the RHF reference determinant and the na parameter of the initiator algorithm was set at 3. In all post-Hartree-Fock calculations, the lowest core orbital was kept frozen and the spherical components of d orbitals were employed throughout.  III. Convergence of the CC(P ) and CC(P ;Q) energies of the A 3 Σ + u state of (HFH) − , as described by the 6-31G(d,p) basis set, at selected H-F distances RH-F toward their parent CCSDT values. The P spaces used in the CC(P ) and CC(P ;Q) calculations were defined as all singly and doubly excited determinants and subsets of triply excited determinants extracted from the i-FCIQMC propagations with δτ = 0.0001 a.u. The Q spaces used to determine the CC(P ;Q) corrections consisted of the triply excited determinants not captured by the corresponding i-FCIQMC runs. The i-FCIQMC calculations preceding the CC(P ) and CC(P ;Q) steps were initiated by placing 1500 walkers on the ROHF reference determinant and the na parameter of the initiator algorithm was set at 3. In all post-Hartree-Fock calculations, the lowest core orbital was kept frozen and the spherical components of d orbitals were employed throughout.  V. Convergence of the CC(P ) and CC(P ;Q) energies of the X 1 B1g and A 3 A2g states of cyclobutadiene, as described by the cc-pVDZ basis set, and of the corresponding vertical singlet-triplet gaps toward their parent CCSDT values. All calculations were performed at the D 4h -symmetric transition-state geometry of the X 1 B1g state optimized in the MR-AQCC calculations in Ref. 167. The P spaces used in the CC(P ) and CC(P ;Q) calculations were defined as all singly and doubly excited determinants and subsets of triply excited determinants extracted from the i-FCIQMC propagations with δτ = 0.0001 a.u. The Q spaces used to determine the CC(P ;Q) corrections consisted of the triply excited determinants not captured by the corresponding i-FCIQMC runs. The i-FCIQMC calculations preceding the CC(P ) and CC(P ;Q) steps were initiated by placing 1500 walkers on the RHF (X 1 B1g state) and ROHF (A 3 A2g state) reference determinants and the na parameter of the initiator algorithm was set at 3. In all post-Hartree-Fock calculations, the four lowest core orbitals were kept frozen and the spherical components of d orbitals were employed throughout.  VI. Convergence of the CC(P ) and CC(P ;Q) energies of the X 3 A ′ 2 and A 1 E ′ 2 states of cyclopentadienyl cation, as described by the cc-pVDZ basis set, and of the corresponding vertical singlet-triplet gaps toward their parent CCSDT values. All calculations were performed at the D 5h -symmetric geometry of the X 3 A ′ 2 state optimized using the unrestricted CCSD/cc-pVDZ approach in Ref. 118. The P spaces used in the CC(P ) and CC(P ;Q) calculations were defined as all singly and doubly excited determinants and subsets of triply excited determinants extracted from the i-CISDTQ-MC propagations with δτ = 0.0001 a.u. The Q spaces used to determine the CC(P ;Q) corrections consisted of the triply excited determinants not captured by the corresponding i-CISDTQ-MC runs. The i-CISDTQ-MC calculations preceding the CC(P ) and CC(P ;Q) steps were initiated by placing 1500 walkers on the ROHF (X 3 A ′ 2 state) and RHF (A 1 E ′ 2 state) reference determinants and the na parameter of the initiator algorithm was set at 3. In all post-Hartree-Fock calculations, the five lowest core orbitals were kept frozen and the spherical components of d orbitals were employed throughout. The %T values are the percentages of triples captured during the i-CISDTQ-MC propagations [the Sz = 1 triply excited determinants of the B 2 (C 2v ) symmetry in the case of the X 3 A ′ 2 state and the Sz = 0 triply excited determinants of the A 1 (C 2v ) symmetry in the case of the A 1 E ′ 2 state]. c Unless otherwise specified, the values of the singlet-triplet gaps are reported as errors relative to CCSDT in kcal/mol.  VII. Convergence of the CC(P ) and CC(P ;Q) energies of the X 3 A ′ 2 and B 1 A1 states of trimethylenemethane, as described by the cc-pVDZ basis set, and of the corresponding adiabatic singlet-triplet gaps toward their parent CCSDT values. The D 3h -and C2v-symmetric geometries of the X 3 A ′ 2 and B 1 A1 states, respectively, optimized in the SF-DFT/6-31G(d) calculations, were taken from Ref. 115. The P spaces used in the CC(P ) and CC(P ;Q) calculations were defined as all singly and doubly excited determinants and subsets of triply excited determinants extracted from the i-CISDTQ-MC propagations with δτ = 0.0001 a.u. The Q spaces used to determine the CC(P ;Q) corrections consisted of the triply excited determinants not captured by the corresponding i-CISDTQ-MC runs. The i-CISDTQ-MC calculations preceding the CC(P ) and CC(P ;Q) steps were initiated by placing 1500 walkers on the ROHF (X 3 A ′ 2 state) and RHF (B 1 A1 state) reference determinants and the na parameter of the initiator algorithm was set at 3. In all post-Hartree-Fock calculations, the four lowest core orbitals were kept frozen and the spherical components of d orbitals were employed throughout. The %T values are the percentages of triples captured during the i-CISDTQ-MC propagations [the Sz = 1 triply excited determinants of the B 2 (C 2v ) symmetry in the case of the X 3 A ′ 2 state and the Sz = 0 triply excited determinants of the A 1 symmetry in the case of the B 1 A 1 state]. c Unless otherwise specified, the values of the singlet-triplet gaps are reported as errors relative to CCSDT in kcal/mol.