Perturbative universal state-selective correction for state-specific multi-reference coupled cluster methods

The single reference coupled cluster (SRCC) method^2,68 has emerged as one of the most successful post-Hartree-Fock methods for the inclusion of electron correlation effects. While the SRCC formalism is efficient in describing closed-shell systems, significant worsening in the quality of calculated energies can be observed when applying approximate SRCC methods to quasi-degenerate systems. One often encounters these problems in studies of homolytic bond breaking, low-spin open-shell and excited electronic states, diradicals, and transition states. The common feature of the corresponding quasi-degenerate electronic states is the fact that they can no longer be approximated by a single electronic state and multireference methods need to be invoked in order to properly describe complicated dynamic and nondynamic correlation effects in a balanced way. Among several classes of multi-reference coupled-cluster (MRCC) formalism (for a review of various MRCC formulations see Ref. 42), the Hilbert space MRCC formulations based on the Jeziorski-Monkhorst (JM) ansatz^{29,34–36,52,54–56,58} have assumed special position for its efficacy in encapsulating various types of correlation effects. The Hilbert space MRCC theory contains two major classes of methodologies: state-universal (SU)^{4,5,29,30,35,36,52,55,56,58} and state-specific (SS)^{7–9,13,16,19,20,26,43,45,46,48,57,60} MRCC formulations. Although driven by different design principles, both SU- and SS-MRCC formulations draw on the notion of model space (⁠$M 0$⁠) which contains Slater determinants necessary to provide a zero-th order description for electronic state/states of interest. The SU formalism, which furnishes as many electronic states as the dimensionality (M) of the model space, in many cases suffers from the intruder-state problem,^64,65 which significantly deteriorates the qualities of SU-MRCC energies and hence drastically limits the application area of the SU formulation based on the so-called complete model space. In addition to standard strategies based on the increasing excitation-rank in the cluster operators,^{3,7,8,19,21,38,39} several techniques have been introduced to address the intruder-state problem by extending SU-MRCC formulations to incomplete (or somewhere referred as general) model spaces.^{37,40,47,49,50,53}

In contrast to the SU-MRCC formulation, the SS-MRCC methods^26,45,48 target only one electronic state, which ameliorates the detrimental effect of intruder states. Over the last two decades, several forms of the equations for cluster amplitudes (or equivalently sufficiency conditions) have been discussed in the literature. The flexibility in defining sufficiency conditions stems from the known issue of the overcompleteness of the SS-MRCC theories. In the Brillouin-Wigner MRCC formulation,^{10,11,14,15,26,27,48,59–62} which is historically the first SS-MRCC approach, the JM ansatz is introduced directly into the Schrödinger equation, which results in a simple form of the equations but also leads to the appearance of the disconnected terms in cluster amplitudes and corresponding energies. The foundations of the Mk-MRCCSD formalism were first discussed in Refs. 44 and 45. However, in many instances one may stumble into severe numerical problems associated with numerical solving of the Mk-MRCC equations. Possible redundancies of the SU and SS formulations can be eliminated in the Hanrath’s MRexpT approach^22,23 by modifying the Jeziorski-Monkhorst ansatz and form of the equations for cluster amplitudes. A wide spectrum of possible choices for the sufficiency conditions has been extensively discussed by Kong in Ref. 31. Based on the numerical results obtained with various SS-MRCC formulations, one can notice that the SS-MRCC energy accuracies strongly depend on the particular choice of the sufficiency conditions. Moreover, it is a highly non-trivial task to infer the optimal form of the sufficiency conditions that would guarantee the most accurate energies.

The universal state selective (USS) approach^12,33 has been proposed to address this issue and possibly unify various SS-MRCC formulations driven by different forms of sufficiency conditions. By providing a simple mean to calculate energy corrections which, in principle, can recover exact full configuration interaction (FCI) energies, the USS approach can be applied to any type of the MRCC formulations employing the JM ansatz. Another attractive feature is that the USS energies are size extensive provided that the MRCC cluster amplitudes are expressed in terms of connected diagrams only. The recent numerical studies on several systems clearly showed that the USS correction is also capable of considerable reducing the size-inextensivity errors of the BW-MRCC energies.^6,12 In our recent paper,¹² we discussed the performance of USS and a posteriori corrections. The Mk-MRCCSD method, which is size-extensive, could still benefit from the inclusion of the higher-order contributions into the effective Hamiltonian. In particular, for each pair of references |Φ_μ〉 and |Φ_ν〉, which are at most doubly excited with respect to each other, the USS formalism at singles and doubles level includes also triply or quadruply excited terms with respect to the reference |Φ_μ〉. These contributions carry a part of the higher-order correlations effects and improve the accuracy of BW- and Mk-MRCCSD methods.

Along with non-iterative USS corrections, we have also developed and implemented an iterative version of the USS correction,⁶ which allows to obtain MRCC amplitudes consistent with the USS-MRCC Hamiltonian.

In this paper, we report on the parallel implementation of the perturbative non-iterative USS model employing singles and doubles (USS(2)).¹² This implementation was developed using the Tensor Contraction Engine (TCE)^24,25 available in NWChem.⁶⁷ As it will be discussed in Sec. II B, the USS(2) formalism requires an inclusion of the subset of triple and quadruple excitations with respect to each reference, which need to be evaluated and stored. Due to storage requirements of the full USS(2) formalism, we also developed a perturbative version of the USS(2) formalism where triply/quadruply excited contributions are evaluated in on-the-fly manner. The efficiency of the full USS(2) and perturbative USS(2) approaches in describing quasi-degenerate electronic states is evaluated by applying USS(2) corrections to the BW-MRCC and Mk-MRCC approaches with singles and doubles (BW-MRCCSD and Mk-MRCCSD methods, respectively). Using the BW-MRCCSD and Mk-MRCCSD formulations, we demonstrate the “universal” character of the USS(2) corrections in improving the quality of the BW-MRCCSD and Mk-MRCCSD correlation energies. We also demonstrate that the problems associated with the lack of invariance of the BW-MRCCSD and Mk-MRCCSD energies upon the rotations of active orbitals are significantly ameliorated by the USS(2) approach.

As benchmark systems/problems, we chose: CH₂ diradical, description of torsional barrier in ethylene, the estimation of the singlet-triplet separation in twisted ethylene, description of singlet and triplet states in tetramethyleneethane (TME), and singlet-triplet gap in α,3-dehydrotoluene.

The asymmetric USS functional³³ in its original form is given by the equation

where the $〈 Ψ ̄ ( μ ) |$’s are the reference specific trial wave functions and the c_μ’s are the normalized right eigenvectors of the approximate MRCC effective Hamiltonian. We note that in the limit of exact trial functions $〈 Ψ ̄ ( μ ) |$⁠, the above functional yields the exact energy.

In our earlier work,⁶ we found that it is advantageous to use a modified form of the USS functional

where the coefficients c_μ and $c μ ′$ obey a normalization condition

It is shown that this expression can be written in a compact form as

where

and the similarity-transformed Hamiltonian

has been introduced. Notice that the denominator term accounts for an arbitrary normalization of the $〈 Ψ ̃ ( μ ) |$ wave functions.

As discussed in our earlier works,^6,12 a reference specific ansatz of the trial wave function

or

is posit to simplify this USS functional. Here, $Λ ν ( μ )$ are the de-excitation operators usually truncated at the same excitation rank as T^(μ) operators. With this parametrization of the trial wave function, the USS functional takes the form

There are two major ways for obtaining the trial wave functions. One is to solve the left Bloch equations to get the $Λ ν ( μ )$ amplitudes. The second way is to approximate the Λ^(μ) operator simply as

Note that in both cases the $Λ ν ( μ )$ are tied to a specific form of the sufficiency conditions and depends on particular form of the T(μ) operators. A detailed analysis of both approaches can be found in Refs. 6 and 12.

In our implementation, we take $c μ ′$ to be the left eigenvectors of the effective Hamiltonian

and for the $d ν ( μ )$ coefficients we assume

Assuming these normalization conditions, the USS functional takes a very simple form

Note that the first term in the square bracket is the matrix element of the MRCC effective Hamiltonian

The diagonalization of the $H μ ν eff$ matrix generates the MRCC energy as an eigenvalue

Thus, we get the equation for the a posteriori USS correction

Alternatively, one can get an iterative form of the USS correction⁶ by augmenting the MRCC effective Hamiltonian matrix with the USS correction term

This augmented effective Hamiltonian is diagonalized in each iteration until convergence is achieved.

In our earlier work,^6,12 we developed implementations of the non-iterative as well as iterative USS corrections in the NWChem code. These implementations employ restricted representation of the $H ̄ μ$ operator in the effective Hamiltonian operator (only scalar, one-, and two-body terms of $H ̄ μ$ are included), which can produce single and double excitations from the corresponding reference function. Thus, one can write, for example, the double excited references as $〈 Φ ν |=〈 ( Φ μ ) k l c d |$⁠, where k, l, c, d are the active indices corresponding to the internal excitation. For mono- and bi-excitations in T(ν), the required matrix element then corresponds to projections of the main term of MRCC amplitude equations to triply and quadruply excited bras, i.e.,

Therefore, in our implementation, we compute a subset of triple and quadruple projections of $H ̄ μ | Φ μ 〉$ and take care of re-indexing them with respect to |Φ_ν〉 Fermi vacuum before contracting them with amplitudes $t i j a b ( ν )$⁠. In this work, we restrict our complete model space (CMS) only to CMS(2,2), i.e., the model space corresponding to two active electrons in two active orbitals. Note that in this case, the required subset of triple projections $〈 ( Φ μ ) i k l a c d | H ̄ μ | Φ μ 〉$ contains at least one active hole and one active particle. Analogously, the quadruple projections $〈 ( Φ μ ) i j k l a b c d | H ̄ μ | Φ μ 〉$ must contain two active holes and two active particles.

The calculation of such triple and quadruple projections remains the main bottleneck for the efficient implementation of the USS correction. When working at the singles/doubles (SD) level, the triple projections can be written as

where F_N(μ) and V_N(μ) are the one- and two-electron parts of the Hamiltonian H in the normal product form (H_N(μ)) with respect to the |Φ_μ〉 reference, i.e.,

The quadruple projections can be obtained in a similar way from the single and double excitation operators of the corresponding reference function. In this paper, we employ perturbative arguments to approximate the triple and quadruple projections at the lowest order,⁶⁶ which are for the triple projections

and for the quadruple projections

The triple and quadruple projections from the μth reference are calculated, then re-indexed with respect to the νth reference function and contracted on-the-fly to avoid their storage. We will refer the above procedure to as perturbative variant of the USS(2) approach (or USS(pt) for short).

The computational cost of the USS(pt) algorithm is mainly determined by the computation of the subset of triple and quadruple projections of the similarity Hamiltonian $H ̄ μ$ acting on the μth reference. Because the subset of triple projections involves at least one pair of active indices and two pairs of active indices for quadruple projections, the USS(2) algorithm scales as O(N⁶) with a prefactor depending on the active space size and mutual excitation level of reference functions.

The full USS (further denoted as USS(full)) and perturbative USS(pt) corrections have been implemented in the development version of NWChem and we have performed several benchmark calculations on following systems: (1) CH₂ diradical at linear and bent geometry, (2) ethylene with an arbitrary torsional angle, (3) the tetramethyleneethane (TME) system, (4) H₂O molecule, and (5) α,3-dehydrotoluene biradical. In this work, we study the performance of the USS(pt) correction and compare USS(pt) results with the ones obtained with the USS(full) or uncorrected MRCCSD methods. We study the invariance of the USS energies with respect to rotations in the active space using the example of twisted ethylene. We also investigated the differences in the description of the triplet state using USS-corrected MRCCSD methods (M_s = 0 component) and by the single-reference CCSD method (M_s = 1 component).

For the CH₂ diradical, the geometries were taken from Ref. 16. In the above reference, the geometry optimizations were performed at the MkCCSDT level using the cc-pVTZ basis set.¹⁸ In the CH₂ diradical coupled cluster calculations, we employed cc-pVTZ, cc-pVQZ and cc-pV5Z basis sets¹⁸ and kept core 1s orbital of the carbon atom frozen in all CC calculations. For the twisted ethylene model, the geometry was optimized at the CASPT2(2,2) level. In studies of active space rotation invariance, we used the cc-pVTZ basis set. In the calculations of the triplet-state energy difference between MRCCSD energy (M_s = 0 component) and single-reference CCSD energy (M_s = 1 component), we employed cc-pVDZ and cc-pVTZ basis sets. The following calculations of the torsion barrier and singlet-triplet gap were performed using four consecutive correlation consistent basis sets (from cc-pVDZ to cc-pV5Z). In this case, the geometry was optimized at the CASPT2 level for each torsional angle and basis set. In all calculations of ethylene, we kept core orbitals frozen. The TME system was optimized for each torsional angle using CASPT2(6,6) and cc-pVTZ basis set. The SRCC and MRCC/CMS(2,2) calculations were performed using CASSCF(6,6) orbitals. The α,3-dehydrotoluene geometry was optimized using CASPT2 in the cc-pVTZ basis set. We used 8 active electrons in 8 active orbitals to define the model space for the CASPT2 geometry optimization. For the study of singlet-triplet gap in the TME system, we used a spherical cc-pVTZ basis set, while for the α,3-dehydrotoluene system we employed cc-pVDZ and cc-pVTZ basis sets. The spherical representation of the correlation consistent basis sets was used. All CASPT2 calculations were carried out using MOLPRO package.¹ Optimized geometries of all systems can be found in the supplementary material.

We now illustrate the accuracy of the USS(pt) correction with respect to the original USS(full) formulation or FCI, if available. In Table I, we present the comparison between the USS-corrected BW-MRCCSD and Mk-MRCCSD results against FCI,¹⁶ MRCCSD(T), CASPT2, CCSD(T), and CCSDT values for both linear and bent CH₂ structures using the spherical cc-pVTZ basis set. The difference between BW-MRCCSD+USS(pt) and FCI energy is 4.1 mhartree for the bent geometry, whereas the difference for the linear structure is 2.2 mhartree. We have observed similar performance also for the USS-corrected Mk-MRCCSD methods. The best agreement with the FCI energy gives the USSD correction, which might be due to a fortuitous error cancellation. The difference between the approximate USS(pt) corrections and full USS corrections is negligible in all cases. In terms of accuracy, the USS-corrected methods can be placed between uncorrected MRCCSD and MRCCSD(T) approaches. In Table II, we present the energy difference between BW- and Mk-methods for the CH₂ system computed in different basis sets. Clearly, the addition of USS corrections decreases differences between BW-MRCC and Mk-MRCC results. The maximum USS(pt) deviation is only 0.019 mhartree for the linear geometry in cc-pVTZ basis, while the difference between BW-MRCCSD and Mk-MRCCSD energies is 0.445 mhartree. The difference is larger for the less multireference system with the bent geometry, but still significantly smaller than the difference between uncorrected MRCCSD methods.

While single-reference CC methods provide a good performance for bent geometry, for the multireference case they require perturbative/iterative inclusion of full triples in order to achieve similar level of accuracy. For linear and bent geometries, the CASPT2 method shows large deviation from the FCI results even when larger model spaces are employed.

The ethylene system has been extensively used for testing of multireference methods due to the bi-radical nature of the twisted conformation. The torsion barrier is defined as the energy difference between twisted and planar geometries. The experimental value reported from the kinetics study of the thermal cis-trans isomerization is 65 kcal/mol,¹⁷ whereas parametrization from the resonance Raman spectra predicts 60 kcal/mol.⁶⁹ In Table III, we summarized energies obtained using various basis sets with Mk-MRCCSD and BW-MRCCSD methods. The planar ethylene is described using the single-reference CCSD method, while in parentheses we provide the barrier energy when MRCCSD methods with the CMS(2,2) model space were used. The twisted ethylene requires multireference approach with the CMS(2,2). Clearly, BW-MRCCSD and Mk-MRCCSD methods without the USS correction predict higher torsion barrier than the experiments. The USS-corrected energies in Table III are lower, between 61-64 kcal/mol, which is in a good agreement with both experiments. Using MRCC methods in calculating energies of the planar system leads to a small increase in the values of calculated barrier energy (by around 1-2 kcal/mol) and effectively to a good agreement with the experimental data of Ref. 17. Also, Mk-MRCCSD and BW-MRCCSD USS-corrected energies for the same basis set are almost identical. The barrier energy does not change much with the basis set size, USS(pt) and USS(full) method with cc-pVTZ, and higher basis set predicts the energy around 63 kcal/mol, while USSD energies are about 1 kcal/mol larger.

The torsion energy profile is shown in Fig. 1, where cc-pVTZ basis set was utilized. The uncorrected MRCC energy is taken as the zero reference, so we can directly compare the size of corrections (and SR-CCSD energy) with respect to the torsional angle. The weight of the second (doubly excited) reference varies between 1.5% for the planar geometry and 50% for orthogonally twisted ethylene. While the BW-MRCCSD energy almost blend into the CCSD energy in the single-reference region, the Mk-MRCCSD energy is invariably about 2 mhartree lower. This difference is compensated by the USS correction, which is about 2 mhartree smaller than for the BW-MRCCSD method, so both USS-corrected BW- and Mk-MRCCSD energies are almost equal. For rotational angle greater than 50°, where the multireference effects start to prevail, the USS corrections grow to approximately 10 mhartree for the BW-MRCCSD method and 7 mhartree for the Mk-MRCCSD approach, respectively. The a posteriori correction to BW-MRCCSD performs similarly to USSD, but for the strongly correlated case (90°) the a posteriori correction is for about 4 mhartree smaller than USS(pt) and USS(full) corrections. Similarly to the torsion barrier study, we also studied the singlet-triplet gap for the twisted geometry of ethylene using the same methods and basis sets. In Table IV, we present S-T energies in eV. While the uncorrected Mk-MRCCSD and BW-MRCCSD methods predict triplet as the ground state, USS corrections push the singlet state about 0.1 eV below the triplet state. As observed previously, there is almost no difference between USS corrected energies of both methods. The largest discrepancy 0.004 eV was found for the cc-pVDZ basis set and USS(pt) or USS(full) correction, which is a negligible value. In Table V, we report the triplet state energies for twisted ethylene calculated using single-reference CCSD approach and MRCCSD methods. In contrast to the triplet state, which is well described by using ROHF M_s = 1 reference determinant in the CCSD calculations, singlet state requires multireference treatment. For uncorrected MRCCSD methods, the difference between CCSD and MRCC energies is approximately 6 mhartree for the Mk-MRCCSD approach and almost 10 mhartree for the BW-MRCCSD method. The energy difference significantly decreases when USS corrections are included and falls to less than 3 mhartree for the USS(full) and USS(pt) formulations. For the USSD formulations, this difference is even less than 0.5 mhartree.

In our next example, we used USS(pt) and USS(full) approaches to evaluate the singlet-triplet energy gap of the TME molecule at different torsional angles. This system poses a significant challenge for electronic structure methodologies due to a combination of strong static and dynamical correlation effects;⁶³ in this paper, we focused only on the comparison of the BW-MRCC and Mk-MRCC methods with or without USS corrections. In our calculations, we utilized spherical cc-pVTZ basis set. To obtain triplet state energies, we employed the single-reference CCSD method. As seen from Fig. 2(a), which compares the twisting potentials described by the BW- and Mk-MRCCSD methods, the uncorrected Mk-MRCCSD method fails to predict the singlet state to be the lowest one. This is a good illustration of the fact that for this system even multireference CC methods can stumble into serious problems with proper inclusion of correlation effects. Adding USS correction (USSD or USS(pt)) to the Mk-MRCCSD energy results in a lowering singlet state energy below the energy of the triplet state. Also, irrespective of the fact how different are the uncorrected MRCCSD energies, the USS-corrected energies do not differ by more than 0.07 eV. This is a quite remarkable feature of the USS approach. For the S-T gap, which is shown in Fig. 2(b), the USS methods correctly predict the minimum gap at the angle 45°.

For the H₂O benchmark system, we tested the accuracy of the USS approach for larger CAS(4,4) model space. In Table VI, we report BW- and Mk-MRCCSD energies against the FCI results for the water molecule at various geometries (corresponding to symmetric stretching of the H–O bonds). For comparison, we also provide CASPT2 energies. The equilibrium geometry of H₂O was defined by the H–O distance R_e = 1.843 45 bohrs and the H–O–H angle α = 110.6°. Moreover, the spherical representation of the cc-pVDZ basis set was employed in calculations. In all cases considered here, USSD and USS(pt) corrections improve the accuracy of uncorrected MRCCSD methods. For the equilibrium geometry (R = R_e), the absolute USS(pt) and USSD errors with respect to the FCI energy are only slightly smaller than those obtained with the uncorrected methods (BW-MRCCSD and Mk-MRCCSD). For larger distances (R = 1.5R_e and R = 2R_e), USS errors are significantly smaller for both CMS(2,2) and CMS(4,4). Despite the fact that the USS(pt)-corrected BW-MRCCSD energy for 2R_e for CMS(2,2) model space shows a quite large deviation from the FCI energy (12.861 mhartree) one should stress that this error is still much smaller than the error corresponding to the uncorrected BW-MRCCSD method (26 mhartree). This big discrepancy is eliminated when using a larger CAS(4,4) model space, where USS-corrected energies are much closer to the FCI energy than the uncorrected ones. The USSD/USS(pt) Mk-MRCCSD energies are consistently better than the corresponding Mk-MRCCSD energies. The CASPT2 method employing CAS(2,2) and CAS(4,4) yields energies that are characterized by significantly larger deviations from FCI results than the corresponding USS energies. To illustrate this fact, one should collate 42 mhartree of CASPT2/CAS(2,2) energy error and 6 mhartree error of the Mk-MRCCSD USS(pt)/CMS(2,2)approach. The CASPT2 method requires much larger active space (CAS(10,12)) in order to achieve the accuracy of the USS-corrected MRCCSD methods obtained with smaller model spaces.

The approaches based on the Jeziorski-Monkhorst ansatz are not invariant with respect to rotations of active orbitals.³² In particular, it is very important to qualitatively estimate how much active orbital rotations impact the USS energies. To address this question we have performed two benchmark calculations for the twisted ethylene using the cc-pVTZ basis set, where two active orbitals were continuously rotated/mixed using rotation angles falling into the [0°,45°] interval. Results for the singlet state are shown in Fig. 3, while the energies for the triplet state are shown in Fig. 4. In all cases, the uncorrected MRCCSD energy at rotational angle 45° is taken as the zero-energy reference. As seen from Fig. 3(a), the USS correction significantly improves the invariance of the BW-MRCCSD method, where 9.2 mhartree deviation (calculated as a difference between the maximum and minimum values of the energy as a function of rotation angle α) is reduced to only 2.4 mhartree when using USS(full), or to 2.8 mhartree when using USS(pt) correction. Smaller discrepancies are the consequence of the fact that the USS(full) and USS(pt) curves are significantly “flatter” which indicates that these methods are less sensitive to the orbital rotation than the original BW-MRCCSD formulation. For comparison, the USSD energy shows a minimum around 22° and the largest discrepancy is 4 mhartree. The a posteriori correction gives a similar non-monotonic shape to USSD correction with the largest difference 3.4 mhartree. As shown in Fig. 3(b), the shape of Mk-MRCCSD is very similar to the BW-MRCCSD case and demonstrates a significant effect of mixing active orbitals on the resulting energies. The Mk-MRCCSD energy discrepancy is 5.4 mhartree. As in the BW-MRCCSD formalism, the USS(full) and USS(pt) decrease the value approximately to 2.5 eV and 2.9 eV, respectively. The triplet state energies at the MRCCSD level are shown in Fig. 4. In analogy to the singlet energies, the USS energies are virtually flat. It is worth mentioning that for the triplet state the USSD is flat almost to the plateau. The largest discrepancy for BW-MRCCSD USS(full) is 2.2 eV, 2.5 eV for USS(pt), and 0.2 eV for USSD. The shape of USS-corrected Mk-MRCCSD energy curves is practically identical to the BW-MRCCSD case.

Our next benchmark system is the α,3-dehydrotoluene diradical. Experimentally, the ground state of this diradical was found to be the open-shell singlet,^41,51,70 which is in contrast to α,2- or α,4-isomers with the triplet ground state. The experimental data indicate that the S-T gap is smaller than 5 kcal/mol. In Table VII we compare the USS-corrected and uncorrected BW-MRCCSD and Mk-MRCCSD energies for the singlet ground state of this diradical. For comparison, we also provide CASSCF and CASPT2 results for different active space sizes. In the first benchmark, we used canonical RHF orbitals, while in the second set, which is much more challenging, we rotated two active orbitals for 45° in order to show the non-invariance to active orbital rotations. Due to the numerical cost of the full USS corrections, we calculated the USS(full) energies only for the cc-pVDZ basis set.

As shown in Table VII, we find that the differences between the USS(pt) and USS(full) corrected MRCCSD energies (BW and Mk) for HF orbitals in the cc-pVDZ basis set are very small and amount to 0.1 kcal/mol, while for rotated active orbitals the analogous differences reach 0.7 kcal/mol (Mk) and 0.5 kcal/mol (BW). For both (cc-pVDZ and cc-pVTZ) basis sets, the differences between HF orbitals based BW-MRCCSD USS(pt) and Mk-MRCCSD USS(pt) results are negligible (less than 0.1 kcal/mol), whereas for rotated orbitals these differences grow to 0.7 kcal/mol and 0.8 kcal/mol for the Mk-MRCCSD and BW-MRCCSD USS(pt) formulations, respectively. However, this is much less than 14.1 kcal/mol (cc-pVDZ) or 11.9 kcal/mol (cc-pVTZ) difference between uncorrected MRCCSD methods (Mk-MRCCSD-BW-MRCCSD) for α = 45. The effect of orbitals rotation on the singlet-triple separation is especially strong for the BW-MRCCSD approach. While for α = 0 the singlet-triplet splitting obtained with uncorrected BW-MRCCSD approach in cc-pVDZ basis set corresponds to 2.8 kcal/mol, for α = 45 the analogous singlet-triplet value reaches −15.5 kcal/mol. It is worth mentioning that despite of these huge differences, the USS(pt) approach can “stabilize” the BW-MRCCSD formalism and yield reasonable results of 2.8 and 4.9 kcal/mol for α = 0 and α = 45, respectively. These results provide a very good illustration of the universality of the USS corrections. The singlet-triplet MRCCSD/cc-pVTZ energy with the USS(pt) correction is approximately 2.6 kcal/mol, which is in agreement with the experimental observations and also with the values reported in Refs. 71 and 28, where calculations were performed at the CI or MCSCF(8,8) level using 6-31G** or cc-pVDZ basis set. The CASSCF(2,2) S-T energy corresponding to our CASPT2/cc-pVTZ optimized geometry is −0.5 kcal/mol, using the larger model space (4,4) or (8,8) the energy increases to 3 kcal/mol. The CASPT2 energies behave more consistently with the size of the model space compared to the CASSCF method. The CASPT2(8,8) gap is approximately by 1 kcal/mol smaller than the corresponding USS-MRCCSD results. In contrast to BW- and Mk-MRCCSD approaches, for rotated orbitals these methods take advantage of the invariance to active orbital rotations.

Concerning the computational cost, the perturbative USS approach has two advantages. First, the required computational time for USS(pt) is much smaller than for the full USS calculation. Second, since the USS(pt) approach calculates the triple and quadruple projections in on-the-fly manner it eliminates the global memory bottleneck of the USS(full) associated with the need of storing expensive recursive intermediates.

In Fig. 5, we report timings of the USS algorithms for α,3-dehydrotoluene. All calculations were performed using NWChem on the Cascade system, which is equipped with many Xeon E5-2670 8C 2.6 GHz 16-core CPUs and a Infiniband FDR network, and maintained at the Pacific Northwest National Laboratory’s EMSL facility. For example, calculation of the BW-MRCCSD USS(pt) energy on 1024 cores takes 441 s compared to 6762 s for the BW-MRCCSD USS(full) calculation, which means that the speed up factor is bigger than 15 times. It is interesting to notice that the USS(full) correction requires more than 6400 cores in order to achieve the same timing as USS(pt) calculation on 512 cores. The scalability of the USS(pt) code on more than 1000 cores is still good even the number of tasks is much smaller than in the USS(full) calculation.

In this work, we proposed and implemented an efficient perturbative approximation of the USS correction for the BW-MRCCSD and Mk-MRCCSD methods. We demonstrated on several numerical examples that the perturbative USS approximation substantially reduces the computational cost of the full USS correction without any significant effect on the accuracy of resulting energies, compared to the full USS algorithm. For example, calculations of the S-T energy gap for twisted ethylene using the cc-pVTZ basis set demonstrated that the energy difference between USS(pt) and USS(full) energies is as small as 0.3 mhartree. For larger system, such as the α,3-dehydrotoluene, the analogous difference corresponds to only 0.16 mhartree. For all small benchmark systems studied in this paper, the discrepancies between USS-corrected BW-MRCCSD and Mk-MRCCSD energies are very small. We have observed that these differences for the α,3-dehydrotoluene are much smaller compared to uncorrected MRCCSD energies. It has also been shown that the USS approach provides more accurate results than standard a posteriori correction added to the uncorrected BW-MRCCSD and Mk-MRCCSD energies. Our studies provide a strong argument in favor of the universality of the USS formalism. We have also shown that for both BW- and Mk-MRCC methods, the USS corrections significantly suppress the lack of invariance to rotations in the active space. On several examples, we demonstrated that the USS energies are significantly less sensitive to the rotations of active orbitals than the iterative BW-/Mk-MRCCSD formulations. The performance tests show a good scalability of the USS(pt) implementation and a significant speed up compared to the original USS(full) implementation, while preserving almost full accuracy.

See the supplementary material for geometries and energies.

We acknowledge the support by the Ministry of education, youth and sports of the Czech Republic (Project No. LH13117) and by the Czech Science Foundation (Project No. 208/11/2222). JB acknowledges the support by the Czech Science Foundation (Project No. 15-00058Y). This work was supported by the Mi nistry of Education, Youth and Sports from the Large Infrastructures for Research, Experimental Development and Innovations project “IT4Innovations National Supercomputing Center - LM2015070.” A large portion of calculations has been performed using EMSL, a national scientific user facility sponsored by the Department of Energy’s Office of Biological and Environmental Research and located at Pacific Northwest National Laboratory. The Pacific Northwest National Laboratory is operated for the U.S. Department of Energy by the Battelle Memorial Institute under Contract No. DE-AC06-76RLO-1830.