Revealing the nature of electron correlation in transition metal complexes with symmetry breaking and chemical intuition

The ab initio modeling of transition metals is a long-sought goal, and to date, the widespread use of Density Functional Theory (DFT) in chemical and biological realms has led to many important discoveries.^1–3 However, the robust accuracy that modern density functionals⁴ and single-reference (SR) wavefunction methods such as coupled cluster with singles, doubles, and perturbative triples (CCSD(T)) have shown for closed-shell organic systems may not necessarily carry over to transition metal compounds. Indeed, there are a number of aspects relevant to transition metals that are either less or not at all relevant for typical organic molecules, such as relativistic effects, a balanced treatment of solvation for redox species, and the likely possibility of converging to extrema other than the global minimum depending on the choice of Self-Consistent Field (SCF) algorithm and/or initial guess. Yet, as these are, in principle, well-defined problems with well-defined (albeit often not perfect) solutions,^5–8 it is often the case that terms such as “strong correlation” and “multi-reference (MR) character” are used as generic explanations in the face of unsystematically erroneous (and thus, apparently mysterious) predictions.

Let us start with some accepted working definitions: By multi-reference (MR) character,⁹ we mean that the wavefunction of a chemical system contains more than one determinant with significant weight. Strong correlation is not necessarily implied by MR character as these determinants may only interact very weakly [e.g., unrestricted Hartree–Fock (UHF) solutions for H₂ in the dissociation limit] but can be said to occur when the MR nature of a system leads to the breakdown or lack of convergence of SR perturbation theory (PT). Herein, we focus only on MR character, also known as “static correlation.” In the cases of stretched H₂ and biradicaloids such as O₂ or benzyne isomers, static correlation can be encountered when HOMO–LUMO and singlet–triplet gaps approach (but do not exactly reach) zero. In these situations, the lowest energy singlet wavefunctions have substantial open-shell character and thus require two determinants for a proper qualitative description. Within DFT, Yang and co-workers have likened such a situation to a fractional spin error.^10,11 Variational single-determinant methods will, in these cases, exhibit spin symmetry breaking (SSB) with the expectation value of the spin operator, ⟨S²⟩, between 0 and 1 (exactly 1 for dissociated H₂ and a perfect biradical). The second aspect that leads to MR character arises from the need of a wavefunction to transform as an irreducible representation of the molecular point group. Variational optimization of a single-determinant wavefunction, here, can only yield one of, at times, multiple configurations that when superposed yield the correct symmetry, and spatial symmetry breaking (spatial SB) then occurs.¹² Simple examples can be found in HF calculations of stretched diatomics such as $F2+$⁠.¹³ While we note that independent-particle theories can also break other intrinsic symmetries of the electronic Hamiltonian,^14,15 e.g., those related to complex conjugation and time-reversal, in this work, we will focus on spin and, to a lesser extent, on spatial SB and their relationship to MR character.

Transition metal atoms with partially filled valence d shells can exhibit substantial MR character due to the close energetic spacing of many low-lying electronic states in contrast to the spectra of typical closed-shell organic compounds.¹⁶ Transition metal diatomics have been studied extensively in recent years,^16–25 and recent precise experiments from Morse²⁶ have enabled meaningful comparisons between theory and experiment. Errors in the computed bond dissociation energies vs experiment as large as 30 kcal/mol were reported from DFT, and neither the “gold standard” CCSD(T)²⁷ nor a MR variant²⁸ showed reliable accuracy for all species.²² It appears that some diatomics were “easier” to treat than others; however, common MR diagnostics did not correlate with the accuracy of any method and were shown to point to inconsistent conclusions depending on the particular diagnostics and furthermore typically lack physically interpretive values.^22,29 Many of these systems were later investigated with quantum Monte Carlo methods,³⁰ which yielded very accurate results vs experiment when consistently employing multi-determinant trial wavefunctions. This level of accuracy, but with respect to near-exact benchmark calculations, could only be attained via high orders of coupled cluster (CC) theory, which are not scalable to larger systems.²⁴ Reference 24 implies, as might be expected, that higher orders of CC theory are needed to describe the increasing number of bonds; however, there were many exceptions to this trend. Considering even the simplest bonding motif (i.e., a single bond) as found in metal hydride diatomics, many glaring irregularities could be found [convergence issues or ∼23 kJ/mol inaccuracies at the CCSD(T) level], and no explanation on simple physical grounds has been provided.

Metal complexes with a higher coordination number are of relatively greater interest to chemists as the coordination environment resembles that of realistic transition metal catalysts or active sites in biology. 3d-containing complexes with partially filled d shells are of particular interest due to their earth-abundance³¹ and, from the perspectives of electronic structure and reactivity, due to the competing accessibility of low spin (LS) and high spin (HS) configurations. As they are non-linear, Jahn–Teller distortions³² occur to avoid the unequal occupancy of degenerate orbitals, i.e., to break orbital degeneracies due to spatial symmetry. However, the presence of MR character has been suggested in Ref. 33 as a way of rationalizing large errors vs experiment of SR methods. Indeed, within a set of apparently simple and similar metal cations from Ref. 34, B3LYP, B97, and DLPNO-CCSD(T) performed very well for a subset, but very poorly for another. The elucidation of underlying reasons for such aberrant behavior, most notably for complexes such as Fe$(NH3)4+$ and Mn$(NH3)4+$⁠, is a primary motivation for the present study. We then proceed to highlight the MR character in states that are LS due to antiferromagnetic coupling, as found in mono-metal complexes with redox-noninnocent ligands and oxygen-bridged Mn(III)/Mn(IV) dimers.

To provide more details on the suitability of commonly used MR diagnostics for transition metal systems, Wilson and co-workers have explored metrics derived from CC and Configuration Interaction (CI), e.g., T1 and D1 (Frobenius and matrix two-norm of the singles amplitude t₁ in CCSD, respectively), $C02$ (square of the HF configuration in the CISD wavefunction or of the leading configuration state function in active-space methods), and %TAE (triples contribution to the atomization energy). Those authors concluded that CC diagnostics are insufficient, $C02$ from the Complete Active Space Self-Consistent Field (CASSCF³⁵) method is only reliable when a full active space can be considered, and %TAE fails for weakly bound systems such as Zn₂. No linear correlation between any investigated diagnostics and the accuracy of the correlation consistent composite approach³⁶ was found.²⁹ We note that even if they were reliable, diagnostics based on CCSD amplitudes and the relative contribution of the (T) component would require calculations that scale as the sixth and seventh power of the system size, respectively. Another quantitative approach that can shed light on the presence of MR character essentially aims to quantify the number of unpaired electrons in a molecule using the eigenvectors and eigenvalues of the one-particle density matrix, known as natural orbitals (NOs) and NO occupation numbers (NOONs), respectively.^37,38 More recent extensions supplement a scalar quantity with real-space descriptors identifying the local regions of the molecule responsible for static correlation,^39,40 and in some cases, canonical Kohn–Sham orbitals are used in combination with a finite-temperature DFT formalism.⁴¹ 

The applicability of CASSCF (e.g., in Ref. 42) is limited because exact CI solvers are feasible only for systems with less than ∼24 active orbitals.⁴³ Approximate solvers, e.g., DMRG⁴⁴ or selected CI,^45–48 enable the use of a larger number of orbitals,^49,50 yet they are still sensitive to active space specification and nonetheless still scale exponentially. This active-space dependence renders CASSCF calculations, including RASSCF methods,⁵¹ rather difficult to properly converge. More generally, misleading conclusions can be drawn when active spaces are not sufficiently large.⁴⁴ References 52 and 53 illustrate this difficulty with regard to the MoFe₇S₉C catalytic center of nitrogenase. Furthermore, the delicate interplay between MR character and dynamic correlation can sometimes make the choice of sufficiently large active spaces quite unintuitive (e.g., active spaces with double d shells).

For the ground-state (GS) of non-linear molecules (with geometries optimized without symmetry constraints), we will explore, in this work, the use of SSB in unrestricted single-determinant wavefunctions formed from orbitals from Kohn–Sham DFT or κ-regularized orbital-optimized MP2 (κ-UOOMP2)⁵⁴ to detect the presence of MR character in transition metal complexes. It is known⁵⁵ that DFT orbitals, which are optimized in the presence of mean-field electron correlations, do not break spin or spatial symmetries as easily as HF orbitals. With regard to κ-UOOMP2, the inclusion of pairwise additive electron correlation in the orbital optimization has a similar effect and has been used to diagnose MR character in organic molecules and fullerenes^56,57 and in one transition metal containing system (neutral iron porphyrin).⁵⁸ Moreover, given the plethora of MR diagnostics, we take a more appealing route using single-particle orbital theory, drawing on chemical concepts from ligand-field and molecular orbital theories.^59–61 We note that these models have recently been connected with sophisticated ab initio quantum-chemical methods and successfully applied to corroborate a variety of experimental observations.^62,63 We take one step further now and seek to explore how ligand-field theoretical arguments can be used to understand and predict the presence of MR character in conjunction with the quantitative use of SSB in SR theories that include some electron correlations. Thus, one main goal of this work is to provide, through a series of illustrative and chemically relevant examples, a link between well-established qualitative concepts and the occurrence of a MR wavefunction.

Correlation is often partitioned into “static” and “dynamic” contributions. If we assume, as most do, that the former is synonymous with MR character, there remains the question of which type of correlation is responsible for the bulk of the errors of the commonly used SR quantum chemical methods. Admittedly, the extent to which static correlation is relevant for transition metal complexes is, at present, unclear. It is often assumed that there is a large degree of dynamic correlation, but can this quantity be connected with the physical properties of the bonding exhibited by transition metal compounds? If so, what order of PT or CC theory is needed? We endeavor to propose a more nuanced description of electron correlation than the strong vs weak distinction, which, at present, permeates the field, and we hope this will guide practitioners in their selection of appropriate quantum-chemical methods. As an example, we focus on decidedly SR metal-carbonyl complexes and suggest that the pairwise additive correlation energy expression of MP2 is inadequate for the quantitative description of dative interactions consisting of σ-donation and π-backbonding, and thus caution against the use of MP2-based double-hybrid density functionals (DHDFs) for organometallic complexes.

A variational theory is the one that minimizes

with respect to one or more parameters. In HF theory, Ψ is a Slater determinant (antisymmetrized occupied orbitals) and the energy is minimized with respect to the orbital coefficients. In the CASSCF method, Ψ is a linear combination of Slater determinants constructed from all possible electron configurations within a specified set of orbitals called the active space. Both the CI and orbital coefficients are treated as variational parameters. A CASSCF calculation with the full orbital space as the active space is equivalent to full CI (FCI) and is exact within the basis set employed. Approximate FCI calculations are performed with the Adaptive Sampling Configuration Interaction (ASCI) method,^47,64 a selected CI approach that iteratively improves a fixed-size CI wavefunction by selecting the most important configurations in the Hilbert space (via a first-order PT based selection rule). The convergence of observables such as energy or molecular properties can be gauged by comparing results from ASCI wavefunctions of increasing size or by considering the second-order perturbative correction⁶⁵ to the ASCI energy.

κ-OOMP2 introduces a regularizer to OOMP2 (which optimizes orbitals to minimize E₀ + E⁽²⁾, where E⁽²⁾ is not variational) to prevent the MP2 energy from diverging and to ensure a continuous restricted (R) to unrestricted (U) transition during bond breaking and is presented and discussed in detail in Refs. 54 and 56. The κ-MP2 total energy, which is minimized with respect to orbital rotations in κ-OOMP2, is

where i, j and a, b represent the occupied and virtual orbitals, respectively. $Δijab=ϵa+ϵb−ϵi−ϵj$⁠, and the antisymmetrized two-electron integrals ⟨ij‖ab⟩ are defined as ⟨ij|ab⟩ − ⟨ij|ba⟩. Equation (2) reveals that E_κ-MP2(κ → 0) = E_HF. SSB in the HF solution can occur even for stable molecules at equilibrium bond lengths, which should be well-described by a closed-shell configuration. This situation has been termed “artificial” SSB.⁵⁶ In the limit of κ → ∞, unregularized OOMP2 is obtained, which does not have Coulson–Fischer (CF) points⁶⁶ due to divergences accompanying small HOMO–LUMO gaps.⁶⁷ Therefore, in this limit, the theory is strongly though artificially biased in favor of spin symmetry restoration. As the tendency to break spin symmetry clearly depends on the choice of regularization parameter, an optimal κ has been chosen in light of two criteria. The first is physically motivated, requiring CF points for single, double, and triple carbon–carbon bonds to occur at increasing bond distances. The second is empirically motivated, selecting κ such that errors for reference reaction datasets are minimized: For thermochemistry, this leads to κ = 1.45 $Eh−1$⁠.⁵⁴ 

We note that a situation analogous to the latter is relevant to the specification of global hybrid DFT functionals with the exchange part of the correlation energy given as

where $ExHF$ is the exact HF exchange (EXX) functional and $ExSL$ is the semilocal exchange functional. The SSB behavior of a global hybrid can be modulated by the fraction of EXX employed, i.e., the parameter a. One advantage of κ-UOOMP2 vs DFT is that the former is free of self-interaction error, in which an electron can spuriously interact with itself, and its many-body generalization known as delocalization error (DE).^68,69 It is well-known that DE tends to (artificially) reduce the extent of SSB (and spatial SB) with pure and global hybrid (i.e., those with a modest amount of EXX, e.g., B3LYP) functionals in cases such as stretched hydrogen fluoride, and is also relevant to the description of ligand-to-metal dative bonding.

For some of the smaller molecules in this study, we will utilize a variant of CC with double excitations that variationally optimizes a set of orbitals rather than including a singles term in the cluster operator (OOCCD).⁷⁰ 

In 1955, Löwdin derived⁷¹ that

where the two-particle density matrix, Γ, is normalized to $N(N−1)2$⁠. The expressions of the two-particle density matrix are known for ROHF- and UHF-based theories along with simple density functional theories.^72,73 The UHF expression is of particular relevance to the present work and can be written as

where $S=Nα−Nβ2$ and $Sij̄$ is the αβ overlap integral between the i-th spin-up orbital and the j-th spin-down orbital. We emphasize that the ⟨S²⟩ value reported for DFT calculations is exact only for the fictitious non-interacting Kohn–Sham system, as the occupied Kohn–Sham orbitals are used instead of those from UHF to compute Eq. (5).

For κ-UOOMP2 and UOOCCD, unless otherwise mentioned, Eq. (5) is used with the orbitals resulting from the minimization procedure. A rigorous expression for ⟨S²⟩_UOOCCD can be found in Ref. 74. Regarding κ-UOOMP2, with the first-order MP wavefunction ψ₁ in hand, the expectation value of the spin operator can be rigorously obtained as

assuming real amplitudes and orbitals. NOONs can be obtained by diagonalizing the one-particle density matrix formed from ψ₁. As can be seen from Eq. (2), the amplitudes of the doubly excited states in the κ-OOMP2 wavefunction are, by construction, encouraged to be small (i.e., when $Δijab$ is small, so is the regularizing term in parentheses and thus the regularized amplitude). Thus, we expect that in most cases, ⟨S²⟩ will be similar when computed either by Eq. (5) from κ-UOOMP2 orbitals or by Eq. (6) with the perturbed wavefunction. Similar considerations apply to the NOONs.

All calculations were performed with Q-Chem 5.2,⁷⁵ except for CASSCF calculations and geometry optimizations for Fe(II) X₆ and the metal carbonyls not in Ref. 33, which were performed with Orca.⁷⁶ These geometry optimizations utilized the DKH Hamiltonian and DK basis sets, for consistency with previous work.³³ CASSCF calculations were initialized from DFT orbitals (in most cases, obtained with the B3LYP functional). For the ASCI calculations, we use approximate NOs^47,64 or MCSCF-like orbital optimization within the active space,⁷⁷ which leads to more compact CI wavefunctions and thus lower variational energies for a given wavefunction size. Further details about ASCI can be found in Refs. 47, 64, and 65. ASCI NOONs were computed from the variational CI wavefunction alone without any perturbative corrections. For ASCI calculations of the hydrogen fluoride molecule and transition metal hydrides, the second-order PT correction to the energy (E_PT2) is smaller than 10⁻³ Ha (i.e., 0.6 kcal/mol) in magnitude. This high level of convergence was not possible for para-benzyne, and so we used the full-valence active space of 28 orbitals. The NOONs are converged by 4 × 10⁶ ASCI determinants. ASCI calculations were also performed for Fe(H₂$O)62+$ and Fe(⁠$CO)62+$⁠, with convergence shown in the supplementary material.

To obtain ⟨S²⟩ values and NOONs, the def2-SV(P) basis set is used⁷⁸(unless mentioned otherwise). UOOCCD calculations use the frozen-core approximation. While we are well-aware that for ground-state properties, the quantitative accuracy of energetic quantities vs experimental measurements requires much larger basis sets, the def2-SV(P) basis set is found to produce the qualitative descriptions required, and at a much-reduced computational cost that would be expected of a practical diagnostic tool.

The GDM algorithm⁷⁹ is used as the default SCF solver for HF, DFT, κ-UOOMP2, and UOOCCD calculations. We confirmed the stability, with respect to orbital rotations, of HF and DFT solutions. For the κ-UOOMP2 calculations, the resolution of the identity approximation was employed with the corresponding auxiliary basis set. The integration grid for DFT calculations consists of 99 radial spheres each with 590 points.⁸⁰ 

We employed an energy decomposition analysis (EDA) based on absolutely localized molecular orbitals to decompose the binding energies. The original scheme decomposes the interaction energy into frozen, polarization, and charge transfer contributions.^81,82 In order to reveal more insight into the nature of the bidirectional charge transfer of some metal–ligand bonds, we augmented the analysis with a further decomposition of the charge transfer energy using the newly developed variational forward–backward decomposition.⁸³ 

In this section, we investigate three simple systems that can be made to exhibit static correlation and compare the frontier NOONs computed from various approximate unrestricted methods with exact or near-exact reference values from ASCI calculations. Analytical expressions for H₂⁸⁴ show that the lowest unoccupied natural orbital (LUNO) occupation number for a single Slater determinant is related to ⟨S²⟩ by

While this relation is exact only for H₂, the ability of a SR method to produce, via spin polarization, (fractional) NOONs in agreement with exact values indicates that the presence of SSB reliably reflects the genuine MR character.

In Fig. 1, we plot n_LUNO vs bond length for three DFT functionals with varying amounts of EXX (0%, 20%, and 50% for BLYP, B3LYP, B5050LYP, respectively), HF, and κ-UOOMP2 vs exact results (obtained via CISD). The exact results show the fractional occupation of the antibonding orbital even at short bond lengths, whereas the other methods have zero occupation of the LUNO until beyond the CF point. The CF point is pushed to longer bond lengths going from HF, κ-UOOMP2, to more pure DFT functionals. While all methods provide a satisfactory description of H₂ in the dissociation limit, we note that in the region beyond ∼1.5 Å, n_LUNO as approximated using DFT orbitals appears to be closer to the exact value than when (rigorously) derived from HF and κ-UOOMP2 theories.

As a second test, we consider para-benzyne, a prototypical biradicaloid.^85,86 Using the cc-pVDZ basis and the geometry from Ref. 87, we compute reference highest occupied natural orbital (HONO) and LUNO populations with full valence CASSCF (28 electrons in 28 orbitals, using the ASCI solver⁷⁷), which are compared with values from other methods in Fig. 2. As found above for H₂, adding exact HF exchange systematically shifts the DFT predictions toward the HF result. We note that both B3LYP and κ-UOOMP2 produce NOONs in good agreement with the ASCI estimates.

While the above two examples would suggest that B3LYP is an excellent choice to describe SSB and fractional NOONs in these MR situations, its susceptibility to DE warrants caution in polar systems such as stretched hydrogen fluoride. In HF theory, the one-particle self-interaction error is canceled exactly—this is also the case in MP2 theories. Figure 3 shows that the LUNO occupation at the B3LYP level at long bond lengths is significantly lower than the ASCI value because DE leads to contamination of the unrestricted solution with closed-shell, ionic contributions. Adding a higher percentage of EXX into the functional form provides the nonlocal orbital-dependent exchange necessary to describe the derivative discontinuity, and DE can thus be reduced by global hybrids with a larger amount of EXX (e.g., B5050LYP) or range-separated hybrid functionals (e.g., CAM-B3LYP⁸⁸).^89,90 B5050LYP provides a substantial improvement over B3LYP, and CAM-B3LYP further improves the predicted NOONs.

In the context of transition metal complexes, DE can have a dramatic effect on the predicted covalency or dative nature of metal–ligand interactions. It has been shown that calculated spin-densities and derived electron-nuclear hyperfine coupling constants (which can be measured by EPR spectroscopy) are sensitive to DE and can lead to extreme variations in predicted paramagnetic nuclear magnetic resonance shifts in the range of $O$(1000) ppm.^91,92 Fortunately, the presence of DE can be diagnosed in straightforward ways, implied, e.g., by deviations from straight-line behavior in plots of energy vs fractional occupation number^69,92 (see Ref. 93 for the analysis of the fluoride anion relevant to dissociated hydrogen-fluoride). In such cases, methods such as κ-UOOMP2, B5050LYP, and range-separated hybrid functionals are to be preferred over pure or typical global hybrids such as BLYP or B3LYP, respectively.

Finally, regarding Fig. 3, we note that using the one-particle density matrix from the κ-UOOMP2 wavefunction to compute n_LUNO tracked the value obtained from using the κ-UOOMP2 orbitals in Eq. (5) nearly indistinguishably, except at bond lengths before the CF point. In this region, n_LUNO is necessarily zero for a single-determinant wavefunction and can be non-zero due to the doubly excited configurations in the MP1 wavefunction.

In this section, we investigate the 3d transition metal hydrides and will glean insights into the connection between static correlation and symmetry breaking. These systems are small enough such that near-exact wavefunctions can be obtained, and we use ASCI to converge FCI-quality wavefunctions (with a Ne frozen-core for the metal atom). As mentioned in the Introduction, 3d metal diatomics can possess surprisingly complicated electronic structures, resulting in a host of literature showing that SR methods (and some MR ones) have pronounced difficulty in predicting experimental thermochemistry. As will be discussed below, metal hydrides are linear molecules, point group C_∞v, for which the Jahn–Teller theorem does not apply. Therefore, MR character can arise in the GS due to both spatial and spin symmetry requirements.

Figure 4 shows a qualitative schematic of an expected molecular orbital (MO) diagram for metal hydride compounds. The shell of formally non-bonding orbitals (consisting of one σ, two π, and two δ orbitals) are, for now, treated as a single degenerate shell on the grounds of very small energetic splittings between the five orbitals. For our analysis of SSB, this will be adequate; however, a closer look is required for our analysis of spatial SB in TiH and CoH.

Table I shows the ⟨S²⟩ values from UHF, UDFT, κ-UOOMP2, and UOOCCD, compared with the exact spin quantum number. Comparing with the MO diagrams in Fig. 4, the first thing to notice is that a number of diatomics that do not show SSB are HS states (VH, CrH, and MnH). The HS states are known to be of SR nature.⁹⁴ SSB can only occur when there are paired valence electrons, as in the Fe–Zn hydrides. In addition, a higher spin state must be sufficiently close in energy to mix into the unrestricted SR wavefunction of the LS state. In the single-particle/MO picture, this means that the LUMO in the dominant electronic configuration must be low-lying enough to be occupied (via spin-flip excitation) in the HS state. As will be discussed further in Sec. III C, tetrahedral and octahedral metal complexes exhibit a splitting between t_2g and e_g levels known as 10Dq. Below, we discuss the analogs of this for linear 3d metal hydrides.

As shown in the red box of Fig. 4, the analog of 10Dq, which can modulate SSB for FeH, CoH, NiH, and CuH, is the gap between the non-bonding 4s and 3d metal orbitals with the antibonding σ*. The lowest excitation energies, as (roughly) estimated via time-dependent density functional theory (TDDFT) within the Tamm–Dancoff approximation (TDA) with the B3LYP functional and def2-SV(P) basis, for CuH and ZnH are 2.2 and 4.4 eV, respectively, suggesting that for these diatomics (though especially the latter), this gap is too large for HS states to compete. In contrast, the spin-contamination found in FeH, CoH, and NiH with all dynamically correlated SR theories considered suggests that, for these three species, this gap is small enough such that SSB can occur. Going from CuH to FeH, the fractional population of the σ* NO, as given by ASCI calculations and shown in Table II, monotonically increases from 0.06 to 0.35. Correspondingly, the gap between the non-bonding manifold and the σ* orbital can be expected to decrease in this direction due to (i) the decrease in metal atomic electronegativity, which, in the MO picture, increases the energy of the metal AOs and, as a result, the metal non-bonding MOs relative to the σ* orbital (which has partial ligand character), and (ii) the increase in metal atomic radius, which, according to ligand-field theory, should decrease the splitting of bonding and antibonding orbitals upon complexation with the hydride. As a result, the HS state formed by populating σ* becomes increasingly more energetically competitive with the LS state, resulting in proportional SSB.

For ScH, which is the only diatomic in the left-half of the row with significant post-HF SSB, the analog of the 10Dq parameter corresponds to the gap between the 4s and 3d MOs (0.4 eV with TDDFT/TDA). In other words, the intruding HS state responsible for spin-contamination is formed when the spin-down electron in the doubly occupied 4s orbital undergoes a spin-flip excitation into the energetically proximate non-bonding 3d shell. For ScH at the UOOCCD level of theory, we note that computing ⟨S²⟩ via Eq. (5) with the optimized orbitals led to a value of 1, whereas using the full UOOCCD wavefunction yielded ⟨S²⟩ = 1.75. The RooCCD energy is found to be lower than that of UOOCCD (but slightly higher than the ASCI energy), which reflects the fact that CC wavefunctions are highly sensitive to (and can be adversely affected by) spin-contamination in the reference state.

In light of the four half-occupied ASCI-computed NOONs shown in Table III, the lack of SSB in TiH is rather provocative. Indeed, due to the quartet multiplicity, there are no paired electrons in the 4s-3d MO manifold, which would imply no SSB, which until now has implied a SR wavefunction. TiH’s four NOONs of 0.5 can be explained in the context of spatial, rather than spin, SB as follows. The ⁴Φ term symbol denotes a doubly degenerate (E) irreducible representation (⁴Φ_x and ⁴Φ_y) where, e.g., $Φx4=12|σδ+πx〉−|σδ−πy〉$⁠.¹⁶ The corresponding electron configurations are illustrated in Fig. 5, where the wavefunction is a linear combination of two configurations varying in their occupation of the π and δ orbitals. The term symbol of the wavefunction for CoH leads to an analogous situation but for two holes rather than two electrons, e.g., $Φx3=12|σ2πy2δ−2πxδ+〉−|σ2πx2δ+2πyδ−〉$⁠,¹⁶ also illustrated in Fig. 5. Thus, for TiH and CoH, the four 0.5 and 1.5 NOONs, respectively, imply the MR character due to spatial symmetry. This can occur in the absence of SSB (TiH) or in addition to it (CoH).

For the rest of this paper, we turn to 3d metal complexes with higher coordination numbers, for which Jahn–Teller distortions readily lift any degenerate electron configurations due to spatial symmetry. In this section, we introduce the correspondence between SSB and the magnitude of the ligand-field paramater, 10Dq, which denotes the splitting between t_2g and e_g orbitals for tetrahedral (T_d) and octahedral (O_h) complexes. Whether or not this can be interpreted as a marker for MR character will be postponed to Sec. III D. Our discussion will center around O_h complexes, and in particular Fe(II) L₆, which implies a d⁶ configuration. In the O_h field, which yields a three below two d orbital ligand-field splitting, the LS state is a singlet with each of the three t_2g orbitals doubly occupied, whereas the HS state is a quintet formed by unpairing and spin-flipping two electrons from the t_2g to the e_g orbitals. These states along with the definition of 10Dq are shown in Fig. 6. A small (large) 10Dq value generally results in the GS being HS (LS), respectively, though a complete analysis involves examining the delicate balance between (i) the energetic cost to promote from t_2g to e_g (10Dq), (ii) the energetic cost to pair two opposite-spin electrons in the same t_2g orbital, and (iii) the stabilizing exchange interaction between the same-spin electrons.

It is well-known that the 10Dq parameter can be modulated by a ligand’s position on the spectrochemical series, which reflects whether the metal (M)–ligand (L) interaction is characterized by L-to-M σ donation only or, in addition, L-to-M π donation or M-to-L π backbonding. MO theory provides an intuitive model, which corroborates the experimentally determined trend that 10Dq gets smaller going from π-accepting to σ-donating-only to π-donating ligands.⁹⁵ 

Table IV shows the calculated ⟨S²⟩ values at the UHF, UB3LYP, UCAM-B3LYP, UB5050LYP, and κ-UOOMP2 levels of theory for O_h Fe(II) complexes in the LS state representing a range of ligand-field strengths. The CO ligand is perhaps the strongest-field π-acceptor in the spectrochemical series, and there is no SSB at any level of theory (not even UHF). We attribute this to 10Dq being sufficiently large such that the HS state is energetically inaccessible and therefore unable to mix into the LS wavefunction. This implies that the LS state is SR with very strong-field/π-accepting ligands (indeed, as will be argued in Sec. III F, the inaccuracy of MP2-based methods for organometallic thermochemistry is due to other reasons). As the ligand-field strength is attenuated (NH₃ and H₂O are σ-donation only ligands and the halides are π-donor ligands), the magnitude of 10Dq decreases, and the increasing accessibility of the HS configuration manifests as deviations from exact ⟨S²⟩ values in the LS state.

The trend of increasing SSB going toward weak-field ligands along the spectrochemical series is present in all methods investigated. In fact, all DFT and κ-UOOMP2 methods yield qualitatively similar values of ⟨S²⟩ [with the exception of the artificial SSB from the B5050LYP functional for Fe(⁠$NH3)62+$⁠, which, as might be expected, is nearly half the value found from UHF]. We have used a polarizable continuum model (ϵ = 78.4) for the DFT calculations of all anionic species, as motivated by the interesting case of the complex with fluoride (F⁻) ligands: In the gas phase, the ⟨S²⟩ value from B3LYP (1.74) exceeds that from UHF (1.49). The fluorine anion is known to have a positive HOMO eigenvalue when employing functionals such as B3LYP in all but exceedingly large basis sets,⁹⁷ and the curvature in a plot of the energy as a function of fractional occupation number is a hallmark of DE, which further encourages long tails in the radial charge density.^93,98 Indeed, for Fe$F64−$ in the gas phase with B3LYP, 21 of the 42 occupied Kohn–Sham orbital eigenvalues were greater than zero. When improving (slightly) the description of a continuum orbital via the def2-SVPD basis set, ⟨S²⟩ is further increased to 1.94, which reflects the expected narrowing of the singlet–triplet gap as continuum orbitals become more involved. Using the dielectric characteristic of water solvent, only the HOMO eigenvalue remained positive with a value of 0.01, and the resulting ⟨S²⟩ decreased from 1.74 to 1.03, now below the UHF value and in agreement with functionals with more EXX (and thus less DE) and κ-UOOMP2, which is free of DE. We note that increasing ϵ and the flexibility of the basis set can make all occupied B3LYP eigenvalues negative, but we did not find any large effect on the SSB behavior (e.g., a HOMO eigenvalue of −0.11 is obtained with ϵ = 1000 and the def2-SVPD basis, and ⟨S²⟩ = 1.20).

Many of these O_h Fe(II) compounds are, in fact, prototypical and well-studied spin-crossover complexes. Under appropriate external conditions, e.g., pressure or protein environment, spin-crossover complexes can exhibit transitions between LS and HS states, enabling the precise control of interesting magnetic phenomena.^99–102 Four of the Fe(II) complexes—[Fe(H₂O)₆]²⁺, [Fe(NH₃)₆]²⁺, [Fe(NCH)₆]²⁺, and [Fe(CO)₆]²⁺—have received significant attention from theoreticians utilizing an array of sophisticated ab initio methods (in the absence of gas-phase experimental measurements).¹⁰³ Spin gaps pose a difficult problem for DFT methods since a range of splittings can be obtained depending on the functional employed (and in particular, for global hybrids, the amount of EXX incorporated).^103,104

A notable study of these four Fe(II) complexes compared the results of a host of DFT functionals and wavefunction methods with those from diffusion Monte Carlo (DMC) calculations (employing pseudopotentials and within the fixed-node approximation).¹⁰⁵ It is claimed therein that the CO species is markedly multi-configurational. In another study, all-electron DMC calculations have been performed. For the complex with the CO ligand, the calculated spin gaps with single-determinant and multi-determinant trial wavefunctions agreed to within the 0.005 Ha error bar, suggesting that this complex is not strongly correlated.¹⁰⁴ Neese and co-workers did not consider the CO species but carried out a detailed investigation of the remaining three Fe(II) complexes with DLPNO-CCSD(T) methods employing large basis sets.¹⁰⁶ 

Now, we will add to the data in Table IV the finding that all HS states of the four Fe(II) complexes presently under consideration did not exhibit any SSB. This observation, found also for the diatomic species mentioned above, is quite general and provides the point of departure for the development of spin-flip approaches to, e.g., TDDFT and wavefunction methods (CI, CC, and more). This can be understood in the context of the quantum theory of angular momentum, where a high spin quantum number can have numerous, e.g., m_s states, at least one of which can typically be well-described by a single determinant. For the LS singlet species, going from CO, NCH, NH₃, to H₂O (i.e., going from strong-field/π-accepting toward weak-field/σ-donation-only), we find that the deviation of ⟨S²⟩_UHF increases consistently. UB3LYP completely restores the SSB found in UHF, except in the case of Fe(H₂$O)62+$⁠. κ-UOOMP2 yields similar conclusions. In disagreement with the aforementioned claims from other groups in the literature, our results strongly suggest that the CO species is not MR, given that no SSB occurs even at the UHF level. As detailed in the supplementary material, the NOONs from an 18e32o ASCI calculation are all either greater than 1.93 or less than 0.06, clearly indicating a closed-shell, SR singlet state. On the other hand, all dynamically correlated methods here suggest that the LS H₂O species has HS character mixed in. This SSB behavior is consistent with calculated 10Dq values from Ref. 96, reproduced in Table IV.

The adiabatic LS/HS gaps, i.e., derived from separately optimized LS and HS geometries, have been calculated in Refs. 106 and 96. Both CASPT2 and CCSD(T) methods predict that only Fe(⁠$CO)62+$ has a LS GS (due to the large 10Dq of 4.9 eV), while the GSs of the three other molecules with weaker-field ligands are HS. Evidently, the NCH ligand reduces 10Dq such that the cost of promoting two electrons to the e_g orbitals is more than compensated by the stabilization provided by exchange (with four spin-up unpaired electrons) and by unpairing two pairs that had been in the t_2g manifold. NH₃ and H₂O ligands continue the trend of reducing 10Dq, and of the three factors mentioned previously, the HS state drops lower in energy due to the increasing ease of promoting from t_2g to e_g.

Yet while the adiabatic energy difference is relevant for determining the proper GS multiplicity when calculating thermochemical properties, it is the vertical gap—i.e., the difference between the LS state and the HS state in the LS geometry—which is directly related to the SSB behavior: At a fixed geometry, when the HS state is the GS, SSB in the LS state (assuming it is a saddle point in the energy surface in orbital space) is to be expected from unrestricted variational methods, wherein constraining a molecule’s multiplicity is done by constraining the projection S_z (equal to N_α - N_β) rather than the value of S². As a consequence, such methods will include suitable HS contributions into the single-determinant reference in order to minimize its energy. However, the presence of HS contributions at the single determinant level for such species does not necessarily imply that the LS state has significant MR character. Indeed, with a large-enough LS–HS gap, it is quite possible that the lowest LS state is dominated by a single determinant, while being above the HS state in energy. UHF/UKS calculations for such LS states can still exhibit SSB due to inclusion of HS character, despite the actual state being fairly SR (as indicated by NOONs). A classic main-group example of this is CH₂, where the triplet state is the GS, while the lowest singlet state is predominantly closed-shell (HONO occupation of 1.89 and LUNO occupation of 0.09). UHF/UKS calculations on S_z = 0 CH₂ would generate SSB, as it is energetically favorable to contaminate the closed-shell singlet with the lower energy triplet. Several species in Table V provide additional examples, as described later.

Table V reveals that, at the LS geometry, UHF puts the HS state energetically below the LS state in all cases except for the carbonyl complex, which may explain why only the carbonyl complex was spin-pure at this level of theory (Table IV). UB3LYP and κ-UOOMP2 predictions of the relative spin-state ordering for the NCH complex are in agreement with the CCSD(T) reference, which suggests that UHF has overstabilized the HS state (a consequence of HF’s neglect of all correlation except for same-spin exchange stabilization) such that its predicted state ordering is incorrect. Going toward weak-field ligands along the spectrochemical series, UB3LYP predicts a LS-below-HS ordering for the CO, NCH, and NH₃ complexes, but the ordering switches for Fe(H₂$O)62+$⁠. This can explain the onset of SSB in the UB3LYP level of theory at this molecule. Interestingly, the spin-state ordering predicted by κ-UOOMP2 switches from a LS to HS GS at the NH₃ complex; however there is no SSB in Table IV; furthermore, we note that the vertical spin gap at this level of theory deviates by 3 eV from the CCSD(T) benchmark for the CO species—a clue that will become relevant in our later discussion of the possible inappropriateness of MP2 in describing interactions such as metal-carbonyl bonds. On the whole, it appears that with a HS GS, SSB in the LS excited state need not imply MR character, but rather what we will refer to as “variational collapse.” Indeed, this is the reason that some sort of spin-projection^57,107–109 is mandatory in such cases.

Consider Fe(H₂$O)62+$⁠, for which SSB persists not only at the UHF level but at all DFT and κ-UOOMP2 theories. A CASSCF calculation with six electrons in five metal d orbitals yields NOONs of 1.960, 1.954, 1.954, 0.066, and 0.066. As shown in the supplementary material, both the 12e14o active space selected following the protocol from Ref. 110 and an 18e20o calculation with the ASCI solver yielded NOONs either greater than 1.95 or less than 0.06. The implied SR character suggests that, in this case, SSB simply reflects the HS-below-LS relative energetics rather than the presence of static correlation. For Fe$X64−$⁠, with X = F, Cl, and Br proceeding toward weaker-field/π-donating ligands, the expected 10Dq decrease is a small effect, with LUNO populations resulting from minimal (e.g., 6e5o) CASSCF calculations of 0.076, 0.100, and 0.105, respectively. For Fe$F64−$⁠, we verified that using a larger active space of 6e15o to include a second d shell yielded a similar LUNO of 0.070 (vs 0.076 from 6e5o). The LS state of Fe(II)$Br64−$ is less MR than NiH, CoH, and FeH (with LUNO occupations of 0.14, 0.22, and 0.35, respectively), and its LUNO value is strikingly similar to that of the lowest singlet state of CH₂ (0.09), which is also predominately closed shell. This analysis suggests that neither the CO nor H₂O LS Fe(II) complexes are MR, in agreement with previously reported D1 diagnostic values of 0.14 and 0.06, which are below the 0.15 threshold suggested by Jiang and co-workers for transition metals.²⁹ Rather, the SSB observed from theories that include dynamic correlation is a manifestation of variational collapse. We thus reiterate that for excited states, SSB should be used together with NOONs from a MR theory in order to probe for the presence of MR character.

To summarize, this set of Fe(II) complexes suggests that strong-field ligands such as CO (and, e.g., CN) yield single-configurational, LS GSs (which is also the case for 4d and 5d metal complexes, for the same reason—i.e., due to the large 10Dq that results from the strong M–L interaction). Weaker-field ligands on the spectrochemical series (e.g., NH₃ and H₂O) favor HS GSs due to the attenuation of 10Dq and the stabilizing exchange interactions among same-spin electrons (indeed, as the 10Dq parameter of T_d compounds is roughly half of that of O_h compounds, the former, in general, have HS GSs). Practically speaking, for the calculation of thermochemical properties (for which only GSs are relevant) of these types of coordination complexes, SSB will therefore be the exception rather than the rule, and encountering MR character need not be much of a concern (though appropriate SR methods must still be chosen with care, vide infra).

We now seek to uncover chemical themes or circumstances, in addition to the diatomics analyzed above, in which static correlation can arise in the GS of transition metal complexes. Such conditions, while admittedly quite rare, are found in LS metal compounds, which have very weak M–L bonding (encountered, e.g., in the gas phase in particular for σ-donating only ligands with low metal oxidation states), or antiferromagnetically coupled spins occupying either separate metal centers or a metal center and a low-lying π^* orbital of a redox-noninnocent ligand.

In Ref. 33, the ligand dissociation energies of 34 metal complexes, formed via ligand coordination of neutral or cationic 3d metals in the gas phase, were investigated. We have computed ⟨S²⟩ with respect to UHF, UB3LYP, and κ-UOOMP2 orbitals for all GSs involved. Table VI shows the cases for which SSB at the UHF level was restored in both UB3LYP and κ-UOOMP2 cases; Table VII shows those cases for which SSB persisted at the UB3LYP and κ-UOOMP2 levels.

We observe that the molecules in the former group, for which UHF SSB is restored with UB3LYP and κ-UOOMP2, consist of metals with CO ligands, except in one case (which will be discussed after). CO ligands can participate in both σ-donation and π-backbonding. The latter substantially strengthens the metal–ligand bond, increasing 10Dq which decreases the chances of finding a truly strongly correlated state; hence, the spin symmetry restoration when dynamic correlation is included via DFT or κ-UOOMP2.

The molecules in the second group (Table VII), for which SSB persists with UB3LYP and κ-UOOMP2, all have weak-field ligands, which predominately participate in σ-donation, i.e., without π back-donation. This is obvious for H₂O and NH₃ ligands. The CH₂O ligand has a double bond between C and O atoms and coordinates at the O end, with sp² hybridization that is, in principle, capable of π interactions with the metal. We performed an EDA calculation, which shows that while the formaldehyde ligand can π-accept, the backbonding charge-transfer component is less than half that found for the CO complex (see Table S2). The weaker π-backbonding to CH₂O vs carbonyl appears to create a 10Dq comparable to σ-only ligands, which enables SSB.

We also note that this second group of complexes contains only Fe and Mn ions. Mn⁺ and Fe⁺ are unique, in that they have 4s¹ 3dⁿ⁻¹ electronic configurations. When a Lewis base donates two electrons, one must go into a bonding orbital, while the other must go into an antibonding orbital, which lowers the bond order and weakens the covalent bond more so than if the metal atom was 3dⁿ. As implied by bond strength trends with H₂ and CO ligands, Fe⁺ and Mn⁺ require 24 and 113 kJ/mol, respectively, to promote the 4s electron to attain the 3dⁿ configuration.³⁴ This trend is consistent with our finding that the SSB is more severe in the Mn vs the Fe ammonia species; indeed, the Mn molecule is the only species, which has SSB with UB3LYP orbitals in excess of 10%.

A parallel can be drawn between the chemical situation that occurs in these very weakly bound metal complexes and that which occurs while stretching bonds. Taking H₂ as an example, the gap between the bonding and antibonding MOs decreases as the bond length is stretched; after a certain distance, the near-degeneracy of the orbitals (or equivalently, of the many-body singlet and triplet states) results in an unpairing of the closed-shell singlet state into a two-configurational, biradicaloid singlet wavefunction. When the wavefunction is constrained to a single-determinant with unrestricted orbitals, beyond the CF point, the broken-symmetry determinant, $|↑↓〉$⁠, acquires partial triplet character since $|S2=2,ms=0〉=12(|↑↓〉+|↓↑〉)$⁠. In these transition metal complexes, we generalize the two involved multiplicities to LS and HS (quartet and sextet for the Fe complexes and quintet and septet for the Mn complex), and the accessibility of the HS state due to the factors mentioned above can be visualized as a vertical compression of the entire MO diagram.

For septet Mn(⁠$NH3)4+$⁠, the spin-up HOMO has predominately s character (as given by partial Mulliken populations from the B3LYP calculation). This suggests that, at least for this case, the fractional closed-shell unpairing corresponds to population of an orbital not in the t_2g or e_g manifolds, but rather one derived from the 4s metal AO. Indeed, the NOONs from a CASSCF(6e6o) calculation of the quintet have occupation numbers 1.06, 1, 1, 1, 1, and 0.94; the orbitals corresponding to the 1.06 and 0.94 NOONs have equal contributions from Mn s and d orbitals. The weak bonding apparently compresses the MOs such that the 4s orbital, which is in typical cases well-separated and energetically above the 3d MOs, is close enough to the 3d levels to stabilize the septet (via extra exchange). Indeed, CASSCF(6e6o) (and also κ-UOOMP2) predicts a septet GS with the quintet ∼5 mHa higher in energy, and thus, the quintet NOONs reflecting six unpaired electrons suggest that, rather than diagnosing an MR quintet GS, what we have witnessed instead is a variationally driven collapse toward the HS state. This underscores the danger of using small active spaces in CASSCF methods, namely, that the spin-state ordering can be in error without dynamic correlation, which can lead to a false diagnosis of MR character. However, even theories that do formally include dynamic correlation lead to inconclusive results: In the def2-TZVPP basis and with the DKH Hamiltonian, BP86 (pure) and B3LYP (global hybrid with 20% EXX) functionals put the quintet below the septet state, while PBE0 (global hybrid with 25% EXX) and HF reverse the ordering. These multiplicities are clearly very close in energy, and the possibility of a septet GS—at odds with conventional chemical intuition—is currently being re-investigated carefully with a number of ab initio methods.

Finally, we comment that in contrast to the SSB found in Mn and Fe complexes with ammonia ligands discussed above, in Co(⁠$NH3)3+$⁠, UHF SSB is artificial (i.e., restored via UB3LYP and κ-UOOMP2) presumably because the 3d⁸ configuration of Co, without s occupancy, enables relatively stronger bonding/ligand-field splitting even in the absence of π-backdonation.

While we were unable to perfectly correlate deviations in the calculated ligand-dissociation energies using SR methods such as DLPNO-CCSD(T) or DFT vs experiments with SSB at the UB3LYP or κ-UOOMP2 levels of theory, we note that the largest error from DLPNO-CCSD(T) vs experiment among the 34 molecules investigated in Ref. 33 was for the Fe(⁠$NH3)4+$ species (9.15 kcal/mol). For the Mn(⁠$NH3)4+$ complex, the B3LYP, B97, M06, PBE0, and ωB97X-V functionals, along with DLPNO-CCSD(T), consistently underestimated the experimentally measured ligand-dissociation energy in the range of 4.74–6.88 kcal/mol (though the use of the septet multiplicity is being explored). We must point out, here, that such weakly bound ions are unlikely to be found in solution. The binding energies are very weak, e.g., 10.0 ± 1.7 and 8.6 ± 1.4 kcal/mol for Fe and Mn tetra-ammonium species, respectively, and in polar solvents such as water, these complexes may readily dissociate. It is also likely that the Mn(I) and Fe(I) states will undergo redox events with the solvent to become the more stable Mn(II) and Fe(II) oxidation states.

The ab initio prediction of electrochemical redox potentials of transition metal catalysts is a long-sought goal. Encouragingly, previous analyses of the deviations of DFT predictions from experimental measurements have revealed some systematic trends,¹¹¹ and there are indications that computing the potential of the reference redox couple can encourage favorable error cancellation.¹¹² However, one can easily find a sizable number of large outliers, and as the presence of MR character in the highly reduced species is one factor that would lead to erroneous predictions, the ability to quickly diagnose these situations is a prerequisite if computational predictions are to be reliable and predictive.

When an electrode reduces a homogeneous metal catalyst one or more times, making it redox-active toward substrates such as CO₂ and O₂, the reduction is most often metal-centered.^113–117 Indeed, a unique property of transition metals is that they can accommodate multiple oxidation states, which can enable remarkable reactivity. However, large delocalized ligands can have low-lying π^* orbitals, which can be preferentially reduced vs a virtual metal d orbital.^118–127 Such reduced species can have a LS GS due to the antiferromagnetic (AFM) interaction which can arise between opposite spins localized separately in metal and ligand orbitals. When this AFM stabilization outweighs the potential exchange stabilization if all unpaired spins were oriented in the same direction, the GS will be LS with a wavefunction requiring more than one determinant.^{51,122,125,126,128}

From the perspective of SSB, we first investigate a set of metal complexes that was the focus of a recent work by Batista et al., in which redox potentials in non-aqueous solution were computed and compared to experimental measurements.⁷ This set included MCp₂ and $MCp2*$ where M = Fe, Co, Ni; and M′bpy₃ where M′ = Fe, Co. Cp, Cp*, and bpy denote cyclopentadienyl, pentamethylCp, and bipyridine ligands, respectively.¹²⁹ As shown in Table VIII, SSB in the UHF wavefunction occurs in nearly all species (with the exception of $NiCp2*$⁠). For the Cp and Cp* complexes, UB3LYP and κ-UOOMP2 restore spin-purity in all cases, implying SR character.

For the tri-bpy complexes, the notably large UHF SSB is approximately restored in all cases with UB3LYP and κ-UOOMP2. Looking more closely at Co(⁠$bpy)32+$⁠, for which ⟨S²⟩_κUOOMP2 calculated from the Slater determinant of optimized orbitals deviates by 6.6% from the exact value, we find that this deviation increases to 10.4% when evaluated with respect to the first-order wavefunction associated with κ-UOOMP2 [i.e., using Eq. (6)]. To investigate the possibility that DE in the B3LYP functional has led to a bias toward spin symmetry restoration for this Co(II) species, as seen previously for the stretched hydrogen-fluoride molecule in Fig. 3, we find ⟨S²⟩_CAM-B3LYP = 0.83, identical to the value from the κ-UOOMP2 wavefunction. Figure 7 also shows that SSB can be modulated by the amount of EXX included in global hybrid functionals, revealing that PBE with 20% EXX yields SSB comparable to CAM-B3LYP, which appears sensible given that CAM-B3LYP has 19% short-range EXX (and 65% long-range EXX). Increasing %EXX still further increases ⟨S²⟩ toward the HF value. Recall that this arbitrariness is not unique to hybrid DFT functionals as the tendency toward SSB in κ-UOOMP2 is dependent on the value of the κ regularizer: Scanning κ gives a full view of the symmetry breaking landscape.⁵⁶ 

Admittedly, it is difficult to draw conclusions regarding the MR character of the Co(⁠$bpy)32+$ species. On the one hand, x-ray absorption spectra have been well-reproduced by linear combinations of simulated spectra from the LS (doublet) and HS (quartet) states, leading to the claim that Co(⁠$bpy)32+$ is 57% HS and 43% LS, whereas Co(⁠$bpy)33+$ is relatively more monoconfigurational (∼ 80% LS).¹³⁰ Taking these distributions and assuming the spectra were taken at room temperature imply spin gaps of roughly 0.01 and 0.03 eV for the Co(II) and Co(III) species, respectively, which indeed would imply SC. On the other hand, the bpy ligand is between NH₃ and CO on the spectrochemical series suggests that 10Dq should be relatively large due to π-backbonding, and therefore, SSB is likely to be restored with a suitable level of dynamic correlation (and indeed, it is with B3LYP). It is also true that UB3LYP was used in Ref. 7 and the calculated Co(III) → Co(II) reduction potential matched the experimental value to a very high degree of accuracy (within a tenth of an eV, when triple-ζ basis sets were used). Thus, negligible SSB from B3LYP and worst case deviations of around 10% from CAM-B3LYP and the κ-UOOMP1 wavefunction suggest that this complex is predominately of SR character with redox-innocent bpy ligands.

We now turn to Fe complexes with terpyridine (tpy) or porphyrin ligands, which due to the more delocalized ligand frameworks are expected to have a relatively lower-lying π^* orbital than bpy. Both complexes are efficient electrocatalysts for the CO₂ reduction reaction (CO2RR) to CO. An important feature of CO2RR electrocatalysts is the substrate selectivity of CO₂ against the hydrogen evolution reaction because proton-coupled CO₂ reduction and direct proton reduction occur at similar potentials. The incorporation of a redox-active ligand, which when reduced yields a relatively Lewis-acidic metal center, results in the metal favoring CO₂ binding over protonation (the formation of a metal hydride is the first step in the hydrogen evolution mechanism). This origin of product-selectivity has been established for the Re(bpy)(CO)₃Cl complex^123,131 and also for iron polypyridine^125,126,132 and porphyrin¹²² catalysts. Furthermore, the placement of the extra electron in a π^* orbital that is lower-lying than a virtual d orbital results in milder reduction potentials—another desirable feature of an electrocatalyst—that are made milder still due to the additional stabilization from the M–L AFM coupling.

Combined spectroscopic and computational work conducted by Neese and co-workers¹²² suggests that the active species of the popular FeTPP catalyst (TPP = tetraphenylporphyrin) is an intermediate-spin Fe(II) center that is antiferromagnetically coupled to a porphyrin diradical anion. The FeTPP catalyst and its derivatives are among the most active CO₂ reduction catalysts.^118,121,133 Recent work by Derrick et al. shows that a tpy-based ligand framework (tpyPY2Me) in combination with an iron center is an efficient CO₂ reduction catalyst at a low overpotential.¹²⁶ The doubly reduced active form of the catalyst is a singlet with a doubly reduced tpy ligand strongly coupled to an intermediate spin Fe(II) center. This electronic structure was established based on both computational and spectroscopical evidence. Both complexes and their reduction reactions are depicted in Fig. 8.

Table IX shows the SSB behavior for the Fe(II)-tpy species with net charge n = 2+, 1+, and 0. The first reduction is known to occupy the noninnocent tpy π^* orbital, leaving the d⁶ Fe(II) center closed-shell and thus forming an overall doublet. Consistent with this picture, both n = 2+ and singly reduced complexes, while severely spin-contaminated at the UHF level, do not exhibit significant SSB when accounting for dynamic correlation via B3LYP, CAM-B3LYP, and B5050LYP functionals. In contrast, the second reduction is accompanied by the loss of a ligand and a spin transition of the iron center from LS to intermediate spin (i.e., S_Fe = 1), while the twice-reduced tpy (S_tpy = 1) couples to the metal center to form an overall singlet. This AFM coupling of opposite spins separately localized on metal and ligand requires a MR description (for fundamentally the same reason that open-shell singlet biradicals require two determinants) and manifests as SSB, which persists upon inclusion of dynamic correlation with the three DFT functionals investigated. As expected, the calculated ⟨S²⟩ value increases with %EXX in the hybrid functional.

In the case of the iron porphyrin, here modeled without the phenyl groups (denoted FeP), the neutral complex has a triplet GS.¹³⁴ As indicated by the complete restoration of the UHF SSB with all DFT functionals and κ-UOOMP2, shown in Table X, this species is predicted to have a SR electronic structure. Both the first and second reductions are ligand-centered and result in M–L AFM coupling, which has been observed experimentally. The SSB that persists from HF through all DFT methods (increasing, as in the Fe tpy systems above, with %EXX), for both FeP⁻ and FeP²⁻ corroborate the presence of M–L AFM coupling involving the S_Fe = 1 center and the reduced non-innocent porphyrin.

As first encountered in Ref. 58, κ-UOOMP2 predicts unphysical GSs for the FeP species. In particular, while the neutral species exhibits very minor spin-contamination, the Mulliken spin density on the Fe atom is found to be 1, as opposed to the expected value of 2 for an Fe-centered triplet. In addition, we found that ⟨S²⟩_κUOOMP2 for FeP¹⁻ is 0.87, roughly twice as small as the values found with the three hybrid DFT functionals used, and the spin-density on the Fe atom is 0.4 (vs the expected value of 2 for a ligand-centered reduction stabilized by AFM coupling to the metal; UB3LYP gives 1.96). Similarly, for the twice-reduced Fe tpy and FeP species, κ-UOOMP2 predicts spin-pure, closed-shell GSs with a net spin of 0 on the metal, inconsistent with the AFM-coupled states deduced from experimental measurements and predicted from B3LYP, CAM-B3LYP, and B5050LYP functionals.

Finally, we remark that MR character arises also when catalysts such as metalloporphyrins bind small-molecule substrates such as O₂ and NO.^135–138 These are of great biological relevance and deserve further investigation alongside catalytic systems for CO2RR and other desirable reactions. The present results show that SSB from these hybrid DFT functionals is a computationally inexpensive hallmark, which can be used in future in silico catalyst design projects and can inform, if not replace, chemical intuition regarding the possible non-innocence of novel ligand frameworks.

The quintessential example of AFM coupling is multi-metal compounds.^139–141 As a simple example system, we consider the four smallest Mn(III)–Mn(IV) compounds from Ref. 142, which have two bridging oxygens directly connecting the Mn centers. These molecules are shown in Fig. 9. The ⟨S²⟩_UB3LYP values for compounds numbered 5, 6, 7, and 11 are all 3.8, which is far from the exact value of 0.75 for these experimentally assigned doublet species. We confirmed that the spin-densities on the Mn atoms are 4 and −3 in all cases (we note that a HS, octet calculation needed to be performed first, and the resulting density was used to initialize a subsequent calculation of the doublet state), indicating strong AFM coupling. As expected, for these types of multi-metal states, MR methods with large active spaces are required to obtain exchange-coupling parameters that agree with experiment.⁴⁴ Rather unexpectedly, in many cases, broken-symmetry DFT (i.e., utilizing some form of spin projection) appears capable of producing accurate coupling constants as well, though relative spin-state energetics are nevertheless sensitive to the functional employed.^143–146 

In this section, we emphasize that even in the absence of persistent SSB (i.e., for species well-described by a single determinant), computational methods should be chosen with care. In other words, in our view, categorizing a system as either MR or not is only the first step to quantitatively accurate predictions. As an illustrative example, we focus on metal complexes with strong-field CO ligands, which are of primary importance in organometallic chemistry. As discussed above, the MO diagrams of these complexes are characterized by large 10Dq values, which result in spin-pure, LS GSs. However, Fig. 10 shows that DHDFs, which incorporate MP2 correlation energies, consistently overestimate experimental ligand dissociation energies^34,147 for M–CO complexes. While one might expect that the highest rung of Jacob’s ladder^4,148 would provide the most accurate results for these organometallic compounds, DHDFs drastically underperform simpler functional forms, e.g., the global hybrids B3LYP and B97 and the range-separated hybrid ωB97X-V.

Next, we consider the 3d metal complexes from Ref. 33 with gas-phase experimental ligand-dissociation measurements and select only those for which DLPNO-CCDS(T)/CBS yields accurate results. The six suitable compounds are V$(H2)4+$⁠, Co$(H2)4+$⁠, Ti(H₂$O)4+$⁠, Cu(⁠$NH3)4+$⁠, Cu(⁠$CO)4+$⁠, and Fe(⁠$N2)4+$⁠. H₂O is a π-donor, and NH₃ ligands can only engage in σ-bonding with the metal. H₂, N₂, and CO ligands are π-acceptors, in the order of increasing backbonding strength (see Table S2). Figure 11 shows, in the cc-pVTZ basis, the deviations of the ligand-dissociation energies calculated with UMP2, κ-UOOMP2, and UCCSD from UCCSD(T). Evidently, while MP2 methods yield reasonable accuracy for transition metal chemistry (∼2–3 kcal/mol³⁶) for the complexes with ligands that can only σ-donate or weakly π-accept, large errors are found for Cu(⁠$CO)4+$ (interestingly orbital optimization makes the accuracy worse). These results imply that the overestimation of ligand-dissociation energies with DHDFs vs experiment (Fig. 10) stems from an overestimation at the MP2 level. Neese et al. also found that MP2 and OOMP2 produce large (4.8–60 kcal/mol) errors vs CCSD(T) for CO dissociation of four Cr and Ni complexes;¹⁵⁰ similarly, Hyla-Kryspin and Grimme found for MP2-calculated CO-dissociation enthalpies an average overestimation of 19.3 kcal/mol over a set of seven 3d metal carbonyl species (interestingly, while MP3 led to underestimations of even larger average magnitude, the performance of spin-component-scaled MP2 and MP3 was more promising, with MAEs of 9.6 and 3.7 kcal/mol, respectively).¹⁵¹ 

We hypothesize that while the pairwise additive approximation for the correlation energy of MP2 theory can adequately describe σ-bond only ligands (in the absence of static correlation), it cannot be expected to produce quantitative results for four- and six-electron-like interactions involving both σ-donation and π-backbonding, regardless of the orbital set employed. σ-donation involves one electron pair, while π-backbonding can additionally involve either one or two more electron pairs, depending on the number of available donor–acceptor orbitals. In the case of a metal-carbonyl bond, π-backbonding can involve two pairs of electrons backdonated into the degenerate $πx*$ and $πy*$ orbitals of CO, resulting in a bond that effectively involves six correlated electrons in total.

This reasoning begs the following question: For metal carbonyl complexes, or more generally for metal complexes with strong-field ligands which can significantly π-accept, will the SSB behavior of κ-UOOMP2 be compromised due to its inability to describe the inter-pair correlations involved in dative bonds characterized by simultaneous σ-donation and π-backdonation? If so, how can SSB due to this be distinguished from SSB due to genuine MR character? The first point is that these inter-pair effects are, in general, smaller than single pair correlation energies. Moreover, when inter-pair correlations are physically significant (as we argue to be the case in organometallic compounds), the gap between LS and HS states is typically already large, as in the Fe(II)(⁠$CO)62+$ system discussed earlier, rendering spin-state mixing an irrelevant concern in most situations. However, for smaller-gap systems such as complexes with redox-noninnocent ligands, e.g., twice-reduced Fe tpy and porphyrin, we have shown above that κ-UOOMP2 yields qualitatively incorrect predictions regarding the presence of MR character, and thus, the use of SSB from DFT orbitals is recommended. Work is ongoing to disentangle errors due to κ-regularization and the MP2 ansatz itself.

While transition metal chemistry has a reputation for being challenging for computational quantum chemistry models, one has to look hard for static correlation in the GSs of mono-metal transition metal complexes. This is because for small 10Dq values, the HS state (which is, in general, SR) is favored due to the exchange stabilization of same-spin electrons and the energy “saving” from not having to pair electrons; for large 10Dq values, the LS is energetically far-below the HS state, which prevents the latter from contaminating the total spin of the former, leading again to a SR state. The latter situation occurs when 4d and 5d transition metals are involved. One of the goals of this work, though, is to find relevant situations in which static correlation is present in molecules. We have encountered MR character in the following:

1. Metal hydride diatomics courtesy of the inapplicability of the Jahn–Teller theorem to linear systems, in which spatial symmetry mandates two important determinants for both TiH and CoH. For less electronegative 3d metals, i.e., going from Ni to Co to, mostly strikingly, Fe, we find increasingly non-negligible LUNOs due to the narrowing gap between the non-bonding metal 4s-3d shell and the antibonding orbital (involving metal $dz2$ and hydrogen 1s), which opens the door to SSB due to the intrusion of higher-spin-state character.

2. Metal complexes with higher coordination numbers and low [i.e., Fe(I) or Mn(I)] oxidation states, which exhibit very weak bond energies (tens of kJ/mol) due to weak-field (σ-donation only) ligands and metals with singly occupied 4s orbitals preventing favorable dative bonding. While less likely to exist in solution, such weak bonding can be encountered in the gas phase (in, e.g., metal-organic frameworks or atmospheric chemistry).

3. Molecules that exhibit AFM coupling resulting in LS GSs. We have demonstrated that this can occur in reduced states of homogeneous catalysts containing redox-noninnocent ligands such as terpyridine and porphyrin and also in oxygen-bridged bimetallic Mn(III)–Mn(IV) dimers.

For para-benzyne along with stretched H₂ and hydrogen fluoride, the SSB behavior and NOON predictions from DFT orbitals agree quite well with those of ASCI (a near-exact FCI approximation). As expected, increasing the fraction of EXX in hybrid functionals pushed CF points to shorter bond lengths and increased the amount of artificial SSB (i.e., LUNOs in excess of the exact value) approaching that of UHF (i.e., 100% EXX). It is remarkable that computing ⟨S²⟩ from the KS orbitals (which represent the solution of a non-interacting problem) yields results similar to those obtained from approximate (e.g., κ-UOOMP2) wavefunction methods and NOONs in close agreement with those from ASCI. However, as the case of stretched hydrogen fluoride reveals, when pure and common hybrid DFT orbitals are used to compute ⟨S²⟩ and NOONs, the static correlation error is intertwined with DE, which leads to erroneous results. In such cases, κ-UOOMP2, which is free of DE, is formally more reliable (in addition to the fact that ⟨S²⟩ and NOONs are well-defined). We note that range-separated hybrids and functionals with a high percentage of EXX are lower-scaling alternatives in which DE is much-reduced, and the use of orbitals from optimally tuned range-separated hybrid functionals may be a potentially interesting way forward in certain cases.

Regarding the diagnosis of static correlation, SSB at the UB3LYP and κ-UOOMP2 levels is shown to be insufficient in case 1 from above as it does not account for a wavefunction that must be multi-configurational to respect spatial symmetry, and generally in LS states which are not the GS due to variational collapse (practical implications of the latter for the calculation of HS–LS spin splittings will be discussed below). In these special cases, SSB should be used in conjunction with NOONs from a MR method, ideally one with an appropriate treatment of dynamic correlation. Methods that use a HS reference and apply spin-flip excitations to generate the Hilbert space of lower spin multiplicities appear to be good candidates to, e.g., probe the possible multi-configurational nature of a LS excited state and have the advantage of a black-box selection of the active space orbitals involved in spin-unpairing from the target LS state.

Another important takeaway suggested by this work is that strong correlation is a term that demands clarification, especially in the context of transition metals. The term is frequently used as a synonym of static correlation and MR character, and we have presented SSB as an intuitive and meaningful diagnostic for the GS when spatial symmetry is not an issue. Having established a connection with chemically revealing models such as ligand-field and molecular orbital theories, a picture is painted which suggests that, fundamentally, static correlation in transition metal systems involves the same phenomenon as in organic molecules, e.g., biradicaloids.^56,152

In contrast, evidently, it is the dynamic correlation that presents additional difficulties compared to the case of typical organic molecules. For organometallics, four- and six-electron-like interactions (or at least inter-electron pair correlations) as relevant to π-donation or π-backdonation on top of σ-donation appear to be important, as illustrated by the failure of MP2 and MP2-based DHDFs in predicting experimental ligand dissociation energies for metal carbonyl complexes. The second idea that supports this conclusion, while admittedly less concrete, is that large errors in DLPNO-CCSD(T) calculations were found when the UHF solution exhibited significant spin-contamination. It is our expectation that when orbitals obtained from theories that include (even approximately) dynamic correlation, such as DFT and κ-UOOMP2, are employed, subsequent CC predictions will improve in accuracy, as has already been seen in main group molecules.^153,154 Thus, in our view, the common notion that transition metals are difficult for traditional electronic structure methods is due to both static correlation and dynamic correlation, yet we stress our finding that truly MR situations are encountered only in the special cases enumerated above (and perhaps a few others), making the proper treatment of dynamic correlation more important in most commonly encountered cases.

There is also an important connection between a theoretical method’s ability to capture dynamic correlation and the ability of SSB to imply MR character, which depends on the accurate prediction of relative energies between two (or possibly more) spin states. In the limit of no correlation other than the exchange interaction between same-spin electrons (required by Fermi statistics/Pauli-exclusion), as is the case in HF theory, HS states are artificially favored relative to LS states. For a system with a LS GS (determined experimentally or by an exact theoretical method), HF (and often CASSCF with small active spaces) will significantly underestimate the LS–HS gap or even incorrectly predict a HS GS. In such cases, the SSB implied by the HF LS wavefunction will be artificially large due to the unphysical inclusion of the HS state. The idea behind using orbitals optimized from Kohn–Sham DFT or κ-UOOMP2 is essentially an attempt to more accurately describe relative spin-state energetics via the (approximate) inclusion of dynamic correlation, which enables true MR character to be reliably diagnosed.

Yet the difficulty involved in accurately predicting LS–HS gaps (in a reasonable amount of time) cannot be overestimated. With regard to DFT, the issue stems from the systematic overstabilization of HS states as the fraction of EXX is increased. Pure DFT functionals, e.g., PBE, BLYP, or BP86, are justifiably advantageous when investigating transition metal systems in the sense that the resulting orbitals are likely to be free of “artificial” SSB; subsequent PT or CC will converge more quickly, and there is no arbitrariness in, e.g., how much EXX to include. However, pure functionals are plagued by high levels of DE and, as a result, are known to yield unphysical charge and spin densities, underestimate barrier heights, and so on. With regard to κ-UOOMP2, we recall (i) the large 3 eV error in the vertical spin gap of the iron hexacarbonyl species [vs CCSD(T) in the same basis] and (ii) the spurious closed-shell GSs implied for the twice-reduced Fe tpy and porphyrin species, which are known experimentally to be open-shell singlet states with AFM coupling between the metal and the ligand. Both (i) and (ii) involve the overstabilization of closed-shell LS states in complexes with strong-field ligands, which form bonds to metal ions involving simultaneous σ-donation and π-backdonation. Evidently, for these types of organometallic complexes, the fact that dynamic correlations among multiple electron pairs are entirely missing in MP2 approaches (and thus to some extent in DHDFs incorporating MP2 correlation) has a detrimental impact on the accuracy of predictions regarding thermochemistry (Fig. 10), relative spin-state energetics, and the presence of static correlation via SSB.

Further development and assessment of SR methods, in addition to local correlation functionals, capable of describing this type of dynamic correlation relevant to transition metal bonding is clearly necessary. CCSD and CCSD(T) are limited by their high scaling. The direct variant of the Random Phase Approximation (dRPA),¹⁵⁵ which can be implemented to scale as the fourth power of system size,¹⁵⁶ has been shown to be equivalent to ring-CCD¹⁵⁷ and therefore includes terms approximately correlating excitations involving more than one electron pair. dRPA has shown promising accuracy for bond energies of metal carbonyl complexes (among others).¹⁵⁸ On the other hand, the improved description of the dissociation limit within the RPA formalism comes at the expense of DE.¹⁵⁹ OORPA approaches have been developed,¹⁶⁰ which do not seem to require regularization, and may be promising in the context of a DHDF.^161,162 OOCCD and MP3 also approximately include inter-pair excitations, and the scaling of these methods might be reduced by localized approximations or tensor decompositions of the two-electron integrals, enabling their use in DHDFs. Our data suggest, though it remains to be seen, that the performance of DHDFs incorporating such correlation energies^163,164 may yield improved accuracy for organometallic thermochemistry.

The implications of our findings for practical quantum-chemical calculations on TM complexes and catalysis modeling are broadly as follows:

• The MR character in the GS, while relatively rare, can be diagnosed in a black-box and active-space-free manner through SSB in hybrid DFT calculations (even employing small basis sets). We suggest B3LYP or B5050LYP and range-separated hybrid functionals if substantial DE is detected. Such SSB is typically found to coincide with fractional NOONs from multi-configurational wavefunctions, which together signal that the relative energetics of frontier orbitals are such that multiple spin states are nearly degenerate.

• In the absence of MR character, we can state the following:

• (i)

  When the LS state of interest is the GS, then SR wavefunction methods (with an appropriate treatment of dynamic correlation) and DFT are, in principle, capable of providing robust predictions. While the use of hybrid functionals is necessary to diagnose MR character, the choice of functional used to model, e.g., catalytic cycles should be based on a judicious mixture of system-specific criteria (e.g., agreement with available experimental measurements) and general benchmark performance across TM compounds.^{30,33,165–167} The performance of hybrids is typically superior to that of semi-local functionals. For organometallic complexes with π donating or accepting ligands, such as metal carbonyls, the use of DHDFs based on MP2 correlation energies may not yield improved results beyond hybrids such as B97 or ωB97X-V, as demonstrated here.

• (ii)

  When the LS state of interest is not the GS, HF and hybrid functionals can still exhibit artificial SSB. Care must be taken when evaluating spin gaps and related properties, and we recommend the use of either spin-restricted orbitals or spin-projection of spin-unrestricted solutions. Such approaches may also reduce (though not completely eliminate) the sensitivity of DFT-predicted spin-splittings on the choice of functional.

• c.

  In the presence of MR character, methods based on a single-determinant reference (e.g., MPn and CC theories) should generally be avoided. DFT, especially when using hybrid functionals with a high fraction of EXX, also appears to be unsuitable (though it is SR only in the context of the fictitious KS system). One could attempt to proceed (semi-quantitatively) using broken-symmetry DFT, where spin-projection can help if there is one leading contaminant¹⁰⁷ or if Heisenberg physics is operative.^109,168 When possible (in light of relatively high computational costs), electronic structure approaches that explicitly account for a multi-configurational wavefunction should be used.

See the supplementary material for NOONs and convergence details for ASCI-SCF calculations of larger molecules (para-benzyne, [Fe(CO)₆]²⁺, and [Fe(H₂O)₆]²⁺), NOONs and active space selection for a selected CASSCF calculation ([Fe(H₂O)₆]²⁺), and complete ALMO-EDA results.

J.S. thanks Romit Chakraborty for insightful discussions. This work was supported by the Director, Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

The data that support the findings of this study are available within the article and its supplementary material.