Single-Root Multireference Brillouin-Wigner Perturbative Approach to Excitation Energies.

The state-specific Brillouin-Wigner multireference perturbation theory [which employs Jeziorski-Monkhorst parametrization of the wave function] using improved virtual orbitals, denoted as IVO-BWMRPT, is applied to calculate excitation energies (EEs) for methylene, ethylene, trimethylenemethane, and benzyne systems exhibiting various degrees of diradical character. In IVO-BWMRPT, all of the parameters appearing in the wave function ansatz are optimized for a specific electronic state. For these systems, the IVO-BWMRPT method provides EEs that are in close agreement with the benchmark results and experiments, where available, indicating that the method does not introduce imbalance in the target-specific treatment of closed- and open-shell states involved. The good performance of the present methodology is primarily related to structural compactness of the formalism. Overall, present findings are encouraging for both further development of the approach and chemical applications on the energy differences of strongly correlated systems.

Introduction

People working in the field of chemistry and spectroscopy are often interested in energy differences (say, excitation energies, EEs) that are usually orders of magnitude smaller than total state energies. Despite various methodological developments, the proper understanding and accurate interpretation of electronic EEs remains one of the most compelling challenges for the modern electronic structure theory because a correct and balanced treatment is needed for the electronic states of interest, which are often plagued by quasidegeneracy.[1] Nonetheless, electronic quasidegeneracy, which is ubiquitous in chemistry (such as electronically excited and open-shell species and molecules with elongated bonds), poses a challenge to the ab initio method as it leads to multiconfigurational (MC) or multireference (MR) wave functions. The methodologies employed for the estimations of EEs can be broadly classified under two categories: (i) approaches formulated for direct energy differences rather than for individual state energies and (ii) methods that entail the explicit evaluation of individual state energies followed by the estimation of the EE as an energy difference of the states involved. Within the coupled cluster (CC) theory, the equation-of-motion (EOM)[2,3] and linear response theory (LRT)[4,5] [and symmetry-adapted cluster-configuration interaction (SAC-CI)[6]] methods, as well as the multireference CC scheme defined in the Fock space (FSCC),[7,8] are widely used computational protocols aimed at the direct calculations of energy differences. The second class of schemes is exemplified by MR configuration interaction (MRCI)[9] and MR perturbation theory (MRPT)[10−12] methods where small errors in total state energies may invite very large errors in the target EEs. However, it is often stated that methods tailored to supply energy differences directly are potentially more accurate than (approximate) schemes of similar complexity designed for the states’ total energies. However, it should be noted that this potential can be satisfied only if the realization assures a balanced description of the states of interest. Although each of the above mentioned methods offers certain merits, the most appealing feature of MRPT methods is their computational cost effectiveness (yet capturing the essential effects of electronic correlation), and hence these protocols are popular quantum chemistry approaches for studying quasidegenerate electronic states.

Over the decades, a large number of different variants of MRPT methods have been proposed and used in the literature,[13−61] among which the Hilbert-space state-specific MRPT (SSMRPT) method (where only one state is studied at a time) introduced more recently by Mukherjee and co-workers[37,42] and further developed by Chattopadhyay and co-workers[39,44,46,47,49] is growing in popularity because of its rigorous size extensivity and explicit insensitivity toward intruder states. The most unsettling issues of several extensively practiced MRPT methods are the well-known intruder-state problem (states within the framework of the perturbative protocol corresponding to excited configurations with vanishing energy denominators) and lack of size extensivity. One can use level shifting of the energy denominators (that adds a small constant to the zeroth-order Hamiltonian in an ad hoc manner)[62−64] to remove or at least to attenuate the intruder problems; however, level shifting can affect the state ordering.[64]

In MRPT, the perturbation correction is applied after constructing the reference space wave function. Regarding the construction of the reference function, the complete active space self-consistent field (CASSCF) protocol is the most widely used, which requires the simultaneous optimization of the CI and orbital coefficients to self-consistency. This nonlinear optimization process often invites multiple solutions, convergence difficulties, and unphysical spatial symmetry breaking of the wave function. The potential difficulties of the CASSCF method can be reduced by employing the complete active space configuration interaction (CASCI) approach in conjunction with useful expansions of the low-lying excited-state wave functions of molecules based on improved virtual orbitals (IVOs).[65] The resulting method, IVO-CASCI,[65,66] is cost-effective since it does not require additional costly iterations beyond an initial SCF calculation. Recently, a new version of the SSMRPT method with the IVO-CASCI reference function within the framework of Brillouin–Wigner (BW) perturbation (called IVO-BWMRPT)[50] has been developed by our group. The IVO-BWMRPT is endowed with all of the redeeming formal properties of the parent method.[37] The elimination of the nonlinear parametrization in the construction of the reference wave function leads to a reduction in computational overhead. The IVO-BWMRPT ansatz is sufficiently flexible to describe changes in target-state wave functions along the bond-breaking coordinates in the ground and excited states.[50,51] It also furnishes competitive estimates, compared to various alternative MRPT approaches.[50,51] The focus of the present work is to investigate this MRPT scheme further to establish its robustness and efficacy for providing excitation energies [calculated as the difference between the energies of specific states] of methylene, ethylene, trimethylenemethane (TMM), and benzyne isomers. As the IVO-BWMRPT method treats one state at a time through a wave operator that lifts only one state, description of excited states using this approach to compute the energy difference is an intriguing challenge. To study the EEs employing the state-specific/single-root MR approaches, the perturbation corrections are separately estimated for each reference wave function corresponding to the state of focus. Apart from the ground state, the target excited state can be accessed through identification of the appropriate eigenvalue of the effective Hamiltonian (which provides several roots in a single calculation) using a suitable root homing scheme in IVO-BWMRPT calculations. Note that the ground-state IVO-BWMRPT shows a faster convergence in comparison to the excited one. Krylov et al.[67−71] developed spin–flip (SF)-based schemes, designed particularly for radical systems, of which the above systems are prominent representatives.

The structure of the remaining part of the paper is as follows. The IVO-BWMRPT method is briefly reviewed in Section . Section presents several test cases and illustrates that the method has a broad range of potential applications.

Synopsis of the IVO-BWMRPT Formalism: A Hilbert-Space MRPT Method

Apart from the convergence and the existence of multiple solution problems, there is no straightforward way to truncate virtual orbitals in the CASSCF procedure, as the ordering of the orbitals is not necessarily linked with their importance that is carried in describing the system. The IVO-CASCI scheme of Freed and co-workers[65] has the ability to tackle these difficulties. Therefore, the main question is how to obtain IVO-CASCI orbitals suitable for ground and excited states without doing an expensive computation. To construct IVO-CASCI wave functions, one can start with (i) a single restricted Hartree–Fock (RHF) or restricted open-shell HF (ROHF) computation for the ground state (or some other reference state), followed by (ii) the diagonalization of a single Fock operator to generate IVOs, and then (iii) a CASCI to capture the MR/MC character of wave functions. One portion of both the IVO-CASCI and CASSCF approaches effectively associates a CASCI calculation, where the CAS comprises all near-degenerate configurations. No iterations are needed in either step (ii) or (iii) [thus, it does not possess properties that invite convergence problems with increasing size of the CAS], and different CASs can be treated without repeating steps (i) and (ii). As the IVO-CASCI wave function consists of both singly and doubly excited configurations in the reference space on the top of higher-order excitations, the method has the ability of treating both singly and doubly excited states on an equal footing, advocating that the IVO-CASCI method is a viable scheme for providing consistent description of multiple electronic states simultaneously as that from the computationally expensive CASSCF one. The IVO-CASCI wave functions satisfy Brillouin’s condition.[65] IVO-CASCI (CASSCF) often yields errors in various processes or situations because of the inadequate description of the dynamic correlation effects. The dynamical correlation effects can be incorporated to a good extent using the BW version of the SSMRPT method (BWMRPT). Note that here orbitals from a state-specific IVO-CASCI are used for each state separately in BWMRPT calculations.

In the BWMRPT method, one can solve perturbatively the Schrödinger equation for a single target state, rather than the Bloch or Lippmann–Schwinger-like equation. The IVO-BWMRPT formalism is a single-root MR method based on the Jeziorski–Monkhorst (JM) ansatz[72] where the exact wave function is written as |Ψ⟩ = ∑μexp(Tμ)|ϕμcμ⟩ (for the sake of simplicity, one can omit the target-state indices). Here, each reference function ϕμ yields the various virtual functions χ due the action of reference-specific cluster operator Tμ (i.e., each reference configuration has its own set of independent amplitudes) and cμ describes the reference space coefficients (not known a priori). Here, the mismatching between the number of wave function parameters and conditions that can emerge from projections of the Schrödinger equation (often termed as a problem of redundant parametrization) has been resolved (to ensure a unique solution) using suitable additional sufficiency conditions.[37,42] It is important to stress that redundancy in the number of cluster amplitudes emerges from the fact that a given excited function χ is created by Tμ acting on ϕμ and by Tν acting on ϕν, ϕμ ≠ ϕν. Within the active or model space, the equation giving the desired energy can be written in a matrix form of the effective Hamiltonian operator (constructed up to order 2) as (by relying on the intermediate normalization)where H̃μν(2) = Hμν + ∑Hμtν(ν). Here, ⟨χμ|Tν|ϕμ⟩ = tμ(ν), ⟨ϕμ|H|ϕν⟩ = Hμν, and ⟨χμ|H|ϕμ⟩ = H. The dynamical correlation effects are “wrapped-in” through the similarity transformation. The eigenvalue equation of H̃μν(2) provides the coefficients and energy of the target at order 2.

The working equation for computing the cluster amplitudes required to solve eq can be given bywhere the superscript “0” represents that the quantities are unperturbed. It should be underlined that the cluster amplitudes, tμ(μ) = ⟨χμ|Tμ|ϕμ⟩, depend on the unknown second-order perturbed state energy, which is characteristic of the BW theory. The intermediate normalization in the theory is maintained by enforcing the conditions (for a complete model space) ⟨ϕν|Tμ|ϕμ⟩ = 0 and ∑μ=1cμ2 = 1, where M describes the dimension of the model space. One can also see that eq involves the effective Hamiltonian matrix elements (H̃μν(2)), target-state energy, and the reference space coefficients for a desired state, indicating the state specificity of the resulting method. It is important to underline that only those cluster amplitudes are coupled together that bear a common set of orbitals, implying computational advantage of this approach. The denominators of the above equation depend on the unknown exact energy eigenvalue E, which is one of the distinguishing features of the BW scheme, in contrast to the Rayleigh–Schrödinger version.[37]Equations and 2 are the principle working equations and need to be solved self-consistently. The effective Hamiltonian must be diagonalized in each iteration and, therefore, the reference space is modified in each iteration. The iteration is initiated by diagonalizing the Hamiltonian matrix in the reference space. This feature of the method allows the updating of the CAS coefficients as a result of mixing of the virtual functions. As a result, IVO-BWMRPT follows a “diagonalize–perturb–diagonalize” philosophy, which permits coupling of the dynamical and nondynamical correlation effects that can be crucial for accurate treatment of EEs, description of energy surfaces in close proximity, the curve crossings, and electronically excited states of mixed valence with a Rydberg character. As the method is single root, only the selected eigenvalue of the effective Hamiltonian, for which we converge the equations in a self-consistent manner, has a physical meaning. Other roots, along with perturbing the model space Hamiltonian in the manner mentioned above, serve for “buffering” and getting rid of the intruder states by shifting them away and, finally, preventing the occurrence of divergence due to intruder states (even when considering a wide range of geometries of the nuclear framework). One can view this issue from another angle. The uncorrelated [E0,μμ – E0,] part is shifted by the correlated one [E(2) – H̃μμ(2)] in the energy denominator of the cluster finding IVO-BWMRPT equation (eq ), a reason why the theory is manifestly robust against intruders (without any need of resorting to the use of unphysical level shifts common to most usually used MRPT methods) irrespective of whether the target state is ground or excited as long as the target-state energy is well separated from the energy of virtual functions. Here, the zero-order energies, E0, are the expectation values of H0 with respect to a reference space configuration. Although it is not restricted to the use of a specific unperturbed Hamiltonian, the multipartitioning Möller–Plesset (MP) technique of Zaitsevskii and Malrieu[73] is adopted here to define the zero-order Hamiltonian, H0 (where a separate H0 is used for each reference space state), in the IVO-BWMRPT method. Note that the multipartitioning protocol is usually more flexible and efficient than MRPT schemes based on the generalized Bloch-like equation.[56] The unitary group approach (UGA) has been used to arrive at a spin-pure formulation, where the cluster operators are spin-free via the use of unitary generators. The details of discussion about H0 and the spin adaptation of the SSMRPT scheme can be found in previous publications[39,42] and are not discussed here. In the context of spin adaptation, the recent proposal of Jeszenszki et al.[43] is also worth mentioning. The present method, IVO-BWMRPT, is applicable to electronic states of closed-shell or open-shell configurations. It treats all configurations in a democratic manner, i.e., without overemphasizing the importance of any configuration, which is very crucial to estimate EEs. The connected nature of the IVO-BWMRPT working equations leads to a strict size-extensive nature of its energy corrections. The IVO-BWMRPT wave functions are size-consistent with localized orbitals. On the other hand, the state-of-the-art internally contracted MRPT methods are not strictly size-extensive.[13,15] Highly popular and successful MRPT methods (such as MRMP2, CASPT2, and NEVPT2[13,15,21]) belong to the family of internally contracted (unrelaxed) protocol (also called the “diagonalize–perturb” scheme). Note that at the CASSCF/CASCI level, if the mixing between valence and Rydberg states occurs, application of an internally contracted single-state perturbative correction, leaving the component of the reference wave function unchanged, is unaccepted. Moreover, it often occurs that two or more electronic states, which are not interacting at the zero-order level, become very close in energy after the incorporation of the perturbative correction in such a manner that a mixing is liable to happen. Note that both these issues play an important role in the accurate estimation of excitation energies using single-root MRPT methods.[74] It is, however, well-accepted that the treatment via multistate (MS) versions of these MRPTs[56−59] is adequate for studying near-degenerate states, in which a strong mixing emerges between the reference and other secondary states via the second-order perturbative interactions. The MS-MRPT builds an effective Hamiltonian matrix in the reference state basis under the influence of the perturbation, and its diagonalization furnishes the mingling of the perturbed states. In this context, development of generalized Van Vleck perturbation theory (GVVPT2) is worth mentioning[27] (a subspace-specific MRPT), in which the perturbed primary functions and the unperturbed secondary functions form a basis of the effective Hamiltonian.

The weight factor in the coupling term in eq merits some discussion. Although one can notice that the IVO-BWMRPT ansatz satisfies a number of important requirements, a very small cμ compared to cν may result in rather large values for cluster amplitudes tμ(μ), which in turn give rise to numerical instability (that is common to state-specific JM-based methods) to H̃μν(2). Then, the solution of IVO-BWMRPT equations appears to achieve convergence to the correct root with difficulty. This phenomenon is different from the intruder-state effects as it is not linked with a vanishing energy denominator in the perturbative expansion rather with a vanishing CAS coefficient(s). This unwanted issue becomes more acute with enhancing dimensions of the reference space because many cμ’s are close to zero at one zone of the energy surface or another. At this point, it should be noted that the contribution from a configuration state function (CSF) with small reference space coefficients is usually notably small. Dropping the cluster amplitudes stemming from a ϕμ with coefficients cμ below a given threshold (combining coefficients less than 10–6) is the simplest solution to get rid of this problem.[75] Since the weight of the excitations out of the dropped configurations is negligibly small, there is a minimal loss of accuracy in the computed energy of the target state. A more sophisticated solution to avoid or at least to attenuate this problem is to use the Tikhonov method of regularization (which preserves all amplitudes).[42] To cure this numerical instability of the IVO-BWMRPT scheme, one may organize the working equations in such a manner that the division by cμ’s is also minimized to the extent possible. Use of an incomplete model space (i.e., shrinking the size of the model space by removing all configurations of low importance) in conjunction with the so-called constraining (C) conditions [on cluster amplitudes that are associated with the internal excitations transforming one reference configuration into another one] originally proposed by Li and Paldus[76] might be a useful means of avoiding the numerical instability.

Results and Discussion

Computations of the EEs of molecular systems of various sizes and characters, namely, the methylene radical, ethylene, TMM, and benzynes, have often been used to calibrate new theoretical methods tailored to compute energies of excited and ground states as the wave function changes character if one moves from the ground to excited states and different types of correlation effects can dominate across the low-lying states. Thus, several considerations should be taken into account when calibrating ab initio strategies. It is useful to gauge the performance of IVO-BWMRPT for evaluating energy differences such as vertical EEs (VEEs) as well as adiabatic EEs (AEEs). It is not the AEEs but the VEEs that are generally reported. The dependence of adiabatic excitation energies on equilibrium geometries of states involved invites additional uncertainties. It is often stated that the VEEs are more useful for benchmarking, since they do not depend on the geometrical parameters of excited states; rather, they are estimated at the ground-state equilibrium geometry. It should be noted that the selection of a common-to-all-states’ and yet compact active space is not a trivial task to treat ground and excited states of different characters for an MR method designed to handle the electronic state in a target-specific manner. The choice of the CAS depends on the system and the states we want to explore. Here, we use optimal CAS to construct the reference function. The IVO-CASCI computations have been set up with the same CAS as the subsequent BWMRPT counterparts. Note that an IVO-CASCI wave function can handle nondynamic correlation effectively and some dynamic correlation also, but it is not practical to incorporate enough dynamic correlation effect in this step (one can add the remaining dynamic correlation in a more systematic and effective way through the MRPT scheme) exploiting larger CAS to gain desired chemical accuracy due to the computational cost. As the cost of a CAS-based ab initio calculation increases steeply (exponentially) with the size of the CAS, for large systems, we require a more affordable protocol.[52−55] Note that physical insight is very important to identify the excited states in SSMR calculations and to choose the CASs correctly. As the development of the gradient code for the IVO-BWMRPT approach has not been done yet in our group, in the present calculations, the geometries of various states of the systems considered here are taken from previously published works. Investigation of the basis set effect on the findings are also addressed here. The present MRPT computations have been performed using our in-house computer codes, which are interfaced with GAMESS.[77] The basis sets used here are obtained from the EMSL website.[78]

The choice of orbitals, especially active orbitals, is very crucial in illustrating the performance of MRPT methods, so it is also important to calibrate the accuracy of BWMRPT excitation energies using different orbitals. Thus, we have also depicted here how well the IVO-BWMRPT scheme compares with the previous method based on CASSCF orbitals (CASSCF-BWMRPT). For a given system, the same CAS was employed in the IVO-CASCI and IVO-BWMRPT procedures. Moreover, it is important to compare the BWMRPT results to the standard and state-of-the-art ab initio methods. However, comparisons with literature data are hindered by the differences in the computational schemes, especially with regard to the basis sets and geometries used in the computations, which are expected to affect the computed EEs.

Carbene Analogues: Methylene (CH2) and Silylene (SiH2)

Tunable diradical nature of the isovalent CH2 and SiH2 (analogues of carbenes) makes them very attractive model molecules for both experimental[79] and theoretical[21,25,68,80−96] investigations. The variation of the diradical nature of these systems, the spin multiplicities of the ground state, and equilibrium geometries can be explained with the help of the Walsh diagrams and Hund’s rules (for details, see ref (68)) exploiting the mutual interplay of correlations between the σ and π valence orbitals (leading electronic configurations of σ2π0, σ1π1, and σ0π2). Hund’s rule suggests that the ground state of CH2 is the triplet X̃3B1 (with one electron occupying the σ orbital and another one occupying the π orbital), whereas the Walsh correlation diagrams favor the closed-shell singlet X̃1A1 as the ground state for SiH2 (mostly the σ orbital doubly occupied and a small contribution from π orbital doubly occupied). This reversal in the spin multiplicity of the ground state on moving from CH2 to SiH2 has a marked effect on their chemistry.[97] CH2 and SiH2 are known for their very small adiabatic ST energy gap and thus are very sensitive to any errors made in the treatment of the differential correlation between the singlet and triplet states.

For CH2, a heterosymmetric diradical (as two orbitals occupied by the unpaired electrons have different symmetries), the ground X̃3B1 and second-excited b̃1B1 (open-shell in nature with singly occupied σ and π is the third state) states, can be described well with the single configuration state function (CSF) with the dominant configuration (1a1)2 (2a1)2 (1b2)2 (1b1) (3a1), whereas the lowest-lying singlet excited ã1A1 (zwitterion) state is a closed-shell MC state with considerable admixture of the configurations (1a1)2 (2a1)2 (1b2)2 (3a1)2 and (1a1)2 (2a1)2 (1b2)2 (1b1)2 and is therefore more challenging to describe. As experimental geometrical parameters are not reported for all of the diradical states of CH2 [and also SiH2] treated here, previously published widely used geometries [where the authors determined the geometries and energies of the states of CH2 using FCI] have been employed in the present work.[80−84] While in the FCI calculations, one core and one virtual orbitals were kept frozen, all electrons/orbitals are active in the present work. The minimum CAS consisting of two active electrons and two active orbitals, corresponding to highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO), CAS(2,2) is used for IVO-CASCI/CASSCF and subsequent BWMRPT calculations.

The vertical ã1A1–X̃3B1 splitting for CH2 calculated at the IVO-BWMRPT(2,2) level using geometrical parameters for the triplet ground state has been reported by Pittner et al.[98] Here, the basis set and scheme were adopted from ref (98). Results for VEEs are given in Table . Note that all approaches mentioned in the table yield acceptable (overestimated) VE gaps relative to the FCI value. The vertical ST gaps due to the present BW expansion-based SSMRPT calculations with both CASSCF (here termed as CASSCF-BWMRPT) and IVO-CASCCI (i.e., IVO-BWMRPT) orbitals are within 0.05 eV of the FCI data. The MR-AQCC and MR-CISD+Q results match well with IVO-BWMRPT and provide deviation from FCI by at most 0.009 eV. The performance of the MR-BWCCSD(RHF) method is much better, deviating from FCI by 0.005 eV. The MR-RSPT2 method overestimated the FCI result by 0.165 eV, which is larger than that from the other methods considered here.

Table 1

Vertical ã1A1–X̃3B1 Excitation Separation (eV) in CH2 (Re = 1.0895 Å, < = 131.48°)[98] Computed Using the cc-pVDZ Basis Set and CAS(2,2)

reference	method	ã1A1
present	IVO-BWMRPT(2,2)	1.052
	CASSCF-BWMRPT(2,2)	1.043
ref (98)	MR-AQCC	1.015
	MR-CISD+Q	1.010
	MR-RSPT2	1.171
	MR-BWCCSD	1.015
	FCI	1.006

The assessment of the accuracy of the IVO-BWMRPT for CH2 can also be calibrated by computing the adiabatic energy splitting. In Table , IVO-BWMRPT and CASSCF-BWMRPT single–triplet (ST) excitations are summarized together with the ones obtained using other commonly used approaches and experimental results. The optimized geometries obtained by FCI/TZ2P[83] (see ref (68) for details) and FCI/DZP[82] calculations have been used for such calculations. The “closed-shell” singlet, ã1A1, obtained using BWMRPT with different orbitals clearly shows that it is much lower in energy than the “open-shell” b̃1B1 ones. The high-level theoretical calculations provide AEEs that agree well with the experimentally derived estimates, yielding an energy range, depending on the specific scheme used, from 0.371 eV for CMRCI+Q[91] and complete basis set (CBS) extrapolation to 0.488 eV for CCSDT.[96] We also observed that established correlated excited-state schemes may remain sensitive to the zeroth-order/reference wave function.[82,83] As expected, the CASSCF-BWMRPT results for AEE are very close to those of IVO-BWMRPT in terms of energy differences. Benchmark investigations of CH2[80,82,83] clearly reflect that even FCI does not achieve very good agreement with experimental values when the FCI method is used with a small orbital basis set. Note that FCI ST energy gaps for ã1A1X̃3B1 are at least 0.077 eV above the best (nonrelativistic, Born–Oppenheimer) experimental estimate of 0.406 eV[79] [which is obtained from the Jensen and Bunker fit of experimental data using the MORBID Hamiltonian (see ref (82))]. The corresponding value for b̃1B1X̃3B1 is 0.117 eV. Apart from 1.425 eV,[79] other reliable estimates are 1.411 and 1.465 eV for the [b̃1B1X̃3B1] energy gap.[82]

Table 2

Adiabatic Excitation Energies (eV) Relative to the Ground State for Two Lowest-Lying States of Methylene (CH2)

reference	method	ã1A1	b̃1B1
present	IVO-BWMRPT(2,2)/TZ2P	0.604	1.646
	CASSCF-BWMRPT(2,2)/TZ2P	0.587	1.607
	IVO-BWMRPT(2,2)/DZP	0.681	1.624
	CASSCF-BWMRPT(2,2)/DZP	0.693	1.613
	IVO-BWMRPT(2,2)/cc-pVQZ	0.570	1.582
ref (68)	SF-CIS(D)/TZ2P	0.613	1.646
	SF-OD/TZ2P	0.514	1.564
	SF-CIS(D)/cc-pVTZ(cc-pVQZ)	0.578 (0.541)	1.579 (1.527)
	SF-OD/cc-pVTZ (cc-pVQZ)	0.483 (0.451)	1.495 (1.447)
	SF-OD/CBS-3pb [SF-CIS(D)/CBS-3pb]	0.430 (0.517)	1.420 (1.498)
ref (83)	FCI/TZ2P	0.483	1.542
ref (82)	FCI/DZP	0.549	1.667
ref (99)	SHCI/aug-cc-pVQZ	0.393	1.398
	CR-EOM-CC(2,3)/aug-cc-pVQZ	0.412	1.460
ref (91)	CMRCI+Q(6,6)/DZP	0.371
refs (86) and (87)	CASPT2(6,6)/DZP	0.670
	CASPT3(6,6)/DZP	0.555
	MR-CISD(6,6)/DZP	0.517
ref (88)	GVVPT2(6,6)/DZP	0.670
	SDS-MSMRPT2(6,6)/DZP	0.542
	MR-CISD(6,6)/DZP	0.420
ref (25)	SC-NEVPT2(6,6)/DZP	0.593
	PC-NEVPT2(6,6)/DZP	0.586
	SC-NEVPT3(6,6)/DZP	0.564
ref (90)	CASPT2(6,6)/DZP	0.677
	SC-NEVPT2(6,6)/DZP	0.670
	cu-SC-NEVPT2(6,6)/DZP	1.062
	cud-SC-NEVPT2(6,6)/DZP	0.620
ref (89)	MRMP2(SA)(6,6)/DZP(cc-pVTZ)	0.689 (0.629)
	MRMP2(SD)(6,6)/DZP(cc-pVTZ)	0.546 (0.438)
ref (71)	EOM-EE-CC(2,3̃)/TZ2P	0.488	1.545
	EOM-SF-CC(2,3̃)/TZ2P	0.497	1.554
	EOM-EE-CC(2,3)/TZ2P	0.454	1.539
	EOM-SF-CC(2,3)/TZ2P	0.477	1.542
ref (79)	experiment	0.406 (0.393)	1.425

As Table demonstrates, the IVO-BWMRPT/DZP overestimates the FCI adiabatic excitation energy for the ã1A1 state by 0.132 eV, and for the b̃1B1 state, by 0.043 eV. The CASPT2 and CASPT3 results with DZP for the ã1A1 excited state deviate from the corresponding FCI values[82] by −0.121 and −0.006 eV, respectively (the negative sign indicates the overestimation). The GVVPT2 and SDS-MSMRPT2 (static–dynamic–static multistate multireference second-order perturbation theory) results for DZP basis exhibit deviations[88] of −0.121 and 0.007 eV for the ã1A1 state. Note that SDS-MSMRPT2 is a CI-like scheme rather than being a perturbation approach and it does not require any level shifts of either orbital or state energies to avoid intruder-state problems. For methylene, the NEVPT2 results[25] are within 0.044 eV of those of FCI/DZP. The X̃3B1 and ã1A1 splitting of CH2 obtained in the case of both CR-CCSD(T)/DZP and CR-CC(2,3)/DZP calculations is 0.538 eV,[85] differing (underestimating) from the corresponding FCI calculations by 0.011 eV. The adiabatic ST splitting provided by IVO-BWMRPT(2,2)/DZP gives very good agreement with MRMP2(SA)/DZP and MRMP2(SD)/DZP[89] calculations with the full-valence reference space, CAS(6,6) [here, SA and SD stand for spin-averaged (SA) and spin-dependent (SD) orbital energies, respectively]. The deviation of MRMP2(SA) from FCI/DZP is −0.14 eV, which is improved by the MRMP2(SD) scheme, the discrepancy being only 0.03 eV. IVO-BWMRPT/DZP with CAS(2,2) successfully yields consistent results, which agree well with the MR-CISD ones with CAS(6,6).[86,87,91] All MRPT methods including the present one generally overestimate the adiabatic transition energy for the ã1A1 state.

The TZ2P basis set is more flexible than the DZP basis set. It is important to note that the estimates with TZ2P are in better harmony with experiment than those obtained in a similar benchmark study using the smaller DZP basis set.[83] The discrepancies between IVO-BWMRPT/TZ2P and FCI/TZ2P for the ã1A1 and b̃1B1 states are −0.121 and −0.104 eV, respectively. The corresponding errors provided by SF-OD/TZ2P are −0.031 and −0.022 eV, respectively. The error for the AEEs introduced by SF-CIS(D)/TZ2P is −0.13 eV for the ã1A1 state and −0.104 eV for the b̃1B1 state. As shown in Table , if more correlated SF strategies are employed, the errors become smaller. The errors against FCI for EOM-SF-CC schemes with TZ2P basis for both states are less than 0.04 eV.[71] Note that EOM-CC(2,3) provides an improvement over both EE-EOM and SF-EOM approaches since the errors against FCI are 0.03 and 0.01 eV, respectively.[71] The accuracy of IVO-BWMRPT/TZ2P also closely follows that of EOM-CC(2,3̃), as EOM-CC(2,3̃) errors against FCI are less than 0.1 and 0.03 eV for EE and SF variants, respectively.[71] From Table , one can observe that the IVO-BWMRPT(2,2)/TZ2P adiabatic transition energies for the ã1A1 and X̃3B1 states are also in agreement with high-level CC calculations, CCSDt/TZ2P (0.506 eV) and CC(t;3)/TZ2P (0.490 eV).[96] In CCSDt and CC(t;3) computations, the minimum active space, CAS(2,2), has also been employed.[96] The SF-DFT/6-31G(d) geometries reported in ref (68) are also used by Piecuch and co-workers[96] in their single-point energy calculations as in the present work. In passing, we note that the AEEs of the ã1A1 and X̃3B1 states of CH2 obtained in the CCSDt/TZ2P and CC(t;3)/TZ2P calculations differ from those resulting from the corresponding FCI calculations by 0.1 eV, which is a major improvement over the basic CCSD approximation that gives an error of 0.24 eV.[96]

Overall, the accordance between the IVO-BWMRPT/cc-pVQZ method and the experiment for AEEs is acceptably good as the results for the ã1A1 and b̃1B1 states deviate from the experimental values[79] by −0.164 and −0.157 eV, respectively. The SF-CIS(D)/cc-pVQZ values for AEE are off by −0.14 and −0.10 eV, respectively,[68] which are close to the errors provided by the present estimates. The corresponding errors for the SF-OD/cc-pVQZ scheme are −0.05 and −0.02 eV. Note that the SF-CIS(D) (SF-OD) results, provided using the CBS-3pb extrapolation [3pb describes three-point extrapolations based on cc-pVDZ, cc-pVTZ, and cc-pVQZ energies],[68] deviate by −0.111 (−0.024) and −0.073 (−0.005) eV from the experimental values. Therefore, SF-based calculations with the CBS-3pb scheme yield predictions in good agreement with IVO-BWMRPT/cc-pVQZ since the IVO-BWMRPT with the cc-pVQZ basis generates AEEs within 0.08 (0.16) eV of SF-CIS(D)/CBS-3pb (SF-OD/CBS-3pb).[68] The MRMP2(6,6)/cc-pVTZ values, obtained using the spin-averaged (SA) orbital, deviate by −0.223 eV from the experimental value for the lowest excited state, and when including the spin-dependent (SD) orbital, it is within −0.032 eV from the experiment. The CAS-BCCC4 AEEs with cc-pVDZ (0.493 eV) and cc-pVQZ (0.395)[93] are marginally more accurate than those resulting from the analogous ΛCCSD(T)/cc-pVDZ (0.525 eV) and CCSDT/cc-pVDZ (0.512) calculations[95] (by 0.02 eV). The CAS-BCCC4/cc-pVQZ method yields an error of 0.013 eV (with respect to experiment) for ã1A1, which is lower than the errors provided by MRPTs considered here. Table shows that IVO-BWMRPT/cc-pVQZ is also in reasonably good agreement with highly expensive CR-EOM-CC(2,3)/aug-cc-pVQZ with differences of up to 0.16 eV. IVO-BWMRPT/cc-pVQZ relative energies are generally within 0.18 eV of the benchmark semistochastic heat-bath configuration interaction (SHCI) values,[99] indicating that the present method is able to correlate electrons sufficiently to provide FCI-quality transition energies. As one can see from Table , AEEs provided by the accurate SHCI algorithm differ from experimental data by about 0.03 eV similar to the CAS-BCCC4 method. One can notice that diffusion Monte Carlo (DMC)[92] is in close agreement with IVO-BWMRPT(2,2)/cc-pVQZ calculations as ã1A1–X̃3B1 and b̃1B1–X̃3B1 gaps obtained from IVO-BWMRPT are 0.16 eV away from the DMC values (0.406 and 1.416 eV). The IVO-BWMRPT(2,2)/cc-pVQZ AEEs are also in reasonable agreement with the computationally expensive EOM-SF-CCSD/cc-pVQZ [0.447 and 1.432 eV for ã1A1 and b̃1B1, respectively] and EOM-SF-CCSD(dT)/cc-pVQZ [0.420 and 1.410 eV for ã1A1 and b̃1B1, respectively] methods.[100] One can observe deviations of IVO-BWMRPT/cc-pVQZ relative to EOM-SF-CC as large as 0.15 eV for ã1A1 and 0.172 eV for b̃1B1. The change of energy gaps is systematic with the size of basis sets used, and the correlation between IVO-BWMRPT/cc-pVQZ-calculated AEE and the available databased transition energies is quite reasonable despite being evaluated at the second-order perturbative level.

Unlike CH2, the first excited state of SiH2 is triplet, ã3B1. For SiH2, two quasidegenerate orbitals are 5a1 and 2b1. While both first [ã3B1] and second [Ã1B1, open-shell in nature] excited states in SiH2 are well characterized by the [core](4a1)2 (2b2)2 (5a1) (2b1) configuration, the ground state X̃1A1 has a moderately strong MR character due to mixing of [core](4a1)2 (2b2)2 (2b1)2 and [core](4a1)2 (2b2)2 (5a1)2. For SiH2, the geometry and basis set have been taken from the work of Bauschlicher and Taylor[81] for DZP. For another set of calculations, the CISD/TZ2P(f,d) optimized geometries and basis set reported in ref (84) have been used.[68] For SiH2, present calculations were performed with CAS(2,2) at the IVO-CASCI and IVO-BWMRPT levels. The results for AEEs for SiH2 have been reported in Table . Note that IVO-BWMRPT/TZ2P(f,d) AEE values deviate from CASSCF(6,6)-SOCI/TZ2P(f,d)[84] by no more than 0.08 eV for both states. Errors in AEEs for ã3B1 and Ã1B1 states between SF-OD/TZ2P(f,d) results and CASSCF(6,6)-SOCI/TZ2P(f,d)[84] data are at most 0.005 eV only. As in the case of CH2, there is a good overall agreement between the IVO-BWMRPT and SF-based methods for the TZ2P(f,d) basis. The IVO-BWMRPT AEEs, obtained using TZ2P(f,d), deviate by 0.125 and −0.113 eV from the experimental values, whereas SF-OD/TZ2P(f,d) yields values at most 0.06 eV away from the experiment. The errors between the SF-CIS(D) results and the experimental data for ã3B1 and Ã1B1 states are 0.140 and 0.194 eV, respectively. The corresponding errors (eV) provided by SF-CIS are 0.413 and −0.271. The AEEs obtained by EOM-SF-CCSD/cc-pVQZ [0.874 and 1.950 for ã3B1 and Ã1B1 states, respectively] and EOM-SF-CCSD(dT)/cc-pVQZ [0.892 and 1.937 eV for ã3B1 and Ã1B1 states, respectively][100] appear to be in reasonable agreement with the present estimates keeping in mind the second-order perturbative nature of the IVO-BWMRPT method.

Table 3

Adiabatic Excitation Energies (eV) Relative to the Ground State for Two Lowest-Lying States of SiH2

reference	method	ã3B1	Ã1B1
present	IVO-BWMRPT(2,2)/TZ2P(f,d)	0.791	2.041
	IVO-BWMRPT(2,2)/DZP	0.642	1.848
ref (68)	SF-CIS/TZ2P(f,d)	0.503	2.199
	SF-CIS(D)/TZ2P(f,d)	0.776	2.122
	SF-OD/TZ2P(f,d)	0.866	1.994
ref (84)	CASSCF(6,6)-SOCI/TZ2P(f,d)	0.871	1.992
ref (90)	FCI/DZP	0.759
	CASPT2/DZP	0.654
	SC-NEVPT2/DZP	0.724
	cu-SC-NEVPT2/DZP	0.299
	cud-SC-NEVPT2/DZP	0.702
ref (101) and (102)	experiment	0.916 ± 0.03	1.928

The X̃1A1 → ã3B1 gap, calculated from the IVO-BWMRPT/DZP method, is in agreement with the FCI/DZP estimate with a discrepancy (underestimated) of 0.117 eV. It can be seen that CASPT2/DZP yields a value 0.105 eV away from the FCI/DZP value for the adiabatic X̃1A1 → ã3B1 transition. The IVO-BWMRPT AEE ST gap agrees with the SC-NEVPT2 one with a difference of 0.082 eV. Tables and 3 clearly indicate that the cud-SC-NEVPT2 (cumulant with diagonal-based SC-NEVPT2) approximated scheme provides notably good estimates. It should be mentioned that the ST gaps for the ã3B1 state calculated with the SC-NEVPT2 and cud-SC-NEVPT2 methods deviate from the FCI data[90] for the DZP basis by only 0.035 and 0.057 eV, respectively, which is in line with the IVO-BWMRPT error. The IVO-BWMRPT/DZP energy gaps deviate from the experimental values by 0.274 and 0.08 eV for ã3B1 and Ã1B1 states, respectively. The CASPT2/DZP and SC-NEVPT2/DZP numbers are off by 0.26 and 0.19 eV, respectively, for the ã3B1 state. Similar to the CH2 system, IVO-BWMRPT also leads to a notable deviation from CAS-BCCC4 (0.904 eV)[93] for X̃1A1 → ã3B1 adiabatic excitation. The ST splitting energies for the ã3B1 state provided by these MRPT methods are very close to each other and relatively far away from experimental results.[101,102] This fact indicates that the second-order MRPT may lead to an unbalanced recovery of correlation for the ã3B1 and Ã1B1 states. Moreover, the ST separation of XH2 [X = C and Si] is a well-known property that has a large differential correlation effect and different basis set requirements for both ã3B1 and Ã1B1 states.[82−84,88] Thus, it is a challenging job to treat both states with a democratic level of accuracy and requires further study.

The vertical energy gaps have also been calculated for SiH2 as for CH2. Calculations were performed using the cc-pVDZ basis set and geometries described in ref (98). BWCCSD, MR-RSPT2, MR-CISD+Q, and MR-AQCC results of Pittner and co-workers[98] have been used here to calibrate the new IVO-BWMRPT vertical excitation energies. The IVO-BWMRPT value of the TS gap (0.395 eV) agrees well with the BWCCSD one (0.35 eV). Note that IVO-BWMRPT and BWCCSD yield TS gaps exhibiting 0.04 and 0.001 eV difference from FCI (0.351 eV), respectively. The MR-RSPT2 value is off by 0.10 eV from FCI.[98] It can be seen that MR-AQCC and MR-CISD+Q methods predict the gaps showing 0.004 eV deviation from FCI.[98]

Ethene (Ethylene)

The next example is ethylene, a homosymmetric diradical in a twisted form, belonging to a class of molecules that is fairly difficult to describe due to the ubiquitous doubly excited state that it possesses. Evaluation of the VEEs from the ground-state singlet 1Ag (N-state) to the first excited triplet 3B1u (T-state) or singlet 1B1u (V state) for ethene has always received great interest of the quantum chemistry community, and proper description of states of different characters, ionic and covalent in this case, is a computational challenge.[57,58,103−112] Note that the accurate description of singlet diradicals is difficult owing to their explicit MR nature. The nature of the V (1B1u) state of ethene especially presents various difficulties owing to its unexpectedly diffuse nature. The interplay of various complex effects [such as dynamic σ-polarization,[103−105] the contraction of the p-orbitals[113,114] as a consequence of the dynamic σ-polarization, the differential dynamical correlation of the electrons (much more localized in the ionic structures than in the neutral ones), and the possible pronounced valence–Rydberg mixing between the qualitatively different configurations arising from π → π* and π → 3d π excitations] must be taken into consideration for the correct description of the ionic π → π* [corresponding in orbital notation to 1b3u2 → 1b3u1 1b2g1] excited states. In 2009, Rydberg and the valence nature of the V state were intensively investigated by Angeli, and it was reported for the first time that the significant charge fluctuation within the π-system needs a dynamic response of the σ frame already at the orbital optimization step,[113,114] which is of central importance for correct evaluation of AEE. It explains the emergence of the poor convergence of the conventional methods where dynamic correlation is computed with a preoptimized set of molecular orbitals (MOs), which does not incorporate such σ–π dynamic mixing. Various established methods are used to cope with this difficulty. Such an effect cannot be captured at the second order in MRPT, even for PT2 methods, allowing the modification of the wave function within the reference space.[113] At this juncture, it should be restated that the IVO-BWMRPT allows revision of the reference coefficients in the first-order wave function through diagonalization of the perturbatively constituted effective Hamiltonian. Nevertheless, the CASCI is a full CI in a given reference space and hence the CASCI wave function has been constructed via the diagonalization of the Hamiltonian in a CASCI space. Therefore, the effects due to orbital optimization can be encompassed to some extent in IVO-BWMRP if the CAS is appropriately selected (although the IVO-CASCI scheme does not involve explicit orbital optimization), a feat that cannot be accomplished by other successful unrelaxed/contracted correlated methods. Consequently, the IVO-BWMRPT wave function is less sensitive to the orbital relaxation/optimization effects compared to the contracted MR methods. The 1B1u state is markedly more diffuse than its high-spin counterpart 1 3B1u [which possesses a clearly valence nature] as described in terms of valence bond theory and hence the (1B1u, 3B1u) states differ in energy by almost a few eV’s. The calculations at SCF and π-CASSCF levels showed the π* orbital to have significant Rydberg property. Due to the lack of dynamical correlation, the CASSCF puts the ionic 1B1u state very close to a Rydberg state (3dπ state 2 1B1u), inviting unwanted valence–Rydberg mixing. The 3B1u state of ethylene is dominated by a single excitation (3π → π) at the ground-state equilibrium geometry and is thus easier to treat theoretically.

Present computations are carried out at the experimental gas-phase equilibrium geometry of the ground state[115] using the ANO[57] and cc-pVTZ[78] basis sets. For ethene, our MR-calculations were performed using CAS(2,2) and CAS(6,6). It should be mentioned that Schreiber et al.[107] and Schapiro et al.[108] have found that geometrical constants obtained by the MP2/6-31G* scheme are very close to the experimental values. The molecule lies in the yz plane with the two carbon atoms on the z axis (z axis extends along the C=C bond). For ethylene, CAS(2,2) has been used as the reference space in the IVO-CASCI and subsequent BWMRPT calculations. Table compiles various theoretical and experimental results for VEEs for low-lying excitation states of ethene. Experimentally and theoretically estimated energies for the N–T splitting are low as the molecule twists immediately upon excitation to the triplet state. Electron energy loss spectroscopy revealed that N–T vertical transition occurs at 4.21 eV, whereas electron impact spectroscopy has predicted the value to lie in the range 4.21–4.36 eV.[116,117] Assembled values in Table reflect that IVO-SSMRPT VEEs depend rather weakly on the nature of the basis set used, except for the 1 1B1u state. Owing to the ionic character of the 1 1B1u state, diffuse functions are essential for proper description of the state. On the other hand, the incorporation of diffuse functions in the basis set also invites an artificial mixing of the 1 1B1u state with Rydberg states, which brings the computed VEE very close to the experimental estimate for the wrong reason(s). Schreiber et al.[107] found that at the CC level, the basis set dependence is smaller but still substantial for the 1 1B1u state (such as CC3 with TZVP provides the value 8.37 eV, whereas the converged CC3 value is 7.88 eV). The BW version of the SSMRPT method with different orbitals, say IVO-CASCI and CASSCF, is able to correctly predict the energy ordering of the states treated here and reproduce the experimental data rather well. The difference between the IVO-BWMRPT and the CASSCF-BWMRPT strategy is not more than 0.23 eV, indicating very similar performance of BWMRPT with both IVO-CASCI and CASSCF orbitals. The BWMRPT schemes provide small deviation from the experimental data. For the 1 1B1u–1 1Ag and 1 3B1u–1 1Ag transitions, IVO-BWMRPT(2,2)/cc-pVTZ (CASSCF(2,2)-BWMRPT/cc-pVZT) theoretical values (eV) are 0.48 (0.37) and 0.18 (0.35) higher than the experimental band maximum.[118,119] The corresponding deviations (eV) for IVO-BWMRPT(2,2)/ANO are 0.60 and 0.25. The IVO-BWMRPT values for the VEEs lie in the same range as other MRPT estimates. The IVO-BWMRPT(2,2)/cc-pVTZ transition energies for N–V and N–T are very close to the MR-CISD+Q(2,2)/cc-PVTZ and MRMP2(2,2)/cc-pVTZ values. Errors (eV) relative to experimental values are 0.29 for the N to V transition and 0.09 for the N to T transition at the MRMP2 level and 0.42 for the N to V transition and 0.06 for the N to T transition at the MR-CISD+Q level of calculations. It should be mentioned that MROPT(2)/cc-pVTZ values match the available experimental N–V and N–T excitation energies quite well with the deviations (eV) 0.04 and 0.17, respectively. For N–V and N–T transitions, MS-CASPT2(2,2) [which diagonalizes an effective Hamiltonian obtained from the single-root CASPT2 wave functions] values with the TZVP basis of Thiel and co-workers[107] are 0.62 and 0.24 eV higher than the experimental dada similar to those of the present estimates. Additionally, one can note that the values of QD-NEVPT2(2,2)/TZVP[108] with various contraction schemes for the N–V and N–T transitions are also in good agreement with those computed by IVO-BWMRPT. The TZVP basis set does not contain diffuse functions and hence may not be suitable for the description of those excited states that are spatially extended and have Rydberg nature. The present IVO-BWMRPT(2,2)/cc-pVTZ result for the N–V transition energy shows acceptable agreement with those of the earlier MR-AQCC(2,2)/MR-CISD+Q(2,2) in the CBS limit (8.06/8.05 eV).[104] One can also find that the IVO-BWMR(2,2)/cc-pVTZ scheme yields excitation energy for the N to T transition that is in good agreement with those generated with the MR-AQCC(2,2)/MR-CISD+Q(2,2) calculations at the CBS limit (4.48/4.47 eV).[104] The N–V energy gap obtained from iFCI/cc-pVTZ (8.06 eV)[111] and the FCIQMC/ANO-L-pVTZ N–T excitation energy (4.64 eV)[110] are also in reasonably good harmony with the present IVO-BWMRPT/cc-pVTZ estimate with errors of 0.4 and 0.1 eV, respectively. Table shows that the semistochastic heat-bath configuration interaction (SHCI) algorithm with ANO-L-pVTZ and CR-EOM-CC(2,3)D with ANO-L-pVT estimates[99] for N–V and N–T transitions exhibit good accordance with IVO-BWMRPT(2,2). Previous work[99] indicates that CC methods must include more than triple excitations to obtain acceptable relative energies of ethylene. Therefore, the IVO-BWMRPT method can also run in situations when valence–Rydberg mixing often happens at the IVO-CASCI/CASSCF level.

Table 4

Vertical Excitation Energies (eV) Relative to the Ground State for Three Low-Lying States of Ethylenea

reference	method	1 1B1u	2 1B1u	3B1u
present	IVO-BWMRPT(2,2)/ANO	8.60	9.52	4.61
	IVO-BWMRPT(6,6)/ANO	8.32	9.18	4.42
	IVO-BWMRPT(2,2)/cc-pVTZ	8.48	9.39	4.54
	CASSCF-BWMRPT(2,2)/cc-pVTZ	8.37	9.45	4.71
ref (106)	MR-CISD(2,2)+Q/cc-pVTZ	8.42		4.53
	MRMP2(2,2)/cc-pVTZ	8.29		4.45
	MROPT(2)(2,2)/cc-pVTZ	8.04		4.53
ref (120)	CASPT2(2,11)/ANO	8.40	9.31	4.39
ref (57)	CASPT2(6,6)/ANO	8.45	8.93
	MS-CASPT2(6,6)/ANO	7.98	9.40
ref (58)	CASSCF(6,6)/ANO	8.32	9.10
	QD-PC-NEVPT2(6,6)/ANO	8.10 (8.06)	9.35 (9.31)
	PC-NEVPT2(6,6)/ANO	8.17 (8.18)	9.18 (9.19)
	MS-CASPT2(6,6)/ANO	8.03	9.42
	CASPT2(6,6)/ANO	8.25	9.20
ref (107)	MS-CASPT2(2,2)/TZVP	8.62		4.60
ref (108)	QD-SC-NEVPT2(2,2)/TZVP	8.69		4.60
	QD-PC-NEVPT2(2,2)/TZVP	8.64		4.60
	CC3/TZVP	8.37		4.48
ref (104)	MR-CISD+Q(2,2)/CBS	8.05		4.48
	MR-AQCC(2,2)/CBS	8.06		4.47
	MR-CISD+Q(6,6)/CBS	7.81		4.54
	MR-AQCC(6,6)/CBS	7.88		4.53
ref (121)	Hv2nd (Hv3rd)	8.00 (8.15)	9.34 (9.39)	5.09 (4.22)
ref (99)	SHCI/ANO-L-pVTZ	8.05		4.59
	CR-EOM-CC(2,3)D/ANO-L-pVTZ	8.25		4.76
refs (118) and (119)	experiment	7.66	9.33
refs (116) and (117)	experiment		9.32	4.36
ref (148)	experiment	8.00

For ref (58) values of SC-NEVPT2(6,6)/ANO are given in parentheses.

To investigate how the size of CAS affects the VEEs in ethene, the CAS has been changed by inclusion of two σCH orbitals (such as 3ag and 1b3g) and two σCH* orbitals (2b2u and 3b1u) on the top of the π and π* active orbitals in the next set of calculations. The same reference space, CAS(6,6), has also been used by various groups in their MR-based calculations.[58,104] The effect of CAS size is substantial, which is reflected in a notable perturbation correction of the transition excitation energies 8.60 and 4.61 eV at the IVO-BWMRPT(2,2) level to 8.32 and 4.42 eV at the IVO-BWMRPT(6,6) level for 1 1B1u and 3B1u and brings them close to experimental estimates. The BWMRPT calculations with IVO-CASCI(6,6)/ANO show errors of 0.32 and 0.06 eV for 1 1B1u and 3B1u states, respectively. One can also note that an overall improvement of the IVO-BWMRPT estimates relative to high-quality estimates such as unc-MR-CISD+Q/aug-cc-pVDZ,[112] EOM-CCSD(T),[109] and MR-AQCC(6,6)/CBS as well as MR-CISD+Q(6,6)/CBS results[104] occurs due to the enhancement of CAS size. Due to enhancing the CAS, the IVO-BWMRPT/ANO VEEs decrease by about 0.3 eV. With an extension of the basis set (d-aug-cc-pV5Z) and the CAS(4,8), MS-CASPT2 calculations lead to lowering of N–V and N–T vertical excitation energies to 7.83 and 4.49 eV, respectively[107] [relative to 8.62 and 4.60 eV for MS-CASPT2/TZ2P]. Thiel and co-workers[107] showed that enhancing the CAS lowers the MS-CASPT2(2,2)/TZVP values by about 0.2 eV and brings them to the databased estimates (say CC3) for both transitions. For the N–V and N–T transitions, CASPT2(2,11)/ANO[120] values exhibit deviations of 0.40 and 0.03 eV from the corresponding experimental values, reflecting that the performance of the IVO-BWMRPT(6,6) scheme is better than that of CASPT2. The N–V vertical transition energies provided by the MS-CASPT2(6,6) level of calculations[57] with the ANO basis set lead to an error of 0.02 eV, which is smaller compared to the single-root CASPT2 version (providing an error of 0.45 eV).[57] This fact indicates that the relaxation of reference coefficients is very crucial to achieve correct results in the case of the CASPT2 protocol. Similar to the CASPT2 calculations, the change of the reference space from CAS(2,11) to CAS(6,6) leads to a decrease of errors of the computed VEEs.[57,120] Angeli et al.[58] also observed the modification in the excitation energies toward the experimental values for both single- and multi-root versions of CASPT2 and NEVPT2 methods. It should be mentioned that the results published so far suggested that the VEE of the 1 1B1u state can still not be identified accurately but should lie between 7.7 and 8.0 eV. The best estimate of the present work is 8.32 eV, which yields an error of 0.32 eV. It is important to note that IVO-BWMRPT and MRMP2 overestimate the vertical excitations, while CASPT2 and NEVPT2 underestimate them. Although IVO-BWMRPT(6,6) exhibits a notable deviation for 1 1B1u excitation, it shows close accordance with the extensive MR-CISD+Q/CBS and MR-AQCC/CBS calculations of Müller et al.[104] that provide values in a narrow range between 4.47 and 4.53 eV. The differences between the MS-CASPT2(4,8)/d-aug-cc-pV5Z[107] [and MR-CISD(6,6)+Q/CBS[104]] and IVO-BWMRPT(6,6) energy gaps for N–V excitation are much larger; the IVO-BWMRPT(6,6) energies correspond to the experiment in a marginally better way. The excitation energies provided by Hv2nd (Hv3rd) computations with eight valence orbital reference space [which includes both Rydberg and valence p-electron improved virtual orbitals][121] also compare favorably with IVO-BWMRPT(6,6)/ANO and match quite well with experiment with deviations within 0.2 and 0.7 eV for N–V and N–T transitions, respectively. The point to be noted is that a single set of orbitals is used in Hv calculations to describe all excited states. However, researchers often argue that a correct ab initio theory needs introducing different orbitals for different electronic states[103,121] (as done in SSMR calculations). Moreover, the Hv method is numerically less robust because of numerical instabilities due to the intruder-state problems.

Comparing with the findings of Roos and co-workers,[57] one can observe a notable difference of CASPT2 calculations of Angeli et al.[58] in the treatment of the 2 1B1u state. This is due to the contraction scheme of the basis set used in ref (57) that leads the 3dπ Rydberg 2 1B1u state to a level of high energy that may affect the correct treatment of the valence–Rydberg mixing in the 1 1B1u state (vide supra). Finley et al.[57] have exploited the contraction of the diffuse functions to provide a more concentrated representation of the Rydberg states as well as to avoid intruder-state problems that are often found in CASPT2 calculations with basis sets with many diffuse functions. This is due to the fact that the low-energy diffuse functions may lead to a generation of small denominators in the perturbation expansion when one-electron zero-order Hamiltonians are used in the MRPT scheme (say CASPT2). Description of the 3d Rydberg state 2 1B1u (π → 3dπ, Rydberg) merits some discussion as characterization of Rydberg states is considerably more difficult to describe compared to other states considered here. The restricted open-shell HF (ROHF) level of calculations on the 1 1B1u state leads to artificially close energy to the 2 1B1u state. This yields strong mixing between the two states, 1 1B1u and 2 1B1u, resulting in the ROHF 1B1u wave function to become spatially diffuse. In general, the comparison of VEEs between the BWMRPT(6,6) energies and reference values indicates that IVO-BWMRPT(6,6) does very well in reproducing 2 1B1u VEE. The deviation of IVO-BWMRPT(6,6)/ANO value from the corresponding experimental data is 0.15 eV, whereas IVO-BWMRPT(2,2)/ANO provides error of 0.2 eV. Here, it should be mentioned that the perturbation correction of the 2 1B1u VEE due to the change of CAS is not very significant for both BWMRPT and CASPT2 schemes. Angeli et al.[58] also observed that the variation in the excitation energies due to change of reference space from the CAS(2,11) to the CAS(6,6) is not significant for 2 1B1u. IVO-BWMRPT(2,2) results indicate that the VEE for the 1 1B1u state changes by 0.34 eV while moving from ANO to cc-pVTZ basis sets. Agreement of CASPT2 (MS-CASPT2)[57] with experimental VEE to within 0.4 (0.07) eV has been found for 2 1B1u. The reason for the large error of CASPT2 is the mixing at the CASSCF level, between 1 1B1u and 2 1B1u states. Roos and co-workers[57] demonstrated that treating both the states democratically using the relaxed version of CASPT2, the MS-CASPT2 scheme can remove this unwanted mixing and provide VEE of 7.98 eV for the 1 1B1u state, which is more close to the experimental data. It should be mentioned here that modification of reference space coefficients (as occurring in IVO-BWMRPT and MS-CASPT2 calculations) permits the 1 1B1u and 2 1B1u states to interact at the dynamical correlated level at which the energy of the 1 1B1u state is much lower. In fact, CASPT2 with the standard Fock operator faces intruder-state problems (partially alleviated by the use of level shifts) that affect the correct description of both 1 1B1u and 2 1B1u states and hence their mutual interaction (mixing). It should be reemphasized that IVO-BWMRPT is free from such a problem. Further improvements of IVO-BWMRPT values for both 1B1u states need the addition of diffuse functions. The 2 1B1u–1Ag splitting yielded by expensive Hv2nd (9.34 eV) and Hv3rd (9.39 eV) calculations are marginally (0.02 and 0.07 eV) higher than experimental estimates. For the sake of completeness, it is useful to address the performance of IVO-BWMRPT values vis-à-vis the results obtained by Mk-MRCCSD and SU-MRCCSD using a mixed basis set (aug-cc-pVTZ for C and aug-cc-pVDZ for H).[122] Although IVO-BWMRPT(6,6) AEE for the 1 1B1u state shows notable deviation from Mk-MRCCSD (7.67 eV) and SU-MRCCSD (7.94 eV), the IVO-BWMRPT(6,6) energies agree very well with the observed SU-MRCCSD values for other two states (9.27 and 4.67 for 2 1B1u and 3B1u states, respectively). Moreover, Mk-MRCCSD AEE (4.80 eV) for the 3B1u state agrees very well with the IVO-BWMRPT(6,6) value. It should be mentioned that the Mk-MRCCSD method generates the same AEE values as does SU-MRCCSD for all of the three target states. The Mk-MRCCSD and SU-MRCCSD transition energies for the 3B1u state are marginally higher than the IVO-BWMRPT one, indicating the usefulness of the BWMRPT method.

The adiabatic N–T transition energy has also been estimated using MP2/cc-pVTZ optimized geometries of ethylene using GAMESS.[77] The adiabatic singlet–triplet energy gap of IVO-BWMRPT/cc-pVTZ (3.01 eV) is found to differ by 0.5 eV from the value reported in a recent photodissociation experiment (2.52 eV), suggesting that the IVO-SSMRPT value is in good accord with the photodissociation finding of Qi et al.[123] The CASSCF(4,8)-based all-electron DMC calculations using effective core potential (ECP)[124] provide a value of 2.88 eV, which is close to the IVO-BWMRPT/cc-pVTZ value. The MRCI adiabatic N–T energy gap (2.78 eV)[124] is also found to be in close agreement with the present IVO-BWMRPT values, indicating the reliability of the present method for estimation of transition energies.

Trimethylenemethane (TMM) Diradical

It is important to ascertain the usefulness of the IVO-BWMRPT on more complicated systems, so TMM has been chosen as a next test molecule. According to Hund’s rule, the ground state of TMM is a triplet (3A′2, D3 symmetry) as four π electrons are distributed over the four π-type molecular orbitals (MOs), two (e′) of which are exactly degenerate in the D3 framework, which is also supported by both electron paramagnetic resonance (EPR) and IR experimental results[125−128] and previous established theoretical investigations.[68−70,94−96,129−132] The first closed-shell 1 1A1 and the open-shell 1 1B2 singlets (similar to the N and V states of ethylene, respectively) are exactly degenerate in D3 geometry. However, owing to the Jahn–Teller effect, their degeneracy at the D3 geometry is lifted in lower C2 symmetry and both states have a manifestly MR biradical character. Despite numerous attempts at the state-of-the-art level, reliable multiplet splitting in a non-Kekulé TMM system is a challenging task for electronic structure theories[68−70,94−96,129−132] because all four of its π orbitals are very close in energy. We should also mention that the use of multidimensional reference spaces along with the state-of-the-art MRCC-like approaches does not necessarily promise good-quality results.

The geometries of various states of TMM are taken as those optimized at the SF-DFT/6-31G(d) level of Slipchenko and Krylov.[68] These geometries have also been used by other workers.[69,94,131] It has been found that geometries obtained with the various methods such as CCSD(T)/cc-pVTZ,[129] SF-DFT/6-31G(d),[68] CASSCF,[130] IVO-CASCI,[133] and IVO-RSMRPT[132] are close to each other because the energy surface is very flat around the equilibrium or optimized structure. We have used the DZP and cc-pVTZ basis sets[78] and CAS(4,4). The four lowest occupied orbitals have been kept frozen in the correlated BWMRPT calculations. Computed adiabatic energy separations between the states are summarized in Table . The 1 1B1 state has not been identified in the experiment. All methods reported here including the present IVO-BWMRPT one overestimate the energy splittings between the ground and the excited 1 1A1 state. On analyzing the values in Table , one can find a strong influence of the scheme used for computing the splitting, especially for the ΔEAEE gap, and see that a high level of correlation is crucial for this system. As those of other previously published well-accepted findings,[68−70,94,95,129−132] the present IVO-BWMRPT calculations correctly predict a triplet ground state, 3B2 (3A′2, D3), for TMM and also reveal that TMM in the open-shell singlet 1 1B1 state is lower in energy than in the closed-shell singlet 1A1 state by a few eV’s as predicted by other established methods. As shown in the table, the BWMRPT estimated ΔEAEE ST splitting for different basis sets ranging from 0.678 to 0.721 eV. The corresponding variation for ΔE1AEE values ranges from 0.895 to 0.837 eV. This fact advocates that IVO-BWMRPT excitation energies are not very sensitive to the basis sets used here. Previous SF-based calculations of Krylov and co-workers[68,69] have also demonstrated this viewpoint in their MCQDPT and SF-based calculations. Comparing results in BWMRCC with cc-pVDZ and cc-pVTZ basis sets, one can also see that the influence of the basis set size on computed ΔEAEE [130] is not significant. Cramer and Smith[129] found that increasing the basis set size from cc-pVDZ to cc-pVTZ has a negligible effect on the relative CASPT2 TMM energies. In the case of IVO-BWMRPT, there is a 0.06 eV drop in the excitation energy of the 1 1A1 state if one moves from the DZP to the cc-pVTZ basis set, whereas AEEs increase by 0.04 eV for the 1 1B1 state. For the DZP basis, there is good agreement of the IVO-BWMRPT findings with the recent CASPT2, CAS-BCCC4, and MR-CISD results with CAS(4,4).[131] Li and co-workers[131] also use the SF-DFT/6-31G(d) optimized geometries of the ground state and two lowest excited states of Krylov and co-workers.[68] The gap between the 1 1B1 and 1 1A1 states predicted by IVO-BWMRPT is within 0.22 eV. The corresponding values due to CASPT2, CAS-BCCC4, and MR-CISD calculations are 0.20, 0.15, and 0.19 eV, respectively. The small adiabatic energy separation between the states considered here permits us to control the electronic structural properties of the system. The IVO-BWMRPT(4,4)/DZP and IVO-BWMRPT(4,4)/cc-pVTZ results for the ΔEAEE transition are within −0.196 and −0.138 eV, respectively, from experiment, whereas for CASPT2(4,4)/DZP,[131] the deviation is about −0.157 eV. The CAS-BCCC4/DZP and MR-CISD(4,4)/DZP numbers lead, respectively, to the deviations of −0.127 and −0.145 eV from the experimental data. The IVO-BWMRPT(4,4) results are in good agreement with not only MCQDPT2(4,4) AEEs but also the expensive MCQDPT2(10,10) estimates. At this juncture, it should be noted that the CAS(4,4) includes four π orbitals and four π electrons. The larger CAS(10,10) includes six additional carbon–carbon σ and σ* orbitals and six additional electrons participating in the C–C bonding. The IVO-BWMRPT(4,4)/DZP [IVO-BRMRPT(4,4)/cc-pVTZ] results differ from the MCQDPT2(10,10)/DZP [MCQDPT2(10,10)/cc-pVTZ] by 0.002 [0.011] eV for the 3A′2 → 1B1 adiabatic transition and by 0.071 [0.009] eV for the 3A′2 → 1A1 one. The performances of IVO-BWMRPT and SF-OD methods appear to be of similar quality. For the cc-pVTZ basis, the SF-OD values deviate from MCQDPT2(10,10) by 0.034 eV for the 3A′2 → 1B1 adiabatic excitation and by 0.113 eV for the vertical excitation from 3A′2 to 1A1. The MCQDPT2(4,4) and MCQDPT(10,10) results with the cc-pVTZ basis for the 3A′2 → 1A1 splitting deviate from the experiment by at most 0.17 eV.

Table 5

Adiabatic Splitting (eV) of Multiplets (1A1 and 1B1) Relative to the Ground State (3A2′) of TMM

references	methods	ΔE1B1	ΔE1A1
present work	IVO-BWMRPT(4,4)/DZP	0.678	0.895
	IVO-BWMRPT(4,4)/cc-pVTZ	0.721	0.837
ref (131)	CAS-BCCC4(4,4)/DZP	0.676	0.826
	CASPT2(4,4)/DZP	0.667	0.856
	MR-CISD(4,4)/DZP	0.659	0.844
ref (68)	SF-CIS(D)/DZP	0.987	0.885
	SF-OD/DZP	0.697	0.936
	MCSCF(4,4)/DZP	0.643	0.843
	MCQDPT2(4,4)/DZP	0.676	0.863
	MCSCF(10,10)/DZP	0.704	0.834
	MCQDPT2(10,10)/DZP	0.674	0.824
	SF-CIS(D)/cc-pVTZ	1.025	0.893
	SF-OD/cc-pVTZ	0.744	0.941
	MCQDPT2(4,4)/cc-pVTZ	0.711	0.862
	MCQDPT2(10,10)/cc-pVTZ	0.710	0.828
ref (129)	CASPT2N(4,4)/cc-pVDZ	0.685	0.872
	CASPT2N(10,10)/cc-pVDZ	0.680	0.828
ref (71)	EOM-SF-CCSD/cc-pVTZ	0.554	0.933
	EOM-SF-CC(2,3)/cc-pVTZ	0.697	0.788
ref (94)	2R-MRCCSD(T)/cc-pVTZ		0.980
	2R-MRCCSD(RHF)/cc-pVTZ		1.327
	2R-MRCCSD(MCSCF)/cc-pVTZ		1.032
	2R-SU-MRCCSD/cc-pVTZ	0.629	0.698
	2R-SU-MRCCSD(T)/cc-pVTZ	0.685	0.763
	CCSD(T)/cc-pVTZ		1.028
	CCSD/cc-pVTZ		2.021
ref (130)	4R-BWCCSD/cc-pVTZ	0.634	0.779
	4R-BWCCSDit/cc-pVDZ	0.631	0.791
	4R-BWCCSDit/cc-pVTZ	0.633	0.771
ref (96)	CCSD/cc-pVTZ		2.033
	CCSDT/cc-pVDZ		0.941
	CR-CC(2,3)/cc-pVTZ		1.383
	CCSD(T)-h/cc-pVTZ		0.837
	CC(t;3)/cc-pVTZ		0.950
ref (125) and (126)	experiment		0.699

We now focus on the performance of IVO-BWMRPT vis-á-vis CC methods. The CCSD level of calculations provides certainly very unsatisfactory estimates for adiabatic 3A′2 → 1A1 (recall that the experimental value is 0.699 eV), as it should be. This is basically due to a poor efficacy of the SD-truncation scheme for the singlet state, which has a pronounced MR nature. EOM-CC values reported in ref (71) reveal that triple excitations are very crucial for getting correct values for excitation energy of this system. In fact, the incorporation of triples in the EOM-SF-CCSD (SF-CC) scheme alters the energies of the excited states by around 0.2 eV in favor of experimental estimates. However, the IVO-BWMRPT(4,4)/cc-pVTZ ΔEAEE value is as impressive as the extrapolated estimates, which emerge for full active space EOM-SF-CC(2,3)/cc-pVTZ calculations.[71] The corresponding CCSDT/cc-pVDZ[96] result is not as good as that of EOM-SF-CC(2,3)/cc-pVTZ, which may be a consequence of using the RHF orbitals in the computations for the closed-shell singlet 1A1 plagued by near-degeneracy. Clearly, when treating 1A1, it is the near-degeneracy of the leading configuration or, equivalently, the MR nature of the state of interest that is decisive for the success or failure of the method employed. One can notice that CCSDt, CCSD(T), and its CR-CC(2,3) variants although improve the poor CCSD gaps, neither CCSD(T) nor CCSD(2)T estimations, and not even CR-CC(2,3), offer an acceptable performance. It should be noted that the λCCSD(T) and DEA-CCSD methods also appear to have larger errors.[95] This is due to the fact that the MR effect in the singlet state is notable. CCSD(T)-h/cc-pVTZ yields [noticeably more accurate than CR-CC(2,3)] consistent results, which agree well with the IVO-BWMRPT ones in describing ΔEAEE.[96][96] Although the best values for ΔEAEE obtained from the sophisticated MRCC and RMRCC computations fall into the 0.694–0.997 eV wide range,[94] the CC(t;3)/cc-pVTZ value of the ΔEAEE of 0.941 eV in TMM[94] is in acceptable accordance with the IVO-BWMRPT/cc-pVTZ result of 0.837 eV. None of the other methods that incorporate noniterative triples’ corrections to the CCSD or CCSDt energies can provide highly impressive estimates. One can notice that the energy separations obtained by the two-reference MRCCSD (2R-MRCCSD) method with RHF and MCSSCF orbitals are also too large, amounting to 1.327 and 1.032 eV with RHF and MCSSCF orbitals, respectively. The SUMRCC results with both SD and SD(T) truncations[94] for ΔEAEE and ΔEAEE reported in ref (94) show very good agreement with the corresponding IVO-BWMRPT ST splittings. The performance of IVO-BWMRPT(4,4)/cc-pVTZ is also close to that of the computationally expensive EOM-SF-CCSD/cc-pVQZ [0.788 and 0.918 eV for ΔE and ΔE, respectively] and EOM-SF-CCSD(dT)/cc-pVQZ [0.696 and 0.846 eV for ΔE and ΔE, respectively] methods.[100] One can notice that the multireference BWCC estimates get disconcertingly worse with the addition of triples, although the four-reference BWCCSD (4R-BWCC) scheme[130] provides very good results as does the present IVO-BWMRPT and EOM-SF-CC(2,3)[71] methods.

At this juncture, it should be emphasized that CASSCF(2,2)/DZP and SF-CIS(D) calculations using SF-DFT/6-31G* geometries[68] provide reverse state ordering of the 1A1 and 1B1 states. However, the SF-OD method adiabatically correctly identifies the location of 1A1 and 1B1 states since a more dynamic correlation has been recovered by the SF-OD approach. We should also note that MR-CISD calculations with CASSCF(2,2) orbitals fail to rectify this error. However, when the larger CAS(4,4) is used, both CASSCF and MR-CISD methodologies put the 1A1 state above the 1B1 one. This fact reveals that the use of the CAS(4,4) is necessary for proper descriptions of the low-lying states and their energy gaps of TMM. The estimated energy gaps of multiplets depicted in Table 2 in ref (129) indicate that the minimum CAS size necessary for accuracy with TMM is (4,4). However, as shown earlier, the RS-based SSMRPT method with the IVO-CASC(2,2) reference function correctly identifies the location of both the singlet states relative to the 3B2 TMM.[132] To illustrate the sensitivity of the estimates for VEE to the methodologies employed for geometry optimization, IVO-BWMRPT values at SF-DFT/6-31G* geometries have been compared with those at IVO-CASCI[133] geometries. The vertical energy gaps so obtained by the IVO-BWMRPT method with different geometries are in proximity because the differences yielded are at most 0.04 eV. Perera et al.[95] also demonstrated that the effect of geometry is rather small on the computed ST splittings in TMM.

Aromatic σ, σ-Biradicals: ortho-, meta-, and para-Benzynes

In contrast to the TMM system, benzyne (didehydrobenzene) diradicals possess a closed-shell singlet [1A1 in C2ortho- and meta-benzynes and 1Ag in D2para-benzyne symmetries] rather than a triplet ground state.[134−136] These benzyne isomers have been scrutinized theoretically using various established ab initio methods.[68,93,94,132,137−147] The lowest excited state of the benzynes is a triplet state [3B2 in ortho- and meta-benzynes and 3B1u in p-benzyne]. As the energy gap between the frontier orbitals decreases in the ortho → meta → para sequence, the MR nature of the closed-shell ground state increases in the same direction. The distance between the unpaired electrons increases from the ortho- to para-isomer, resulting in a decrease in ST energy splitting in the same direction as determined in ultraviolet photoelectron spectroscopy experiments of Wenthold et al.[136] Therefore, the para-benzyne has the largest diradical (and MR) character. Consequently, the SR-based calculations cannot be relied upon in the case of para-benzyne.[94] Being on adjacent carbon atoms in ortho-benzyne, the interaction between the orbitals containing lone pair electrons is quite significant and leads to a generation of well-separated pairs of bonding and antibonding MOs (as that of strong covalent interactions).

The BWMRPT estimations of benzynes were carried out employing the cc-pVTZ basis and the optimized geometries (obtained at the SF-TDDFT/50-50/6-31G* level) as reported in ref (68). Here, a CAS of two electrons in two orbitals [CAS(2,2)] is employed. Computed adiabatic transition energies are presented in Tables , 7, and 8 for the ST splittings in ortho-, meta-, and para-benzyne, respectively. It should be mentioned that the knowledge of the magnitude of the ST splittings is important for the reactivity and specificity of benzynes as reaction intermediates. Here, it should be mentioned that SF-based values mentioned here have been estimated by Krylov and co-workers using the same geometry as employed by present calculations. The present work indicates that the ordering of ST splittings (eV) of the benzynes as provided by IVO-BWMRPT is ortho (1.701) > meta (0.819) > para (0.201). The corresponding values (eV) for CASSCF-BWMRPT are ortho (1.690) > meta (0.831) > para (0.224). Therefore, the BWMRPT AEEs computed with different reference functions are very close to each other for all of the three benzynes and BWMRPT calculations conclusively predict the adiabatic ST splitting to shrink when moving in the ortho → meta → para sequence as that of the well-accepted current generation MRPT methods.[143,144] The findings furnished by BWMRPT methods are also endorsed by the experimentally observed ordering of the diradical character of the benzynes. Wenthold et al.[136] experimentally found that the ordering for ΔEST (eV) was ortho (1.628) > meta (0.911) > para (0.165). However, Wenthold et al. also suggested (but less likely) an alternative value for the ST gap of the para-isomer, 0.091 eV. The parent method of BWMRPT, at the Mk-CCSD(2,2)/cc-pVDZ level of calculations,[140] predicted ordering of the ΔEST gap (eV) for the benzynes as ortho (1.522) > meta (0.811) > para (0.195). Therefore, IVO-BWMRPT provides larger ST splittings than the values obtained from the Mk-CCSD calculations for all three isomers, and BWMRPT calculations are also found to overestimate the experimental values. The use of CASSCF orbitals instead of IVO-CASCI marginally improves the agreement of BWMRPT scheme with the experimental values for o- and m-benzynes. However, in the case of p-isomer, use of CASSCF orbitals slightly increases the errors.

Table 6

Adiabatic Energy Gaps, ΔEST = ETriplet – ESinglet (eV), for ortho-Benzyne

reference	methods	ΔESTbenzyne
present work	IVO-BWMRPT(2,2)/cc-pVTZ	1.701
	CASSCF-BWMRPT(2,2)/cc-pVTZ	1.690
ref (68)	SF-CIS(D)/cc-pVTZ	0.842
	SF-OD/cc-pVTZ	1.632
ref (139)	SF-CCSD/cc-pVTZ	1.578
	SF-CCSD(fT)/cc-pVTZ	1.615
	SF-CCSD(dT)/cc-pVTZ	1.619
ref (100)	EOM-SF-CCSD/cc-pVTZ	1.578
	EOM-SF-CCSD(dT)/cc-pVTZ	1.619
ref (93)	CAS-BCCC4(2,2)/cc-pVDZ	1.466
ref (143)	CASPT2(8,8)/ANO	1.474
ref (144)	CASPT2(8,8)/ANO-L	1.531
	CASPT2(12,12)/ANO-L	1.592
ref (146)	CASPT2(8,8)/aANO	1.413
	CCSD(T)/cc-pVTZ	1.531
ref (141)	ic-MRCCSD(2,2)/cc-pVTZ	1.588
	ic-MRCCSD(T)(2,2)/cc-pVTZ	1.598
ref (94)	CR-CC(2,3)/cc-pVTZ	1.520
	RMRCCSD/cc-pVTZ	1.496
	RMRCCSD(T)/cc-pVTZ	1.539
	CCSD(T)/cc-pVTZ	1.522
ref (140)	Mk-MRCCSD(2,2)/cc-pVDZ	1.522
ref (136)	experiment	1.628 ± 0.013

Table 7

Adiabatic Energy Gaps, ΔEST = ETriplet – ESinglet (eV), for m-Benzyne

reference	methods	ΔESTbenzyne
present work	IVO-BWMRPT(2,2)/cc-pVTZ	0.819
	CASSCF-BWMRPT(2,2)/cc-pVTZ	0.831
ref (68)	SF-CIS(D)/cc-pVTZ	0.842
	SF-OD/cc-pVTZ	0.837
ref (139)	SF-CCSD/cc-pVTZ	0.782
	SF-CCSD(fT)/cc-pVTZ	0.875
	SF-CCSD(dT)/cc-pVTZ	0.892
ref (100)	EOM-SF-CCSD/cc-pVTZ	0.782
	EOM-SF-CCSD(dT)/cc-pVTZ	0.892
ref (93)	CAS-BCCC4(2,2)	0.750
ref (143)	CASPT2(8,8)/ANO	0.824
ref (144)	CASPT2(8,8)/ANO-L	0.880
	CASPT2(12,12)/ANO-L	0.880
ref (146)	CASPT2(8,8)/aANO	0.824
	CCSD(T)/cc-pVTZ	0.898
ref (94)	CCSD/cc-pVTZ	0.312
	CR-CC(2,3)/cc-pVTZ	0.711
	RMRCCSD/cc-pVTZ	0.702
	RMRCCSD(T)/cc-pVTZ	0.767
	CCSD(T)/cc-pVTZ	0.742
ref (140)	Mk-MRCCSD(2,2)/cc-pVDZ	0.811
ref (134)	experiment	0.911 ± 0.014

Table 8

Adiabatic Energy Gaps, ΔEST = ETriplet – ESinglet (eV), for p-Benzyne

reference	methods	ΔESTbenzyne
present work	IVO-BWMRPT(2,2)/cc-pVTZ	0.201
	CASSCF-BWMRPT(2,2)/cc-pVTZ	0.224
ref (68)	SF-CIS(D)/cc-pVTZ	0.014
	SF-OD/cc-pVTZ	0.092
ref (139)	SF-CCSD/cc-pVTZ	0.147
	SF-CCSD(fT)/cc-pVTZ	0.169
	SF-CCSD(dT)/cc-pVTZ	0.172
ref (100)	EOM-SF-CCSD/cc-pVTZ	0.147
	EOM-SF-CCSD(dT)/cc-pVTZ	0.172
ref (149)	CAS-BCCC4(2,2)/cc-pVDZ	0.130
ref (142)	CASPT2(2,2)/cc-pVTZ	0.195
	PC-NEVPT2(2,2)/cc-pVTZ	0.159
	CASPT2(8,8)/cc-pVTZ	0.248
	PC-NEVPT2(8,8)/cc-pVTZ	0.207
	DSRG-MRPT2(2,2)/cc-pVTZ	0.128
	DSRG-MRPT2(8,8)/cc-pVTZ	0.189
ref (143)	CASPT2(8,8)/ANO	0.143
ref (144)	CASPT2(8,8)/ANO-L	0.230
	CASPT2(12,12)/ANO-L	0.234
ref (146)	CASPT2(8,8)/aANO	0.252
	CCSD(T)/cc-pVTZ	0.009
ref (140)	Mk-MRCCSD(2,2)/cc-pVTZ	0.195
ref (141)	ic-MRCCSD(2,2)/cc-pVTZ	0.1616
	ic-MRCCSD(T)(2,2)/cc-pVTZ	0.2115
	MR-CISD+Q(2,2)/cc-pVDZ	0.116
	Mk-MRCCSD(T)(2,2)/cc-pVDZ	0.195
	DSRG-MRPT2(2.2)/cc-pvDZ	0.111
ref (140)	Mk-MRCCSD(2,2)/cc-pVDZ	0.195
ref (94)	CCSD/cc-pVTZ*	–0.830
	2R-MRCCSD/cc-pVTZ*	–0.091
	CR-CC(2,3)/cc-pVTZ*	–0.078
	2R-MRCCSD(SU-MRCCSD)/cc-pVTZ*	0.030
	2R-MRCCSD(T)/cc-pVTZ*	0.139
	CCSD(T)/cc-pVTZ*	0.134
ref (134)	experiment	0.165 ± 0.016

In spite of possessing MR nature, especially for meta- and para-benzynes, CCSD(T) performs rather well in the present situations. For other schemes that handle triples effectively, one can observe that CR-CC(2,3) also operates within “a safe zone” of quasidegeneracy in the case of o- and m-benzynes. It is important to note here that the CCSD-, CR-CC(2,3)-, and CCSD-based RMR separations have wrong signs (incorrectly predicting the triplet to be the ground state) for para-benzyne. However, the spin–flip version of the CCSD (SF-CCSD) scheme is free from such an objection.[139] In the para-isomer, SF-CCSD/cc-pVTZ provided the ST gap that is within 0.003 eV from the experimental data. ortho-Benzyne is the only isomer for which the CCSD ST energy gap is qualitatively correct. It should be mentioned that the CCSD scheme overestimates the bonding interactions between the two centers and provides an incorrect bicyclic structure for m-benzyne; however, incorporation of higher excitations, such as triples in CCSD, can, of course, solve the problem. The RMR protocol combined with state universal CCSD (SU-CCSD) with a modified cc-pVTZ basis set [where the f functions on carbons and the d functions on hydrogens have been deleted from the cc-pVTZ basis] underestimates the ST gaps significantly.[94] The Mk-MRCCSD(2,2)/cc-pVDZ ST gap corrected with (CBS + ZPVE) corrections for para-benzyne, 0.104 eV, is in closer agreement with the alternative value of 0.091 rather than the well-accepted one (0.165 eV).[140]

In o-, m-, and p-benzynes, IVO-BWMRPT (CASSCF-BWMRPT) calculations lead to overestimations (relative to the experimental data) of the adiabatic ST gaps by 0.073 (0.062), 0.092 (0.080), and 0.036 (0.059) eV. The absolute errors of the SF-OD/cc-pVTZ method against experimental values for ortho-, meta-, and para-benzynes are 0.004, 0.074, and 0.073 eV, respectively.[68] In all three isomers, the SF-CCSD(fT) and SF-CCSD(dT) adiabatic splitting gaps are within 0.04 and 0.02 eV, respectively, from the experimental estimates. The errors indicate that IVO-BWMRPT yields adiabatic EEs that are in line with the ones obtained using SF-CCSD with triples’ corrections and SF-OD methods. It should be noted that the connected triples’ correction enhances the adiabatic ST transition energies in all of the isomers, improving the accordance with the corresponding experimental findings for both ortho- and meta-benzynes. The performance of IVO-BWMRPT(2,2)/cc-pVTZ is more close to the values of the computationally demanding EOM-SF-CCSD(dT)/cc-pVTZ method compared to the EOM-SF-CCSD/cc-pVTZ one.[100] At this point, it is important to mention that the performance of SF protocols considered here can be plagued by spin contamination of the reference state, which can be removed or at least attenuated using restricted open-shell references. The present IVO-BWMRPT is a JM-ansatz-based explicit spin-adapted MR scheme. The quantitative agreement of IVO-BWMRPT performance with SA and SD variants of the MRMP2 method of Hirao and co-workers[89] is worth mentioning for ortho-benzyne. The ST energy gaps provided by MRMP(SA)/cc-pVTZ and MRMP(SD)/cc-pVTZ with CAS(8,8) are 1.47 and 1.78 eV, respectively, exhibiting deviations of around 0.16 eV for both methods.[89] The adiabatic ST splitting obtained with the previously reported CASPT2[143,144,146] that additionally included the six π-orbitals’ approach agrees well with the computationally less expensive IVO-BWMRPT(2,2) one. The CASPT2(8,8)/ANO[144] ST estimates deviate from the experimental values by 0.097, 0.031, and 0.065 eV for ortho-, meta-, and para-benzyne, respectively, which are slightly lower than the present results. For para-benzyne, the IVO-BWMRPT result is also in line with the recent predictions of DSRG-MRPT2(8,8) [driven similarity renormalization group-based second-order MRPT] and PC-NEVPT2(8,8)/cc-pVTZ calculations with the cc-pVTZ basis of Li and Evangelista.[142] For the adiabatic ST gap of p-benzyne, the predictions of DSRG-MRPT2(8,8) and PC-NEVPT2(8,8) are within 0.024 and 0.042 eV, respectively, from the experimental value. Note that their calculations suggest that both CASPT2 and PC-NEVPT2 methods with CAS(2,2) show satisfactory agreement with an alternative value of 0.091 rather than the well-accepted one (0.165 eV) for the p-isomer. Findings of Evangelista and co-workers[142] have shown that the selection of the active space plays a crucial role in perturbative methods. The present IVO-BWMRPT(2,2) calculations yield results that are also in line with the state-of-the-art MRCC estimates,[140,141] despite the second-order perturbative nature of the method. Overall, the IVO-BWMRPT results agree with the general trends predicted by the ic-MRCC(2,2)/cc-pVTZ scheme with both SD and SD(T) truncation strategies,[141] which demonstrates the suitability of the IVO-BWMRPT methodology. Li and Evangelista[142] demonstrated that the effect of CAS on the ic-MRCC theories is not significant. For example, the ST gap of p-benzyne suggested by the ic-MRCCSD(T)/cc-pVDZ scheme only differs by less than 0.009 eV between CAS(2,2) and CAS(8,8). Overall, one can find that the IVO-BWMRPT with minimal CAS can describe the adiabatic ST gaps of benzynes, with an accuracy similar to that of other high-level MR perturbative and nonperturbative theories.

The relative energies of the m- and p-benzynes relative to the o-isomer calculated using the IVO-BWMRP(2,2)/cc-pVTZ scheme along with other previously reported values obtained by various levels of the method are described in Table . One can find that the IVO-BWMRPT method underestimates the o/m and o/p energy splittings. As shown in Table , the IVO-BWMRPT method clearly predicts an energy ordering of o- < m- < p- <, which is in agreement with the findings of preexisting theoretical methods[145,147] and with Wenthold’s experimental relative energy values.[136] Although the CCSD/cc-pVTZ o/m energy gap is consistent with the experimentally and theoretically estimated data, the o/p energy splitting poses a greater challenge for the CCSD method. The accuracy of the IVO-BWMRPT scheme’s relative-energy predictions for both o/m and o/p compares to that of the CCSD(T) values[147] in terms of the magnitude of error. The IVO-BWMRPT(2,2)/cc-pVTZ estimates reported here are in good harmony with those obtained from the parametric two-electron reduced-density matrix (p2-RDM) with the cc-pVTZ basis.[147] Our results for relative energies of m- and p-benzynes relative to the o-isomer demonstrate that the performance of the IVO-BWMRPT(2,2) model is better than that of CASPT2. It can be seen that the agreement of relative energies provided by IVO-BWMRPT with CR-CC [0.659 and 1.613 eV for o/m and o/p, respectively], CCSD(T) [0.603 and 1.587 eV for o/m and o/p, respectively], and p2-RDM [0.698 and 1.158 eV for o/m and o/p, respectively] estimates in the CBS limit is good for both m- and p-isomers. Numerical analysis about the EEs of the systems treated here reflects that the IVO-BWMRPT is a potentially useful means to treat a variety of biradical systems, including those with a lesser degree of static correlation.

Table 9

Computed Relative Energies of the m- and p-Benzynes Relative to o-Benzyne, ΔEm/ (eV), Obtained Using a Variety of Methods and Experiment

reference	methods	ΔEm	ΔEp
present work	IVO-BWMRPT(2,2)/cc-pVTZ	0.558	1.123
ref (147)	CCSD/cc-pVTZ	0.598	1.860
	CR-CC/cc-pVTZ	0.650	1.579
	p2-RDM/cc-pVTZ	0.711	1.145
	CCSD(T)/cc-pVTZ	0.581	1.188
ref (145)	CASPT2/cc-pVDZ	0.438	0.941
ref (134)	experiment	0.663 ± 0.13	1.353 ± 0.13

Rather than resorting to a brute-force scheme, the correct representation of the low-lying electronic excited states including the ground state of the systems treated here is achieved when the relevant configurations are democratically treated and dynamical correlation is accounted for that can be accomplished using the present IVO-BWMRPT formalism with adequate model spaces. The main difficulty in modeling excited states is in quantifying varying degrees of dynamical and nondynamical correlations equally well among states of different nature. The IVO-BWMRPT method has the ability to describe not only the vertical excitation spectra but also adiabatic transitions, a very notable property for photochemistry applications because the BWMRPT method furnishes a balanced treatment of all of the low-lying excited states as they include different correlation effects in a balanced way. The change of the degree of MR nature among systems has implications for the congruence of the IVO-BWMRPT approach’s predictions with other theoretical estimates. When the extent of MR correlation is low or moderate, the IVO-BWMRPT energies are comparable to those of the recent-generation SR-based CC calculations. It is encouraging that the BWMRPT results for VEE are not dramatically different for different choices of orbitals. This suggests that the consistent and good performance of the BWMRPT is due to the inherent structural nature of the method and ansatz of the wave function and not the orbital choice. More studies need to be carried out on the zero-order Hamiltonian or on the basis set used to have an overall notable accordance of IVO-BWMRPT findings with high-quality current generation estimates for the systems considered here. As IVO-SSMRPT is capable of incorporating the dynamical correlation at only the second-order perturbative level, it cannot pretend to compete against computationally demanding high-level approaches or to be very close to the experimental data.

Conclusions

The efficacy of the recently suggested Brillouin–Wigner version of the SSMRPT (BWMRPT) method with the IVO-CASCI reference function (termed as IVO-BWMRPT), a Hilbert-space MR approach, has been examined to describe excitation energies for low-lying excited states of different characters. Although the estimation of energy surfaces has been mostly a prerogative of the Hilbert-space-based MR methodologies, albeit such protocols can also be used in computations of the energy differences of the electronic states of very different characters. In the formalism, each target state must be accessed individually, thus requiring a repeated solution of the IVO-BWMRPT equations. The systems chosen to assess the method are methylene, ethylene, TMM, and three isomers of benzyne, which possess strongly interacting electronic states and serve as popular benchmark systems. Due to the “diagonalize-then-perturb-then-diagonalize” nature of the formalism, the IVO-BWMRPT method is very effective to handle quasidegenerate states, in which the mixing emerges between the reference and other secondary CASCI states via the second-order interactions.

The excitation energies obtained for the systems considered here indicate that the use of the minimum CAS instead of optimal CAS may not always be sufficient for very accurate predictions of the energy gaps. Employing adequate model spaces, the IVO-BWMRPT formulation has been capable of furnishing accurate estimates of excitation energies of strongly correlated systems. The results obtained by IVO-BWMRPT are in agreement and consistent with those reported by other MRPT, MRCI, and expensive CC calculations for both the vertical and adiabatic excitations. In many situations, IVO-BWMRPT estimates are competitive with the high-level CC methods. The present work also reveals that the IVO-BWMRPT and high-level SF-based methods (provide efficient and robust computational tools for EEs) furnish very similar results. The IVO-BWMRPT method is useful not only for valence excitations but also for Rydberg excitations that involve transitions from localized to diffuse orbitals. The present estimates clearly advocate that the IVO-BWMRPT describes both open-shell and closed-shell singlet states with equal ease. Keeping in mind the limitations of the basis sets used, the BWMRPT results are in most cases reasonably close to the experimental values. This indicates its capability to quantitatively describe, even at the second-order perturbative level, the electronic properties of molecules. Due to its structural compactness, the IVO-BWMRPT method has the potential to quantify varying degrees of the active–active and active–virtual electron correlation effects, which are equally important among states of different characters in a balanced and accurate manner. This flexibility of the method in the description of electron correlation permits the IVO-BWMRPT to be nearly consistent with experimental findings. One can notice that the BWMRPT scheme with CASSCF orbitals (CASSCF-BWMRPT) is as accurate as IVO-BWMRPT (i.e., BW-SSMRPT with IVO-CASCI orbitals), which is extremely encouraging since replacement of the CASSCF step by IVO-CASCI reduces the computational cost while preserving the inherent quality of the desired results. Thus, one can hope that the present strategy can be exploited to more diverse MR situations/processes with an affordable computational cost. Finally, although MR-based correlated approaches have inherent advantages in various challenging chemical applications, futuristic perspectives of MR ab initio protocols are not easy to assess.