Quantum chemical studies of redox properties and conformational changes of a four-center iron CO₂ reduction electrocatalyst

Abstract

The CO₂ reduction electrocatalyst [Fe₄N(CO)₁₂]⁻ (abbrev. 1⁻) reduces CO₂ to HCO₂⁻ 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³⁻. 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³⁻ with a reaction energy of +10.6 kcal mol⁻¹ and an activation energy of 18.8 kcal mol⁻¹; including harmonic free energy terms, this yields ΔG_sep = 1.4 kcal mol⁻¹ for fully separated species and ΔG^‡ = +17.4 kcal mol⁻¹, indicating CO dissociation is energetically accessible at ambient conditions. The analogous dissociation pathway from 1²⁻ has a reaction energy of 22.1 kcal mol⁻¹ and an activation energy of 22.4 kcal mol⁻¹ (ΔG_sep = 12.8 kcal mol⁻¹, ΔG^‡ = +18.1 kcal mol⁻¹). Our computed harmonic vibrational analysis of [Fe₄N(CO)₁₁]³⁻ or 2³⁻ 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.

---

Introduction

The development of economically viable technologies for reducing the CO₂ concentration in Earth’s atmosphere is one of the global environmental problems that we must solve in the near future. One of the major research fields in modern chemistry is the development of CO₂ capture, utilization, and storage strategies. Electrochemical CO₂ reduction has been studied as a CO₂ utilization technique, which can provide us with the possibility to produce useful products from CO₂.

The discovery of CO₂ reduction electrocatalysts represents a significant advance in CO₂ utilization.¹ 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.² 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₂ reduction.^3–8 In addition to the heterogeneous catalysts, a number of molecular catalysts have also been investigated for CO₂ 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₄N(CO)₁₂]⁻ or 1⁻ in its resting state, can act as a selective electrocatalyst for CO₂ reduction to formate in aqueous solution.¹⁴ The preference of the catalyst for hydrogen evolution vs. CO₂ reduction can be adjusted by tuning the strength of the acid used as a proton donor. An isoelectronic compound, [Fe₄C(CO)₁₂]²⁻, was found to be a catalyst for hydrogen evolution only.¹⁵ In more recent work, Taheri and Berben further characterized the CO₂ reduction mechanism and proposed the reduced hydride H-1⁻ as a key reaction intermediate.¹⁶ The hydricity, or hydride donor ability of H-1⁻ was proposed as a thermodynamic predictor of the selectivity for hydrogen evolution or CO₂ reduction; a free energy window was proposed to explain the activity of 1⁻ for CO₂ reduction, as opposed to its isoelectronic analogues.

[Fe₄N(CO)₁₂]⁻ is experimentally known to undergo two reduction events. When 1⁻ is electrochemically reduced to 1³⁻, slow CO dissociation from the cluster is observed, resulting in the [Fe₄N(CO)₁₁]³⁻ or 2³⁻ 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³⁻ is reported.¹⁷ 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⁰, 1⁻, 1²⁻, and 1³⁻, and provide mechanistic insight into the CO dissociation side reaction: 1³⁻ → 2³⁻ + 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.⁴⁷ The dissociation pathway was found using high temperature ab initio molecular dynamics and relaxed to the minimum energy path to calculate the activation barrier.⁴⁸ The calculations predict a structure of 2³⁻ in remarkable agreement with the X-ray crystal structure that was determined concurrently,¹⁷ 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²⁻ + 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

Computational methods and results

1. Redox potential calculations

We evaluated the relative free energies between the redox intermediates 1⁰, 1⁻, 1²⁻, and 1³⁻ using unrestricted Kohn–Sham DFT⁵¹ with the implicit solvent environment, conductor-like screening model (COSMO)⁵² for comparison to experimentally determined redox potentials. These calculations were carried out using the TeraChem software, which uses graphics processing units to accelerate the computation of the Coulomb and exchange matrices,^53–55 effective core potentials (ECPs),^56,57 and solvent response⁵⁸ that appear in the SCF calculation. A recently developed geometry optimization method using translation-rotation internal coordinates was employed to accelerate the energy minimization calculations.⁵⁹

G = G^solv + H_SCF + ZPE + H_tr,rot,vib − TS_tr,rot,vib

ΔG = G_red. − G_ox. = −FE⁰

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,⁶⁰ although this quantity is difficult to measure and values in the range of 4.2–4.7 V have been reported in the literature.⁶¹

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,⁴⁷ and the PBE0 hybrid functional.⁶² 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,⁷⁰ 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.⁷¹ 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₂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.⁶³ 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⁻¹ (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.

Table 1 Functional dependence of the redox potential calculation and comparison with the experimentally determined redox potentials. The BP86 results are shown in bold as they are judged to be the most reliable from the present data and literature precedent

	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

---

Table 2 Addition of diffuse functions and solvent dependence of the redox potential calculation. The line shown in bold is identical to the bold line in Table 1

	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

---

2. Computational discovery of dissociation pathways by ab initio molecular dynamics

We used ab initio molecular dynamics (AIMD) to explore the chemical and structural rearrangements of Fe₄N(CO)₁₂ in its different electronic states. In the Born–Oppenheimer MD (BOMD) framework, the motion of atoms is simulated by applying the nuclear gradients of the energy as classical forces to the atoms, and then accelerating the atoms along the force vectors using Newton’s second law. In order for the simulations to broadly sample the chemical space and discover many reaction pathways while keeping the computational cost affordable, accelerated sampling techniques are needed to cross over potential barriers more rapidly.^72–77 In this study, we simply ran unbiased AIMD at elevated temperatures to accelerate the sampling. A velocity Verlet integrator was used with a time step of 1.0 fs. A Langevin thermostat was used with the equilibrium temperature set to 1000 K and a collision frequency of 1.0 ps⁻¹. Several simulations were started from the energy minimized structures of 1²⁻ and 1³⁻, as well as from the protonated isoelectronic species, i.e. H-1¹⁻ and H-1²⁻. These simulations used the B3LYP functional and a hybrid basis set combining 6-31G*⁷⁸ for light elements and the LANL2DZ basis set/ECP for Fe,⁷⁹ abbreviated as 6-31G*-LDZ.

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¹⁻. 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 22 500, 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 37 000.

Fig. 1 RMSD time series to the initial optimized structure for a high temperature AIMD simulation of H-1¹⁻. Several trajectory snapshots are shown, along with blue arrows indicating their corresponding simulation time. Fe, orange; C, grey; N, blue; O, red; H, white.

3. Characterization of optimized structures

The AIMD simulation explores the potential landscape very broadly, but a closer examination of the optimized structures and barriers is needed to assess the feasibility of the discovered pathways under experimental conditions. We focused on trajectory frames numbered 22 000–22 390, where the CO is observed to dissociate from the complex, and optimized a total of 40 trajectory frames evenly spaced by 10 frames (i.e. spanning 400 simulation time steps). The proton was deleted from the trajectory frames prior to optimization. The charge and spin multiplicity of the twice-reduced state were set to −3 and 1, respectively, prompted by CASSCF (8,8) calculations that indicated that the lowest-energy state is a closed shell singlet.

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⁸⁰ range from 1.05 to 1.25, indicating single bond order.

Fig. 2 Optimized structures of 1³⁻ (left) and 2³⁻ + CO (right) at the B3LYP/6-31G*-LANL2DZ level of theory, before and after CO dissociation. The structures are characterized by a nearly planar isosceles triangular face (left) and a rectangular face (right) that contain the nitrogen atom. Fe, orange; C, grey; N, blue; O, red; H, white.

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¹⁷ 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³⁻ structure; these were not included in the present calculations.

Fig. 3 Optimized, doubly-occupied molecular orbital of 2³⁻ at the CASSCF(8,8)/6-31G*-LANL2DZ level of theory, indicating a delocalized bond that connects the three Fe centers and two bridging C atoms in the foreground. The orbital is plotted with an isosurface value of 0.07.

To assess the possibility that 2 may be a catalyst for CO₂ reduction, we computed redox potentials of the 2⁰/2⁻, 2⁻/2²⁻, and 2²⁻/2³⁻ 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₂ reduction reaction.

4. Calculation of barrier heights of CO dissociation

The AIMD simulation that discovered the dissociation pathway is a good starting point for estimating the activation barrier separating the initial and final states. An initial reaction pathway is obtained by concatenating the MD trajectory frames with the output frames from the geometry optimization. From these structures, an “initial chain” of 21 equally spaced frames is selected. Because the initial chain may contain kinks that interfere with the convergence of reaction path optimization methods, we performed an initial smoothing by minimizing an elastic band energy function that depends solely on internal coordinate displacements. The resulting “smoothed chain” is free of kinks and has a shorter arc length than the initial chain, and is input into a nudged elastic band (NEB) calculation. The NEB uses a climbing image approach to ensure the highest-energy structure is as close as possible to the true transition state.

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⁸¹, with input parameters that reproduce the TeraChem total energies to within 0.0005 a.u. (<0.5 kcal mol⁻¹). 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.

Fig. 4 Comparison of relative energies along the CO dissociation coordinate from 1³⁻ and 1²⁻, calculated using BP86/6-31G*-LDZ. The Fe–C distance for the dissociating CO ligand is highlighted. Fe, orange; C, grey; N, blue; O, red; H, white.

Table 3 Relative energies and free energies (via harmonic approximation) for CO dissociation from 1³⁻ and 1²⁻. Each group of four rows refers to calculations performed using a different DFT approximation. In the fourth row of each group, energies and free energies are calculated as a sum of the separated species in the product. All energies are reported in kcal mol⁻¹

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³⁻ 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₁ = +14.5 kcal mol⁻¹ and an activation barrier of E_a1 = +15.5 kcal mol⁻¹; 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⁻¹ (ΔG^‡ = 17.4 kcal mol⁻¹) 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⁻¹ (ΔG = +7.2 kcal mol⁻¹). We also computed the reaction energy by treating the products as completely separate species and obtained ΔE_sep = 13.3 kcal mol⁻¹ (ΔG_sep = +1.4 kcal mol⁻¹). 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³⁻ species.

We also investigated CO dissociation from the 1²⁻ electronic state; because dissociation is not observed from 1²⁻ in the experimental studies, we presume that the calculated thermodynamic and/or kinetic parameters should be less favourable compared to 1³⁻. In searching for the reaction energies and activation barriers for the 1²⁻ 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⁻¹ and E_a = 22.5 kcal mol⁻¹ (ΔG = +16.6 kcal mol⁻¹; ΔG^‡ = 18.1 kcal mol⁻¹). The reaction energy calculated using separated species as the products is ΔE = 25.2 kcal mol⁻¹ (ΔG = 12.8 kcal mol⁻¹).

Comparison of the dissociation pathways from 1³⁻vs.1²⁻ gives reaction free energies of ΔG = +7.2 vs. +16.3 kcal mol⁻¹; with separated product species, ΔG_sep = +1.4 vs. +12.8 kcal mol⁻¹. These values indicate that CO dissociation from 1²⁻ is thermodynamically less favourable than from 1³⁻, consistent with the experimental findings. On the other hand, although the energy barrier for 1³⁻ is lower than for 1²⁻ (ΔE = 18.8 vs. 22.4 kcal mol⁻¹), the calculated activation free energies are nearly the same (ΔG^‡ = 17.4 vs. 18.1 kcal mol⁻¹). Comparison of the overall shape of the dissociation curve shows some other important differences: whereas the 1³⁻ pathway has two clearly defined barriers and an intermediate, the 1²⁻ pathway is nearly barrierless, which indicates that tunnelling effects may play a significant role in determining the reaction rate.⁸² In summary, CO dissociation from 1²⁻ 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.

5. Validation of electronic structure method

The veracity of our predictions regarding CO dissociation rests upon the choice of method. In this section we provide some justifications for our use of DFT in general, and the BP86/6-31G*-LDZ level of theory in particular. Our comparison tests include four DFT approximations (BP86, B3LYP, the meta-GGA functional M06-L⁸³, and the hybrid meta-GGA M06 ⁸⁴). Whereas the first two functionals contain minimal empiricism, the latter two functionals contain 30+ parameters fitted to databases of diverse molecular properties. Optimized IRCs from 1³⁻ and 1²⁻ were computed using all four functionals in the 6-31G*-LDZ basis (Table 3). We also tested for basis set effects in the BP86 and B3LYP calculations by comparing a smaller double-zeta basis 6-31G*-LDZ (6-31G* for main group, LANL2DZ for Fe) and a larger triple-zeta basis TZVP-LTZ (TZVP^85,86 for main group, LANL2TZ for Fe). ΔE and E_a in the large basis set were estimated by taking differences in single-point energies along the small basis optimized pathway following the IRCMax approach.⁸⁷ Our results for comparing BP86 vs. B3LYP and the basis set effects in the 1³⁻ dissociation pathway are shown in ESI Table S2 and ESI Fig. S4.†⁸⁸

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³⁻, ΔE is essentially unchanged from the 6-31G*-LDZ result (10.6 kcal mol⁻¹); E_a is slightly lower at 18.3 kcal mol⁻¹. For the 1²⁻ pathway, BP86/TZVP-LTZ predicts a slightly higher value of ΔE = 22.7 kcal mol⁻¹, 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⁻¹ 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⁻¹ and E_a = 25.3 kcal mol⁻¹ for CO dissociation from 1³⁻; notably, E_a is 6 kcal mol⁻¹ higher than in BP86. Despite differences in the barrier height, the structures along the 1³⁻ 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²⁻ dissociation pathway. For the reactant (1²⁻) structure, B3LYP predicts a pyramidal structure with an isosceles triangular base, almost identical to the structure of 1³⁻ 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.

Fig. 5 Optimized structure of 1²⁻ at the BP86/6-31G*-LANL2DZ level of theory, characterized by a “crooked butterfly” structure with a single elongated Fe–Fe distance of 3.0 Å. Fe, orange; C, grey; N, blue; O, red; H, white.

The differences in BP86 vs. B3LYP in the 1²⁻ 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²〉 to measure any deviations of the KS wavefunction from a pure doublet (ESI Fig. S5†). Along the BP86 pathway, the 〈S²〉 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²〉 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²〉 value close to 2.0 at the dissociated state. BP86 also predicts 〈S²〉 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²⁻ 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³⁻ and 1²⁻ 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³⁻ and 1²⁻ 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⁻⁴ of the diagonal elements, and off-diagonal elements are all <10⁻⁴. 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³⁻ 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₂.⁴⁹

Fig. 6 Natural orbital occupation numbers calculated from CASSCF. The input geometries are from the IRCs of CO dissociation from 1³⁻ (left) and 1²⁻ (right) optimized using BP86 (top) and B3LYP (bottom). Active spaces of (4,6) and (3,6) were used for all structures from 1³⁻ and 1²⁻, respectively.

For CO dissociation from 1²⁻, the CASSCF calculations using BP86-optimized structures show a similar pattern to 1³⁻, 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²⁻ 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²⁻ dissociation pathway is low, with 〈S²〉 = 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²⁻ dissociation pathways.

6. Calculated vibrational analyses

Infrared (IR) absorption spectra provide a meaningful connection between theory and experiment; a harmonic vibrational analysis calculation provides a series of frequencies and intensities that may be converted to a simulated spectrum by applying artificial broadening to each absorption peak. The results of two frequency calculations are shown in Fig. 7, where we compare the IR absorption peaks of 1²⁻ and 2³⁻, the presumed initial and final states of the CO dissociation. The approximate spectra cannot accurately reproduce the widths of the experimental peaks, and only the shift in the peak positions, or the appearance of new peaks, is meaningful. The most notable feature in the spectrum of 2³⁻ is a new peak that appears in a region red-shifted from the main CO-stretching band by about 150 cm⁻¹. The vibrational mode of this peak corresponds to CO-stretching of the bridging CO ligands. The reduced frequency indicates a slightly lower force constant in the CO bond of the bridging ligands that donate more electron density to the Fe centers. This red-shifted stretching peak may be used as a vibrational signature of CO dissociation.

Fig. 7 Comparison of vibrational spectra calculated at the BP86/6-31G*-LDZ level for 1²⁻ (blue) and 2³⁻ (orange). The IR spectrum of 2³⁻ has a distinct peak red-shifted from the main band of CO-stretches by about 150 cm⁻¹ (green arrow), corresponding to a symmetric and antisymmetric stretch of the bridging CO ligands (bottom). An artificial Lorentzian broadening of 10 cm⁻¹ is used.

Conclusions

In this study, we calculated the redox properties of the CO₂ reduction electrocatalyst [Fe₄N(CO)₁₂]³⁻ (1⁻) and investigated the possibility of CO dissociation from the twice-reduced state, 1³⁻. Our calculated redox potentials show close agreement with experimentally measured values. The structure of the product of CO dissociation (2³⁻) was predicted and found to be in close agreement with the experimental X-ray crystal structure. The CO dissociation pathway from 1³⁻ is energetically accessible under ambient conditions (in kcal mol⁻¹: ΔE = +10.6, E_a = 18.8; ΔG_sep = +1.4, ΔG^‡ = +17.4). The analogous CO dissociation from 1²⁻ has a higher reaction energy and similar activation free energy (in kcal mol⁻¹: ΔE = 22.1, E_a = 22.4; ΔG_sep = +12.8, ΔG^‡ = +18.1) with a nearly barrierless dissociation curve. Vibrational analysis of 2³⁻ shows a distinct CO-stretching peak red-shifted from the main CO-stretching band, indicating a possible vibrational signature of CO dissociation. Our calculations indicate that the BP86 semi-local functional gives more reliable results than the B3LYP hybrid functional in the study of this system. Future studies will focus on the potentially important role of counterions in stabilizing redox intermediates, as well as the strong solvent dependence in the selectivity of this catalyst for H₂ evolution vs. CO₂ reduction.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

L. A. Berben is grateful to the Department of Energy Office of Science for support from award number DE-SC0016395; L.-P. Wang and Y. Qiu acknowledge a gift from Walt Disney Imagineering. H. Jang acknowledges support from the ACS Petroleum Research Fund, award number 58158-DNI6. We would like to thank Atefeh Taheri, Tobias Sherbow and Natalia Loewen for helpful discussions.

References

1. T. P. Senftle and E. A. Carter, Acc. Chem. Res., 2017, 50, 472–475.
2. Y. Hori, A. Murata and R. Takahashi, J. Chem. Soc., Faraday Trans. 1, 1989, 85, 2309–2326.
3. X. Feng, K. Jiang, S. Fan and M. W. Kanan, J. Am. Chem. Soc., 2015, 137, 4606–4609.
4. X. Min and M. W. Kanan, J. Am. Chem. Soc., 2015, 137, 4701–4708.
5. M. Asadi, K. Kim, C. Liu, A. V. Addepalli, P. Abbasi, P. Yasaei, P. Phillips, A. Behranginia, J. M. Cerrato, R. Haasch, P. Zapol, B. Kumar, R. F. Klie, J. Abiade, L. A. Curtiss and A. Salehi-Khojin, Science, 2016, 353, 467.
6. D. R. Kauffman, D. Alfonso, C. Matranga, H. Qian and R. Jin, J. Am. Chem. Soc., 2012, 134, 10237–10243.
7. S. Zhang, P. Kang and T. J. Meyer, J. Am. Chem. Soc., 2014, 136, 1734–1737.
8. J. Rosen, G. S. Hutchings, Q. Lu, S. Rivera, Y. Zhou, D. G. Vlachos and F. Jiao, ACS Catal., 2015, 5, 4293–4299.
9. J.-M. Savéant, Chem. Rev., 2008, 108, 2348–2378.
10. E. E. Benson, C. P. Kubiak, A. J. Sathrum and J. M. Smieja, Chem. Soc. Rev., 2009, 38, 89–99.
11. C. Costentin, M. Robert and J.-M. Saveant, Chem. Soc. Rev., 2013, 42, 2423–2436.
12. M. Tachikawa and E. L. Muetterties, J. Am. Chem. Soc., 1980, 102, 4541–4542.
13. M. Tachikawa, J. Stein, E. L. Muetterties, R. G. Teller, M. A. Beno, E. Gebert and J. M. Williams, J. Am. Chem. Soc., 1980, 102, 6648–6649.
14. M. D. Rail and L. A. Berben, J. Am. Chem. Soc., 2011, 133, 18577–18579.
15. A. D. Nguyen, M. D. Rail, M. Shanmugam, J. C. Fettinger and L. A. Berben, Inorg. Chem., 2013, 52, 12847–12854.
16. A. Taheri, E. J. Thompson, J. C. Fettinger and L. A. Berben, ACS Catal., 2015, 5, 7140–7151.
17. A. Taheri, D. B. Cluff and L. A. Berben, Organometallics, 2017 DOI:10.1021/acs.organomet.7b00824.
18. H. Wang, Y. Xie, R. B. King and H. F. Schaefer, Inorg. Chem., 2006, 45, 3384–3392.
19. H. Wang, Y. Xie, R. B. King and H. F. Schaefer, Organometallics, 2008, 27, 3113–3123.
20. H. Wang, Y. Xie, R. B. King and H. F. Schaefer, Inorg. Chem., 2008, 47, 3045–3055.
21. Z. Zhang, Q.-s. Li, Y. Xie, R. B. King and H. F. Schaefer, Inorg. Chem., 2009, 48, 6167–6177.
22. Z. Zhang, Q.-s. Li, Y. Xie, R. B. King and H. F. Schaefer, Inorg. Chem., 2009, 48, 1974–1988.
23. H. Wang, Z. Sun, Y. Xie, R. B. King and H. F. Schaefer, Organometallics, 2010, 29, 630–641.
24. L. Xu, Q.-s. Li, Y. Xie, R. B. King and H. F. Schaefer, Inorg. Chem., 2010, 49, 1046–1055.
25. L. Xu, Q.-s. Li, R. B. King and H. F. Schaefer, Organometallics, 2011, 30, 5084–5087.
26. Y. Zeng, S. Wang, H. Feng, Y. Xie, R. B. King and H. F. Schaefer III, New J. Chem., 2011, 35, 920–929.
27. J. Chen, S. Chen, L. Zhong, F. Hao, Y. Xie, R. Bruce King and H. F. Schaefer III, Theor. Chem. Acc., 2012, 131, 1090.
28. Z. Zhang, Q.-s. Li, R. B. King and H. F. Schaefer, Eur. J. Inorg. Chem., 2012, 2012, 1104–1113.
29. Q. Fan, H. Feng, W. Sun, Y. Xie, R. B. King and H. F. Schaefer, Organometallics, 2012, 31, 3610–3619.
30. R. Jin, X. Chen, Q. Du, H. Feng, Y. Xie, R. B. King and H. F. Schaefer, Organometallics, 2012, 31, 5005–5017.
31. R. Jin, X. Chen, Q. Du, H. Feng, Y. Xie, R. B. King and H. F. Schaefer, RSC Adv., 2016, 6, 82661–82668.
32. X. Gong, L. Zhu, J. Yang, X. Gao, Q.-s. Li, Y. Xie, R. B. King and H. F. Schaefer III, New J. Chem., 2014, 38, 3762–3769.
33. H. Li, H. Feng, W. Sun, Y. Xie, R. B. King and H. F. Schaefer, Organometallics, 2013, 32, 88–94.
34. G. Li, L. Zhou, X. Zhai, Q.-s. Li, Y. Xie, R. B. King and H. F. Schaefer III, New J. Chem., 2013, 37, 3294–3302.
35. J. Deng, Q.-s. Li, Y. Xie, R. B. King and I. H. F. Schaefer, New J. Chem., 2013, 37, 2902–2910.
36. T. T. Shi, Q.-S. Li, Y. Xie, R. B. King and H. F. Schaefer III, New J. Chem., 2010, 34, 208–214.
37. H. Wang, Y. Xie, R. B. King and H. F. Schaefer, J. Am. Chem. Soc., 2006, 128, 11376–11384.
38. G. Wang, J. Cui, C. Chi, X. Zhou, Z. H. Li, X. Xing and M. Zhou, Chem. Sci., 2012, 3, 3272–3279.
39. C. Chi, J. Cui, Z. H. Li, X. Xing, G. Wang and M. Zhou, Chem. Sci., 2012, 3, 1698–1706.
40. M. Mirmohades, S. Pullen, M. Stein, S. Maji, S. Ott, L. Hammarström and R. Lomoth, J. Am. Chem. Soc., 2014, 136, 17366–17369.
41. S. Kuppuswamy, J. D. Wofford, C. Joseph, Z.-L. Xie, A. K. Ali, V. M. Lynch, P. A. Lindahl and M. J. Rose, Inorg. Chem., 2017, 56, 5998–6012.
42. M. H. Baik and R. A. Friesner, J. Phys. Chem. A, 2002, 106, 7407–7412.
43. A. V. Marenich, J. M. Ho, M. L. Coote, C. J. Cramer and D. G. Truhlar, Phys. Chem. Chem. Phys., 2014, 16, 15068–15106.
44. L. P. Wang and T. Van Voorhis, J. Chem. Theory Comput., 2012, 8, 610–617.
45. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100.
46. J. P. Perdew, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 33, 8822–8824.
47. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
48. L.-P. Wang, R. T. McGibbon, V. S. Pande and T. J. Martinez, J. Chem. Theory Comput., 2016, 12, 638–649.
49. M. S. Gordon, M. W. Schmidt, G. M. Chaban, K. R. Glaesemann, W. J. Stevens and C. Gonzalez, J. Chem. Phys., 1999, 110, 4199–4207.
50. P. Pulay and T. P. Hamilton, J. Chem. Phys., 1988, 88, 4926–4933.
51. R. G. Parr and W. Yang, Annu. Rev. Phys. Chem., 1995, 46, 701–728.
52. A. Klamt, J. Phys. Chem., 1995, 99, 2224–2235.
53. I. S. Ufimtsev and T. J. Martinez, J. Chem. Theory Comput., 2008, 4, 222–231.
54. I. S. Ufimtsev and T. J. Martinez, J. Chem. Theory Comput., 2009, 5, 2619–2628.
55. I. S. Ufimtsev and T. J. Martinez, J. Chem. Theory Comput., 2009, 5, 1004–1015.
56. C. C. Song, L. P. Wang, T. Sachse, J. Preiss, M. Presselt and T. J. Martinez, J. Chem. Phys., 2015, 143, 014114.
57. C. C. Song, L. P. Wang and T. J. Martinez, J. Chem. Theory Comput., 2016, 12, 92–106.
58. F. Liu, N. Luehr, H. J. Kulik and T. J. Martinez, J. Chem. Theory Comput., 2015, 11, 3131–3144.
59. L. P. Wang and C. C. Song, J. Chem. Phys., 2016, 144, 214108.
60. H. Reiss and A. Heller, J. Phys. Chem., 1985, 89, 4207–4213.
61. W. A. Donald, R. D. Leib, J. T. O’Brien, M. F. Bush and E. R. Williams, J. Am. Chem. Soc., 2008, 130, 3371–3381.
62. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
63. F. Furche and J. P. Perdew, J. Chem. Phys., 2006, 124, 044103.
64. C. J. Cramer and D. G. Truhlar, Phys. Chem. Chem. Phys., 2009, 11, 10757–10816.
65. J. M. Galbraith and H. F. Schaefer, J. Chem. Phys., 1996, 105, 862–864.
66. A. A. Jarcki and E. R. Davidson, Chem. Phys. Lett., 1999, 300, 44–52.
67. F. Jensen, J. Chem. Theory Comput., 2010, 6, 2726–2735.
68. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
69. F. Weigend, Phys. Chem. Chem. Phys., 2006, 8, 1057–1065.
70. L. E. Roy, P. J. Hay and R. L. Martin, J. Chem. Theory Comput., 2008, 4, 1029–1031.
71. J. Zheng, X. Xu and D. G. Truhlar, Theor. Chem. Acc., 2011, 128, 295–305.
72. A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 12562–12566.
73. J. Pfaendtner and M. Bonomi, J. Chem. Theory Comput., 2015, 11, 5062–5067.
74. F. Pietrucci and W. Andreoni, Phys. Rev. Lett., 2011, 107, 085504.
75. L.-P. Wang, A. Titov, R. McGibbon, F. Liu, V. S. Pande and T. J. Martínez, Nat. Chem., 2014, 6, 1044–1048.
76. L. Xie, Q. Zhao, K. F. Jensen and H. J. Kulik, J. Phys. Chem. C, 2016, 120, 2472–2483.
77. N. Goldman, E. J. Reed, L. E. Fried, I. F. William Kuo and A. Maiti, Nat. Chem., 2010, 2, 949–954.
78. W. J. Hehre, R. Ditchfield and J. A. Pople, J. Chem. Phys., 1972, 56, 2257–2261.
79. P. J. Hay and W. R. Wadt, J. Chem. Phys., 1985, 82, 299–310.
80. I. Mayer, J. Comput. Chem., 2007, 28, 204–221.
81. Y. Shao, Z. Gan, E. Epifanovsky, A. T. B. Gilbert, M. Wormit, J. Kussmann, A. W. Lange, A. Behn, J. Deng, X. Feng, D. Ghosh, M. Goldey, P. R. Horn, L. D. Jacobson, I. Kaliman, R. Z. Khaliullin, T. Kuś, A. Landau, J. Liu, E. I. Proynov, Y. M. Rhee, R. M. Richard, M. A. Rohrdanz, R. P. Steele, E. J. Sundstrom, H. L. Woodcock, P. M. Zimmerman, D. Zuev, B. Albrecht, E. Alguire, B. Austin, G. J. O. Beran, Y. A. Bernard, E. Berquist, K. Brandhorst, K. B. Bravaya, S. T. Brown, D. Casanova, C.-M. Chang, Y. Chen, S. H. Chien, K. D. Closser, D. L. Crittenden, M. Diedenhofen, R. A. DiStasio, H. Do, A. D. Dutoi, R. G. Edgar, S. Fatehi, L. Fusti-Molnar, A. Ghysels, A. Golubeva-Zadorozhnaya, J. Gomes, M. W. D. Hanson-Heine, P. H. P. Harbach, A. W. Hauser, E. G. Hohenstein, Z. C. Holden, T.-C. Jagau, H. Ji, B. Kaduk, K. Khistyaev, J. Kim, J. Kim, R. A. King, P. Klunzinger, D. Kosenkov, T. Kowalczyk, C. M. Krauter, K. U. Lao, A. D. Laurent, K. V. Lawler, S. V. Levchenko, C. Y. Lin, F. Liu, E. Livshits, R. C. Lochan, A. Luenser, P. Manohar, S. F. Manzer, S.-P. Mao, N. Mardirossian, A. V. Marenich, S. A. Maurer, N. J. Mayhall, E. Neuscamman, C. M. Oana, R. Olivares-Amaya, D. P. O’Neill, J. A. Parkhill, T. M. Perrine, R. Peverati, A. Prociuk, D. R. Rehn, E. Rosta, N. J. Russ, S. M. Sharada, S. Sharma, D. W. Small, A. Sodt, T. Stein, D. Stück, Y.-C. Su, A. J. W. Thom, T. Tsuchimochi, V. Vanovschi, L. Vogt, O. Vydrov, T. Wang, M. A. Watson, J. Wenzel, A. White, C. F. Williams, J. Yang, S. Yeganeh, S. R. Yost, Z.-Q. You, I. Y. Zhang, X. Zhang, Y. Zhao, B. R. Brooks, G. K. L. Chan, D. M. Chipman, C. J. Cramer, W. A. Goddard, M. S. Gordon, W. J. Hehre, A. Klamt, H. F. Schaefer, M. W. Schmidt, C. D. Sherrill, D. G. Truhlar, A. Warshel, X. Xu, A. Aspuru-Guzik, R. Baer, A. T. Bell, N. A. Besley, J.-D. Chai, A. Dreuw, B. D. Dunietz, T. R. Furlani, S. R. Gwaltney, C.-P. Hsu, Y. Jung, J. Kong, D. S. Lambrecht, W. Liang, C. Ochsenfeld, V. A. Rassolov, L. V. Slipchenko, J. E. Subotnik, T. Van Voorhis, J. M. Herbert, A. I. Krylov, P. M. W. Gill and M. Head-Gordon, Mol. Phys., 2015, 113, 184–215.
82. U. Bozkaya, J. M. Turney, Y. Yamaguchi and H. F. Schaefer, J. Chem. Phys., 2012, 136(16), 164303.
83. Y. Zhao and D. G. Truhlar, J. Chem. Phys., 2006, 125, 194101.
84. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
85. A. Schäfer, C. Huber and R. Ahlrichs, J. Chem. Phys., 1994, 100, 5829–5835.
86. A. Schäfer, H. Horn and R. Ahlrichs, J. Chem. Phys., 1992, 97, 2571–2577.
87. D. K. Malick, G. A. Petersson and J. A. Montgomery, J. Chem. Phys., 1998, 108, 5704–5713.
88. D. K. Malick, G. A. Petersson and J. A. Montgomery, J. Chem. Phys., 1998, 108, 5704–5713.
89. D. Schweinfurth, M. G. Sommer, M. Atanasov, S. Demeshko, S. Hohloch, F. Meyer, F. Neese and B. Sarkar, J. Am. Chem. Soc., 2015, 137, 1993–2005.
90. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 73–78.

---

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²⁻. See DOI: 10.1039/c7sc04342b

---