Multireference Configuration Interaction Methods – An Application to the Valence Isomerism in Cyclobutadieno-p-benzoquinone and its Diprotonated Form

Multireference averaged quadratic coupled cluster (MR-AQCC) calculations for cyclobutadieno-p-benzoquinone indicate that valence bond isomers 1a and 1b can exist as distinct species. The energy barrier height for their interconversion are 4.6 and 4.5 kcal mol−1, respectively, what is by ca. 2 kcal mol−1 lower than in the parent cyclobutadiene, implying that they could perhaps exist only under extreme conditions, namely at very low temperatures. For double protonated cyclobutadieno-p-benzoquinone, the CASSCF calculations erroneously predict existence of two valence isomers, 2a and 2b, whereas the MR-AQCC calculations reveal that geometry of the double protonated species could be best described by structure 2b. This nicely illustrates the crucial role of dynamic correlation and the need for using a highly-correlated theoretical method including geometry optimization in studied molecules.


INTRODUCTION
[3][4] A paradigmatic example of the bond stretch isomerism which includes the localized π-double bond shift mechanism as an underlying electronic effect is provided by cyclobutadiene (CBD), the smallest neutral antiaromatic organic compound, with a high angular strain which is reflected in its highly pronounced reactivity. [5][8][9] In a seminal experiment, Cram and coworkers have succeeded to isolate CBD in a hemicarcerand matrix but the crystal structure of cyclobutadiene could not been determined. [10]More recently, Barboiu and coworkers have claimed to identify two kinetically stable isomeric forms (rectangular and square) of 1,3-dimethylcyclobutadiene in the matrix at 175 K by X-ray spectroscopy, [11] which according to numerous calculations correspond to ground state and transition structure, respectively. [12,13]However, this finding initiated heated debate in literature on whether the measured crystal structure corresponds to 1,3-dimethylcyclobutadiene at al. [14][15][16][17] Valence bond isomerism of cyclobutadiene has also aroused considerable interest of theoretical chemists, to mention only the most recent high level computational studies by Balkova and Bartlett, [12] Sancho-Garcia et al., [18] Shen et al., [19] Lefrancois et al., [20] and our research group. [13]e have shown that the barrier for automerization reaction of the parent CBD molecule obtained at the MR-AQCC level of theory including extrapolations to the basis set limit is 6.3 kcal mol −1 . [13]This value is 2.2 kcal mol −1 above the two determinant TDCCSD(T) result of Balková and Bartlett. [12]The state specific multireference Brillouin-Wigner CC (MR BWCC) theory gives 4.5 kcal mol −1 , [18] B Croat.Chem.Acta 2015, 88 (4), 495-503 DOI: 10.5562/cca2770 whereas the recent block correlated coupled cluster CAS-BCCC4 method yields 6.2 kcal mol −1 . [19]Except for cyclobutadiene MR-AQCC calculations were also performed on derivatives of CBD, such as benzo [1,2:4,5]dicyclobutadiene, [21,22] and cyano substituted CBDs. [23]For all these molecules isomerization barrier compared to parent CBD was found to be somewhat lower, being around 5 kcal mol −1 .These indicates that the existence of valence bond isomers in these molecules could be possible only at very low temperatures due to the low energy barrier interconnecting two minima on the potential energy surface.[26][27] Similarly to cyclobutadiene, cyclobutadieno-p-benzoquinone has not been synthesized as yet, but a synthesis of naphthoquinone-fused cyclobutadiene ruthenium complexes [28] and dibenzoquinone cyclobutadiene derivative [29] has been reported.The first computational study of the π-bond shifting in cyclobutadieno-p-benzoquinone using HF/6-31G(d) and B3LYP/6-31G(d) method was reported by McKee et al. [24] They found that valence isomers 1a and 1b (Figure 1) are true minima at the potential energy surface, but due to single-determinant approach used in calculations, they were not able to locate the transition structure for their interconversion.
The same holds for recent MP2 and B3LYP calculations of Golas et al. [30] The only study carried at higher CASPT2//GVB2 level of theory which adequately included electron correlation, has revealed that structure 1b is not kinetically stable at the potential energy surface with isomerization barrier of only 0.3 kcal mol −1 . [25]It should be mention, however, that these results have been obtained by single point CASPT2 calculations which did not provide efficient geometry optimization of minima and transition structures.In order to remedy for this drawback we deemed it worthwile to reexamine possibility of existence of two valence bond isomers 1a and 1b by using more advanced multireference average quadratic coupled cluster method (MR-AQCC) [31][32][33] which was proved to provide reliable description of valence bond isomerism in our previous study of cyclobutadiene. [13]This will give us opportunity to compare the results with those of cylobutadiene obtained at the same level of theory and thus allow reliable estimate of annelation effect on the electronic structure and π-bond shifting in the parent molecule.We should mention at this point that practically identical results for automerization of cyclobutadiene were obtained using MR-CISD+Q method, [13] but for the sake of comparison with previously reported MR-AQCC results for variety of cyclobutadiene derivatives studied in our group, [21,23] only MR-AQCC calculations will be performed here.In addition to the study of π-bond shifting in the neutral cyclobutadieno-p-benzoquinone, we shall briefly address effect of diprotonation of the quinone ring on the overall geometry of molecule.A key question of interest here is whether Kekulé structures related to isomers 1a and 1b in the neutral molecule have their counterparts (2a and 2b) in the double protonated molecule (Figure 1).Special attention will be paid to impact of redistribution of electron density induced by protonation of the carbonyl groups on location of the double bonds in the cyclobutadiene moiety.

METHODOLOGY
Multiconfigurational and multireference computational methods are widely used in cases where nondynamical electron correlation effects, together with the ever-present dynamic electron correlation, represent an important contribution in total electronic wavefunction. [34]In case of standard single reference methods, such as Hartree-Fock, post-Hartree-Fock, and density functional based methods, a total electronic wavefunction is determined by a single Slater determinant and nondynamical electron correlation does not play an important role.In most of the routine quantum chemical calculations, this approximation is more than sufficient for an adequate description of electronic properties if dynamic electronic correlation effects are properly taken into account.Today, this is readily achieved by use of standard density functional theory (DFT) or post-Hartree-Fock methods, such as Møller-Plesset perturbation theory (MP2) or coupled cluster (CC) methods provided that sufficiently flexible basis sets of at least double-zeta quality are used.However, in special cases, such as bond creating/breaking, [35] excited states calculation [36] or quasidegenerative ground states, [37] the single reference approximation breaks down resulting in a qualitatively wrong description of a given electronic system.This occurs due to the inability of single reference methods to properly account for the nondynamic electron correlation arising from degenerate molecular orbitals which are always present in given situations.In this respect, a number of multi-configurational methods, based on a linear combination of selected degenerate molecular orbitals participating in a chemical reaction, were developed.One of the commonly used methods is the multi-configurational self-consistent field (MCSCF) method. [38]It uses a linear combination of several electronic configurations, each with its own coefficient, and by using a variational principle the total wavefunction with the lowest energy is obtained by simultaneous optimization of electronic configurations and their coefficients.In order to account for dynamical electron correlation effects, MCSCF wave function is often used as a reference state for the multireference configuration interaction (MRCI) [39] or complete active space perturbation theory (CASPT2) [40] calculations.In case of CASPT2 methods dynamical correlation is described by the perturbation theory of the multiconfigurational MCSCF wavefunction, similar to the MP2 method in single reference methods.On the other hand, in MRCI methods, the total wavefunction is obtained by all single and double electron excitations from the reference state to the virtual space limited by the size of a given basis set.Particularly valuable among them are analytic energy gradient methods, like the multireference averaged quadratic coupled cluster [31] (MR-AQCC) method which enables straightforward computation of analytic energy gradients, thus allowing an efficient optimization of the energy minima and transition-state structures on the PES and straightforward calculation of the ZPVEs.Consequently, it allows detailed inspection of chemical reactions along the reaction paths by treating both dynamic and nondynamic electron correlation on the same footing offering a unified description of widely different electronic structures occurring in the ground, transition and excited states of the same molecule.This feature is of the outmost importance in the study of the bond-stretch isomerization where the transition structure for isomerization cannot be properly described by the single reference methods. [13,21,23]47][48][49][50][51] Finally, it should be stressed that use of MCSCF, and especially MRCI methods, is often computationally too expensive which limits their application to relatively small molecules containing no more than 10−15 heavy atoms.

COMPUTATIONAL DETAILS
Complete active space self-consisted field (CASSCF) [40] and multi reference average-quadratic coupled cluster (MR-AQCC) [31][32][33] calculations were performed for the molecules presented in Figure 1.C2v symmetry restriction was always employed unless stated otherwise.All molecules are planar and lie in the yz plane.In the CASSCF calculations, the complete active space (CAS) wave function was comprised of 10 π electrons distributed over 10 π orbitals (1b1, 2b1, 3b1, 4b1, 5b1, 1a2, 2a2, 3a2, 4a2 and 5a2).The dynamic correlation energy was taken into the account by using MR-AQCC method.It should be pointed out that this modified MR-CISD approach includes size-extensivity corrections in an iterative way.Also, it incorporates analytical gradients [39,41,52] which provide the geometry optimization feasible at very high level of theory.The final expansion space in terms of configuration state functions (CSFs) for MR-AQCC calculations, was obtained by allowing all single and double excitations from all reference configurations into all internal and virtual orbitals.In the geometry optimizations, the symmetry of the reference configurations was restricted to the state symmetry, and the interacting space restriction was applied. [53]The ten core electrons of carbon and oxygen atoms have been kept frozen in all calculations.
It should be stressed that the use of full CAS(10,10) active space as a reference space would be too demanding for the present computational capabilities.Therefore, the construction of the final reference space was modified based on the natural orbital (NO) occupation (nocc) obtained at the CASSCF level of theory. [43]Namely, the 1b1 orbital has the largest orbital occupation (nocc = 1.944, 1.943 and 1.943 in structures 1a, 1b and transition structure 1TS, respectively).Therefore, this orbital was kept doubly occupied in all reference configurations.On the other hand, the orbital occupation of 5a2 orbital has a very small value (nocc = 0.045, 0.046 and 0.047 in structures 1a, 1b and transition structure 1TS, respectively) and was moved to the external space.The remaining eight orbitals (2b1, 3b1, 4b1, 5b1, 1a2, 2a2, 3a2 and 4a2) were also classified based on their NO occupation numbers (Table 1).
Thus, data in Table 1 show that two orbitals (2b1 and 1a2) in the neutral molecule have occupation number close to the value of 1.9 in all structures.Also, two orbitals (5b1 and 4a2) are very weakly occupied, with the occupation number close to 0. This enables grouping of orbitals 1a2 and 2b1 into restricted active space (RAS), orbitals 3b1, 4b1, 2a2 and 3a2 in the complete active space (CAS) and orbitals 5b1 and 4a2 in auxiliary space (AUX).Therefore, in addition to all possible CSFs obtained by distributing four electrons among CAS orbitals, all CSFs obtained by single and double excitations from RAS to CAS, and from CAS to AUX space were included in the reference space.The final expansion space in the MR-AQCC calculations consisted of a total number of 67.5 million CSFs obtained with the 6-31G(d) basis set.This method will be henceforth denoted as MR-AQCC(4*,4*).An analogous procedure was used to construct the active and reference space for molecules 2a and 2b and transition structure 2TS.
Zero point vibrational energy (ZPVE) corrections to the total energies of all studied molecules was computed at the CASSCF(10,10)/6-31G(d) level of theory because MR-AQCC method handles second derivatives calculations only numerically and harmonic force constants were computed by finite differences of the energy gradients.The harmonic vibrational frequencies were obtained by the SUSCAL program developed by Pulay et al. [54] The minima on the potential energy surfaces were verified by frequency analysis.All of the CASSCF and MR-AQCC calculations have been performed employing Pople's 6-31G(d) basis set [55] and using COLUMBUS suite of codes. [39,56,57]The atomic orbital (AO) integrals and AO gradient integrals have been computed with program module taken from DALTON. [58]Geometry optimizations were performed in natural internal coordinates as defined by Fogarasi et al. [59] by using the GDIIS method. [60]
Detailed analysis of bond lengths in isomers 1a and 1b shows that they could be classified as distinct valence bond isomers as predicted in earlier computational studies. [25]Focusing on the cyclobutadiene subunit, the bond length d(C2−C7) obtained at the MR-AQCC level of theory in isomer 1a is 1.546 Å whereas in isomer 1b it equals 1.362 Å.Bond lengths d(C7−C8) and d(C2−C3) in isomer 1a are 1.356 Å and 1.364 Å, respectively.The corresponding bond lengths in isomer 1b are 1.554 Å and 1.538 Å, respectively.Moreover, analysis of the transition structure 1TS reveals that bond lengths d(C2−C7), d(C7−C8) and d(C2−C3) are 1.445 Å, 1.442 Å, and 1.451 Å, respectively.This clearly indicates that cyclobutadiene subunit is localized in minimum structures and delocalized in transition structure (Figure 2), in analogy to trends observed in the previously studied parent cyclobutadiene and its derivatives. [13,23]It is also evident that annelation results in some angular strain in the benzoquinone ring, as indicated by deviation of the C2−C1−C6 angle from 120°.Specifically, its value in isomer 1a is 112.7°,111.1° in isomer 1b and 111.8° in transition structure 1TS.This feature is not surprising due to a spillover of some angular strain caused by fusion of the strained four-membered ring and was also documented in previous calculations. [25]We also notice that geometries obtained by  the CASSCF and MR-AQCC methods are in qualitative agreement with the largest deviations in the transition structure where MR-AQCC method predicts larger delocalization in the cyclobutadiene subunit.
In contrast to the neutral molecule CASSCF(10,10)/6-31G(d) and MR-AQCC(4*,4*)/CASSCF(10,10)/6-31G(d) calculations for the double protonated species yield contradictory results (Table 2).Namely, CASSCF calculations led to location of energy minima 2a and 2b and transition structure for their interconversion 2TS, while MR-AQCC method predicts that valence tautomerism is not operative here and that structure 2b best represents the molecule.Closer analysis of the calculated CASSCF geometries shows that only geometry of dication 2b resembles geometry of the related valence bond isomer (namely 1b) encountered in the neutral molecule.In structure 2a, the length of the bonds forming the cyclobutadiene ring are practically equal, while in the 2TS localization of the double bonds expected for isomer 2a are found.Both of these features strongly indicate that valence bond isomerism found in the parent molecule is not operative in the double protonated species (see also energetic properties reported below).It is of some interest to compare geometries of 1b and 2b.In the latter the C=O bond is elongated due to the π-electron depopulation in this bonding region as result of formation of O−H bond.Concomitantly, the presence of significant portion of the positive charge on the C1 atom results in the polarization of adjacent CC (C1−C2 and C1−C6, resp.)bonds, which as consequence become shortened, and to the opening of the C6−C1−C2 angle.As to the cyclobutadiene subunit, the annelated bond (C2−C3) undergoes only slight shortening, while the C7−C8 bond gets shorter by almost 0.06 Å.On the Table 2. Selected geometrical parameters (a) for molecules 1a, 1b and their corresponding transition structure 1TS obtained by using MR-AQCC(4*,4*)/CASSCF(10,10)/6-31G(d) and CASSCF(10,10)/6-31G(d) and levels of theory, respectively.Distances are given in Å, and angles are given in degrees  a) The numeration of atoms is shown in Figure 1.
other hand, the double bonds (C2−C7 and C3−C8), emanating from the sites of fusion become longer by ca 0.3 Å indicating that the four-membered ring becomes less antiaromatic.This is corroborated by calculating Mulliken partial atomic charges for the two species, which predict a withdrawal of electronic densities upon protonation from all ring carbon atoms, with the largest ones from the C1 and C4 (from 0.39 to 0.29, in 1b and 2b, respectively).Table 3 summarizes the weights of the leading electronic configurations φ1 and φ2 in isomers 1a, 1b and corresponding transition structure 1TS.It appears that the electronic configuration in 1a and 1b is (3b1)2(4b1)2(2a2)0 and (3b1)2(2a2)2(4b1)0, respectively, whereas in the transition state 1TS these two electronic configurations are mixed with a similar weight.Figure 3 schematically shows related orbitals of isomers 1a and 1b.
Analysis of the data in Table 4 shows that isomer 1a is by 4.1 kcal mol −1 more stable than isomer 1b at the CASSCF(10,10)/6-31G(d) level of theory, whereas at the MR-AQCC(4*,4*)/CASSCF(10,10)/6-31G(d) level of theory this difference is merely 0.3 kcal mol −1 , respectively.Zeropoint vibrational energies, calculated at the CASSCF(10,10) /6-31G(d) level of theory, differ only by 0.2 kcal mol −1 , thus affecting the relative stabilities of the isomers only slightly (Table 4).As expected, ZPVE for the transition structures are lower by ca. 2 kcal mol −1 which is a consequence of one inactive vibration in total zero-point vibrational energy.In the case of diprotonated isomers, total energies of minima and transition structure could be calculated only at the CASSCF(10,10)/6-31G(d) level of theory.At this level of theory structures 2a and 2TS were found to be by 11.8 kcal mol −1 less stable than 2b, whereas the barrier for isomerization from isomer 2a and transition structure 2TS drops to only 0.4 kcal mol −1 .Inclusion of ZPVEs leads to disappearance of the energy barrier and isomerization becomes spontaneous.This conclusion is further corroborated by calculations at the MR-AQCC(4*,4*)/CASSCF(10, 10)/6-31G(d) level of theory which predict that only one minimum (i.e.structure 2b) exists at the potential energy surface.Thus, we can safely conclude that CASSCF method fails to give reliable answer to the question whether   valence bond isomerism operates in protonated cyclobutadieno-p-benzoquinone.Inspection of Table 4 shows that the energy barrier for interconversion of isomer 1a to isomer 1b calculated at the CASSCF(10,10)/6-31G(d) level of theory is 6.7 kcal mol −1 , and 2.6 kcal mol −1 for interconversion of isomer 1b into isomer 1a, respectively.Inclusion of ZPVE correction diminishes energy barrier to 4.6 kcal mol −1 for the former, and only to 0.7 kcal mol −1 for the latter process.These results are close to previously reported CASSCF(10, 10)//GVB2 calculated energy barrier values of 5.9 of kcal mol −1 for the interconversion of isomer 1a to isomer 1b and 0.8 kcal mol −1 for the interconversion of isomer 1b to isomer 1a. [25]The situation is quite different at the MR-AQCC(4*,4*)/CASSCF (10,10)/6-31G(d) level of theory which gives an almost symmetric barrier with heights of 6.4 kcal mol −1 , (4.5 kcal mol −1 with ZPVE correction) and 6.7 kcal mol −1 (4.6 kcal mol −1 with ZPVE correction) relative to isomers 1b and 1a, respectively.It is noteworthy that these results are significantly different than those obtained employing CASPT2(10, 10)//GVB2 level of theory, [25] where the ZPVE corrected energy barrier of 5.3 kcal mol −1 for the interconversion of isomer 1a to 1b, and only 0.3 kcal mol −1 for the interconversion of isomer 1b to 1a were found.The main reason for this discrepancy is due to the fact that PT2 procedure describes dynamical correlation only to the second order and that results in Ref. [25] were obtained by single point calculations since geometry optimization procedures were not available for this computational level.Hence, the transition structure in the previous study [25] could not be precisely described.Thus the present results obtained by more advanced MR-AQCC(4*,4*)/CASSCF(10, 10)/6-31G(d) method provide more realistic description of the studied system strongly indicating that isomers 1a and 1b are practically isoenergetic and that they can exist as distinct molecules under appropriate conditions.
Comparison of the calculated interconversion energy barriers to the MR-AQCC benchmark value for cyclobutadiene automerization which is 6.3 kcal mol −1 , [13] reveals that annelation of cyclobutadiene to p-quinone ring lowers energy barriers for the interconversion of isomers 1a and 1b by ca. 2 kcal mol −1 (4.6 and 4.5 kcal mol −1 , respectively).This implies that isomers 1a and 1b are not kinetically stable at room temperature where they freely interconvert between two forms.However, their existence could be possible at very low temperatures as experimentally suggested for 1,3dimethylcyclobutadiene. [11] On the other hand, CASSCF and MR-AQCC calculations give contradictory results for diprotonated species.The CASSCF calculations erroneously predict existence of two valence isomers, 2a and 2b, whereas the MR-AQCC calculations predict that geometry of double protonated species could be best described by structure 2b.This nicely illustrates the crucial role of dynamic correlation and the need for using a highly-correlated theoretical method including geometry optimization in studied molecules.

CONCLUSIONS
Complete active space (CASSCF) and multireference average coupled cluster (MR-AQCC) calculations were performed in order to study valence bond isomerism in neutral and diprotonated cyclobutadieno-p-benzoquinones.Energy barriers with included zero-point energy correction obtained by MR-AQCC(4*,4*)/CASSCF(10,10)/6-31G(d) level of theory for valence bond isomerization of neutral molecules 1a to 1b and from 1b to 1a via transition state 1TS are 4.6 and 4.5 kcal mol −1 , respectively.Notably, energy barriers are lower by ca. 2 kcal mol −1 as compared to the MR-AQCC benchmark cyclobutadiene automerization energy barrier of 6.3 kcal mol −1 . [13]This finding is at variance with previously reported results of CASPT2(10, 10)//GVB2 calculations by Maksic and coworkers [25] who found that isomer 1b is not kinetically stable at the potential energy surface due to the small interconversion barrier of only 0.3 kcal mol −1 .It is strongly pointed out that theoretical treatment of the studied systems requires multireference methods capable of accurate reproduction of the Born-Oppenheimer potential energy surface (PES) along the isomerization pathways.

Table 4 .
CASSCF(a)and MR-AQCC (b) calculated energy barrier heights for the isomerization reaction of molecules 1a, 1b, 2a and 2b with and without ZPVE correction(c), respectively.All values are given in kcal mol−1