Is There a Quadruple Fe-C Bond in FeC(CO) 3 ?

: A recent computational paper (Kalita et al., Phys. Chem. Chem. Phys. 2020 , 22 , 24178–24180) reports the existence of a quadruple bond between a carbon and an iron atom in the FeC(CO) 3 molecule. In this communication, we perform several computations on the same system, using both density functional theory and post-Hartree–Fock methods and ﬁnd that the results, and in particular the Fe-C bond length and stretching frequency depend strongly on the method used. We ascribe this behavior to a strong multireference character of the FeC(CO) 3 ground state, which explains the non-conclusive results obtained with single-reference methods. We therefore conclude that, while the existence of a Fe-C quadruple bond is not disproved, further investigation is required before a conclusion can be drawn.

In a recent paper [1], a computational investigation on FeC(CO) 3 revealed the unprecedented result that in such a molecule, a quadruple Fe-C bond may be present. Unfortunately, FeC(CO) 3 has never been synthetized, and therefore no experimental characterization exists to confirm such a result: a computational investigation is therefore, at the moment, the only tool available to investigate such a system. In their work, Kalita et al. use the Fe-C bond length and the associated stretching force constant to assess the existence of a quadruple bond. They use Density Functional Theory (DFT), and in particular the M062X exchange correlation functional together with the Def2-TZVP basis set to perform such an analysis. As DFT calculations may be unreliable for multireference systems, they perform CASSCF and CCSD(T) single point calculations on the DFT geometry to confirm the single-reference character of the wavefunction, and thus the reliability of the DFT results.
Given the exceptionality of the results found by Kalita et al., we decided to perform further calculations on FeC(CO) 3 to confirm their results. We performed three different analyses. First, we compared the Fe-C bond distance and harmonic stretching frequency computed with different DFT exchange-correlation functionals to assess the robustness of the results with respect to the latter choice. In particular, we compare results obtained using five popular hybrid functionals: M062X [2], B3LYP [3], MN15 [4], PBE0 [5], and ωB97X-D [6]. All the DFT calculations were run using the Gaussian 16 [7] suite of programs, employing Dunning's cc-pVTZ basis set [8,9], using analytical first and second derivatives and a pruned (99,590) grid for the integration of the exchange-correlation potential. The DFT results are reported in Table 1.
The DFT results show a large variability, with bond distances ranging from 1.481 Å to 1.526 Å and harmonic stretching frequencies ranging from 1113 to 1211 wavenumbers, with the M062X functional predicting the shortest and strongest bond and the B3LYP one the longest and weakest. We therefore decided to extend our computational investigation by performing, to the best of our ability, high-level post-HF calculations, increasing the accuracy of the methods systematically. We performed calculations at the Hartree-Fock (HF), second-order Møller-Plesset perturbation theory [10] (MP2), coupled cluster theory with single and double excitations [11] (CCSD), and with single, doubles, and perturbative triples [12] (CCSD(T)), using Dunning's cc-pVTZ basis set and the frozen-core approximation. Analytical gradients [13][14][15] were used for all the geometry optimizations, while harmonic frequencies were computed by numerically differentiating analytical gradients. All the post-HF calculations were performed using the CFOUR [16,17] suite of programs. The results are reported in Table 2. Table 1. Equilibrium Fe-C distance (R e , in Å) and harmonic vibrational Fe-C stretching frequencies (in cm −1 ) computed with various DFT exchange-correlation functionals. The values reported with a "*" are taken from Reference [1] and reported here for comparison. The harmonic frequencies obtained at every level of theory are all real, confirming thus that the structure is a minimum. Unfortunately, it is apparent from the results that the treatment of electronic correlation that we could achieve is not sufficient to obtain converging results. The MP2 method fails to give a correct description of the system, with a very long bond distance and an unphysically very low stretching frequency. Going from MP2 to CCSD and CCSD(T), the bond distance seems to approach convergence, but the same cannot be said for the harmonic frequency. It is interesting to note that the CC bond length is significantly longer than the DFT one. Furthermore, while the fluctuations in the computed harmonic frequencies leave us unable to give a quantitative estimate, the converged value is likely to be in between the CCSD and CCSD(T) result, and in any case lower than all the DFT predictions. Unfortunately, the CCSD(T) calculation is already stretching the limits of computational feasibility: a more accurate treatment of correlation, including full triples excitations or even quadruples ones, seems at the moment out of reach.

Method
The large oscillations observed in the results motivated us to investigate whether the source of the lack of robustness of both the DFT and post-HF results stems from static correlation. We therefore performed further complete active space self-consistent field (CASSCF) calculations [18,19]. Using Dunning's cc-pVDZ basis set, which is sufficient for qualitatively correct results, and the CFOUR suite of programs, we optimized the geometry of FeC(CO) 3 using analytical gradients and computed harmonic frequencies, again by numerically differentiating analytical gradients. CASSCF calculations can be an important source of insight in the electronic structure of a molecule; however, the results depend strongly on the choice of the active space, which therefore introduces some elements of arbitrariness. The choice of the active orbitals can be motivated by chemical considerations, as done by Kalita et al. [1], who include in the active space all the orbitals involved in Fe-C bonds, resulting in a CASSCF calculation with 12 correlated electrons in 12 orbitals. Here, we adopt a different strategy, based on an automatic, a priori strategy to select the active orbitals. In particular, we select the active space according to the Unrestricted-Natural-Orbitals (UNO) criterion [20,21]. In this procedure, symmetry-broken unrestricted HF solutions are first found. An average density matrix is computed using all the independent UHF solution and then diagonalized to obtain the UNO. The active space is built by including all the UNO with occupation between 0.01 and 1.99. The Restricted HF wavefunction of FeC(CO) 3 presents six triplet instabilities; the resulting active space consists in 10 electrons in 10 orbitals, whose representation is given in Figure 1. The active space mostly involves the d orbitals of the metal and the p orbitals of the investigated C atom. Looking at the orbitals and at their occupation numbers, obtained by diagonalizing the final CASSCF one-body density matrix at the CASSCF equilibrium geometry, we note that the system definitely exhibits a relevant static correlation. Looking closer at the CASSCF wavefunction, we find that the contribution of the HF determinant has a weight of just 0.64, confirming the multireference nature of the molecule. This is in contrast to what was found by Kalita et al., who found a CASSCF wavefunction dominated by the HF reference (they report a coefficient of 0.952). We believe that this contradiction could be due to the fact that their calculation was performed on the DFT-optimized geometry, while we relaxed the geometry. Furthermore, our choice of the active space, resulting from a well-established automated procedure, is different, and in particular slightly more compact. While this may accentuate the multireference character of our CASSCF description, we believe that our choice of active space is robust, which is confirmed empirically by the smooth convergence of the CASSCF wavefunction optimization, which is performed using a second-order method. The CASSCF calculations predict a Fe-C distance of 1.573 Å and an associated harmonic stretching frequency of 1016 cm −1 , which is in qualitative agreement with the CC results.
We believe that our results show that it is too soon to conclude that there exists a quadruple Fe-C bond in FeC(CO) 3 . The evanescent nature of the concept of bond order itself has stimulated many discussions in the literature [24,25], which makes the topic particularly interesting. However, assessing something as unprecedented as the existence of a quadruple Fe-C bond requires a very careful validation of the results. Therefore, we conclude that the possibility of such a quadruple bond can, at the moment, be neither confirmed nor denied, and further, more accurate calculations are required.