Hyesu
Jang
,
Yudong
Qiu
,
Marshall E.
Hutchings
,
Minh
Nguyen
,
Louise A.
Berben
and
Lee-Ping
Wang
*
Department of Chemistry, University of California, Davis, CA 95616, USA. E-mail: leeping@ucdavis.edu
First published on 29th January 2018
The CO_{2} reduction electrocatalyst [Fe_{4}N(CO)_{12}]^{−} (abbrev. 1^{−}) reduces CO_{2} to HCO_{2}^{−} in a two-electron, one-proton catalytic cycle. Here, we employ ab initio calculations to estimate the first two redox potentials of 1^{−} and explore the pathway of a side reaction involving CO dissociation from 1^{3−}. Using the BP86 density functional approximation, the redox potentials were computed with a root mean squared error of 0.15 V with respect to experimental data. High temperature Born–Oppenheimer molecular dynamics was employed to discover a reaction pathway of CO dissociation from 1^{3−} with a reaction energy of +10.6 kcal mol^{−1} and an activation energy of 18.8 kcal mol^{−1}; including harmonic free energy terms, this yields ΔG_{sep} = 1.4 kcal mol^{−1} for fully separated species and ΔG^{‡} = +17.4 kcal mol^{−1}, indicating CO dissociation is energetically accessible at ambient conditions. The analogous dissociation pathway from 1^{2−} has a reaction energy of 22.1 kcal mol^{−1} and an activation energy of 22.4 kcal mol^{−1} (ΔG_{sep} = 12.8 kcal mol^{−1}, ΔG^{‡} = +18.1 kcal mol^{−1}). Our computed harmonic vibrational analysis of [Fe_{4}N(CO)_{11}]^{3−} or 2^{3−} reveals a distinct CO-stretching peak red-shifted from the main CO-stretching band, pointing to a possible vibrational signature of dissociation. Multi-reference CASSCF calculations are used to check the assumptions of the density functional approximations that were used to obtain the majority of the results.
The discovery of CO_{2} reduction electrocatalysts represents a significant advance in CO_{2} utilization.^{1} Certain metallic electrodes have been reported to have catalytic activity for carbon dioxide reduction; Hori reported the formation of hydrocarbons and alcohols in the electrochemical reduction of carbon dioxide at copper electrodes in aqueous solution and discussed the reaction mechanism in 1989.^{2} In recent years, several metal and metal dichalcogenide nanostructured catalysts with high surface areas have been proposed as promising candidates for electrocatalysts for the CO_{2} reduction.^{3–8} In addition to the heterogeneous catalysts, a number of molecular catalysts have also been investigated for CO_{2} reduction and reviewed in several papers.^{9–11}
In 2011, Rail and Berben found that an Earth-abundant metal complex, first described by Muetterties and coworkers^{12,13} and denoted as [Fe_{4}N(CO)_{12}]^{−} or 1^{−} in its resting state, can act as a selective electrocatalyst for CO_{2} reduction to formate in aqueous solution.^{14} The preference of the catalyst for hydrogen evolution vs. CO_{2} reduction can be adjusted by tuning the strength of the acid used as a proton donor. An isoelectronic compound, [Fe_{4}C(CO)_{12}]^{2−}, was found to be a catalyst for hydrogen evolution only.^{15} In more recent work, Taheri and Berben further characterized the CO_{2} reduction mechanism and proposed the reduced hydride H-1^{−} as a key reaction intermediate.^{16} The hydricity, or hydride donor ability of H-1^{−} was proposed as a thermodynamic predictor of the selectivity for hydrogen evolution or CO_{2} reduction; a free energy window was proposed to explain the activity of 1^{−} for CO_{2} reduction, as opposed to its isoelectronic analogues.
[Fe_{4}N(CO)_{12}]^{−} is experimentally known to undergo two reduction events. When 1^{−} is electrochemically reduced to 1^{3−}, slow CO dissociation from the cluster is observed, resulting in the [Fe_{4}N(CO)_{11}]^{3−} or 2^{3−} species; the catalytic activity of this species is unknown but presumed to be inactive. In an accompanying experimental work, the X-ray crystal structure of 2^{3−} is reported.^{17} Simulations that uncover the mechanisms of side reactions are important to the overall strategy for designing molecular catalysts which are resistant to them. In this respect, this article describes the redox properties and CO dissociation pathway of this complex using computational quantum chemistry to complement the experimental findings and provide atomic-resolution insights.
The use of density functional theory (DFT) to study the electronic properties of metal carbonyl clusters has precedent in the literature. In particular, Schaefer and coworkers have produced a series of studies on the structures and metal–metal bonding of iron carbonyls and their derivatives.^{18–37} Several other groups have also carried out DFT studies for geometry optimization and vibrational frequency analysis of iron carbonyl complexes.^{38–41} Presumably, the strong fields from the CO ligands promote a low-spin and single-reference electronic state, making DFT a qualitatively appropriate method for studying these otherwise daunting multi-center inorganic clusters. Likewise, the application of DFT and solvent models for calculating redox potentials is well established.^{42–44} On the other hand, we are not aware of any theoretical studies that have investigated the redox properties and reactivity of 1; the significant metal–metal bonding and variation of charge states in this cluster may pose significant challenges for the density functional approximation and solvent model. For this reason, it is vital to compare calculated observables with experimental data where available.
In this theoretical study, we characterize the structures and energetics of the series of redox states: 1^{0}, 1^{−}, 1^{2−}, and 1^{3−}, and provide mechanistic insight into the CO dissociation side reaction: 1^{3−} → 2^{3−} + CO. Our calculations of the one-electron reduction potentials show close agreement with the experimentally measured values and provide some evidence that the BP86 density functional approximation^{45,46} performs more accurately for this system than the hybrid B3LYP functional.^{47} The dissociation pathway was found using high temperature ab initio molecular dynamics and relaxed to the minimum energy path to calculate the activation barrier.^{48} The calculations predict a structure of 2^{3−} in remarkable agreement with the X-ray crystal structure that was determined concurrently,^{17} lending further confidence to the level of theory used in this study. We also compare the CO dissociation barrier height to the analogous reaction after only one reduction event: 1^{2−} → 2^{2−} + CO, and show that dissociation from this electronic state is energetically uphill, though the activation free energy of CO dissociation is similar in both states. Our usage of DFT approximations is checked using natural orbital occupation numbers from multi-reference complete active space self-consistent field (CASSCF) calculations at key geometries.^{49,50}
G = G^{solv} + H_{SCF} + ZPE + H_{tr,rot,vib} − TS_{tr,rot,vib} |
ΔG = G_{red.} − G_{ox.} = −FE^{0} |
Geometry optimization was used to derive the self-consistent field (SCF) electronic energy together with the solvation free energy. Vibrational frequency calculations were used to derive the zero point energy and Gibbs free energy within the harmonic approximation. To calculate the relative redox potential, we took the differences of the free energies of the redox pairs and subtracted the absolute potential of the reference electrode, which is 4.67 V for the saturated calomel electrode (SCE). This value is based on the absolute potential of the NHE, which was determined by Reiss and Heller to be 4.43 V,^{60} although this quantity is difficult to measure and values in the range of 4.2–4.7 V have been reported in the literature.^{61}
We tested the dependence of the results on the choice of DFT approximation by performing calculations using three functionals: the BP86 gradient corrected semi-local functional,^{45,46} the B3LYP hybrid functional,^{47} and the PBE0 hybrid functional.^{62} Previous studies have noted that BP86 may perform more reliably than B3LYP in the study of the compounds in this paper.^{63,64} We also investigated whether adding diffuse basis functions affects the calculation results, because previous gas phase DFT studies suggest that diffuse basis functions are needed for the description of anions.^{65–67} For light elements (H, C, N, and O) we used the def2-TZVP triple-valence Gaussian basis set^{68,69} with f and higher angular momentum functions removed, denoted as def2-TZVP(-f). For the iron atoms we used either the LANL08 or LANL08+ basis set/ECP combination,^{70} which differ by the addition of a diffuse d angular momentum function in the latter. We further tested the effects of adding a minimal set of diffuse functions on light elements.^{71} The combined basis sets are called def2-TZVP(-f)-LTZ, def2-TZVP(-f)-LTZ+, and ma-def2-TZVP(-f)-LTZ+. Finally, since the cyclic voltammetry experiments to measure the redox potentials were carried out in a MeCN/H_{2}O (95:5) solvent, we also conducted the calculations employing the dielectric constants of water (78.4) and MeCN (37.5). From Tables 1 and 2, we concluded that the system has a minor dependence on the choice of basis set and solvent, while the functional dependence is significant. The BP86 functional gave closer agreement with experimental results (root-mean-square error, RMSE < 0.2 V) than the B3LYP and PBE0 hybrid functionals (RMSE > 0.4 V); the improved agreement is not due to shifting the absolute electrode potential, because the BP86 RMSE is still much lower than the other two functionals with the average gap subtracted out. Overall, the combination of the BP86 functional, the def2-TZVP(-f)-LTZ+ basis, and the dielectric constant of water yields good agreement with experimental data, with a RMSE of 0.15 V. The relatively high accuracy of BP86 (compared to hybrid functionals) for this system is consistent with previously published DFT studies of 3d transition metal containing complexes.^{63} To check the possible higher spin multiplicities for each state of the catalyst, we also calculated the energies of the higher spin multiplicities for each redox state (triplet and quintet for even-electron systems, quartet and sextet for odd-electron systems) and found that increasing the spin multiplicity significantly increases the total energy by over 10 kcal mol^{−1} (ESI Table S1†). From these findings we conclude that higher spin multiplicities do not participate in the redox chemistry and reaction pathways in this paper.
def2-TZVP(-f)_LTZ+/water (V) | ||||
---|---|---|---|---|
1 ^{0/1−} | 1 ^{1−/2−} | 1 ^{2−/3−} | RMSE | |
Exp | >0.2 | −1.23 | −1.60 | |
B3LYP | −0.13 | −1.37 | −2.07 | 0.42 |
PBE0 | −0.26 | −1.42 | −2.19 | 0.50 |
BP86 | +0.37 | −1.16 | −1.54 | 0.15 |
BP86/water (V) | ||||
---|---|---|---|---|
1 ^{0/1−} | 1 ^{1−/2−} | 1 ^{2−/3−} | RMSE | |
Exp | >0.20 | −1.23 | −1.60 | |
def2-TZVP(-f)_LTZ | +0.42 | −1.24 | −1.54 | 0.15 |
def2-TZVP(-f)_LTZ+ | +0.37 | −1.16 | −1.54 | 0.15 |
Ma-def2-TZVP(-f)_LTZ+ | +0.38 | −1.18 | −1.43 | 0.19 |
BP86/acetonitrile (V) | ||||
---|---|---|---|---|
1 ^{0/1−} | 1 ^{1−/2−} | 1 ^{2−/3−} | RMSE | |
Exp | >0.20 | −1.23 | −1.60 | |
def2-TZVP(-f)_LTZ | +0.41 | −1.30 | −1.62 | 0.14 |
def2-TZVP(-f)_LTZ+ | +0.34 | −1.22 | −1.64 | 0.13 |
Ma-def2-TZVP(-f)_LTZ+ | +0.39 | −1.27 | −1.55 | 0.15 |
The AIMD trajectories at elevated temperatures feature highly fluxional behavior of the CO ligands. Fig. 1 shows the all-atom root-mean-square deviation (RMSD) of the trajectory frames to the initial structure, and several trajectory snapshots of the simulation of H-1^{1−}. The RMSD rapidly reaches 1 Å after 1000 simulation steps (1 ps) and increases steadily over the course of ∼15 ps to almost 3 Å as larger geometric rearrangements take place. The conformational changes include the concerted rotation of multiple CO groups bonded to the same iron atom (analogous to torsion about a single bond), as well as the exchange of CO ligands on different iron atoms. At frame 22500, we observed a significant increase of the RMSD to >6 Å, where a CO ligand dissociated from the cluster. The distance between the dissociated CO and the catalyst molecule continued to increase until the simulation was terminated at frame 37000.
Fig. 2 summarizes the main results when the cluster is optimized in the −3 charge, singlet electronic state. The lowest energy structure (Fig. 2, left) is close to Cs-symmetric with a single mirror plane; the Fe atoms surround the central N in an isosceles trigonal pyramidal arrangement. The central N is nearly in the plane made by three iron atoms, with three in-plane Fe–N–Fe angles of 137, 137, and 84 degrees, summing up to 358 degrees. The other three Fe–N–Fe angles are between 85 and 90 degrees. Each Fe atom has three CO ligands with a tight distribution of Fe–C distances ranging from 1.74–1.77 Å; the ab initio bond order (BO) indices computed using Mayer’s method^{80} range from 1.05 to 1.25, indicating single bond order.
The lowest-energy structure with a dissociated CO ligand (Fig. 2, right) features two CO ligands bridging a pair of Fe atoms. The three Fe atoms, central nitrogen, and two bridging carbons form nearly a planar rectangle, with Fe–N–Fe and C–Fe–C angles of 168 and 174 degrees, respectively. The cluster is also nearly Cs-symmetric with a single mirror plane. Moreover, the bridging CO ligands have significantly larger Fe–C distances of 1.83 Å (left and right edges of rectangle) and 2.08 Å (bottom edge). The increased lengths of the Fe–C bonds along the bottom edge of the rectangle suggest that they possess a different electronic character; indeed, these two bonds have ab initio bond orders of 0.55, which are almost exactly half of the others. Our interpretation is that the C–Fe–C is a three-center two-electron bond, which compensates for the two σ-electrons that are lost in the dissociation process. To support this interpretation, Fig. 3 shows a doubly-occupied CASSCF (8,8) optimized molecular orbital that shows significant electron delocalization across the C–Fe–C bond; this is the only orbital we observed that possesses bonding character for these atoms. A comparison of the optimized structure with the experimentally determined X-ray crystal structure^{17} revealed an excellent agreement of 0.13 Å RMSD, lending confidence to the accuracy of the theoretical methods used; the calculations were performed without knowledge of the crystal structure, and the comparison was only performed later. The experimental crystal structure also contains three Na^{+} counterions that further stabilize the 2^{3−} structure; these were not included in the present calculations.
To assess the possibility that 2 may be a catalyst for CO_{2} reduction, we computed redox potentials of the 2^{0}/2^{−}, 2^{−}/2^{2−}, and 2^{2−}/2^{3−} couples in analogy to 1. Our computed potentials are +0.30, −0.45, and −1.05 V vs. SCE, respectively. Because all of these potentials are more positive than the applied potential for electrocatalysis, we do not think these species are participating redox intermediates in the main CO_{2} reduction reaction.
The transition state estimate from the NEB is input into a calculation of the ab initio Hessian, followed by a geometry optimization to precisely locate the transition state structure; we then verify, using a second Hessian calculation, that the optimized structure has only one imaginary vibrational frequency. Finally, an intrinsic reaction coordinate (IRC) calculation follows the energy downhill in the forward and backward directions of the imaginary vibrational mode to provide a continuous path connecting the TS and energy minima. The TS optimization and IRC calculations were performed in the Q-Chem software package^{81}, with input parameters that reproduce the TeraChem total energies to within 0.0005 a.u. (<0.5 kcal mol^{−1}). We provide harmonic free energy corrections at the transition state and final geometries, as well as reaction energies and free energies with the dissociated species completely separated. Our main results in this section are summarized in Fig. 4 and the first four rows of Table 3.
Structure | 1 ^{3−} → 2^{3−} + CO | 1 ^{2−} → 2^{2−} + CO | ||
---|---|---|---|---|
ΔE | ΔG | ΔE | ΔG | |
Optimized IRC at BP86/6-31G*-LDZ | ||||
Initial | 0.0 | 0.0 | 0.0 | 0.0 |
TS | 18.8 | 17.4 | 22.4 | 18.1 |
Final | 10.6 | 7.2 | 22.1 | 16.3 |
Separated | 13.3 | 1.4 | 25.2 | 12.8 |
Optimized IRC at M06-L/6-31G*-LDZ | ||||
Initial | 0.0 | 0.0 | 0.0 | 0.0 |
TS | 20.1 | 18.3 | — | — |
Final | 11.0 | 7.6 | — | — |
Separated | 15.9 | 3.6 | 22.9 | 12.9 |
Optimized IRC at B3LYP/6-31G*-LDZ | ||||
Initial | 0.0 | 0.0 | 0.0 | 0.0 |
TS | 25.3 | 23.1 | 25.2 | 16.5 |
Final | 12.2 | 7.8 | 24.9 | 14.7 |
Separated | 14.8 | 3.3 | 27.2 | 11.5 |
Optimized IRC at M06/6-31G*-LDZ | ||||
Initial | 0.0 | 0.0 | 0.0 | 0.0 |
TS | 23.5 | 22.5 | — | — |
Final | 11.1 | 9.9 | — | — |
Separated | 15.6 | 4.2 | 21.9 | 10.0 |
The blue curve in Fig. 4 shows the total energy for CO dissociation from 1^{3−} along the BP86/6-31G*-LDZ optimized reaction coordinate. The first part of the path involves a torsional motion of six CO ligands, allowing the two highlighted Fe–C distances to come into closer contact. An intermediate is found with a relative energy of ΔE_{1} = +14.5 kcal mol^{−1} and an activation barrier of E_{a1} = +15.5 kcal mol^{−1}; the structure contains an additional Fe–C bond (distance = 2.09 Å; BO = 0.63). The second transition state has energy E_{a} = +18.8 kcal mol^{−1} (ΔG^{‡} = 17.4 kcal mol^{−1}) and has the CO ligand beginning to dissociate from the cluster; this is followed by a relatively flat energy basin where the two newly formed Fe–C bonds (the three-center bond) become equal in length. The final CO-dissociated structure gives a reaction energy ΔE = +10.6 kcal mol^{−1} (ΔG = +7.2 kcal mol^{−1}). We also computed the reaction energy by treating the products as completely separate species and obtained ΔE_{sep} = 13.3 kcal mol^{−1} (ΔG_{sep} = +1.4 kcal mol^{−1}). The higher value of ΔE_{sep} is attributed to dissociating intramolecular interactions, and the lower value of ΔG_{sep} to the translational and rotational entropy of separated dissociation products. The slightly uphill ΔG and moderate ΔG^{‡} values indicate that this mechanism may be operative for forming the experimentally observed 2^{3−} species.
We also investigated CO dissociation from the 1^{2−} electronic state; because dissociation is not observed from 1^{2−} in the experimental studies, we presume that the calculated thermodynamic and/or kinetic parameters should be less favourable compared to 1^{3−}. In searching for the reaction energies and activation barriers for the 1^{2−} state, we proceeded from the same initial structures from the AIMD trajectory; the charge and spin multiplicity were set to −2 and 2, respectively. Our BP86 calculations found an uphill and nearly barrierless dissociation pathway (orange curve in Fig. 4) with ΔE = 22.1 kcal mol^{−1} and E_{a} = 22.5 kcal mol^{−1} (ΔG = +16.6 kcal mol^{−1}; ΔG^{‡} = 18.1 kcal mol^{−1}). The reaction energy calculated using separated species as the products is ΔE = 25.2 kcal mol^{−1} (ΔG = 12.8 kcal mol^{−1}).
Comparison of the dissociation pathways from 1^{3−}vs.1^{2−} gives reaction free energies of ΔG = +7.2 vs. +16.3 kcal mol^{−1}; with separated product species, ΔG_{sep} = +1.4 vs. +12.8 kcal mol^{−1}. These values indicate that CO dissociation from 1^{2−} is thermodynamically less favourable than from 1^{3−}, consistent with the experimental findings. On the other hand, although the energy barrier for 1^{3−} is lower than for 1^{2−} (ΔE = 18.8 vs. 22.4 kcal mol^{−1}), the calculated activation free energies are nearly the same (ΔG^{‡} = 17.4 vs. 18.1 kcal mol^{−1}). Comparison of the overall shape of the dissociation curve shows some other important differences: whereas the 1^{3−} pathway has two clearly defined barriers and an intermediate, the 1^{2−} pathway is nearly barrierless, which indicates that tunnelling effects may play a significant role in determining the reaction rate.^{82} In summary, CO dissociation from 1^{2−} is found to be thermodynamically less favourable, but more detailed reaction rate and free energy calculations may be needed to accurately compare the kinetics of these two pathways.
In all of our results, we found that increasing the basis set size has a relatively small effect. In the BP86/TZVP-LTZ calculations of CO dissociation from 1^{3−}, ΔE is essentially unchanged from the 6-31G*-LDZ result (10.6 kcal mol^{−1}); E_{a} is slightly lower at 18.3 kcal mol^{−1}. For the 1^{2−} pathway, BP86/TZVP-LTZ predicts a slightly higher value of ΔE = 22.7 kcal mol^{−1}, and there is no energy maximum on the pathway; this is perhaps not surprising, given the nearly barrierless dissociation curve. In the B3LYP/TZVP-LTZ calculations, the ΔE and E_{a} values changed by <1 kcal mol^{−1} compared to the corresponding B3LYP/6-31G*-LDZ values. The choice of DFT functional has a more significant impact. B3LYP/6-31G*-LDZ predicts ΔE = 12.2 kcal mol^{−1} and E_{a} = 25.3 kcal mol^{−1} for CO dissociation from 1^{3−}; notably, E_{a} is 6 kcal mol^{−1} higher than in BP86. Despite differences in the barrier height, the structures along the 1^{3−} IRCs are highly similar for both functionals, as evidenced by the B3LYP single-point calculations along the BP86 optimized pathway and vice versa (ESI Table S2†).
The most significant DFT functional dependence is seen in the 1^{2−} dissociation pathway. For the reactant (1^{2−}) structure, B3LYP predicts a pyramidal structure with an isosceles triangular base, almost identical to the structure of 1^{3−} in Fig. 2, left. On the other hand, BP86 predicts a “crooked butterfly” structure (Fig. 5) that is closer to the 1^{−} resting state; the largest Fe–N–Fe angle is 165 degrees, and one of the Fe–Fe distances is elongated to 3.01 Å (the other Fe–Fe distances are between 2.55 and 2.65 Å). These structures are only stable on the potential surfaces of their respective functionals, as a BP86 optimization started from the B3LYP-optimized structure leads to the BP86 minimum and vice versa. Clearly, a more objective measure is needed to determine which DFT approximation is more appropriate for this system.
The differences in BP86 vs. B3LYP in the 1^{2−} state originate from the electronic character of the ground state Kohn–Sham (KS) wavefunction. We computed the expectation value of the squared total spin operator 〈S^{2}〉 to measure any deviations of the KS wavefunction from a pure doublet (ESI Fig. S5†). Along the BP86 pathway, the 〈S^{2}〉 value of the BP86 KS wavefunction is stable around 0.77, close to the value of 0.75 for a pure doublet; on the other hand, the B3LYP wavefunction has higher 〈S^{2}〉 values, ranging from 0.84 to 1.08, indicating a higher degree of spin contamination. The spin contamination is even greater along the B3LYP IRC, where the B3LYP wavefunction has a 〈S^{2}〉 value close to 2.0 at the dissociated state. BP86 also predicts 〈S^{2}〉 values around 1.6–1.7 for these structures, indicating a broken symmetry KS wavefunction, containing more than one unpaired electron.
The significant spin contamination along the B3LYP IRC for 1^{2−} points to a multi-reference ground state that is not well described by a KS determinant. To investigate this further, we carried out single-point CASSCF calculations at the initial, TS, and final geometries along the CO dissociation pathway for both the 1^{3−} and 1^{2−} IRCs calculated using BP86 and B3LYP. These calculations employ the same 6-31G*-LDZ basis set as the DFT calculations, and active spaces of (4,6) and (3,6) were used for all states from the 1^{3−} and 1^{2−} state pathways, respectively. These calculations were carried out in the ORCA software package.^{89,90} CASSCF-optimized molecular orbitals for all 12 critical points are provided in the ESI.†
The optimized CASSCF molecular orbitals are very close to the natural orbitals that diagonalize the density matrix; the eigenvalues are within 10^{−4} of the diagonal elements, and off-diagonal elements are all <10^{−4}. The natural orbital occupation numbers for initial, TS, and final structures optimized using B3LYP and BP86 are plotted in Fig. 6; the more the occupation numbers deviate from 2.0 and 0.0 (for occupied and virtual orbitals), the greater the multi-reference character. Our analysis for 1^{3−} shows that the natural orbitals at the “frontier” have occupation numbers in the range of 1.9–1.7 and 0.1–0.2. The variations in these values are small when comparing the initial, TS, and final structures, indicating that there is no qualitative change in the electronic character along the reaction pathway. Moreover, none of the natural orbitals have occupation numbers near 1.0, which is a hallmark of wavefunctions that display strong multi-reference character; this is the case for diradicals and the homolytic dissociation of N_{2}.^{49}
For CO dissociation from 1^{2−}, the CASSCF calculations using BP86-optimized structures show a similar pattern to 1^{3−}, except a singly-occupied molecular orbital is present. On the other hand, a major change in the electronic character is seen for the B3LYP-optimized structures. The TS and final structures have occupation numbers close to 1.0 in three orbitals, indicating strong ground state multi-reference character; this result agrees with the spin contamination observed in the DFT wavefunctions for the same structures. When comparing the BP86 and B3LYP functionals, only the BP86-optimized structures have CASSCF ground states with consistent electronic character; we thus conclude that BP86 gives the more reliable result overall.
We also calculated reaction energies and activation energies of the reactions using the M06 and M06-L functionals to confirm the accuracy of BP86 for this system (Table 3). These calculations were performed in Q-Chem 5.0. We could not find a TS structure for CO dissociation from 1^{2−} using these functionals, again possibly owing to the nearly barrierless dissociation curve. The M06-L results are in close agreement with BP86, which is reasonable given that both functionals contain no Hartree–Fock (HF) exchange; spin contamination along the BP86-optimized 1^{2−} dissociation pathway is low, with 〈S^{2}〉 = 0.79–0.80. The M06 results are closer to B3LYP, perhaps because both functionals contain a similar amount of HF exchange (28% vs. 20%). M06 also shows similar amounts of spin contamination to B3LYP along both the BP86-optimized and B3LYP-optimized 1^{2−} dissociation pathways.
Footnote |
† Electronic supplementary information (ESI) available: XYZ coordinates of energy minimized and transition state structures reported in this paper, as well as wavefunction analyses of CO dissociation from 1^{2−}. See DOI: 10.1039/c7sc04342b |
This journal is © The Royal Society of Chemistry 2018 |