Accuracy of the distinguishable cluster approximation for triple excitations for open-shell molecules and excited states

Open-shell molecular systems play an important role in numerous biological and non-biological processes. Ionization energies (IEs), electron affinities (EAs), excitation energies, or spin gaps are relevant in many fields including material design and in the explanation of biological processes. Coupled-cluster (CC) methods¹ are often used for high accuracy ab initio benchmark calculations of the thermochemical properties.² Thus, it is highly desirable to develop and benchmark new CC based methods with an improved accuracy to cost ratio.

The CC hierarchy of methods systematically improves the description of the wavefunction. This improvement is typically associated with a high computational cost and an increased scaling of the computational resources with the molecular size, which often impedes the use of arbitrarily elaborate CC methods. Applications of CC methods, therefore, usually have to make an inevitable compromise between accuracy and computational cost.

However, it has been shown in various studies that often one can gain accuracy without increasing the cost,^3–10 which has been termed “internal correction”¹¹ or “addition by subtraction.”^8,12 These methods work by emulating the effects of left-out clusters through modification of the amplitude equations. The philosophy and motivation behind the various modifications differ, but commonly the methods leave out some terms in the amplitude equations while retaining important features of the full configuration interaction (FCI) solution. The features are exactness for N electrons if the hierarchy is truncated at the Nth excitation, size-extensivity, and invariance toward virtual–virtual and occupied–occupied orbital rotations. The distinguishable-cluster (DC) approximation¹⁰ enforces additionally the particle–hole symmetry, the importance of which was recently emphasized in Ref. 13, and focuses primarily on the quadratic exchange terms in the amplitude equations.

The DC approximation was first applied to the doubles amplitude equations of Brueckner CC doubles (BCCD),^14,15 which resulted in the Brueckner DC doubles (BDCD) method.¹⁰ Later, it was combined with single excitations and a complete orbital optimization.¹⁶ The DC methods have been shown to yield a qualitatively correct description of strongly correlated systems;^10,16 it was also recognized that the DC approximation improves the CC accuracy in the weakly correlated regime.^16–18 Recently, the DC approximation was applied to the doubles–triples quadratic part of the triples amplitude equations in coupled-cluster with singles doubles and triples (CCSDT),¹⁹ which was termed DC-CCSDT.²⁰ Opposed to CCSDT, DC-CCSDT can be formulated to scale with $O(N7)$ with the molecular size N in the real space representation, whereas for CCSDT, the initial $O(N8)$ scaling remains.²⁰ Since all the quadratic doubles terms are reintroduced into the equations, this method shows similarly poor results in the strongly correlated regime as CCSDT, but was shown to improve considerably upon CCSDT and coupled-cluster with singles and doubles (CCSD) and perturbative triples [CCSD(T)]²¹ in calculating the total and reaction energies of small closed-shell molecules as well as breaking single bonds.^20,22

Open-shell molecules possess several difficulties to quantum chemical methods,²³ which are absent in closed-shell systems, and therefore, a method working well for closed-shell systems might perform worse for systems with unpaired electrons. In this paper, we present for the first time the DC-CCSDT calculations involving open-shell species and excited states. Excitation energies have been calculated using the ΔCC approach,^24–28 i.e., by converging the coupled-cluster and the distinguishable cluster methods to ground and excited states and by computing energy differences. For open-shell excited states, we have introduced a tailoring technique to fix the reference configuration state function (CSF) to the desired spin multiplicity, termed the fixed-reference approach.

In Sec. II, the theory of DC-CCSDT will be revisited and the excited state optimization strategy will be introduced. In Sec. III, various calculations using our new DC-CCSDT implementation, applicable to open-shell molecules, are presented.

The DC-CCSDT method is obtained by excluding exchange terms of the $T̂2T3̂$ part of CCSDT triples amplitude equations that scale as $O(N8)$ in the real space representation and modifying the rest—in a particle–hole symmetric fashion—to ensure exactness for three electrons. The theory of DC-CCSDT has already been presented in previous publications.^20,22 In the following, we will show explicitly how the correct cancellations can be found. Alternatively, one can use a symbolic algebra program.^22,29 Explicit amplitude equations of the CCSDT method in the spin-orbital basis can be found, e.g., in Ref. 30. The mixed doubles–triples quadratic contribution to the triples residual, $Rabcijk$⁠, is as follows:

where i, j, k, l, m and a, b, c, d, e denote the occupied and virtual spin-orbitals, respectively, and $P̂$ is a permutation operator, e.g.,

Here and in the following, the implicit summation over repeated indices is assumed. Following the nomenclature from Ref. 10, we denote the terms in Eq. (1) as A, A′, B₁, B₂, C₁, C₂, D₁, and D₂, respectively. In the distinguishable cluster approximation, the Coulomb term A remains unmodified and the exchange term A′ is omitted. The terms B₁ and B₂ are eliminated according to the local particle–hole symmetry.^13,31 The remaining terms have to be modified to ensure exactness for three electrons. Since the indices i, j, and k in the residual $Rabcijk$ have to be different, and for three electrons there are only three occupied spin-orbitals, we assign index i with the occupied spin-orbital 1, j with 2, and k with 3. Writing out the permutation operators and renaming indices to ensure that they all correspond to the same integral lead to

Considering different occupied spin-orbitals for l and m leads to

which can also be extended to the third occupied spin-orbital. Summing up the equations corresponding to each type, the following equalities can be recognized:

We note that because of Eq. (27), in principle, any linear combination could have been used to cancel term A′. The one chosen in Eq. (26) corresponds to the particle–hole symmetric variant. The above-mentioned equations define the DC-CCSDT method with the doubles–triples contribution to the residual

In the ΔCC approaches,^24–28 first a reference configuration state function (CSF) has to be optimized for the desired state. For the closed-shell or high-spin CSFs, i.e., single Slater determinants, the optimization can be achieved by fixing the occupation pattern in the (restricted open-shell) Hartree–Fock (ROHF) calculations, e.g., by ensuring the maximum overlap of the orbitals with the starting set (nitord option in molpro). The starting set of orbitals and the corresponding occupation patterns can be obtained from state-averaged complete active-space self-consistent field (CASSCF) or low-level equation-of-motion (EOM) CC calculations. Here, the former approach is employed.

For states dominated by a single excitation from the ground state, one can use the multiconfiguration self-consistent field (MCSCF) procedure, after a CASSCF preoptimization (or any other procedure to define the pair of orbitals mainly involved in the excitation), to perform the single-CSF optimization by restricting the orbital occupation according to the reference CSF, which can be done in molpro using select and con keywords. For the M_s = 0 triplet states, we simply employ the corresponding M_s = 1 triplet orbitals from the high-spin restricted open-shell Hartree–Fock (ROHF) calculations. Subsequently, one of the determinants of the open-shell M_s = 0 CSF is chosen to be the reference determinant, while the second determinant can be reached by a double excitation from the reference. In order to restrict the DC-CCSDT optimizations to the singlet or (M_s = 0)-triplet excited states, the doubles amplitude corresponding to the second determinant is fixed to 1.0 or −1.0, which equals to the weight of the reference determinant in the intermediately normalized wavefunction. All other amplitudes are relaxed according to the DC-CCSDT amplitude equations.

This optimization strategy bears resemblance to the recent work on ΔCC for doubly excited states and core-ionization potentials,^24,25 to optimizing the mean-field solutions for the excited states,^32,33 to the two-determinant coupled-cluster methods,³⁴ and to the tailored coupled-cluster methods.³⁵ 

Compared to the usual EOM approach, the ΔCC technique treats the ground and excited states at the same level, and the excited states for which EOM requires high-level excitations, e.g., doubly excited states, are much more accurate with ΔCC.²⁴ The fixed-reference technique extends this approach to open-shell excited states. The main advantage of this approach is the simplicity of the implementation and potentially higher accuracy compared to the EOM methods. However, one should also be aware of some potential drawbacks. The accuracy of the single-reference methods for excited states with multiple important configurations will usually be low, and so this technique is not well suited for such states. The excited state calculations can be hard to converge, and the quality of the starting guess is very important. In addition, the calculations are much less of the black-box type, and therefore, prior investigations using CASSCF or EOM-CCSD are necessary. Additionally, the fixed-reference approach is not spin adapted, and since the determinants of the reference CSF are treated differently, some spin contamination will occur. In order to investigate the severeness of the spin contamination, we calculate the S² expectation value numerically by adding the (scaled) S² integrals to the one- and two-body parts of the Hamiltonian and using a three-point forward finite difference quotient.

The DC-CCSDT method in the spin-orbital formalism was implemented using the Integrated Tensor Framework (ITF)^36,37 into the development version of Molpro.^38,39

The accuracy of DC-CCSDT for calculating the thermochemical properties and excited states of small molecules, the 3d transition metal atoms, the ozone molecule, and the active region of an iron–porphyrin model system is compared to various coupled-cluster methods, such as CCSD(T), CCSDT, and CCSDTQ.⁴⁰ 

The CCSDT and CCSDTQ calculations have been performed with the mrcc^41,42 program interfaced to Molpro. The CCSDT calculations with the fixed-reference approach have been done using an ITF implementation. For the CCSD(T) calculations, we used the unrestricted implementation available in Molpro.

Unless specified otherwise, restricted open-shell Hartree–Fock (ROHF) orbitals and the frozen core approximation have been utilized.

Ionization energies (IEs), electron affinities (EAs), and atomization energies (AEs) of several molecules have been calculated. We have used the aug-cc-pVTZ basis set⁴³ for first-row elements and the aug-cc-pV(T+d)Z basis set⁴⁴ for second-row elements, abbreviated as AVTZ and AV(T+d)Z, respectively. For EAs, to describe the extra electron, a more diffuse A2VTZ basis set with additional diffuse functions as described in Ref. 45 has been used. The thermochemical properties are benchmarked against CCSDTQ values. Calculations with CCSD, DCSD, CCSD(T), and CCSDT are also presented for perspective. Calculations of thermochemical properties have been performed on a subset of the original G3 benchmark set.⁴⁶ The geometries were taken from Ref. 45.

In Fig. 1, the deviations of adiabatic IEs calculated using DC-CCSDT, CCSDT, CCSD(T), DCSD, and CCSD for C, N, O, F, Si, P, S, Cl, OH, HF, N₂ [$→N2+$(²Σ_g)], N₂ [$→N2+$(²Π_u)], CO, O₂, PH, SH, HCl, CS, H₂O, PH₂, SH₂ [$→SH2+$(²B₁)], and SH₂ [$→SH2+$(²A₁)] from CCSDTQ values are plotted. As shown previously,¹⁸ the accuracy of DCSD is much higher compared to CCSD. Methods with explicit triples are even more accurate. DC-CCSDT is consistently closer to the CCSDTQ values than CCSDT and CCSD(T). The DC-CCSDT maximal (MaxD) and root-mean squared (RMSD) deviations for this set are 30 and 9 meV, respectively, which can be compared to 50 and 13 meV, and 38 and 13 meV of CCSDT and CCSD(T), respectively (cf. Table I).

The performance of DC-CCSDT and other methods for calculating adiabatic electron affinities of C, O, F, Si, P, S, Cl, CH, NH, OH, SiH, PH, SH, CN, NO, O₂, CH₂, NH₂, SiH₂, PH₂, and CH₃ systems is shown in Fig. 2. As before, DCSD is much more accurate than CCSD; however, explicit triples are crucial to achieve the chemical accuracy below 0.04 eV, and iterative triples in CCSDT and DC-CCSDT perform better than the perturbative triples correction in CCSD(T), although the gain in the case of the conventional coupled cluster is rather low. The accuracy of DC-CCSDT in these cases is remarkably good: for all studied systems, the difference in EAs of DC-CCSDT compared to CCSDTQ is less than 7 meV, and on average, it is just 2 meV. The maximal deviations of CCSD(T) and CCSDT are four and two times larger, respectively. RMSD of DC-CCSDT is just 3 meV, which is almost three times smaller compared to CCSDT.

Additionally, atomization energies of CH, NH, OH, HF, SiH, PH, SH, HCl, CN, N₂, CO, NO, O₂, F₂, Si₂, SiO, CS, P₂, CH₂ (¹A₁ and ³B₁), NH₂, H₂O, SiH₂ (¹A₁ and 3B₁), PH₂, SH₂, HCN, and CH₃ systems have been calculated. In Fig. 3, it can be seen that DC-CCSDT outperforms CCSDT for the calculation of atomization energies in almost every system tested. For some of the systems, CCSD(T) results are more accurate than DC-CCSDT ones; however, for most of the systems, DC-CCSDT AEs are much closer to the CCSDTQ reference values.

The mean of the absolute deviations (MADs), the mean of deviations (MDs), RMSD, and MaxD of the thermochemical properties compared to CCSDTQ reference values are provided in Table I. The statistical measures underline that DC-CCSDT is consistently closer to the CCSDTQ benchmark results than the other tested methods in calculating IEs, EAs, and AEs.

Excited states of a few small molecules from Ref. 47 have been calculated using the fixed-reference technique outlined in Sec. II B. The AVTZ basis set has been utilized, and the excitation energies are compared to the extrapolated FCI values given by Loos et al.⁴⁷ denoted as exFCI. As explained in Sec. II B, in the calculations of excited states using DC-CCSDT, only one single element of the doubles amplitudes, which corresponds to the second determinant of the open-shell singlet or M_s = 0 triplet “reference” CSF, has been fixed to ensure convergence to the corresponding states. For the M_s = 0 triplet states, M_s = 1 triplet ROHF orbitals have been utilized, and for the singlet states, the MCSCF single-CSF procedure outlined in Sec. II B has been employed. The vertical excitation energies for several states and molecules are shown in Table II. DC-CCSDT clearly outperforms EOM-CCSDT for the calculated excitation energies and rivals EOM-CCSDTQ in accuracy. The ³$Σu+$ state of the nitrogen molecule has a large multi-configurational character, since the four determinants are equally important to describe the M_s = 0 triplet state, and a large error is to be expected. Note that the corresponding M_s = 1 triplet state does not suffer from this problem, and the excitation energy 7.711 eV is close to the reference value.

For comparison, in addition to DC-CCSDT, fixed-reference CCSDT excitation energies have been computed. The accuracy of CCSDT is better or comparable to EOM-CCSDT in most of the cases. For example, for the water molecule, all excitation energies from the fixed-reference CCSDT approach match exFCI and/or EOM-CCSDTQ numbers, and the EOM-CCSDT ones differ by up to 40 meV. However, in nearly all cases, the distinguishable cluster approximation yields even closer (or similarly close) results to the reference values.

The S² expectation value of the excited states from the fixed-reference approach demonstrates that the states are nearly pure spin states. Specifically, states calculated using DC-CCSDT exhibit a very low amount of spin contamination, usually well below 0.01. This shows that the high excitations can adapt to the fixed reference well, and the triples can describe many of the single and double excitations from the second Slater determinant (in particular, double excitation involving at least one of the open-shell electrons). The doubles-only theories are obviously more affected by the unbalanced treatment of the two determinants of the reference CSF. In Table III, the corresponding data for CCSD and DCSD are shown in comparison to DC-CCSDT and EOM-CCSD. As expected, the spin contamination is much larger for the methods without triples, with DCSD again having an edge over CCSD, and also the accuracy of the excitation energies is lower compared to EOM-CCSD. It can be concluded that the fixed-reference approach can be recommended only for methods with explicit triple excitations.

The accuracy of the fixed-reference technique can additionally be estimated by comparing the M_s = 0 and M_s = 1 triplet energies with the high-spin triplet energies calculated in the common way without the fixed reference approach. The corresponding energy differences are shown in Table IV (excluding the problematic ³$Σu+$ nitrogen state). The largest discrepancies for DC-CCSDT are 18 and 11 meV in the ³Π_g state of N₂ and the ³A₁ state of formaldehyde, respectively. In all other cases, the energies differ by only 0–3 meV. The largest discrepancy for CCSDT is 27 meV.

Ionization energies and excitation energies have been calculated for 3d transition metal atoms and compared to CCSDT, extrapolated FCI, and experimental values of Balabanov and Peterson⁴⁸ in Tables V and VI. The DC-CCSDT values have been calculated using the same composite approach as in Ref. 48 according to

where $EAVTZ-DKDC-CCSDT$ are valence-only DC-CCSDT energies computed using the second-order Douglas–Kroll–Hess Hamiltonians and the aug-cc-pVTZ-DK basis set, ΔE_CV are semi-core energy contributions computed at the CCSD(T)/CBS level as (E_3s3p3d4s–E_3d4s), and ΔE_CBS are the basis-set corrections computed at the CCSD(T)/CBS level of theory. The two corrections ΔE_CV and ΔE_CBS are taken from the publications of Balabanov and Peterson.^48,49 Here, only high-spin states have been calculated, and therefore, the pure ΔCC approach without the fixed-reference technique has been used.

In most of the cases, the DC-CCSDT results are closer to the extrapolated FCI results than CCSDT ones. However, the improvement is less pronounced than in the case of molecular systems (cf. Secs. III A and III B), and in the case of ionization energies of copper and zinc atoms, the DC-CCSDT results are noticeably worse with discrepancies of 0.3 and 0.13 kcal/mol, respectively. One should mention that the extrapolated FCI and the experimental ionization energies differ in some cases by 0.9 kcal/mol. The performance of DC-CCSDT for the atomic excitation energies is better. For most of the atoms, the distinguishable cluster approximation reduces the CCSDT error compared to the extrapolated FCI values by 1.5–2 times.

The ozone molecule is challenging from a computational point of view because of its notorious multiconfigurationality, and high excitations or the multireference treatment are required for an accurate description. In 1972, Hay and Goodard predicted a metastable ring geometry of ozone,⁵⁰ which is theoretically well established after numerous studies, but still lacks experimental evidence. We have calculated the energetics of the ring minimum (RM) and the transition state (TS) leading to it from the ground-state geometry called open-ring minimum (OM) with DC-CCSDT using the cc-pVTZ basis set. Furthermore, the 2¹A₁–1¹A₁ vertical excitation energies of the RM, TS, and OM geometries have been calculated as well, which correspond to doubly excited states and for which the simple ΔCC without the fixed-reference approach has been used. The geometries have been taken from Ref. 51. The order of states corresponds to Ref. 52. The reference determinants have been obtained by running state-averaged CASSCF calculations for two lowest states using the full valence active space (18, 12) and reoptimizing the orbitals of the individual closed-shell reference determinants at the HF level by fixing the occupation pattern. In the OM geometry, the ground state reference determinant corresponds to $(a1)12(b1)2(b2)8(a2)2$ occupation and the first excited state to $(a1)10(b1)4(b2)8(a2)2$⁠; in the TS geometry, to $(a1)12(b1)4(b2)6(a2)2$ and to $(a1)12(b1)2(b2)8(a2)2$⁠; and in the RM geometry, to $(a1)12(b1)4(b2)6(a2)2$ and to $(a1)12(b1)4(b2)8$⁠, respectively. In Table VII, the DC-CCSDT results are compared to CCSDT, CCSDT(Q), and tailored DCSD (TDCSD) using a large active space of (18, 39),⁵³ adaptive shift FCI Quantum Monte Carlo (AS-FCIQMC),⁵⁴ and semistochastic heat-bath CI (SHCI)⁵¹ values.

The DC-CCSDT absolute energy for the OM ground state differs from the AS-FCIQMC result by 2 mE_h, which is similar to the AS-FCIQMC–CCSDT(Q) discrepancy and nearly three times smaller than that of CCSDT. The relative energies with respect to the OM geometry and the excitation energies from DC-CCSDT agree very well with the reference values from SHCI, AS-FCIQMC, or TDCSD. In most of the cases, DC-CCSDT is much closer than CCSDT or even CCSDT(Q) to the reference values. Only the TS-OM energy difference is somewhat less accurate at the DC-CCSDT level, which can be explained by the fact that the TS states are nearly degenerate for DC-CCSDT (as well as for TDCSD and SHCI), and the degeneracy is less pronounced at the CCSDT and CCSDT(Q) level of theory.

Transition metal complexes are another example of challenging molecular systems, which often require multi-reference or high-order methods for qualitatively or quantitatively accurate description. Although single-center transition metal compounds usually can be treated using single-reference methods, high excitations are required for sub-kJ/mol accuracy. One possibility to treat such systems is by optimizing a large active space using, e.g., the stochastic CASSCF⁵⁵ and evaluating the spin-state energies by performing explicitly correlated CCSD(T) calculations on the full molecular system and correcting the results using a subtractive embedding approach and high-level coupled-cluster calculations inside the active space.⁵⁶ However, high-level coupled-cluster methods quickly become expensive, even for active-space only calculations. For example, one iteration of CCSDTQ in (40, 38) active spaces employed in Ref. 56 required 5 h wall time using eight central processing unit (CPU) cores. At the same time, one iteration of CCSDT for the same system can be run in 1 min wall time. However, the spin gap calculated at the CCSDT level is 1.7 kJ/mol smaller than the CCSDTQ gap. Thus, it would be extremely beneficial to approach the CCSDTQ accuracy at the CCSDT cost. In Table VIII, the DC-CCSDT results are compared to the CCSD(T), CCSDT, and CCSDTQ values. The calculations have been performed using the active space (40, 38), as optimized in Ref. 56 at the stochastic-CASSCF level and the ANO-RCC-VTZP basis set. The active space includes four doubly occupied σ-orbitals pointing at the metal center, the full π system, five 3d orbitals of the iron atom together with five empty correlating d′, four semicore 3s3p, and four empty 4s4p orbitals.

Already the absolute energies from DC-CCSDT are much closer to CCSDTQ results than the CCSDT energies. For both spin states, the difference is just a few tenths of millihartree. As in our previous benchmarks, the high accuracy of absolute energies carries over to the relative energies. The DC-CCSDT spin gap differs from the CCSDTQ one by 0.6 kJ/mol, which is nearly three times smaller than the CCSDT–CCSDTQ discrepancy. Thus, DC-CCSDT can be employed as a replacement for expensive CCSDTQ calculations without large sacrifice of the accuracy.

The DC-CCSDT method has been implemented using spin-orbitals and can be applied now to open-shell systems. The thermochemical properties and excited states of different systems involving closed-shell and open-shell species have been calculated using the new implementation. Computation of open-shell excited states has been performed using a fixed-reference technique, which ensures that the coefficients of determinants within the corresponding reference CSF are fixed according to the correct spin.

Ionization energies, electron affinities, and atomization energies of various molecules and atoms calculated with DC-CCSDT have been consistently closer to the CCSDTQ values than the CCSDT and CCSD(T) results, especially large is the improvement in the electron affinities, which are on average four times more accurate using DC-CCSDT than with CCSDT. The fixed-reference technique combined with DC-CCSDT yields extremely accurate excitation energies for excited states dominated by a single CSF, outperforming EOM-CCSDT and reaching the EOM-CCSDTQ accuracy, with only a negligible amount of the spin contamination. Comparison of M_s = 0 and M_s = 1 triplet states confirms the accuracy of the method. The results of ionization energies of the 3d transition metal atoms are less unequivocal. For some of the atoms, the distinguishable cluster approximation improves the results, while for others, CCSDT turned out to be more accurate. On the contrary, the DC-CCSDT excitation energies of the 3d transition metal atoms are consistently more accurate than with the usual CCSDT method.

The distinguishable cluster approach also improves the description of the ground and excited states of the ozone molecule at various geometries. In most of the cases, the DC-CCSDT absolute and relative energies are closer to the best theoretical values available in the literature than CCSDT and CCSDT(Q) results. The enhanced accuracy of the distinguishable cluster method can be utilized to replace costly high-level coupled-cluster calculations without significant loss of accuracy, as demonstrated in the calculation of the triplet–quintet spin gap of the iron–porphyrin model system. In this case, the distinguishable cluster approximation reduces the error of CCSDT by a factor of 3.

It can be concluded that in majority of the studied systems and properties, DC-CCSDT showed large improvements compared to CCSDT. Thus, it can be used to benchmark the accuracy of other methods⁵⁷ or as part of composite methods for the calculation of ground and excited states.

See the supplementary material for absolute energies of the systems.

Financial support from the Max-Planck Society is acknowledged.

The data that support the findings of this study are available within the article and its supplementary material.