Insights into dioxygen binding on metal centers: an ab initio multireference electronic structure analysis

Abstract

Why does binding of dioxygen (O₂) to metal centers, the initial step of O₂ storage, transportation, and activation, almost inevitably induce metal-to-O₂ single-electron transfer and generate superoxo (O₂⁻˙) species, instead of genuine O⁰₂ adducts? To address this question, this study describes highly correlated wavefunction-based ab initio calculations using CASSCF/NEVPT2 (CASSCF = complete active space self-consistent field, and NEVPT2 = N-electron valence state second-order perturbation theory) approaches to explore the electronic-structure evolution of O₂ association on Fe(II)(BDPP) (H₂BDPP = 2,6-bis((2-(S)-diphenylhydroxylmethyl-1-pyrrolidinyl)methyl)pyridine) and Co(II)(BDPP) to produce S = 3 Fe(III)(BDPP)(O₂⁻˙) (1) and Co(III)(BDPP)(O₂⁻˙) (2). CASSCF/NEVPT2 calculations suggest that the processes furnishing 1 and 2 feature an avoided crossing resulting from interactions of two diabatic curves, of which one is characterized as Co(II) and Fe(II) centers interacting with a triplet O₂ ligand and the other as Co(III) and Fe(III) centers bound to a superoxo ligand. In both cases, the avoided crossing induces a one-electron transfer from the divalent metal center to the incoming O₂ ligand and leads to formation of trivalent metal–O₂⁻˙ complexes. To facilitate the interpretation of complicated multireference wavefunctions, we formulated two-fragment spin eigenfunctions utilizing Clebsch–Gordan coefficients (CGCs) to rationalize computed spin populations on the metal centers and the O₂ ligand and compared these results with usual valence bonding (VB) analyses. It turns out that both methods give the same results and are complementary to each other. Finally, the limitation of DFT approaches in describing complex electronic structures involving metal–ligand magnetic couplings is delineated.

---

Introduction

Binding of dioxygen (O₂) on metal cofactors of metalloproteins initiates a diverse array of physiologically pivotal processes, including O₂ storage, transportation, and activation for substrate functionalization.^1–6 Myoglobin, a heme enzyme containing a high spin ferrous center, reacts with O₂ yielding a singlet O₂ complex,^7–9 while for nonheme enzymes isopenicillin N synthase (IPNS)^10–12 and a variant of homoprotocatechuate 2,3-dioxygenase (Fe-HPCD) with the active site His200 mutated to Asn, labelled as 4-NC,¹³ and utilizing an alternative substrate 4-nitrocatechol, O₂ addition results in a quintet O₂ adduct that is best described as a high spin ferric center antiferromagnetically coupled to a superoxo ligand (O₂⁻˙) yielding an overall S_t = 2 ground state. To mimic enzymatic O₂ binding, numerous model complexes have been synthesized and characterized.^14,15 The reactions of O₂ with a range of ferrous complexes, such as Fe(II)(BNPA^Me₂S)(Br) (BNPA^Me₂SH = (bis((6-(neopentylamino)pyridinyl)methyl)amino)-2-methylpropane-2-thiol,¹⁶ (L^AdH)Fe (L^AdH = tris(phosphinimine) ligand),¹⁷ Fe(S₂^Me₂N₃(Pr,Pr))¹⁸ and Fe(ttppc) (ttppc³⁻ = 5,10,15-tris((2,4,6-triphenyl)phenyl)corrolate),¹⁹ all yield similar quintet Fe–O₂ adducts. In contrast to the examples discussed above, treatment of O₂ with a high spin Fe(II) precursor (S_Fe = 2), Fe(II)(BDPP) (H₂BDPP = 2,6-bis((2-(S)-diphenylhydroxylmethyl-1-pyrrolidinyl)methyl)pyridine), at 193 K generates S_t = 3, instead of S_t = 2, Fe(III)(BDPP)(O₂⁻˙) (1) attained by a ferromagnetic coupling of Fe(III) center and a O₂⁻˙ ligand .²⁰ While the corresponding reaction with high spin Co(II)(BDPP)²¹ at 193 K furnishes a doublet Co(III)–O₂⁻˙ species, Co(III)(BDPP)(O₂⁻˙) (2), which has been characterized as consisting of a low spin Co(III) center (S_Co = 0) interacting with a superoxo ligand .²² Furthermore, the same product was also found for O₂ addition to Co-HPCD,²³ the Co derivative of Fe-HPCD, and a plethora of high spin or low spin Co(II) complexes, including Co-TPP (TPP²⁻ = tetraphenylporphyrinate),²⁴ Co(Me₃TACN)(S₂SiMe₂) (Me₃TACN = 1,4,7-trimethyl-1,4,7-triazacyclononane),²⁵ Co(Tp^Me₂)(^tBu2APH) (Tp^Me₂ = hydrotris(3,5-dimethylpyrazolyl)borate,^{26 tBu2}APH = 2-amino-4,6-di-tert-butylphenol), Co(salen)²⁷ and Co(Por) (Por = porphyrin).²⁸ Over-decade intensive investigations revealed that O₂ association on divalent nonheme metal centers almost inevitably entails a metal-to-O₂ electron transfer and thus produces a trivalent metal–O₂⁻˙ complex.^29,30 To the best of our knowledge, there is only one exception; O₂ binding on a variant of Fe-HPCD labelled as Y257F-HPCA was proposed to form an Fe(II)–O₂ intermediate, instead of Fe(III)–O₂⁻˙.^31–35 An intriguing question thus arises about whether it is viable to generate such an authenticated O₂ adduct from O₂ association on metal centers.

In parallel to intensive experimental research, a great deal of computational studies have been devoted to probing electronic structures of O₂ adducts and evaluating electronic-structure evolution over the course of O₂ binding processes.^36–44 In this regard, density functional theory (DFT) calculations usually cannot yield satisfactory results, especially for those systems featuring antiferromagnetic coupling of two fragments. In principle, such a complicated electronic structure cannot be appropriately described by a single Slater determinant. Consequently, broken symmetry formalism has to be invoked, because of the inherent single determinant restriction of DFT approaches.⁴⁵ The Kohn–Sham solution thus obtained is no longer the eigenfunction of Ŝ²; consequently, the resulting spin density is not correct, and the subsequent spin population analysis is physically unjustified. To obtain more accurate electronic structures, one has to resort to highly correlated wavefunction-based multireference approaches.^{36–38,42,44} To date, these methods have been successfully applied to O₂ adducts of heme systems, but hardly to nonheme ones.³⁰ To the best of our knowledge, only one ab initio study has been reported thus far where potential energy curves of varying spin states were computed for O₂ association on a heme center.⁴⁶ This work does not focus on the elucidation of electronic-structure changes occurring during this process but on the quintet-to-singlet spin crossover; hence, it does not provide sufficient information to address the aforementioned question.

In the present work, we describe a detailed multireference analysis of the electronic-structure evolution as O₂ steadily approaches the Fe(II) and Co(II) centers in Fe(II)(BDPP) and Co(II)(BDPP) affording 1 and 2, respectively. Specifically, highly correlated wavefunction-based complete active space self-consistent field (CASSCF)⁴⁷/N-electron valence state second-order perturbation theory (NEVPT2)^48,49 computations were carried out on a series of geometric structures on the reaction trajectory generated by DFT computations. To interpret complicated multireference wavefunctions of O₂ adducts, we constructed two-fragment spin eigenfunctions using Clebsch–Gordan coefficients (CGCs)⁵⁰ to rationalize computed spin populations. CASSCF/NEVPT2 calculations suggest that the O₂ binding processes furnishing 1 and 2 feature an avoided crossing of two diabatic curves that represent the limiting valence structures of Fe(II)/Co(II)–O⁰₂ and Fe(III)/Co(III)–O₂⁻˙. Their interaction induces a spontaneous single-electron transfer from the divalent metal center to the incoming O₂ ligand, thereby leading to the formation of trivalent metal–O₂⁻˙ complexes as the final product for the O₂ binding processes.

Computation setup

To obtain the reaction trajectory for the formation of complexes 1 and 2, the M–O₂ distance steadily varied from 1.75 to 4.50 Å in an S = 3 and an state, respectively, in relaxed surface scans conducted by using the B3LYP functional^51,52 in conjunction with the def2-TZVP basis set^53–55 for N, O, Fe, and Co, and the def2-SVP basis set^56,57 for the other atoms (C and H atoms). Moreover, RIJCOSX approximations⁵⁸ were employed in combination with the def2/J auxiliary basis set^59,60 to speed up DFT computations. To take solvation effects into consideration, the CPCM model⁶¹ with THF as the solvent was applied, and D3BJ corrections^62,63 were used to account for dispersion corrections. Frequency analysis on a given optimized geometry structure was performed to compute its zero point energy (ZPE) and thermal corrections. To compute reliable electronic structures, CASSCF calculations using the def2-TZVP basis set for all atoms were carried out at each geometric structure of the reaction trajectory. As verified below, the appropriate active spaces of complexes 1 and 2 were chosen to distribute 12 and 13 electrons, respectively, into 9 orbitals including five metal 3d-based orbitals (five spherical 3d orbitals are used here), two O₂ π bonding (π_ip and π_op) and the corresponding π* antibonding orbitals. For complex 1, CASSCF calculations averaged a triplet, a quintet, and a septet state, and for complex 2, a doublet, a quartet, and a sextet state were considered. Based on the converged state-average CASSCF wavefunction, complete active space configuration interaction (CASCI) calculations were employed to calculate spin density and spin populations of a specific spin state. On top of CASSCF wavefunctions, NEVPT2 computations were undertaken to obtain more precise electronic energies. To accelerate ab initio computations, RIJK approximations⁶⁴ in conjunction with the def2/JK auxiliary basis set⁶⁵ were used. Thus, the final potential energy surface was constructed by using the electronic energy derived from CASSCF/NEVPT2 calculations combined with the ZPE contribution and thermal corrections. All computations were performed using the ORCA-5.0.3 program package.⁶⁶

Results and discussion

Choice of active spaces

Unlike DFT, CASSCF is not a black box method. As the first step to initiate CASSCF calculations, one should prudently select a proper active space,⁶⁷ for which some preliminary understanding is required of the electronic structure of the system under investigation. To ensure that the chosen active space is sufficient to provide reliable electronic structures, further calculations are often necessary with enlarged and/or truncated active spaces. Critically, CASSCF computations with suitable active spaces not only yield precise electronic structures but also result in fast convergence. Usually, because d and f orbitals of transition metal compounds are nearly degenerate which leads to strong multiconfiguration character, these orbitals should be included in the active space. Equally important is that to properly describe a metal–ligand covalent bond, the active space ought to contain its bonding and antibonding molecular orbitals (MOs). Otherwise, the CASSCF computations with an unbalanced active space would result in an incorrect electronic structure. Both criteria were used to choose a suitable active space for our CASSCF calculations.

In the present case, in addition to 5 metal 3d-based orbitals, two π bonding (π_ip and π_op) and two π* antibonding orbitals of the O₂ ligand should be included in the active space. As elaborated below, this active space can properly describe the metal–O₂ covalent bonding built from the metal 3d and O₂ π* fragment orbitals of appropriate symmetry. Therefore, the active spaces of complexes 1 and 2 were chosen to distribute 12 and 13 electrons, respectively, in 9 orbitals. Depicted in Fig. 1 is the active space of complex 2 obtained by the CASSCF(13,9) computation at the Co–O₂ bond length of 2.50 Å. An analogous active space was employed for complex 1.

Fig. 1 Different choices of active spaces (m,n) (m electrons in n orbitals) for complex 2 at the Co–O₂ bond length of 2.50 Å. For clarity, hydrogen atoms are omitted for clarity.

To test the appropriateness of the selected active space, we took complex 2 at the Co–O₂ bond length of 2.50 Å as an example and performed a series of CASSCF computations with different active spaces. First, we added BDPP²⁻ O and N 2p donor based σ_eq into the active space that is the σ bonding counterpart with respect to d_x²−y². Accordingly, two more electrons ought to be added into the active space, thereby resulting in an active space of CAS(15,10). Second, we constructed an active space of CAS(15,11) that includes O₂ σ and σ* orbitals to take the entire 2p shells of the O₂ moiety into consideration. Third, an even larger active space of CAS(17,12) was built up to encompass all the aforementioned effects. As summarized in Table 1, the occupation numbers of σ_eq are invariably 2 for all key electron configurations of each spin state derived from CASSCF(15,10) calculations. Similarly, the occupation numbers of O₂ σ and σ* remain constants at 2 and 0, respectively, for all important electron configurations predicted by CASSCF(15,11) computations. Both findings hold true for the CASSCF(17,12) results. As such, for a given spin state, all CASSCF computations give qualitatively the same electronic structures as evidenced by almost identical percentage of each dominant electron configuration and hence spin populations on the Co center and the O₂ moiety. Moreover, we have also performed CASSCF calculations with these three distinct active spaces on complex 2 at a longer (4.50 Å) and a shorter (1.91 Å) Co–O₂ bond length, and again the calculations produced analogous results, as shown in Fig. S1 and S2 (ESI†). In addition, we have also carried out CASSCF(13,9) calculations using different basis sets, such as def2-SV, def2-SVP, def2-TZVP and def2-TZVPP, and the same electronic structures were obtained as listed in Tables S1–S3 (ESI†). In summary, the active space of CAS(13,9) combined with the def2-TZVP basis set is sufficient to yield reliable electronic structures and bonding features for complex 2 at varying Co–O₂ bond distances. Likewise, CASSCF(12,9) with the def2-TZVP basis set should produce satisfactory results for the O₂ association process affording complex 1 as well.

Table 1 Spin populations on the Co center and O₂ moiety and electron configurations for each spin state from CASSCF calculations with different choices of active orbitals for complex 2 at the Co–O₂ bond length of 2.50 Å

	Configuration			Spin populations
	Sextet	Quartet	Doublet	Sextet	Quartet	Doublet
CASSCF(13,9)	87% 221112211	42% 221112211	21% 221112211	Co	Co	Co
		14% 221102212	21% 221102212	2.83 (2.88)	2.08 (2.11)	1.56 (1.59)
		10% 221122210	14% 221122210	O2	O2	O2
		9% 221012212	13% 221012212	1.98 (2.00)	0.79 (0.80)	−0.65 (−0.66)
		6% 221212210	9% 221212210
CASSCF(15,10)	87% 2221112211	43% 2221112211	22% 2221112211	Co	Co	Co
		15% 2221102212	22% 2221102212	2.83 (2.88)	2.08 (2.11)	1.56 (1.59)
		10% 2221122210	15% 2221122210	O2	O2	O2
		8% 2221012212	11% 2221012212	1.98 (1.99)	0.79 (0.80)	−0.65 (−0.66)
		5% 2221212210	8% 2221212210
CASSCF(15,11)	87% 22111221120	46% 22111221120	20% 22111221120	Co	Co	Co
		19% 22110221220	30% 22110221220	2.83 (2.88)	2.07 (2.11)	1.56 (1.59)
		13% 22112221020	21% 22112221020	O2	O2	O2
		3% 22101221220	5% 22101221220	1.98 (1.99)	0.80 (0.80)	−0.65 (−0.66)
		2% 22121221020	4% 22121221020
CASSCF(17,11)	87% 222111221120	47% 222111221120	21% 222111221120	Co	Co	Co
		19% 222110221220	31% 222110221220	2.83 (2.88)	2.07 (2.11)	1.56 (1.59)
		13% 222112221020	21% 222112221020	O2	O2	O2
		3% 222101221220	5% 222101221220	1.98 (1.99)	0.80 (0.80)	−0.65 (−0.66)
		2% 222121221020	3% 222121221020

---

O₂ binding on Fe(II)(BDPP) affording complex 1

Depicted in Fig. 2 is the S_t = 3 potential energy surface of O₂ bonding to high spin Fe(II)(BDPP) calculated using CASSCF(12,9)/NEVPT2 computations. As expected, starting point ⁷A with a Fe–O₂ bond length of 4.00 Å is best interpreted as a high spin Fe(II) ion (S_Fe = 2) ferromagnetically coupled with a triplet O₂ ligand (S_O2 = 1) yielding an overall septet state. Consistent with this assignment, ab initio calculations indicate that this state is distinguished by a leading electron configuration of (95%). As shown in Fig. 2, the computed spin density and spin populations (Löwdin analysis, 3.86 for Fe and 2.00 for O₂; Mulliken analysis, 3.91 for Fe and 2.00 for O₂) further corroborate this electronic-structure formulation. Of note, because spin populations estimated by Löwdin and Mulliken analyses are essentially the same, in the following section, we only discuss the former but list the latter in Fig. 2. As the O₂ ligand steadily approaches the Fe center, the energy of this valence structure invariably rises as shown by the corresponding diabatic potential curve (magenta line). During this process, the nature of the Fe–O₂ bonding remains essentially the same as indicated by the unchanged spin populations on the Fe(II) center and the O₂ ligand, for instance, complex ⁷B. However, starting from the Fe–O₂ distance of 2.30 Å, CASSCF(12,9) calculations could converge to another solution as exemplified by complex ⁷1. It features a sole electron configuration of (100%) and is best described as a high spin Fe(III) center ferromagnetically coupled with a superoxo ligand as suggested by the computed spin density and spin populations (4.74 for Fe and 1.03 for O₂). Its energy first descends for the Fe–O₂ distance lowering from 2.30 Å to 2.06 Å and then ascends as the Fe–O₂ distance further decreases; thus, its diabatic potential curve (crimson curve) displays a local minimal (⁷1) at a Fe–O₂ distance of ∼2.00 Å.

Fig. 2 Active spaces of complexes A (a) and 1 (b); septet potential energy surface for the reaction of Fe(II)(BDPP) with O₂ as a function of the Fe–O₂ bond distance obtained from CASSCF(12,9)/NEVPT2 computations; dot lines denote diabatic potential curves (magenta line denotes an S_t = 3 Fe^II–O₂ complex consisting of as an S_Fe = 2 Fe^II center ferromagnetically coupled (FC) with an S_O2 = 1 O₂ ligand; the crimson line denotes an S_t = 3 Fe^III–O₂⁻˙ complex consisting of an Fe^III center ferromagnetically coupled (FC) with an O₂⁻˙ ligand) and the black line denotes an adiabatic curve for the reaction on the S = 3 potential surface; spin density and Löwdin and Mulliken (in parentheses) spin populations on the Fe center and the O₂ unit calculated for complexes ⁷A (Fe(II)(BDPP) + O₂), ⁷B, and ⁷1 on the reaction trajectory are shown at the bottom, yellow and red denote negative and positive spin density, respectively, and local spin populations are also listed in the ESI† (c); the quintet potential energy surface for the reaction of Fe(II)(BDPP) with O₂ as a function of the Fe–O₂ bond distance obtained from CASSCF(12,9)/NEVPT2 computations; the dot lines denote diabatic potential curves (violet line denotes a quintet state Fe^II–O₂ complex, interpreted as an S_Fe = 2 Fe^II compound magnetically coupled with an S_O2 = 1 O₂ ligand; the blue line denotes a quintet state Fe^III–O₂⁻˙ complex, interpreted as an Fe^III center antiferromagnetically coupled (AFC) with an ligand) and the black lines denote adiabatic curves for the reaction on the S = 2 potential surface; spin density and Löwdin and Mulliken (in parentheses) spin populations on the Fe center and O₂ unit calculated for complexes ⁵A (Fe(II)(BDPP) + O₂), ⁵B, and ⁵1 on the reaction trajectory are shown at the bottom, yellow and red denote negative and positive spin density, respectively, and local spin populations are also listed in the ESI† (d); hydrogen atoms are omitted for clarity.

The interaction of the aforementioned two diabatic potential curves results in an avoided crossing near the Fe–O₂ bond length of 2.30 Å and a barrier of 13.6 kcal mol⁻¹. Moreover, the entire process is moderately endothermic with a driving force of 8.0 kcal mol⁻¹. This value should be somewhat overestimated by our computations, because experimentally ⁷1 is stable at low temperatures (<193 K) but releases O₂ affording Fe(II)(BDPP) upon raising the temperature. The overestimation of the binding barrier can be ascribed to the fact that the active space employed does not suffice to capture escalating electron correlations at the transition state where two diabatic surfaces strongly interact. Nevertheless, the entire process can be interpreted as shifting one β electron from the Fe d_xz orbital to the O₂ orbital, thus leading to the formation of ⁷1. As indicated by the septet potential surface, because in the avoided crossing region an energetically favorable electron transfer from the Fe(II) center to the O₂ ligand takes place, which allows the system to ultimately convert into ⁷1, it is unlikely to trap the proposed Fe(II)–O₂ complex during this O₂ binding process. In line with this prediction, UV-vis measurements at 193 K showed that no intervening intermediate could be detected en route to ⁷1.²⁰

Experimentally, the O₂ adducts of IPNS^10–12 and the 4-NC variant of Fe-HPCD¹³ have been spectroscopically identified to be quintet Fe(III)-superoxo species. However, our CASSCF(12,9)/NEVPT2 computations could not allow to determine the exact ground state of complex 1 because the estimated septet–quintet energy gap of 2.4 kcal mol⁻¹ falls within the range of uncertainty for ab initio calculations. Thus, we also investigated the reaction of Fe(II)(BDPP) with O₂ in an S_t = 2 state. As shown in Fig. 2, the quintet potential energy surface computed for the O₂ binding process is essentially parallel to the septet one in the sense that there exists an analogous avoided crossing around the Fe–O₂ distance of 2.30 Å derived from interactions of two diabatic potential curves, yielding a barrier of 11.3 kcal mol⁻¹. Furthermore, the quintet reaction also entails a comparable driving force of 5.6 kcal mol⁻¹. As elaborated below, the diabatic potential curve starting with complex ⁵A describes the energy variation of an S_t = 2 state attained by a magnetic coupling of a high spin Fe(II) center (S_Fe = 2) and a triplet O₂ ligand (S_O2 = 1) as a function of the Fe–O₂ distance; the other one represents a high spin Fe(III) center antiferromagnetically coupled with a doublet superoxo ligand .

Ab initio calculations reveal that complex ⁵A is isoenergetic to ⁷A and, importantly, the leading electron configuration (94%) of the former species is also the same as that of the latter. This can be traced back to the large Fe–O₂ separation that prevents any non-negligible metal–ligand interactions from occurring. Following this line of reasoning, ⁵A should contain a quintet Fe(II) center (S_Fe = 2) that is exchange coupled to a triplet O₂ ligand (S_O2 = 1). However, the spin populations calculated for the Fe center (3.21) and the O₂ ligand (0.67) differ significantly from those expected for non-interacting S_Fe = 2 Fe(II) and S_O2 = 1 O₂ fragments. To solve this problem, the vector coupling model was invoked to construct the spin eigenfunction with a well-defined total spin S of the resulting state derived from the magnetic coupling of a high spin Fe(II) center and a triplet O₂ ligand using CGCs. In the present case, the former entity can be characterized by spin eigenfunction |S₁,M₁〉 = |2,M₁〉 (M₁ = ±2, ±1, 0), and the latter by |S₂,M₂〉 = |1,M₂〉 (M₂ = ±1, 0). The quintet state originating from the interaction of these two fragments should be distinguished by spin eigenfunction |S,M〉 = |2,M〉 (M = ±2, ±1, 0), and its standard component reads , where the functions on the right hand of the equal sign are written as |S₁,M₁|S₂,M₂〉. This spin eigenfunction features Ŝ_z expectation values (the spin projection along the z axis) of and for the Fe(II) center and the O₂ unit, respectively. Indeed, both values are in line with the computed spin populations (3.21 for Fe and 0.67 for O₂), because ideally spin populations are twice the Ŝ_z expectation values, thus affirming the proposed electronic-structure assignment.

Regarding the O₂ binding process furnishing ⁵1, as the Fe–O₂ distance decreases from 4.00 Å to 2.30 Å (violet curve), the wavefunction steadily acquires multireference character having three dominant electron configurations of , , and . These three electron configurations differ in their varying occupation numbers in and orbitals, both of which, in fact, describe the covalent interaction between the Fe(II) center and the O₂ ligand (Fig. 1). During this process, the weight of the first one plummets from 94% to 57%, while those of the last two rise from 0% to 21% and 12%, respectively, likely suggesting the progressively strengthened covalent Fe–O₂ interaction. However, the spin density and spin populations computed for ⁵B are analogous to those for ⁵A. This finding demonstrates that the electronic structure of ⁵B is nearly the same as that of ⁵A, indicative of weak covalent Fe–O₂ interaction of ⁵B. As the Fe–O₂ bond length further declines, one electron is shifted from the Fe(II) center to the O₂ ligand near the avoided crossing region accompanied by a rapid adjustment of the electronic structure of the nascent Fe(III) center to afford ⁵1 in the end.

CASSCF(12,9) calculations reveal that the quintet state of complex ⁵1 features strong multireference character and consists of three key electron configurations of (39%), (30%), and (27%). Apparently, this wavefunction is too complex to extract any chemically sensible bonding information from its predominant electron configurations. Specifically, the first one could be interpreted as a high spin Fe(III) center antiferromagnetically coupled to a doublet superoxo ligand , the second one as a quintet Fe(IV) center (S_Fe = 2) bound to a singlet peroxo ligand (S_O2²⁻ = 0), and the last one as a high spin Fe(II) center (S_Fe = 2) coordinated by a singlet O₂ ligand (S_O2 = 0). However, as demonstrated in Fig. 2, the spin density calculated for complex ⁵1 exhibits a spherical shape for the Fe center, reminiscent of the high spin Fe(III) center in complexes ⁷B and ⁷1, but different from the high spin Fe(II) center (S_Fe = 2) in complexes ⁷A and ⁵A. Moreover, the spin density reveals substantial negative spin primarily populating the O₂ orbital, thereby signaling that the O₂ ligand is best described as a doublet superoxo ligand . Based on these findings, we surmised that complex ⁵1 is best described as a high spin Fe(III) center antiferromagnetically coupled with a superoxo ligand affording an overall S_t = 2 state. Following the route used to interpret the electronic structure of ⁵A, such a bonding situation should be described by . This spin eigenfunction furnishes Ŝ_z expectation values of and for the Fe center and the O₂ ligand, respectively, congruent with the computed spin populations (4.41 for Fe and −0.61 for O₂), affirming the bonding interpretation.

Alternatively, to interpret the multireference wavefunction computed for complex ⁵1, valence bonding (VB) reading^68–71 was also employed as exemplified by the VB analysis of the electronic structure of the myoglobin O₂ adduct reported by Sason and coworkers.⁷² As shown in Fig. 3, upon localization of and , representing bonding and antibonding interactions of d_yz and , the ground state wavefunction re-expressed in the basis of the localized orbitals (d_yz and ) contains merely one electron configuration of (99%), whereas the original one in the basis of the natural orbitals exhibits strong multireference character having three dominant electron configurations as discussed above. It should be pointed out that both wavefunctions are equivalent because they are related by a unitary transformation of orbitals in the active space. Moreover, the computed spin density and spin populations reveal that the majority spin is located on the Fe center and the minority one of the opposite sign on the O₂ ligand. As such, the spin-up and spin-down electrons, which are spin paired, are located in the d_yz and MOs, respectively. Both findings therefore approve of the electronic-structure assignment presented above.

Fig. 3 d_yz and MOs in natural and localized forms with the electron configurations and their weights for ⁵1 complex; Fe^III–O₂⁻˙ denotes the electron configuration of , Fe^IV–O₂²⁻ represents the electron configuration of , and Fe^II–O₂ indicates the electron configuration of , and hydrogen atoms are omitted for clarity.

Clearly, in the case of complex ⁵1, VB analyses is as efficient as constructing two-fragment spin eigenfunction in interpreting its multireference wavefunction. However, this approach is not applicable to complex ⁵A where the absence of any covalent interactions between the Fe center and the O₂ ligand renders localization fruitless. It follows that VB analysis cannot be used to characterize multireference wavefunctions without strong electron correlations. In this respect, our simple approach does not have such a limitation.

The reaction of Fe(II)(BDPP) with O₂ in an S_t = 1 state also features two diabatic potential curves, one denotes a high spin Fe(II) center (S_Fe = 2) antiferromagnetically coupled with a triplet O₂ ligand (S_O2 = 1), and the other describes an intermediate spin Fe(III) center antiferromagnetically coupled with a doublet superoxo ligand . Although the triplet transformation is barrierless, it is highly endothermic by 30.2 kcal mol⁻¹. (See the ESI† for more details) and cannot contribute to the actual reaction. Taken together, O₂ association on a high spin ferrous center proceeds either on a septet or on a quintet surface, in line with the experimental findings.²⁰

O₂ binding on Co(II)(BDPP) affording complex 2

As shown in Fig. 4, starting from ²C at a Co–O₂ distance of 4.50 Å, the diabatic potential curve (magenta line) of the doublet state arising from an antiferromagnetic coupling of a high spin Co(II) center and a triplet O₂ ligand (S_O2 = 1) is constantly uphill as the Co–O₂ bond length decreases. Furthermore, starting from the Co–O₂ distance of 2.50 Å, we could find another diabatic potential curve (crimson line) that is best described as a low spin Co(III) ion (S_Co = 0) interacting with a doublet superoxide radical . Its energy first decreases for the Co–O₂ distance descending from 2.50 Å to 1.90 Å and then increases for the even shorter Co–O₂ distances; hence, a local minimal (²2) is located at the Co–O₂ length of ∼1.90 Å. Two diabatic potential curves interact in the vicinity of the Co–O₂ distance being 2.40 Å and generate an avoided crossing with a barrier of 24.7 kcal mol⁻¹ for the O₂ binding reaction affording complex ²2. Moreover, the entire process is slightly endothermic by 8.4 kcal mol⁻¹. In analogy to the O₂ association yielding ⁷1, our ab initio computations likely overestimate the barrier and the driving force, because experimentally complex ²2 is stable at low temperatures (<193 K) but releases O₂ upon elevating the temperature.²²

Fig. 4 Active spaces of complexes C (a) and 2 (b); the doublet potential energy surface for the reaction of Co(II)(BDPP) with O₂ as a function of the Co–O₂ bond distance obtained by CASSCF(13,9)/NEVPT2 computations. The dot lines denote diabatic potential curves, namely, the magenta line denotes a doublet Co^II–O₂ complex, interpreted as an Co^II center antiferromagnetically coupled (FC) with a triplet O₂ ligand (S_O2 = 1); the crimson line denotes a doublet Co^III–O₂⁻˙ complex, interpreted as an S_Co = 0 Co^III ion interacting with a doublet O₂⁻˙ ligand ; the black line denotes the adiabatic curve; spin density and Löwdin and Mulliken (in parentheses) spin populations on the Co center and the O₂ ligand calculated for complexes ²C (Co(II)(BDPP) + O²), ²D and ²2 on the reaction trajectory are shown at the bottom, yellow and red denote negative and positive spin density, respectively, and local spin populations are also listed in the ESI† (c); hydrogen atoms are omitted for clarity.

CASSCF(13,9) calculations of complex ²C produce a predominant electron configuration of accounting for 80% of the wavefunction. As demonstrated in Fig. 4, the computed spin density plot exhibits a double-layer donut shape of the O₂ ligand, which implies that unpaired electrons equally populate its orbitals. Therefore, both findings suggest that complex ²C is best interpreted as a high spin Co(II) center antiferromagnetically coupled to a triplet O₂ ligand (S_O2 = 1) yielding an overall ground state. Again, we constituted two-fragment spin eigenfunctions using CGCs to interpret this wavefunction and it reads . This spin eigenfunction gives Ŝ_z expectation values of and for the Co(II) center and the O₂ ligand, respectively. Indeed, both values are in line with the computed spin populations, which confirms the proposed electronic-structure assignment.

On the other hand, ab initio calculations suggest that complex ²2 features a leading configuration of (95%). Clearly, its electronic structure is best described as a low spin Co(III) center (S_Co = 0) bound to a doublet superoxo ligand , congruent with the computed spin density and spin populations of 0.01 for the Co center and of 0.98 for the O₂ moiety as demonstrated in Fig. 4.

Regarding the O₂ binding process furnishing ²2, as the Co–O₂ distance decreases from 4.50 Å to 2.30 Å (magenta curve), the wavefunction steadily acquires multireference character with three dominant electron configurations of , , and . These three electron configurations differ in their varying occupation numbers in and orbitals that describe the covalent interaction between the Co(II) center and O₂ ligand as depicted in Fig. 1. During this process, the weights of the last two terms rise from 0% to 34% and 17%, respectively, which comes at the price of the percentage of the first one that lowers from 80% to 22%, perhaps signaling the increased Co–O₂ bonding strength. However, the electronic structure of ²D is similar to that of ²C as evidenced by the comparable spin density and spin populations computed for both species, which suggests that the Co–O₂ covalent interaction is rather weak. As the Co–O₂ distance further declines, a one-electron transfer from the Co(II) center to the O₂ ligand occurs around the avoided crossing region concurrent with a rapid adjustment of the electronic structure of the Co(III) center furnishing complex ²2.

The doublet potential energy surface of the O₂ association process yielding ²2 is also constructed by two diabatic potential curves, one represents the reactant complex containing a high spin Co(II) center antiferromagnetically coupled to a triplet O₂ ligand (S_O2 = 1) yielding an overall ground state, and the other describes the product complex consisting of a low spin Co(III) ion (S_Co = 0) interacting with a doublet superoxide radical . For the reaction of Co(II)(BDPP) with O₂, the emergence of the avoided crossing region suggests that no intermediate other than ²2 can be observed and that the Co(II)–O₂ adduct does not exist, consistent with the experimental findings.²²

For O₂ binding on Fe(II)(BDPP) and Co(II)(BDPP), our ab initio calculations show that the former process has a lower barrier of 13.6 kcal mol⁻¹ to ⁷1 than the latter (24.7 kcal mol⁻¹) to ²2. This theoretical prediction is congruent with the experimental observation that the generation of ⁷1 is much faster than that of ²2.^20,22 Furthermore, both transformations feature an avoided crossing that realizes one-electron transfer from a divalent metal center to O₂ and furnishes a trivalent metal center and a superoxo ligand as the final product; consequently, the divalent–O₂ adduct cannot be formed.

As elaborated in the section of Introduction, except one case, reactions of O₂ with a wide variety of Fe(II) and Co(II) precursors reported thus far all generate Fe(III) or Co(III) superoxo products, which intimates that the avoided crossing found for the O₂ addition processes on Fe(II)(BDPP) and Co(II)(BDPP) occurs in those transformations as well. In fact, the O₂ association to related Mn(II)(BDPP) also produces a Mn(III)–O₂⁻˙ complex, instead of Mn(II)–O₂.²⁹ Likewise, the treatment of a series of Cu(I) species including [Cu(NH₂-TMPA)]⁺ (TMPA = tris-methyl pyridine amine),⁷³ [Cu(TMG₃tren)]⁺ (TMG₃tren = tris(tetramethylguanidino)tren)⁷⁴ and [Cu(Ar₃-TMPA)]⁺ (Ar = tpb, dpb and dtbpb)⁷⁵ with O₂ invariably affords Cu(II)–O₂⁻˙, whereas the authenticated O₂ adduct, Cu(I)–O₂, is yet to be disclosed. For more reactions of O₂ and metal complexes affording superoxo adducts, readers may refer to review articles.^76–80 However, the reaction of O₂ with Y257F-HPCA was proposed to furnish Fe(II)–O₂ rather than Fe(III)–O₂⁻˙. This proposition was solely based on ⁵⁷Fe Mössbauer measurements that give an unambiguous assignment of the Fe oxidation state, but the detailed information about the exact valence state of the O₂ ligand remains lacking. Therefore, thorough experimental research is required to verify the advocated electronic-structure assignment and disprove and/or substantiate our theoretical prediction.

Comments on DFT computations

For a given state distinguished by a total spin of S, M_s, its spin projection along the z axis, possesses 2S + 1 values (−S, −S + 1, …, S − 1, S) and constitutes 2S + 1 microstate labelled as |S,M_s〉. Among them, only |S,S〉 and |S, −S〉 could be represented by a single Slater determinant, whereas others have to be constructed as a linear combination of a series of Slater determinants, all having the same eigenvalue of M_s with respect to Ŝ_z, in order to be the eigenfunctions of Ŝ². Therefore, single-Slater-determinant approaches such as DFT could only deal with |S,S〉 and |S, −S〉. Both are equivalent except for the different eigenvalue for Ŝ_z; hence, typically merely |S,S〉 called the standard component is computed. A representative example is complex ²2 in Fig. 5, for which the spin populations of Co and O2 from B3LYP calculations with 〈Ŝ²〉 = 0.76 are −0.02 and 1.00, respectively, congruent with the anticipated spin eigenfunction of .

Fig. 5 Spin density and Löwdin and Mulliken (in parentheses) spin populations for the metal center and the O₂ ligand from DFT calculations; ligand atoms are omitted for clarity.

For systems featuring magnetic couplings of two fragments with their respective spins of S₁ and S₂, the resulting state may possess a total spin S_t = S₁ − S₂, S₁ − S₂ + 1, …, S₁ + S₂, (S₁ > S₂) according to the well-established rule of angular momentum addition. As elaborated above, the spin wavefunction, |S,M〉, could be built up by the spin eigenfunctions of two constituent fragments, |S₁,M₁〉 and |S₂,M₂〉. Of note, as required by the symmetry of CGCs, M = M₁ + M₂. Geometrically, this means that the z component of the vector sum of two angular momenta is equal to the algebraic sum of their z components. For ferromagnetic coupling with S_t = S₁ + S₂, the spin eigenfunction of its standard component, |S₁ + S₂, S₁ + S₂〉, could also be represented by a single Slater determinant constructed by the combination of two single Slater determinants distinguished by |S₁,S₁〉 and |S₂,S₂〉 for the two components into one. Thus, such systems can be properly treated with DFT, for instance, the septet reaction of O₂ with Fe(II)(BDPP). As demonstrated in Fig. 5, the spin populations on the Fe center and O₂ ligand estimated by B3LYP computations for ⁷A (3.62 for Fe(II) and 1.99 for O₂) and ⁷1 (4.10 for Fe(III) and 1.20 for O₂⁻˙) are all in line with those (3.86 for Fe(II) and 2.00 for O₂ in ⁷A, 4.74 for Fe(III) and 1.03 for O₂⁻˙ in ⁷1) derived from CASSCF(12,9) calculations. Moreover, the computed Ŝ² expectation values of 〈Ŝ²〉 = 12.00 for ⁷A and 〈Ŝ²〉 = 12.00 for ⁷1 indicate that both B3LYP solutions are devoid of considerable spin contamination and are eigenfunctions of Ŝ² having well-defined energies.

The situation is distinctly different for systems involving antiferromagnetic couplings where the resulting state has a total spin of S_t < S₁ + S₂, and the wavefunction of its standard component, |S_t,S_t〉, cannot be represented as a single determinant but has to be constructed as a linear combination of a range of Slater determinants with the eigenvalue of Ŝ_z being S_t to fulfill the requirement of spin symmetry. Therefore, in principle, DFT cannot be used to compute the electronic structures of such systems.

For one particular case with S_t = S₁ − S₂, its standard component |S₁ − S₂, S₁ − S₂〉 again should be built up as a linear combination of |S₁,S₁|S₂, −S₂〉, |S₁,S₁ − 1|S₂, −S₂ + 1〉, …, all featuring the same the Ŝ_z eigenvalue of S₁ − S₂. Importantly, the first one could be described as a single Slater determinant formed by the combination of two single Slater determinants characterized by |S₁,S₁〉 and |S₂, −S₂〉 for the two constituent fragments into one. Broken symmetry formalism exploits this feature and approximate |S₁ − S₂,S₁ − S₂〉 as |S₁,S₁|S₂, −S₂〉, which means just choosing the first component from the above series. Consequently, the resulting solution is no longer the eigenfunction of Ŝ² and spin density and spin populations derived from broken symmetry calculations are not correct. Generally speaking, nonnegligible spin contamination should be regarded as a characteristic of broken symmetry solutions that can be used to differentiate it from usual “spin-pure” solution discussed above. A case in point is complex ²C where spin populations of 2.59 for Co(II) and −1.99 for O₂ provided by broken symmetry B3LYP calculations are in line with the expected values for (Fig. 5). These erroneous values differ significantly from those (1.58 for Co(II) and −0.66 for O₂) derived from CASCSF(13,9) calculations, but nevertheless indicate that the local spins of the Co(II) center and the O₂ ligand are and 1, respectively, thus yielding a qualitatively correct bonding description of complex ²C. However, this DFT solution suffers from unacceptable spin contamination with 〈Ŝ²〉 = 2.77. As a consequence, its energy is not reliable because it is not the eigenfunction of the Hamiltonian operator, Ĥ, either. As shown in Fig. 5, complex ⁵1 also belongs to this category.

While for other resulting state with S₁ − S₂ < S_t < S₁ + S₂, its standard component cannot be approximately by a single Slater determinant any more. For example, for S_t = S₁ + S₂ − 1, the standard component characterized by |S₁ + S₂ − 1,S₁ + S₂ − 1〉 should be constructed as a linear combination of the wavefunction distinguished by |S₁,S₁ − 1|S₂,S₂〉 and |S₁,S₁|S₂,S₂ − 1〉. As elaborated above, both |S₁,S₁ − 1〉 and |S₂,S₂ − 1〉 cannot be represented by a single Slater determinant; hence, neither can |S₁,S₁ − 1|S₂,S₂〉 nor |S₁,S₁|S₂,S₂ − 1〉. Thus, DFT completely fails for this kind of systems as exemplified by complex ⁵A. As depicted in Fig. 5, spin populations on Fe(II) of 1.90 and O₂ of 1.99 for O₂ erroneously predicted by DFT calculations indicate that the converged B3LYP solution is best interpreted as a triplet Fe(II) center (S_Fe = 1) ferromagnetically coupled to a triplet O₂ ligand (S_O2 = 1); therefore, the local spins of Fe and O₂ are qualitatively wrong.

In summary, only systems possessing two ferromagnetically coupled fragments, DFT calculations can be used to compute their electronic structures and spin density and spin populations and give reliable energies; otherwise, DFT computations in principle produce incorrect electronic structures and spin density and spin populations. For states with S_t = S₁ − S₂, the electronic structure can be approximately computed by invoking broken symmetry formalism in the framework of DFT calculations. But interpreting the resulting spin density and populations requires particular caution. In particular, for open-shell singlets, the artificial spin density obtained from broken symmetry calculations is qualitatively incorrect, because there does not exist any spin density everywhere for diamagnetic systems irrespective of being closed-shell or open-shell singlets. For comparison, the potential energy surfaces computed by DFT are summarized in the ESI.†

Conclusions

CASSCF/NEVPT2 calculations were employed to explore the electronic-structure evolution over the course of O₂ binding on Co(II) and Fe(II) precursors affording Fe(III)–O₂⁻˙ (1) and Co(III)–O₂⁻˙ (2) complexes. For the former reaction on the septet and quintet potential surfaces, starting complex ^7,5A is best formulated as a high spin Fe(II) center (S_Fe = 2) ferromagnetically or magnetically coupled with a triplet O₂ ligand (S_O2 = 1), while product ^7,51 consists of a high spin Fe(III) center ferromagnetically or antiferromagnetically coupled a superoxo ligand . Regarding the latter reaction in an state, initial complex ²C could be best characterized as a high spin Co(II) center antiferromagnetically coupled to a triplet O₂ ligand (S_O2 = 1), while O₂ adduct ²2 contains a low spin Co(III) center (S_Co = 0) interacting with a superoxo ligand . As O₂ steadily approaches the divalent metal center, the diabatic energy curve of M(II)–O₂ constantly rises in energy and interacts with that of M(III)–O₂⁻˙ at a M–O₂ distance of around ∼2.40 Å. This interaction generates an avoided crossing that effects a spontaneous single-electron transfer from the divalent metal center to the incoming O₂ ligand, thereby leading to M(III)–O₂⁻˙ as the final product and ruling out the possibility of formation of a bona fide M(II)–O₂ complex. This view has been used to explain the reaction outcomes of related O₂ association processes.

To extract electronic-structure information from complicated multireference wavefunctions, we constructed two-fragment spin eigenfunctions with CGCs to rationalize computed spin density and fragment spin populations on metal centers and O₂ moieties and deduced corrected electronic-structure formulations. This approach is complementary to usual VB analyses as demonstrated by applications of both methods to analyze the electronic-structure of O₂ adduct ⁵1. This species features antiferromagnetic coupling of a high spin Fe(III) center and a superoxo ligand ; however, directly reading out electron configurations of the resulting authenticated multireference wavefunctions hardly gains any chemically sensible bonding information.

For systems involving metal–ligand or metal–metal antiferromagnetic couplings with S_t = S₁ − S₂, DFT calculations, despite giving a qualitatively right electronic structure, invariably suffer from spin contamination. Consequently, the resulting spin density and spin populations are not reliable. Furthermore, DFT fails completely for systems with S₁ − S₂ < S_t < S₁ − S₂ in predicting qualitatively correct electronic structures. Only for ferromagnetically coupled systems with S_t = S₁ + S₂, DFT is able to give exact electronic structures and energies.

Author contributions

S. Y. conceived the presented idea. P. Z. performed the computations and wrote the original draft. S. Y. and W. L. reviewed and edited the final draft. All authors discussed the results and contributed to the final manuscript.

Data availability

The data supporting this article have been included as part of the ESI.†

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The work conducted in Dalian was supported by the National Key R&D Program of China (No. 2021YFA1502803), the National Natural Science Foundation of China (92161204) and the DaIian Institute of Chemical Physics, Chinese Academy of Sciences (Grant: DICP I202312). W. L. is grateful for the financial support from the Ministry of Science and Technology of Taiwan (MOST 111-2113-M-003-008-MY3 to W.-Z. L.). We are indebted to Mr Haowei Chen for his helpful discussions.

Notes and references

1. M. Costas, M. P. Mehn, M. P. Jensen and L. Que Jr, Chem. Rev., 2004, 104, 939.
2. E. G. Kovaleva and J. D. Lipscomb, Nat. Chem. Biol., 2008, 4, 186.
3. E. I. Solomon, D. E. Heppner, E. M. Johnston, J. W. Ginsbach, J. Cirera, M. Qayyum, M. T. Kieber-Emmons, C. H. Kjaergaard, R. G. Hadt and L. Tian, Chem. Rev., 2014, 114, 3659.
4. C. E. Elwell, N. L. Gagnon, B. D. Neisen, D. Dhar, A. D. Spaeth, G. M. Yee and W. B. Tolman, Chem. Rev., 2017, 117, 2059.
5. X. Huang and J. T. Groves, Chem. Rev., 2017, 118, 2491.
6. A. J. Jasniewski and L. Que, Jr., Chem. Rev., 2018, 118, 2554.
7. L. Pauling and C. D. Coryell, Proc. Natl. Acad. Sci. U. S. A., 1936, 22, 210.
8. B. A. Springer, S. G. Sligar, J. S. Olson and G. N. Phillips Jr., Chem. Rev., 1994, 94, 699.
9. K. Shikama, Chem. Rev., 1998, 98, 1357.
10. J. E. Baldwin and M. Bradley, Chem. Rev., 1990, 90, 1079.
11. P. L. Roach, I. J. Clifton, V. Fülöp, K. Harlos, G. J. Barton, J. Hajdu, I. Andersson, C. J. Schofield and J. E. Baldwin, Nature, 1995, 375, 700.
12. P. L. Roach, I. J. Clifton, C. M. H. Hensgens, N. Shibata, C. J. Schofield, J. Hajdu and J. E. Baldwin, Nature, 1997, 387, 827.
13. M. M. Mbughuni, M. Chakrabarti, J. A. Hayden, E. L. Bominaar, M. P. Hendrich, E. Münck and J. D. Lipscomb, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 16788.
14. L. Que Jr. and R. Y. N. Ho, Chem. Rev., 1996, 96, 2607.
15. E. I. Solomon, T. C. Brunold, M. I. Davis, J. N. Kemsley, S. Lee, N. Lehnert, F. Neese, A. J. Skulan, Y. Yang and J. Zhou, Chem. Rev., 2000, 100, 235.
16. R. A. Baglia, J. P. T. Zaragoza and D. P. Goldberg, Chem. Rev., 2017, 117, 13320.
17. S. Yadav, V. Yadav, M. A. Siegler, P. Moënne-Loccoz, G. N. L. Jameson and D. P. Goldberg, J. Am. Chem. Soc., 2024, 146(12), 7915.
18. C. Winslow, H. Lee, M. J. Field, S. J. Teat and J. Rittle, J. Am. Chem. Soc., 2021, 143(34), 13686.
19. M. N. Blakely, M. A. Dedushko, P. C. Y. Poon, G. Villar-Acevedo and J. A. Kovacs, J. Am. Chem. Soc., 2019, 141(5), 1867.
20. J. Sacramento, T. Albert, M. Siegler, P. Moënne-Loccoz and D. P. Goldberg, Angew. Chem., Int. Ed., 2022, 61, e202111492.
21. C. Chiang, S. T. Kleespies, H. D. Stout, K. K. Meier, P. Li, E. L. Bominaar, L. Que Jr., E. Münck and W. Lee, J. Am. Chem. Soc., 2014, 136, 10846.
22. P. Zhang, Y. Tian, Z. Wang, W. Lee and S. Ye, Chem. – Eur. J., 2024, 30, e202400336.
23. C. Wang, H. Chang, Y. Lai, H. Fang, C. Li, H. Hsu, Z. Li, T. Lin, T. Kuo, F. Neese, S. Ye, Y. Chiang, M. Tsai, W. Liaw and W. Lee, J. Am. Chem. Soc., 2016, 138, 14186.
24. A. J. Fielding, J. D. Lipscomb and L. Que Jr., J. Am. Chem. Soc., 2012, 134, 796.
25. F. A. Walker and J. Bowen, J. Am. Chem. Soc., 1985, 107, 7632.
26. J. B. Gordon, A. C. Vilbert, M. A. Siegler, K. M. Lancaster, P. Moënne-Loccoz and D. P. Goldberg, J. Am. Chem. Soc., 2019, 141(5), 3641.
27. P. Kumar, S. V. Lindeman and A. T. Fiedler, J. Am. Chem. Soc., 2019, 141, 10984.
28. A. Huber, L. Müller, H. Elias, R. Klement and M. Valko, Eur. J. Inorg. Chem., 2005, 1459.
29. J. Yang, P. Li, X. Li, L. Xie, N. Wang, H. Lei, C. Zhang, W. Zhang, Y.-M. Lee, W. Zhang, R. Cao, S. Fukuzumi and W. Nam, Angew. Chem., Int. Ed., 2022, 61, e202208143.
30. Y. Lin, H. Hinrika Cramer, M. van Gastel, Y. Tsai, C. Chu, T. Kuo, I. Lee, S. Ye, E. Bill and W. Lee, Inorg. Chem., 2019, 58, 9756.
31. H. Pan, H. Chen, Z. Wu, P. Ge, S. Ye, G. Lee and H. Hsu, JACS Au, 2021, 1, 1389.
32. E. G. Kovaleva, M. B. Neibergall, S. Chakrabarty and J. D. Lipscomb, Acc. Chem. Res., 2007, 40, 475.
33. S. M. Adam, G. B. Wijeratne, P. J. Rogler, D. E. Diaz, D. A. Quist, J. J. Liu and K. D. Karlin, Chem. Rev., 2018, 118, 10840.
34. E. G. Kovaleva and J. D. Lipscomb, Science, 2007, 316, 453.201.
35. M. M. Mbughuni, K. K. Meier, E. Munck and J. D. Lipscomb, Biochemistry, 2012, 51, 8743.
36. M. Radon and K. Pierloot, J. Phys. Chem. A, 2008, 112, 11824.
37. H. Chen, M. Ikeda-Saito and S. Shaik, J. Am. Chem. Soc., 2008, 130, 14778.
38. S. Shaik and H. Chen, J. Biol. Inorg. Chem., 2011, 16, 841.
39. H. Chen, K. Cho, W. Lai, W. Nam and S. Shaik, J. Chem. Theory Comput., 2012, 8, 915.
40. G. J. Christian, S. Ye and F. Neese, Chem. Sci., 2012, 3, 1600.
41. S. Ye, C. Riplinger, A. Hansen, C. Krebs, J. M. Bollinger, Jr. and F. Neese, Chem. – Eur. J., 2012, 18, 6555.
42. K. P. Jensen, B. O. Roos and U. Ryde, J. Inorg. Biochem., 2005, 99, 45.
43. K. P. Kepp, ChemPhysChem, 2013, 14, 3551.
44. Q. M. Phung and K. Pierloot, Phys. Chem. Chem. Phys., 2018, 20, 17009.
45. L. Noodleman, J. Chem. Phys., 1981, 74, 5737.
46. J. Ribas-Arino and J. J. Novoa, Chem. Commun., 2007, 3160.
47. P. Malmqvist and B. O. Roos, Chem. Phys. Lett., 1989, 155, 189.
48. C. Angeli, R. Cimiraglia and J.-P. Malrieu, J. Chem. Phys., 2002, 117, 9138.
49. C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger and J. P. Malrieu, J. Chem. Phys., 2001, 114, 10252.
50. M. E. Rose, Elementary theory of angular momentum, Dover Publications Inc., New York, 1957.
51. A. D. Becke, J. Chem. Phys., 1993, 98, 5648.
52. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785.
53. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098.
54. A. Schäfer, C. Huber and R. Ahlrichs, J. Chem. Phys., 1994, 100, 5829.
55. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
56. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297.
57. A. Schäfer, H. Horn and R. Ahlrichs, J. Chem. Phys., 1992, 97, 2571.
58. S. Kossmann and F. Neese, J. Chem. Theory Comput., 2010, 6, 2325.
59. K. Eichkorn, O. Treutler, H. Öhm, M. Häser and R. Ahlrichs, Chem. Phys. Lett., 1995, 242, 652.
60. F. Weigend, Phys. Chem. Chem. Phys., 2006, 8, 1057.
61. V. Barone and M. Cossi, J. Phys. Chem. A, 1998, 102, 1995.
62. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456.
63. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
64. C. Kollmar, K. Sivalingam, B. Helmich-Paris, C. Angeli and F. Neese, J. Comput. Chem., 2019, 40, 1463.
65. F. Weigend, J. Comput. Chem., 2008, 29, 167.
66. F. Neese, F. Wennmohs, U. Becker and C. Riplinger, J. Chem. Phys., 2020, 152, 224108.
67. H. Chen, P. Zhang and S. Ye, Encyclopedia of Inorganic and Bioinorganic Chemistry, 2024 DOI:10.1002/9781119951438.eibc2883.
68. C. Angeli, C. J. Calzado, C. de Graafcd and R. Caballol, Phys. Chem. Chem. Phys., 2011, 13, 14617.
69. C. Angelia, R. Cimiragliaa and J. Malrieu, Mol. Phys., 2013, 111, 1069.
70. L. Blancafort, P. Celani, M. J. Bearpark and M. A. Robb, Theor. Chem. Acc., 2003, 110, 92.
71. J. Malrieu, C. Angeli and R. Cimiraglia, J. Chem. Educ., 2008, 85, 150.
72. H. Chen, M. Ikeda-Saito and S. Shaik, J. Am. Chem. Soc., 2008, 130, 14778.
73. D. Maiti, H. C. Fry, J. S. Woertink, M. A. Vance, E. I. Solomon and K. D. Karlin, J. Am. Chem. Soc., 2007, 129, 26.
74. C. Wrtele, E. Gaoutchenova, K. Harms, M. C. Holthausen, J. Sundermeyer and S. Schindler, Angew. Chem., Int. Ed., 2006, 45, 3867.
75. S. Y. Quek, S. Debnath, S. Laxmi, M. van Gastel, T. Krämer and J. England, J. Am. Chem. Soc., 2021, 143(47), 19731.
76. R. Trammell, K. Rajabimoghadam and I. Garcia-Bosch, Chem. Rev., 2019, 119, 2954.
77. C. E. Elwell, N. L. Gagnon, B. D. Neisen, D. Dhar, A. D. Spaeth, G. M. Yee and W. B. Tolman, Chem. Rev., 2017, 117, 2059.
78. S. Fukuzumi and K. D. Karlin, Coord. Chem. Rev., 2013, 257, 187.
79. S. Hong, Y. Lee, K. Ray and W. Nam, Coord. Chem. Rev., 2017, 334, 25.
80. H. Noh and J. Cho, Coord. Chem. Rev., 2019, 382, 126.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4cp02915a

---