Can single reference density functional theory methods describe spin state crossing for 3d transition metal carbon monoxide association?

Abstract

The bonding between a transition metal (TM) and carbon monoxide is found in many chemical complexes that are important in biological and catalytic processes. Density functional theory (DFT) methods have been widely used for many calculations; however, the validity of such calculations for complex spin state coupling induced by TM–CO bond formation remains unclear. Interestingly, binding between a 3d TM atom and carbon monoxide induces a spin state change, making it a good benchmark system for evaluating various functionals. There have been numerous studies evaluating the TM atom spin state splitting or spin state energy difference of CO-bound TM complexes. However, we did not find benchmark studies evaluating the potential energy surfaces of the TM–CO association for various spin states to determine the spin-state crossing point. In the present study, we calculated the 3d TM + CO association potential energy curve using various DFT functionals. For TM = Sc, Ti, Fe, Co, and Ni, we performed multireference calculations with multiple excited states to clarify the spin state crossing points. We found that hybrid functionals can give spin state crossing points within 0.15 Å of those obtained by multireference methods, and confirmed that accurate atomic spin splitting results in more accurate crossing point geometry. In addition, we found minimal basis set dependence for the association potential energy curve for B3LYP. To our surprise, this study showed that hybrid DFT functionals can describe spin crossing phenomena in the TM + CO association.

---

1 Introduction

The spin state of transition metals (TMs) can affect the selectivity in electrocatalytic O₂ and CO₂ reduction reactions.^1,2 TM carbonyl complexes have attracted great interest from chemists due to their importance in CO poisoning of heme proteins in the human body. In the 1930s, Pauling and Coryell observed a spin state change, or spin crossover, upon CO binding to the heme-related iron–porphyrin complex.³ Recent studies have also evaluated the importance of spin crossover for the transport of CO and O₂ on iron–porphyrin complexes.^4,5 These spin crossover phenomena play a vital role when molecules bind to a TM complex with many low-lying electronic states; however, the effect of spin on reactivity is still elusive.⁶ To evaluate such spin crossover reactions, the two-state reaction model^7,8 has been utilized. In this model, one evaluates the spin crossover probability from the spin–orbit coupling (SOC) interaction energy at the geometry where the potential energy surfaces of the two spin states have the same energy and therefore cross. This idea of the minimum energy crossing point (MECP) greatly reduces computational effort, since one can focus on performing the tedious SOC calculation only near the MECP.^9–11 Studies have noted issues with defining the crossing point when the SOC is not included in the electronic energy.^12–14 However, the MECP concept has been used to evaluate various spin crossover reactions.^15–17 Thus, to determine the catalytic activity of TM complexes for various reactions, obtaining accurate MECP geometries becomes critical. Since it is difficult to obtain such properties experimentally, density functional theory (DFT) and quantum chemistry calculations have been used. There have been many studies evaluating the TM atom spin state splitting energies and the binding energies between TM complexes and simple molecules.^18–24 However, to the best of our knowledge, no systematic benchmark studies have compared DFT predicted spin-state-crossing geometries with those obtained from multireference methods.

Due to the complex nature of the problem, selecting the test system becomes critical for this systematic study. To understand the bonding nature of CO chemisorption, chemists have focused on the electronic interaction between a single TM atom and CO.^25–27 In particular, 3d TM–CO^28–38 adduct complexes have been experimentally synthesized and studied using infrared CO stretching vibrational spectra. Isolated 3d TM atoms (except for Cr and Cu) favor the electronic configuration of [Ar]3dⁿ4s², and the early TM atoms such as Sc, Ti, and V favor low spin states while late TM atoms such as Mn, Fe, Co, and Ni favor high spin states. However, for TM = Mn, Fe, Co, and Ni, the linear TM–CO complexes take one lower (next-highest) spin multiplicity, so the binding of CO inverts the spin multiplicities.^39–45 To physically understand the electronic structure nature behind the spin state change, we have to understand how a TM–CO bond is formed. The “Dewar–Chatt–Duncanson model”⁴⁶ relates the bonding to σ-donation and π-back donation between the adsorbate and the TM atom. Blyholder also introduced a similar concept for metal–inorganic bonds.⁴⁷ Bauschlicher et al.^48,49 showed that while the TM atom in the ground 3dⁿ4s² configuration repels CO due to the strong repulsive interaction between the CO-σ highest occupied molecular orbital and TM-4s electrons, the excited 3dⁿ⁺¹4s¹ configurations form TM–CO bonds by reducing the σ-repulsion. This is schematically shown with orbital energies in Fig. 1 for Sc + CO. Therefore, in addition to evaluating the spin state crossing geometry for the TM + CO association reaction, we should also assess the atomic spin state splitting arising from different TM 4s occupations. In conclusion, the association reaction of a TM atom with a CO molecule provides a simple system for systematically benchmarking quantum chemistry methods and evaluating their ability to describe energy crossings between spin states during TM–CO adduct formation. For this TM atom carbon monoxide system, multireference calculations involving multiple excited states are feasible and will provide a good benchmark for evaluating the validity of different single reference methods.^50–52

Fig. 1 Schematic interaction energy diagram for the Sc + CO reaction for the quartet spin state. Orbital energies, given in eV, are calculated using the B3LYP/aug-cc-pVTZ method. The σ, π, and δ interactions are given by dotted, dashed, and dash-dotted lines, respectively.

Although there are no gas phase experimental studies conducted for the TM + CO reaction, Andrews and coworkers have evaluated the formation of various TM(CO)_n complexes using matrix isolation techniques with laser ablated metal atoms in the presence of CO molecules.^53–60 They performed a systematic study of various 3d TM atoms and concluded that for Sc, Ti, and V, the ground spin state atoms will not react with CO to produce mono carbon monoxide adducts.^53,54 For these systems, their DFT calculations showed that the TMCO adducts and TM atoms have different spin states, and they concluded that excited spin state TM atoms, which are produced during the laser ablation process, likely produce these adducts. For the Co + CO reaction, Tremblay et al. considered a spin crossover reaction induced by photoexcitation and postulated a crossing between the doublet and quartet states during the Co–CO association.³² However, there are no detailed calculations on the crossing point for this reaction. Lastly, for Ni + CO, Andrews and coworkers did not discuss whether the triplet ground state Ni atom can react with CO to produce singlet NiCO, but they were able to observe the NiCO adduct in their matrix isolation experiments.⁵⁵ So there are unanswered questions regarding the MECP and the possibility of spin crossover for the reaction between ground spin state 3d TM atoms and carbon monoxide. Although the binding of a single carbon monoxide molecule and TM atom is not reported, there are various studies concerning carbonyl complexes TM(CO)_n, n = 4–6. Photodissociation studies have highlighted the importance of spin state changes.^9,61,62 For example, a recent study on the photodissociation of Fe(CO)₅ proved that the sequential dissociation of CO slows down at n = 4 due to the singlet–triplet spin state change for Fe(CO)₄ and TM(CO)₃.⁶³

Concerning theoretical calculations on TM + CO association, Head-Gordon and coworkers evaluated the collinear adsorption potential energy curve (PEC) for several 3d, 4d, and 5d TM atoms. Although they focused on electronic structures with a singlet ground electronic state,⁶⁴ their DFT calculation showed a spin state crossing for Ni + CO. Recently, we evaluated the Ni–CO association using a multireference spin–orbit calculation,¹⁴ and showed that the spin crossover Ni–C bond lengths determined using multireference methods and Becke's 3-parameter hybrid DFT functional, B3LYP,^65,66 were the same. In addition, we also showed that B3LYP can give a smooth PEC for 3d TM + CO association.^67,68 Extending these previous studies, we have evaluated the atomic spin state splitting for 3d TM using various DFT functionals and wave function methods. In addition, we calculated the collinear 3d TM + CO association PEC using various quantum chemistry methods to compare the spin state crossing geometries. Since DFT simulations identified spin state crossings for Sc, Ti, Fe, Co, and Ni, we also report multistate multireference potential energy calculations for the CO reaction with these TM atoms. The main goal is to systematically benchmark if single reference DFT methods can be used for spin state crossing for the 3d TM + CO association.

In addition, the CO binding on the TM bulk surface has been studied by DFT methods to evaluate the most stable binding sites: on-top, bridging, or hollow sites.^69–77 In these bulk calculations, spin-polarized DFT is utilized, and the spin state can change during the calculation upon CO adsorption. These studies have shown that modulating the gap between the CO highest occupied molecular orbital and lowest unoccupied molecular orbital energy can control the CO binding energy on different TM surfaces.^71–74,77 The present study of collinear TM + CO association will also provide a better understanding of the validity of such approximations for on-top binding site interactions. To clarify the trends in DFT functionals, we also evaluated the TM atomic orbital energies to elucidate the trend in TMCO binding energies and spin state crossing points.

2 Computational methods

We evaluated various DFT functionals in this study. The generalized gradient approximation functional by Perdew, Burke, and Ernzerhof, PBE⁷⁸ has been used in transition metal bulk calculations.⁷⁹ Hybrid functionals, including exact exchange, such as Becke's hybrid 3-parameter functional, B3LYP,^65,66 have also been utilized in studies involving molecular catalysts. Various hybrid functionals also use the PBE generalized gradient approximation functional, and in the present study, we focus on the full Heyd–Scuseria–Ernzerhof functional, HSE06^80,81 and the PBE0 functional of Adamo and Barone.⁸² A recent study evaluating TM complexes suggested that the MN15 functional⁸³ developed by Truhlar and coworkers gives good binding energies.^23,24 We also evaluated the hybrid functional based on the Tao, Perdew, Staroverov, and Scuseria functional, TPSSh.^84,85 Lastly, given the improvements seen by double hybrid functionals, we also evaluated the B2PLYP⁸⁶ functional by Grimme. In addition to various DFT methods, it is important to evaluate single-reference wave function methods. We also evaluated the validity of coupled cluster singles and doubles with perturbative triples, CCSD(T).^87–90 Unless specified, all the calculations were performed using Dunning's aug-cc-pVTZ basis set^91–94 in the Gaussian 16 program.⁹⁵ We note that we confirmed the atomic occupation for the various DFT results by using the natural atomic orbital analysis using Weinhold's natural bonding orbital analysis.^96,97

For 3d TM atoms from Sc to Cu, we calculated the atom energy for all possible spin states. For some atoms, we altered the initial guess orbital to obtain the correct electron occupation, which comes from the competition of 4s²3dⁿ and 4s¹3dⁿ⁺¹ occupations. For the TM–CO adduct, we focused on the linear geometry and evaluated various spin states: singlet, triplet, quintet, and septet for even numbers of electrons and doublet, quartet, and sextet for odd numbers. To evaluate the spin change for the TM + CO association reaction, we calculated the collinear association potential energy curve (PEC) of the TM–C bond length for the 2 low-lying spin states, while fixing the C–O bond length. For all TMs studied in the present study, the C–O bond length was fixed to the gas phase isolated C–O bond length obtained by each respective method. For a proper evaluation, we require a 3-dimensional potential energy surface that accounts for changes in the TM–C and C–O bond lengths, as well as the TM–C–O angle. However, we have previously found that this effect is likely to be less than 0.05 eV, and in the present study, we focus on the collinear association.⁶⁸ The collinear approach allows us to control the occupation symmetry, which is necessary for efficient multireference calculations. For the various 3d TM atoms, we calculated the potential energy from TM–C equilibrium distance R_eq −0.1 to R_eq +2.0 Å in grid intervals of 0.1 Å and from R_eq +2.0 Å to 8.0 Å in grid intervals of 0.25 Å. Since the electronic occupation during the PEC calculation depends strongly on the initial guess, we performed scan calculations by elongating and shortening the TM–C bond. The lowest energy obtained at each bond length from this dual direction energy scan is used for the PEC presented in this manuscript. Lastly, since we are calculating the PEC, we neglected the zero-point vibration, and the electronic energy is given for the different TM–C bond lengths.

Lastly, it is well-known that multireference treatment is required to properly account for the different spin states, so we also performed calculations using the multireference singles and doubles excitation configuration interaction (MRCI) method.^98,99 The electronic orbitals are obtained by the state-averaged complete active space self-consistent field (SA-CASSCF) calculation,¹⁰⁰ where both electronic orbitals and wave function expansion coefficients are optimized while minimizing the weighted sum over the total energies of multiple spin states. We included the Davidson correction,¹⁰¹ MRCI+Q, to account for the dynamic correlation and partly restore the size consistency of MRCI truncations. We used the atomic natural orbital-relativistic core-correlated triple-zeta basis set, including polarization and diffuse functions, ANO-RCC-VTZP, of the ANO-RCC family.^102,103 All MRCI+Q calculations were performed using the MOLPRO program.^104–106 For all the TM atoms and TMCO adducts, all TM 3d and CO 2p derived orbitals were taken into account for the active space of the SA-CASSCF calculations. The electronic states associated with the spin and orbital angular momentum terms for each system, which are taken for state-averaging, are presented in Table S1 of the SI. We have used more than 10 electronic states in these calculations to ensure that we are obtaining the correct symmetry spin states. We will consider the relative energies between the lowest-energy high-spin and low(next-lowest)-spin states obtained from the results of the SA-CASSCF/MRCI+Q calculations. In the PEC calculated by DFT, consistent with previous studies, we observed spin state crossings for two early TMs (Sc and Ti) and three late TMs (Fe, Co, and Ni).^{32,44,53,54,57,59} So, we calculated the PECs as a function of TM–C bond length for TM = Fe, Co and Ni over the ranges 1.5–2.4 Å (0.05 Å intervals), 2.4–3 Å (0.1 Å intervals), and 3–3.4 Å or beyond (0.2 Å intervals). For Sc and Ti, we also performed MRCI+Q calculations, but only for the short TM–CO bond lengths in the range of 1.9–2.7 Å (0.05 Å intervals). Due to the near 1 eV splitting between the two spin states for the dissociation limit (i.e., Sc and Ti atoms), the convergence of the SA-CASSCF decreases greatly at large TM–CO bond lengths for these two TMs. Since DFT calculations predicted that spin state crossings for these two TM atoms occur at short TM–CO bond lengths, we focused our MRCI+Q calculation on this short TM–CO bond length.

3 Results and discussion

3.1 Atomic spin state energies

As a first step, we evaluate the spin state splitting of the 3d TM atoms. To obtain experimental spin state energies, we follow previous methods by Truhlar and coworkers and take the average of the experimental values to remove the effect of spin–orbit splitting.²³ All calculated energies are presented in the SI. We will present the single-reference results in the order of PBE, TPSSh, PBE0, HSE06, B3LYP, MN15, B2PLYP, and CCSD(T). In Table 1, we summarize the calculated atomic ground state by the various methods for the different atoms. One can see that different DFT functionals and quantum chemistry methods can properly predict the ground electronic spin state for most 3d TM atoms. For Ti and V, some hybrid functionals have problems in predicting the ground spin state. This finding is consistent with the findings of Schultz et al.¹⁹ who found that DFT failed to predict the spin state for V. They found that DFT functionals tend to favor s¹d⁴ occupation over s²d³ for V.

Table 1 The most stable atomic state calculated by various methods in comparison with experimental results. The calculated values that are different from the experimental values are given in bold underline

Atom	PBE	TPSSh	PBE0	HSE06	B3LYP	MN15	B2PLYP	CCSD(T)	MRCI+Q	Exp
Sc	^2D (s^2d^1)	^2D (s^2d^1)	^2D (s^2d^1)	^2D (s^2d^1)	^2D (s^2d^1)	^2D (s^2d^1)	^2D (s^2d^1)	^2D (s^2d^1)	^2D (s^2d^1)	^2D (s^2d^1)
Ti	^5 F̲ (s̲ ^1 d̲ ^3 )	^5 F̲ (s̲ ^1 d̲ ^3 )	^5 F̲ (s̲ ^1 d̲ ^3 )	^5 F̲ (s̲ ^1 d̲ ^3 )	^3F (s^2d^2)	^3F (s^2d^2)	^3F (s^2d^2)	^3F (s^2d^2)	^3F (s^2d^2)	^3F (s^2d^2)
V	^4F (s^2d^3)	^6 D̲ (s̲ ^1 d̲ ^4 )	^6 D̲ (s̲ ^1 d̲ ^4 )	^6 D̲ (s̲ ^1 d̲ ^4 )	^6 D̲ (s̲ ^1 d̲ ^4 )	^6 D̲ (s̲ ^1 d̲ ^4 )	^4F (s^2d^3)	^4F (s^2d^3)	^4F (s^2d^3)	^4F (s^2d^3)
Cr	^7S (s^1d^5)	^7S (s^1d^5)	^7S (s^1d^5)	^7S (s^1d^5)	^7S (s^1d^5)	^7S (s^1d^5)	^7S (s^1d^5)	^7S (s^1d^5)	^7S (s^1d^5)	^7S (s^1d^5)
Mn	^6S (s^2d^5)	^6S (s^2d^5)	^6S (s^2d^5)	^6S (s^2d^5)	^6S (s^2d^5)	^6S (s^2d^5)	^6S (s^2d^5)	^6S (s^2d^5)	^6S (s^2d^5)	^6S (s^2d^5)
Fe	^5D (s^2d^6)	^5D (s^2d^6)	^5D (s^2d^6)	^5D (s^2d^6)	^5D (s^2d^6)	^5D (s^2d^6)	^5D (s^2d^6)	^5D (s^2d^6)	^5D (s^2d^6)	^5D (s^2d^6)
Co	^4F (s^1d^8)	^4F (s^1d^8)	^4F (s^1d^8)	^4F (s^1d^8)	^4F (s^1d^8)	^4F (s^1d^8)	^4F (s^1d^8)	^4F (s^1d^8)	^4F (s^1d^8)	^4F (s^1d^8)
Ni	^3F (s^1d^9)	^3F (s^1d^9)	^3F (s^1d^9)	^3F (s^1d^9)	^3F (s^1d^9)	^3F (s^1d^9)	^3F (s^1d^9)	^3F (s^1d^9)	^3F (s^1d^9)	^3F (s^1d^9)
Cu	^2S (s^1d^10)	^2S (s^1d^10)	^2S (s^1d^10)	^2S (s^1d^10)	^2S (s^1d^10)	^2S (s^1d^10)	^2S (s^1d^10)	^2S (s^1d^10)	^2S (s^1d^10)	^2S (s^1d^10)

---

Next, we evaluated the calculated atomic spin state splitting energies and compared them with the experimental values. Accurate experimental energies have been obtained from spectroscopic studies. We summarize the mean absolute errors (MAEs) of the various DFT functionals and quantum chemistry methods in Table 2. The detailed values are given in the SI, Tables S2–S10. As expected, the two wave function based methods CCSD(T) (0.13 eV) and MRCI+Q (0.25 eV) give the best agreement, with the former being slightly more accurate. When we evaluate the different DFT methods, we find that B3LYP gives the best agreement, with a MAE of 0.40 eV, while PBE gives the worst, with a MAE of 0.73 eV. Interestingly, TPSSh, PBE0, HSE06, and MN15 all gave similar MAE for the atomic spin splitting energies. The double hybrid B2PLYP functional gives MAE close to B3LYP. In conclusion, for atomic spin state splitting energies, the best that DFT can do is about half an eV, and we do not see much difference among the hybrid functionals, although B3LYP performs the best.

Table 2 Mean absolute error, in eV, of the atomic spin state splitting calculated by various density functionals, CCSD(T), and MRCI+Q, in comparison to experimental results

	PBE	TPSSh	PBE0	HSE06	B3LYP	MN15	B2PLYP	CCSD(T)	MRCI+Q
MAE	0.73	0.56	0.52	0.50	0.40	0.52	0.42	0.13	0.25

---

3.2 CO binding energies and geometries

In this section, we evaluate the binding energy and the equilibrium TM–C bond lengths by various methods. We divide the discussion into (1) early TM: Sc, Ti, V; (2) late TM: Fe, Co, Ni; and (3) weak binding TM: Cr, Mn, Cu. As can be seen from the electronic energies in Tables S2–S4, the early TM atom favors low spin states, but when the TM–CO bond is formed, the higher spin states are lower in energy. On the other hand, late TM atoms given in Tables S7–S9 favor high spin states, but the binding with CO results in a decrease in spin state. As for the weak binding TM atoms seen in Tables S5, S6, and S10, the spin states do not change in most methods. This finding is consistent with the previous studies by Fournier.⁴⁴ Looking at the plots in Fig. 2, we can see that the binding energies obtained by different DFT methods can result in a variation of 1 eV for each TM atom. In general, PBE tends to show stronger binding, while PBE0 and HSE06 (two hybrid functionals derived from PBE) have similar values with medium binding among the functionals we have studied. Importantly, B2PLYP and CCSD(T) tend to result in weak binding for most TM–CO adducts.

Fig. 2 TM–CO binding energy, in eV, calculated by various DFT and quantum chemistry methods for (a) TM = Sc, Ti, and V; (b) Fe, Co, and Ni; and (c) Cr, Mn, and Cu.

As shown in Fig. 2(b), experimental values for binding energies are available only for FeCO and NiCO adducts.¹⁰⁷ However, as mentioned in the introduction, there are no gas phase studies for the TM + CO reaction, so it is important to discuss the experimental detection method. Sunderlin et al. utilized collision induced dissociation of TM(CO)_n⁻ to determine the CO dissociation energy of the anion complex: D[(CO)_nTMCO⁻] = E[TM(CO)_n⁻] + E[CO] − E[TM(CO)_n+1⁻]. The electron affinity of the TM(CO)_n, EA[TM(CO)_n] = E[TM(CO)_n] − E[TM(CO)_n⁻], was used to obtain the dissociation energy of neutral D[(CO)_nTMCO] = D[(CO)_nTMCO⁻] + EA[TM(CO)_n] − EA[TM(CO)_n+1]. Therefore, the reported neutral binding energy is not directly derived from the gas phase association of TM + CO. However, these energies correspond to the difference between the most stable spin state TM + CO and the most stable spin state TMCO adduct. In conclusion, we should keep this in mind when evaluating the accuracy of the binding energy calculated by various quantum chemistry methods. As given in Table 3, we find that MRCI+Q gives the best agreement with a MAE of 0.18 eV. For the DFT functionals, all hybrid functionals result in similar MAE values, while PBE shows the largest MAE, which is the result of overbinding (see Fig. 2). In addition, we find that both the double-hybrid B2PLYP and the single reference wave function method, CCSD(T), give much worse binding energy than computationally cheaper hybrid DFT functionals. When comparing the different quantum chemistry methods to the MRCI+Q binding energy, we find that PBE0, HSE06, B3LYP, and MN15 yield the MAE < 0.2 eV. Here, we note that for the MAE comparison with MRCI+Q, we used the binding energy for both the lowest and 2nd lowest spin states, which is not available experimentally (Table 3).

Table 3 Mean absolute error (MAE), in eV, of the TM–CO binding energy compared to the experimental results for Fe and Ni, as well as to MRCI+Q results for Fe, Co, and Ni

	PBE	TPSSh	PBE0	HSE06	B3LYP	MN15	B2PLYP	CCSD(T)	MRCI+Q
MAE (vs. Exp)	1.00	0.33	0.24	0.24	0.22	0.29	0.41	0.67	0.18
MAE (vs. MRCI+Q)	0.67	0.24	0.14	0.15	0.18	0.15	0.34	0.59	—

---

In addition to the energies, it is also important to compare the binding geometries. In Fig. 3, we compare the equilibrium TM–C bond distance for the early and late TM systems. As seen in Fig. 3(a), for the early TM atoms, the general trend is consistent among all three TM atoms, where PBE0, B3LYP and B2PLYP result in longer TM–C bond lengths compared to the other functionals. As for the late TM atoms given in Fig. 3(b), PBE0, HSE06, B3LYP, and B2PLYP also result in a longer TM–C bond length for Fe and Co. On the other hand, MN15 tends to give short TM–CO bond lengths. For CCSD(T), we obtain the longest TM–C bond length among all the methods for both early and late TM systems, and this is consistent with the weak binding energy seen for this method. Compared to the TM–C equilibrium bond length obtained from MRCI+Q, we see that hybrid DFT results in shorter bond lengths for Sc and Ti, while having longer bond lengths for Fe, Co, and Ni. So the trend in the early and late TM atoms is different. However, we note that all hybrid DFT methods yield values within 0.05 Å of each other, indicating that the variation in geometry is very minor.

Fig. 3 TM–CO bond length, in Å, calculated by various DFT and quantum chemistry methods for (a) TM = Sc, Ti, and V; and (b) Fe, Co, and Ni.

3.3 Spin state crossing geometries and energies

After evaluating the atom spin state splitting and the TM–CO adduct binding characteristics, we focus on the spin state crossing geometries. In Fig. 4, we compare the association PECs calculated by B3LYP and MN15 versus those obtained by MRCI+Q. A comparison of other DFT functionals, such as PBE, TPSSh, HSE06, PBE0, and B2PLYP, is provided in Fig. S1–S3 of the SI. MRCI+Q PECs for Sc and Ti are given in Fig. S4 of the SI. For Ni, we see that in the bonding region between Ni–C bond lengths of 1.5 to 2.2 Å, the DFT methods closely match the MRCI+Q singlet PEC. For the triplet Ni + CO PEC, although B3LYP and MN15 overestimate the binding energy, the crossing between the singlet and triplet PEC occurs at a Ni–C bond length of 2.35 Å. This is consistent with 2.3 Å obtained from the MRCI+Q calculations. In the case of Co + CO, we observe that both DFT methods underestimate the binding energy for the quartet PECs relative to MRCI+Q, and the crossing point is at the Co–C bond length of 2.05 and 2.15 Å for B3LYP and MN15, respectively. Although this is slightly longer than the 1.85 Å obtained from MRCI+Q, considering the approximations used in this study, such as the collinear association geometry and the use of a fixed C–O bond length, we can say that this is also a decent agreement. Lastly, for Fe + CO, B3LYP greatly underestimates the binding energies compared to MRCI+Q, while MN15 gives similar values. However, both methods have the crossing at the minima of the two spin state PECs at an Fe–C distance of 1.95 Å, which is 0.2 Å longer than 1.75 Å obtained for MRCI+Q. As given in Fig. S1–S3 in the SI, although the absolute binding energies showed deviations from the MRCI+Q results, all DFT functionals give smooth PECs. In addition, Fig. 5 shows that hybrid functionals, surprisingly, give consistent results compared to MRCI+Q for the TM–C spin state crossing bond length. Regarding Ti + CO, we observed spin state crossings for B3LYP, MN15, B2PLYP, CCSD(T), and MRCI+Q at Ti–C distances of 4.0, 2.75, 2.95, 2.45, and 2.50 Å, respectively. Unfortunately, we do not see a simple relationship in the crossing point compared to the MRCI+Q results. In Fig. 6, we present the PECs of CCSD(T), which gave the best results for the TM atom spin state splitting. The CCSD(T) wave function method results in bumpy PECs. For Fe and Co, they vary greatly from the MRCI+Q PEC and fail to give the general trend. This surprising result shows the robustness of hybrid DFT functionals and is consistent with the finding of Morgante and Peverati, who reported that double hybrids failed to give good binding energies for CO binding to TM porphyrin complexes.²⁴ In conclusion, the present study shows that hybrid functionals are capable of describing spin state coupling for the TM–CO bond formation.

Fig. 4 Potential energy curve for the TM–CO collinear association calculated by the (a) B3LYP and (b) MN15 DFT functionals and the MRCI+Q methods. The solid lines represent the results for MRCI+Q, while the dotted ones represent the DFT results. Red and blue colors represent the low and high spin states for MRCI+Q, while black and green colors represent the low and high spin states for DFT. The crossing point is given by red and black circles for MRCI+Q and DFT, respectively.

Fig. 5 Spin state crossing TM–CO bond length, in Å, calculated by various DFT and quantum chemistry methods for TM = Sc, Fe, Co, and Ni.

Fig. 6 Potential energy curve for the TM–CO collinear association calculated by the CCSD(T) and MRCI+Q methods for (a) Fe, (b), Co, and (c) Ni. The solid lines represent the results for MRCI+Q, while the dotted ones represent the CCSD(T) results. Red and blue colors represent the low and high spin states for MRCI+Q, while black and green colors represent the low and high spin states for CCSD(T).

3.4 Density functional dependence

From the systematic study performed on the simple TM + CO system presented here, it would be ideal to obtain systematic knowledge on which DFT functionals would be reliable for larger TM complexes. In the previous sections, we saw that the TM–CO bond length of the stable adduct shows a small DFT functional dependence. On the other hand, we did not observe a simple trend in the spin state crossing geometries shown in Fig. 5. Since crossing geometries depend on the energies at the asymptotes, we plotted the correlation between the spin state crossing TM–CO bond lengths and the atomic splitting energies in Fig. 7. Although not perfect, the smaller the splitting, the longer the TM–CO bond lengths. Considering that the atomic spin state splitting depends on the 4s–3d energy splitting, one may expect to see a similar correlation from the calculated 4s–3d orbital energy difference. Unfortunately, as shown in Fig. S4, we did not observe a simple correlation. However, this confirms that the accuracy of the spin-state splitting of the bare TM complexes serves as a basis for evaluating the method for MECP calculation. On the other hand, CCSD(T), which gave the best atomic splitting compared to the experimental results in Table 1, fails to yield a smooth PEC and underestimates the binding energy. We will discuss this issue in the next paragraph.

Fig. 7 Relationship between crossing TM–C bond length and atomic splitting energies calculated by various DFT and wave function methods.

If we evaluate the binding energy for the TM–CO in Fig. 2, we observe the general trend E[PBE] < E[TPSSh] < E[PBE0] < E[HSE06] < E[B3LYP] < E[MN15] < E[B2PLYP] < E[CCSD(T)]. As shown in Fig. 1 and discussed in many previous studies, the π-back donation is critical for TM–CO binding, so one could compare the energy difference between the occupied donating 3d orbitals (specifically the 3d_xz and 3d_yz) and the lowest unoccupied orbital (LUMO) of the accepting CO. In Fig. 8, we plot the relationship between the binding energy and orbital energy difference for all DFT functionals. Although the relationship is not perfect, we can see that the smaller the orbital energy difference, the stronger the TM–CO binding (more negative TM–CO binding energy). Since B2PLYP (Fig. 8 triangle symbol) includes energy correction by perturbation terms, the orbital energy is not a direct indicator of the electronic energy. However, it is interesting to see that it follows the general trend. We note that for Hartree–Fock orbitals, the orbital energy difference between the 3d orbital and the CO LUMO is much larger than DFT. Consistently, the TM–CO binding energy is very small in Hartree–Fock calculations, and we believe that the issues with CCSD(T) stem from the failure in the Hartree–Fock method. For atomic spin state splitting, CCSD(T) correlation energies can correct the limitations of Hartree–Fock, but it fails for the TM–CO binding energy. We note that such a relationship between orbital energy differences and binding energies is consistent with comparing the adduct orbital energies with the d-band center of the bulk metal surface.^71–74,77

Fig. 8 Relationship between TM + CO binding energies and orbital energy difference between the CO LUMO and TM occupied d orbitals calculated by various DFT methods for (a) early, (b) later and (c) weak binding TMs. We subtracted the occupied 3d orbital energy from the CO LUMO energy calculated by the respective methods.

3.5 Basis set dependence

Before concluding the paper, it is also important to discuss the basis set dependence of the spin state characteristics. In the present study, we evaluated a small system of a TM and a CO molecule, allowing us to use a fairly large triple-zeta basis set. However, in many studies involving TM complexes, to account for the large ligands that coordinate to the TM atoms, one usually resorts to smaller double-zeta basis sets or effective core potentials (ECPs). Therefore, we think it will also be critical to test the basis set dependence. So we also evaluated the B3LYP spin state splitting and 3d TM–CO binding energy using Pople's 6-31+G(d,p) basis, Los Alamos ECP plus DZ (LANL2DZ),^108,109 and Stuttgart/Dresden ECPs (SDD).^110,111 In Tables S11–S19, we present the energies calculated by various basis sets. We performed the calculation using 6-31+G(d,p) for all TM, C, and O atoms, as well as those that used 6-31+G(d,p) for C and O atoms, while ECP was used for the TM. The Gaussian16 program uses the 6D basis set as the default for the Pople basis sets, so the ECP calculations were also performed using the 6D basis for LANL2DZ and SDD although 5D is the default for these ECPs. For the mean absolute error compared to the aug-cc-pVTZ results, we found that 6-31+G(d,p) gives a very small value of 0.0 eV. However, when we used ECPs, the MAE increased to 0.25 and 0.28 eV for LANL2DZ and SDD, respectively. The large difference comes from the large error in MnCO binding energy and if we remove this data, we obtain MAE of 0.07 and 0.08 for LANL2DZ and SDD, respectively. In Fig. 9, we present the B3LYP TM + CO association PECs for TM = Sc, Ti, V, Mn, Cr, Fe, Co, and Ni calculated by aug-cc-pVTZ and 6-31+G(d,p). As DFT is known to have small basis set dependence, we expected similarity between the two full electron calculations. However, we did not expect to find such a perfect match between the PECs for all 8 TM atoms.

Fig. 9 Potential energy curve for the TM–CO collinear association calculated by B3LYP with aug-cc-pVTZ and 6-31+G(d,p) for TM = (a) Sc, (b) Ti, (c) V, (d) Cr, (e) Mn, (f) Fe, (g) Co, and (h) Ni. The solid lines represent the results for aug-cc-pVTZ, while the dashed ones represent those by 6-31+G(d,p). Blue and red colors represent the low and high spin states for aug-cc-pVTZ, while green and black colors represent the low and high spin states for 6-31+G(d,p).

4 Conclusions

In this study, we evaluated the atomic spin state splitting of 3d TM atoms and the 3d TM–CO binding energy using various DFT functionals and quantum chemistry methods. Furthermore, we also calculated the TM + CO collinear association potential energy curve to determine the spin state crossing the TM–C bond length. We compared our DFT and CCSD(T) potential energy curves with those obtained by multireference methods for Fe, Co, and Ni. For DFT, the PBE functional resulted in the largest error for atomic spin state splitting energies and showed an over-binding characteristic for the TM–CO adduct. The latter over-binding is attributed to the small orbital energy difference between the occupied 3d orbital of the TM and the CO lowest unoccupied orbital. As a general trend, TM–CO equilibrium bond lengths and spin state crossing TM–CO distances, i.e. the MECP, showed small hybrid functional dependence. Among the various single reference DFT functionals, B3LYP and MN15 gave consistent atomic spin state splitting with experiments and showed the best agreement with MRCI+Q for TMCO binding energy and spin crossing characteristics. Our studies confirm that accurate spin state splitting for the bare TM complex is a key quantity to evaluate for obtaining good MECP geometries. We found that CCSD(T) and B2PLYP underestimates the TM–CO binding energy. Accordingly, the potential energy curves for the various spin states of the 3d TM + CO association obtained by these two methods showed trends different from those of MRCI+Q. To our surprise, B3LYP TM + CO association potential energy curves calculated by aug-cc-pVTZ and 6-31+G(d,p) showed a perfect match. This shows that one may use the smaller 6-31+G(d,p) basis set to study CO binding to larger TM complexes.

Lastly, we found that for TM = Sc, Fe, and Co, the crossing between the two lowest spin states occurs at the TM–CO bond length close to the adduct minima. Our previous studies on Ni + CO showed that spin orbit coupling decreases greatly as the Ni–CO bond is formed.¹⁴ Therefore, we believe that this is why Lester and coworkers did not observe the CO association reaction for the ground spin state of these TM atoms.^53,54

Author contributions

Natthakrij Nipanutiyan: software and investigation. Daisuke Yoshida: conceptualization, software, validation, formal analysis, investigation, writing – review and editing, and visualization. Kaito Takahashi: conceptualization, software, validation, investigation, resources, writing – original draft, writing – review and editing, supervision, project administration, and funding acquisition.

Conflicts of interest

The authors declare no competing financial interest.

Data availability

The atomic spin state splitting and binding energies, potential energy curve, density functional dependence, and basis set dependence are included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5cp04203h.

Acknowledgements

We thank the computational resources from NSTDA supercomputer center of Thailand (ThaiSC) and the Institute of Atomic and Molecular Sciences, Academia Sinica Taiwan. KT thanks the Young Researcher Grant from Sirindhorn International Institute of Technology, Thammasat University (SIIT-2024-YRG-KT02). We thank Prof. Jumras Limtrakul for his continued efforts in promoting computational chemistry, catalysis, and materials science in Thailand.

References

1. V.-H. Do and J.-M. Lee, ACS Nano, 2022, 16, 17847–17890.
2. Y. Wang, J. Sun, N. Sun, M. Zhang, X. Liu, A. Zhang and L. Wang, Chem. Commun., 2024, 60, 7397–7413.
3. L. Pauling and C. D. Coryell, Proc. Natl. Acad. Sci. U. S. A., 1936, 22, 210–216.
4. K. Falahati, H. Tamura, I. Burghardt and M. Huix-Rotllant, Nat. Commun., 2018, 9, 4502.
5. M. E. Ali, B. Sanyal and P. M. Oppeneer, J. Phys. Chem. B, 2012, 116, 5849–5859.
6. J.-H. Wu, J.-W. Wang, B. M. Aramburu-Trošelj, F.-J. Niu, L.-J. Guo and G. Ouyang, Nanoscale, 2024, 16, 11496–11512.
7. S. Shaik, Isr. J. Chem., 2020, 60, 938–956.
8. D. Schröder, S. Shaik and H. Schwarz, Acc. Chem. Res., 2000, 33, 139–145.
9. M. Besora, J.-L. Carreón-Macedo, Á. Cimas and J. N. Harvey, in Adv. Inorg. Chem., ed. R. van Eldik and C. D. Hubbard, Academic Press, 2009, vol. 61, pp. 573–623.
10. L. Zhao, K.-J. Watanabe, N. Nakatani, A. Nakayama, X. Xu and J.-Y. Hasegawa, J. Chem. Phys., 2020, 153, 134114.
11. V. D. Dergachev, M. Rooein, I. D. Dergachev, A. O. Lykhin, R. C. Mauban and S. A. Varganov, Top. Curr. Chem., 2022, 380, 15.
12. B. Yang, L. Gagliardi and D. G. Truhlar, Phys. Chem. Chem. Phys., 2018, 20, 4129–4136.
13. L. Zhao and W. Zou, J. Chem. Phys., 2023, 158, 224110.
14. D. Yoshida and K. Takahashi, Inorg. Chem., 2025, 64, 286–294.
15. T. Takayanagi, Comput. Theor. Chem., 2022, 1211, 113682.
16. D. Ricciarelli, Q. M. Phung, L. Belpassi, J. N. Harvey and P. Belanzoni, Inorg. Chem., 2019, 58, 7345–7356.
17. D. Ricciarelli, L. Belpassi, J. N. Harvey and P. Belanzoni, Chem. – Eur. J., 2020, 26, 3080–3089.
18. E. R. Davidson, K. L. Kunze, F. B. C. Machado and S. J. Chakravorty, Acc. Chem. Res., 1993, 26, 628–635.
19. N. E. Schultz, Y. Zhao and D. G. Truhlar, J. Phys. Chem. A, 2005, 109, 4388–4403.
20. S. M. Tekarli, M. L. Drummond, T. G. Williams, T. R. Cundari and A. K. Wilson, J. Phys. Chem. A, 2009, 113, 8607–8614.
21. W. Jiang, N. J. DeYonker and A. K. Wilson, J. Chem. Theory Comput., 2012, 8, 460–468.
22. K. Pierloot, Q. M. Phung and A. Domingo, J. Chem. Theory Comput., 2017, 13, 537–553.
23. D. Zhang and D. G. Truhlar, J. Chem. Theory Comput., 2020, 16, 4416–4428.
24. P. Morgante and R. Peverati, Molecules, 2023, 28, 3487.
25. R. Hoffmann, Angew. Chem., Int. Ed. Engl., 1987, 26, 846–878.
26. R. G. Pearson, Inorg. Chem., 1984, 23, 4675–4679.
27. S. B. H. Bach, C. A. Taylor, R. J. Van Zee, M. T. Vala and W. Weltner, J. Am. Chem. Soc., 1986, 108, 7104–7105.
28. R. L. DeKock, Inorg. Chem., 1971, 10, 1205–1211.
29. E. P. Kündig, D. McIntosh, M. Moskovits and G. A. Ozin, J. Am. Chem. Soc., 1973, 95, 7234–7241.
30. H. Huber, E. P. Kündig, G. A. Ozin and A. J. Poe, J. Am. Chem. Soc., 1975, 97, 308–314.
31. L. A. Hanlan, H. Huber, E. P. Kündig, B. R. McGarvey and G. A. Ozin, J. Am. Chem. Soc., 1975, 97, 7054–7068.
32. B. Tremblay, M. E. Alikhani and L. Manceron, J. Phys. Chem. A, 2001, 105, 11388–11394.
33. H. Huber, E. P. Kündig, M. Moskovits and G. A. Ozin, J. Am. Chem. Soc., 1975, 97, 2097–2106.
34. L. Hanlan, H. Huber and G. A. Ozin, Inorg. Chem., 1976, 15, 2592–2597.
35. P. C. Engelking and W. C. Lineberger, J. Am. Chem. Soc., 1979, 101, 5569–5573.
36. C. H. F. Peden, S. F. Parker, P. H. Barrett and R. G. Pearson, J. Phys. Chem., 1983, 87, 2329–2336.
37. A. E. Stevens, C. S. Feigerle and W. C. Lineberger, J. Am. Chem. Soc., 1982, 104, 5026–5031.
38. K. Tanaka, K. Sakaguchi and T. Tanaka, J. Chem. Phys., 1997, 106, 2118–2128.
39. L. Manceron, M. E. Alikhani and H. A. Joly, Chem. Phys., 1998, 228, 73–80.
40. M. R. A. Blomberg, U. B. Brandemark, P. E. M. Siegbahn, J. Wennerberg and C. W. Bauschlicher, J. Am. Chem. Soc., 1988, 110, 6650–6655.
41. L. A. Barnes and C. W. Bauschlicher, Jr., J. Chem. Phys., 1989, 91, 314–330.
42. M. R. A. Blomberg, P. E. M. Siegbahn, T. J. Lee, A. P. Rendell and J. E. Rice, J. Chem. Phys., 1991, 95, 5898–5905.
43. K. Hirai and N. Kosugi, Can. J. Chem., 1992, 70, 301–308.
44. R. Fournier, J. Chem. Phys., 1993, 99, 1801–1815.
45. M. Castro, D. R. Salahub and R. Fournier, J. Chem. Phys., 1994, 100, 8233–8239.
46. J. Chatt and L. A. Duncanson, J. Chem. Soc., 1953, 2939–2947,  10.1039/JR9530002939.
47. G. Blyholder, J. Phys. Chem., 1964, 68, 2772–2777.
48. C. W. Bauschlicher, S. R. Langhoff and L. A. Barnes, Chem. Phys., 1989, 129, 431–437.
49. C. W. Bauschlicher, Jr., P. S. Bagus, C. J. Nelin and B. O. Roos, J. Chem. Phys., 1986, 85, 354–364.
50. C. W. Bauschlicher, Chem. Phys. Lett., 1996, 249, 244–248.
51. C. W. Bauschlicher, Jr., J. Chem. Phys., 1994, 100, 1215–1218.
52. X. Xu, X. Lü, N. Wang, Q. Zhang, M. Ehara and H. Nakatsuji, Int. J. Quantum Chem., 1999, 72, 221–231.
53. M. Zhou and L. Andrews, J. Phys. Chem. A, 1999, 103, 2964–2971.
54. M. Zhou and L. Andrews, J. Phys. Chem. A, 1999, 103, 5259–5268.
55. B. Liang, M. Zhou and L. Andrews, J. Phys. Chem. A, 2000, 104, 3905–3914.
56. L. Andrews, M. Zhou, X. Wang and C. W. Bauschlicher, J. Phys. Chem. A, 2000, 104, 8887–8897.
57. X. Wang, M. Zhou and L. Andrews, J. Phys. Chem. A, 2000, 104, 10104–10111.
58. M. Zhou, G. V. Chertihin and L. Andrews, J. Chem. Phys., 1998, 109, 10893–10904.
59. M. Zhou and L. Andrews, J. Phys. Chem. A, 1998, 102, 10250–10257.
60. M. Zhou and L. Andrews, J. Chem. Phys., 1999, 111, 4548–4557.
61. J. J. Turner, M. W. George, M. Poliakoff and R. N. Perutz, Chem. Soc. Rev., 2022, 51, 5300–5329.
62. R. Paciotti, C. Coletti, E. Berrino, F. Arrighi, A. Maccelli, A. Lasalvia, M. E. Crestoni, D. Secci, S. Carradori, C. T. Supuran and F. Carta, Int. J. Mol. Sci., 2024, 25, 11644.
63. J. Troß, J. E. Arias-Martinez, K. Carter-Fenk, N. C. Cole-Filipiak, P. Schrader, L. M. McCaslin, M. Head-Gordon and K. Ramasesha, J. Am. Chem. Soc., 2024, 146, 22711–22723.
64. E. Rossomme, C. N. Lininger, A. T. Bell, T. Head-Gordon and M. Head-Gordon, Phys. Chem. Chem. Phys., 2020, 22, 781–798.
65. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
66. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B:Condens. Matter Mater. Phys., 1988, 37, 785–789.
67. K. Takahashi, PACCON2025 e-Proceedings, 2025, vol. 1.
68. N. Natthakrij and T. Kaito, Sci. Technol. Asia, 2025, 30, 120–127.
69. A. Stroppa and G. Kresse, New J. Phys., 2008, 10, 063020.
70. K. M. Gameel, I. M. Sharafeldin, A. U. Abourayya, A. H. Biby and N. K. Allam, Phys. Chem. Chem. Phys., 2018, 20, 25892–25900.
71. A. Patra, H. Peng, J. Sun and J. P. Perdew, Phys. Rev. B, 2019, 100, 035442.
72. M. Gajdoš and J. Hafner, Surf. Sci., 2005, 590, 117–126.
73. A. Föhlisch, M. Nyberg, J. Hasselström, O. Karis, L. G. M. Pettersson and A. Nilsson, Phys. Rev. Lett., 2000, 85, 3309–3312.
74. G. Kresse, A. Gil and P. Sautet, Phys. Rev. B:Condens. Matter Mater. Phys., 2003, 68, 073401.
75. P. Janthon, F. Viñes, J. Sirijaraensre, J. Limtrakul and F. Illas, J. Phys. Chem. C, 2017, 121, 3970–3977.
76. F. Abild-Pedersen and M. P. Andersson, Surf. Sci., 2007, 601, 1747–1753.
77. M. T. M. Koper, R. A. van Santen, S. A. Wasileski and M. J. Weaver, J. Chem. Phys., 2000, 113, 4392–4407.
78. J. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
79. P. Janthon, S. Luo, S. M. Kozlov, F. Viñes, J. Limtrakul, D. G. Truhlar and F. Illas, J. Chem. Theory Comput., 2014, 10, 3832–3839.
80. J. Heyd and G. E. Scuseria, J. Chem. Phys., 2004, 121, 1187–1192.
81. J. Heyd, J. E. Peralta, G. E. Scuseria and R. L. Martin, J. Chem. Phys., 2005, 123, 174101.
82. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
83. H. S. Yu, X. He, S. L. Li and D. G. Truhlar, Chem. Sci., 2016, 7, 5032–5051.
84. J. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria, Phys. Rev. Lett., 2003, 91, 146401.
85. V. N. Staroverov, G. E. Scuseria, J. Tao and J. P. Perdew, J. Chem. Phys., 2003, 119, 12129–12137.
86. S. Grimme, J. Chem. Phys., 2006, 124, 034108.
87. R. J. Bartlett and G. Purvis, Int. J. Quantum Chem., 1978, XIV, 561–581.
88. G. E. Scuseria, C. L. Janssen and H. F. Schaefer, III, J. Chem. Phys., 1988, 89, 7382–7387.
89. J. D. Watts, J. R. Gauss and R. J. Bartlett, J. Chem. Phys., 1993, 98, 8718.
90. J. A. Pople, M. Head-Gordon and K. Raghavachari, J. Chem. Phys., 1987, 87, 5968–5975.
91. T. H. Dunning, J. Chem. Phys., 1989, 90, 1007–1023.
92. R. A. Kendall, T. H. Dunning and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796–6806.
93. D. E. Woon and T. H. Dunning, Jr., J. Chem. Phys., 1993, 98, 1358–1371.
94. K. A. Peterson, D. E. Woon and T. H. Dunning, J. Chem. Phys., 1994, 100, 7410–7415.
95. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian 16, Revision C.02, Gaussian, Inc., Wallingford CT, 2019.
96. J. P. Foster and F. Weinhold, J. Am. Chem. Soc., 1980, 102, 7211–7218.
97. A. E. Reed, R. B. Weinstock and F. Weinhold, J. Chem. Phys., 1985, 83, 735–746.
98. H. J. Werner and P. J. Knowles, J. Chem. Phys., 1988, 89, 5803–5814.
99. P. J. Knowles and H.-J. Werner, Theor. Chim. Acta, 1992, 84, 95–103.
100. D. A. Kreplin, P. J. Knowles and H.-J. Werner, J. Chem. Phys., 2019, 150, 194106.
101. S. R. Langhoff and E. R. Davidson, Int. J. Quantum Chem., 1974, 8, 61–72.
102. B. O. Roos, R. Lindh, P.-Å. Malmqvist, V. Veryazov and P.-O. Widmark, J. Phys. Chem. A, 2004, 108, 2851–2858.
103. B. O. Roos, R. Lindh, P.-Å. Malmqvist, V. Veryazov and P.-O. Widmark, J. Phys. Chem. A, 2005, 109, 6575–6579.
104. H.-J. Werner and P. J. Knowleset al., MOLPRO, version 2021.3, a package of ab initio programs, https://www.molpro.net).
105. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schütz, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2012, 2, 242–253.
106. H.-J. Werner, P. J. Knowles, F. R. Manby, J. A. Black, K. Doll, A. Heßelmann, D. Kats, A. Köhn, T. Korona, D. A. Kreplin, Q. Ma, T. F. Miller, III, A. Mitrushchenkov, K. A. Peterson, I. Polyak, G. Rauhut and M. Sibaev, J. Chem. Phys., 2020, 152, 144107.
107. L. S. Sunderlin, D. Wang and R. R. Squires, J. Am. Chem. Soc., 1992, 114, 2788–2796.
108. P. J. Hay and W. R. Wadt, J. Chem. Phys., 1985, 82, 270–283.
109. W. R. Wadt and P. J. Hay, J. Chem. Phys., 1985, 82, 284–298.
110. H. Stoll, P. Fuentealba, P. Schwerdtfeger, J. Flad, L. V. Szentpály and H. Preuss, J. Chem. Phys., 1984, 81, 2732–2736.
111. M. Dolg, U. Wedig, H. Stoll and H. Preuss, J. Chem. Phys., 1987, 86, 866–872.

---