Gas-phase and solid-state electronic structure analysis and DFT benchmarking of HfCO

Abstract

Ab initio multi-reference configuration interaction (MRCI) and coupled cluster singles doubles and perturbative triples [CCSD(T)] levels of theory were used to study ground and excited electronic states of HfCO. We report potential energy curves, dissociation energies (D_e), excitation energies, harmonic vibrational frequencies, and chemical bonding patterns of HfCO. The ³Σ⁻ ground state of HfCO has an 1σ²2σ²1π² electron configuration and a ∼30 kcal mol⁻¹ dissociation energy with respect to its lowest-energy fragments Hf(³F) + CO(X¹Σ⁺). We further evaluated the D_e of its isovalent HfCX (X = S, Se, Te, Po) series and observed that they increase linearly from the lighter HfCO to the heavier HfCPo with the dipole moment of the CX ligand. The same linear relationship was observed for TiCX and ZrCX. We utilized the CCSD(T) benchmark values of D_e, excitation energy, and ionization energy (IE) values to evaluate density functional theory (DFT) errors with 23 exchange–correlation functionals spanning GGA, meta-GGA, global GGA hybrid, meta-GGA hybrid, range-separated hybrid, and double-hybrid functional families. The global GGA hybrid B3LYP and range-separated hybrid ωB97X performed well at representing the ground state properties of HfCO (i.e., D_e and IE). Finally, we extended our DFT analysis to the interaction of a CO molecule with a Hf surface and observed that the surface chemisorption energy and the gas-phase molecular dissociation energy are very similar for some DFAs but not others, suggesting moderate transferability of the benchmarks on these molecules to the solid state.

---

I. Introduction

Transition metal (TM) carbonyl complexes are frequently encountered in synthetic organometallic chemistry and are cornerstones of modern coordination chemistry. Owing to the diverse physicochemical properties and the tunability of their reactivity through modification of the metal center, they have numerous industrial,¹ biochemical,² and catalytic^3,4 applications. Hence, isolated gas-phase TM(CO)_n models and the chemisorption of CO on metal surfaces are often investigated to shed light on their catalytic potencies.^5–9

The fact that the dominant binding mode of CO with electropositive metals is via the less electronegative C rather than the more electronegative O is at first glance surprising.¹⁰ This unique characteristic of CO arises from the two-electron σ lone pair present on the C, making it the strongest ligand in the spectrochemical series. Often, M⋯CO interactions are further supported by the superior π-acceptor ability of CO that induces metal-to-ligand π back-donation.¹⁰ However, the electron flow through the π-frame (or to π* orbitals of CO) simultaneously weakens the C–O bond.¹⁰

To date, a majority of transition metal CO studies are focused on the first- and second-row TMs, with less attention on the heavier third-row TMs. For example, according to the Cambridge Structural Database (CSD, 2021.1.0) there are 164 first-row (TM = Sc–V) and 185 second-row (TM = Y–Nb) crystal structures available with the TM–C≡O bonding configuration, whereas the number of complexes for the corresponding search for third-row structures (TM = La–Ta) is only 97.¹¹ Among the third-row TMs, one of the least studied systems is Hf, of which only 9 hits are obtained from the CSD for a Hf–C≡O search.¹¹ While this analysis does not include a review of the solution or gas-phase chemistry studies for third-row structures, it nevertheless gives a sense of the proportion of studies that have been dedicated to third-row complexes.

In reviewing the literature, we were only able to locate two experimental gas-phase studies on Hf(CO)_n systems, which are related to infrared spectra analysis.^12,13 Work by Zhou and Andrews provides evidence for the existence of all Hf(CO)_1–4 and Hf(CO)₂⁻ species.¹² They also report the stretching frequency of CO of HfCO to be 1869 cm⁻¹ (obtained from laser-ablated hafnium in solid neon) and a trend of greater d → π* back-donation going from TiCO to HfCO.¹² This trend is consistent with the spectroscopic study of M(CO)_6–8⁺ (M = Ti–Hf) by Brathwaite and Duncan.¹³

Computational modelling is vital in understanding the chemistry of systems that are challenging to handle under laboratory conditions.^14–16 For example, Hf compounds are known to be highly toxic, motivating theoretical studies.^17,18 However, selection of appropriate computational tools is crucial to make accurate predictions. In general, unsaturated TM(CO)_n complexes have several close-lying electronic states, which make them challenging to model accurately with computation.^14,19–23 High-level multi-reference theories provide a platform for predictions of systems with many low-lying electronic states.^24,25 The gold-standard coupled cluster level of theory can also deliver predictions with higher precision relative to most computationally tractable multi-reference (MR) methods as long as the targeted states are sufficiently single-reference (SR) in nature.²⁶ However, both MR and SR coupled cluster levels of theory require a great deal of quantum chemical expertise and computing power, which makes them less widely used, especially in larger systems.^24,25 On the other hand, density functional theory (DFT) is widely popular due to its black box nature and reasonable accuracy at low computational cost.^27–29 However, a practitioner must select an appropriate exchange–correlation functional because the best-performing functional depends on the system and property.³⁰

In this work, we probe the interactions between Hf and CO with high-level multi-reference configuration interaction theory (MRCI), coupled cluster singles doubles and perturbative triples [CCSD(T)], as well as with DFT. Specifically, numerous low-lying electronic states of HfCO are analyzed with MRCI and their potential energy profiles, equilibrium electronic configurations, corresponding chemical bonding patterns, and various energy-related properties are reported. The single-reference electronic states are investigated with CCSD(T) and used as benchmarks for evaluating DFT errors associated with 23 exchange–correlation functionals that span multiple rungs of “Jacob's ladder”.³¹ We assess whether performance trends are transferable to solid-state systems by studying the chemisorption of one CO on a Hf surface with a set of DFT functionals.

II. Computational details

Internally contracted multi-reference configuration interaction (MRCI)^32–34 and coupled cluster singles doubles and perturbative triples [CCSD(T)]²⁶ correlated wavefunction theory (WFT) calculations were carried out with the MOLPRO 2015.1 code.³⁵ These calculations used the largest Abelian subgroup, C_2v, of the parent C_∞v point group associated with the linear species HfCX (X = O, S, Se, Te, Po).

First, at the MRCI level we obtained full potential energy curves (PECs) of HfCO arising from Hf(³F) + CO(X¹Σ⁺), Hf(³P) + CO(X¹Σ⁺), Hf(¹D) + CO(X¹Σ⁺), and Hf(⁵F) + CO(X¹Σ⁺) interactions as a function of Hf⋯C distance. For each potential energy scan, the C–O distance was kept fixed to the experimental bond distance (i.e., 1.128 Å) of the isolated CO molecule.³⁶ This choice was made to simplify the scans, although we had previously noted backbonding could weaken the CO bond. Hence, we do not attempt to compute fundamental frequencies from these potential energy curves since they would be purely approximate. We nevertheless expect the effect on the bond distance to be modest and comparable across all species and electronic states compared. MRCI calculations were initiated from reference complete active space self-consistent field (CASSCF) wavefunctions. Specifically, the CASSCF wavefunctions were constructed by placing 4 electrons in 6 active orbitals [CAS (4,6)]. At long Hf⋯C distances (>5 Å), the selected orbitals are purely 6s and all 5d atomic orbitals of Hf, corresponding to the four active electrons. At the employed C_2v point group, the active orbitals are 3a₁ [5d_z², 5d_x²–y², 6s], 1b₁ [5d_xz], 1b₂ [5d_yz], and 1a₂ [5d_xy] in symmetry. All valence electrons were correlated in the subsequent MRCI calculation. A full MRCI geometry optimization was performed only for the ground triplet and lowest-lying quintet electronic states of HfCO. Davidson-corrected MRCI (MRCI+Q)³⁷ energies obtained at the ground state MRCI geometry were also used to compute the dissociation energy (D_e) and the excitation energies (T_e) of HfCO. At the MRCI level, spin–orbit coupling effects were evaluated by a single-point calculation of the MRCI ground state geometry using the Breit–Pauli Hamiltonian as implemented in MOLPRO.

Based on the CASSCF dominant electronic configurations, Hartree–Fock wavefunctions were constructed for single-reference electronic states of HfCO and used as starting points for CCSD(T) geometry optimizations. At the same level of theory, geometries of isovalent HfCX (X = S, Se, Te, Po) were optimized to evaluate their D_es.

For WFT calculations of HfCO, correlation-consistent cc-pVXZ-PP³⁸ basis sets were chosen for Hf and aug-cc-pVXZ³⁹ for C and O, where X = T, Q or 5. Specifically, only a quadruple-ζ quality (X = Q) basis set was used for all MRCI analysis, but we also compared results from triple-ζ, quadruple-ζ, and quintuple-ζ quality sets for CCSD(T) calculations. The plain cc-pVTZ(-PP) basis set was also tested for all atoms of HfCO and exclusively employed for the isovalent HfCS, HfCSe, HfCTe, and HfCPo species to reduce the computational expense.^40–42 In all cases, the inner 60 electrons of Hf (1s²2s²2p⁶3s²3p⁶4s²3d¹⁰4p⁶4d¹⁰4f¹⁴) were substituted with the Stuttgart relativistic pseudopotential (ECP60).³⁸

HfCO, HfCS, and HfCSe WFT energetics were used to investigate density functional theory (DFT) errors associated with 23 density functional approximations (DFAs) that fall into six rungs of “Jacob's ladder” using the Psi4⁴³ package. The 23 DFAs corresponded to semi-local generalized gradient approximations (GGAs) (BLYP, BP86, and PBE), meta-GGAs (TPSS, SCAN, M06-L, and MN15-L), global GGA hybrids (B3LYP, B3P86, B3PW91, and PBE0), meta-GGA hybrids (TPSSh, SCAN0, M06, M06-2X, and MN15), range-separated hybrids (LRC-ωPBEh and ωB97X), and double hybrids (B2GP-BLYP, PBE0-DH, DSDBLYP-D3BJ, DSD-PBEB95-D3BJ, and DSD-PBEP86-D3BJ) functionals, as implemented in a recently introduced workflow.^31,44 In this workflow, the density of a B3LYP calculation is converged first and used as the starting point for all other DFA calculations. Starting from CCSD(T) geometries, the results from these DFAs were used to evaluate the single-point D_e for HfCO, HfCS, and HfCSe. In the past we have applied def2-XZVP (X = T, Q) basis sets to study TM-based systems,^31,44 and in this work the larger X = Q set was applied for all DFT calculations.⁴⁵

The dissociation energy of CO from a periodic Hf(001) surface of hexagonal Hf was determined using Quantum-ESPRESSO (ESI† Fig. S1).⁴⁶ For all calculations, a norm-conserving pseudopotential⁴⁷ with a kinetic energy cutoff of 70 Ry was employed. The Hf surface was approximated as a 3 × 3 × 1.5 slab and the CO molecule was placed vertically on top of an Hf atom. The positions of the Hf atoms in the top layer and the vertical distance of CO were optimized using the BLYP functional and a 3 × 3 × 1 k-point mesh. Subsequently, the dissociation energy was calculated with BLYP, PBE, BP86, B3LYP, and PBE0 through single-point calculations on the BLYP-optimized geometry at the Γ point of both the slab with CO present and the slab without CO along with an isolated CO molecule.

III. Results and discussion

III A. Ab initio analysis of HfCO

To understand bonding in HfCO, we start by analyzing the electronic structure of the individual Hf atom and CO molecule. The ³F ground electronic state of Hf has a valence 5d²6s² electron configuration. Simple electron promotions within the 5d shell produce the first two excited electronic states, which are experimentally only 16–26 kcal mol⁻¹ higher in energy (i.e., ³P at 15.8–25.7 kcal mol⁻¹, ¹D at 16.1 kcal mol⁻¹).⁴⁸ As a result of the low-lying electronic states of the isolated atom, Hf is expected to form complexes with a variety of chemical bonding configurations. We investigate the reaction between the aforementioned electronic states of Hf with the ground state of the CO ligand (X¹Σ⁺). The first excited state of CO (a³Π) lies well separated from the ground state (by 139.2 kcal mol⁻¹) and hence its interaction with Hf was not pursued in this work.⁴⁸ The combination of Hf(³F) + CO(X¹Σ⁺), Hf(³P) + CO(X¹Σ⁺), Hf(¹D) + CO(X¹Σ⁺), and Hf(⁵F) + CO(X¹Σ⁺) produces 7 triplet, 3 triplet, 5 singlet, and 7 quintet spin electronic states of the HfCO molecule, respectively. In total, we have studied 10 triplet, 5 singlet, and 7 quintet spin PECs of HfCO at the MRCI level of theory to identify its low-lying electronic states (Fig. 1). Note that all quintet states except b⁵Δ are high in energy and are not shown in Fig. 1.

Fig. 1 Full MRCI PECs of HfCO with respect to the Hf⋯C distance [r(Hf⋯C)] in Å. The CO distance is kept fixed to the experimental bond distance of CO (X¹Σ⁺) 1.128 Å. Only the b⁵Δ state originating from the interaction of Hf(⁵F) + CO(X¹Σ⁺) is depicted. Relative energies are referenced to the total energy of the Hf(³F) + CO(X¹Σ⁺) fragments, which are set to 0 kcal mol⁻¹. The PECs with triplet, singlet, and quintet spins are shown in solid, dotted, and dashed lines.

Analyzing the results of the MRCI potential energy scan, we observe that the Hf(³F) + CO(X¹Σ⁺) fragment produces three attractive states, X³Σ⁻, A³Φ, B³Π, and one strongly repulsive ³Δ state. X³Σ⁻, A³Φ, and B³Π are the ground state and the third and fourth lowest-lying electronic states of the molecule, respectively. The first excited state (i.e., a¹Δ) dissociates to the Hf(¹D) + CO(X¹Σ⁺) asymptote. We analyzed the electron configurations of the electronic states and the associated molecular orbitals (Table 1 and Fig. 2). These electron configurations highlight the σ-dative bonding between Hf and CO. Specifically, the dative interaction between these fragments is described by the doubly occupied 1σ orbital, which arises from 5d_z²(Hf) + σ(CO) hybridization (Fig. 2 and ESI† Fig. S2). The occupied 6s atomic orbital of Hf polarizes away from CO to facilitate an efficient σ-dative attack (2σ orbital of Fig. 2). Other than the σ-dative bond, a strong metal-to-ligand π back-donation (d_xz/d_yz of Hf to π_x*/π_y* of CO) was also observed in all electronic states. This π back-donation is evident from the 1π_x and 1π_y molecular orbitals of HfCO, which are occupied in all cases (Fig. 2). Specifically, according to natural bonding orbital analysis, 24% of the two electrons is back-donated from Hf to the two π* orbitals of CO in the ground state. On the other hand, only A³Φ, B³Π, b⁵Δ, and d¹Π excited states occupy non-bonding 1d_xy (1δ_xy) and 1d_x²–y² (1δ_x²–y²) orbitals (Table 1 and Fig. 2).

Table 1 Dominant electronic configurations at the equilibrium distance of the seven lowest-energy electronic states of HfCO as determined from MR calculations

State	% Contribution	Coeff.	1σ	2σ	1πx	1πy	1dxy	1dx^2–y^2
X^3Σ^−	90	0.95	2	2	α	α	0	0
a^1Δ	40	−0.63	2	2	0	2	0	0
	40	0.63	2	2	2	0	0	0
A^3Φ	45	0.67	2	2	α	0	0	α
	45	0.67	2	2	0	α	α	0
B^3Π	44	0.66	2	2	α	0	0	α
	44	−0.66	2	2	0	α	α	0
b^5Δ	98	0.99	2	α	α	α	α	0
c^1Σ^+	38	0.62	2	2	0	2	0	0
	38	0.62	2	2	2	0	0	0
d^1Π	16	0.40	2	2	α	0	0	β
	16	−0.40	2	2	β	0	0	α
	16	−0.40	2	2	0	α	β	0
	16	0.40	2	2	0	β	α	0

---

Fig. 2 Select molecular orbitals of HfCO. Hf, C, and O are shown in yellow, gray, and red colors, respectively. An isovalue of 0.02 e Å⁻³ was used for all contours. Red and blue colors correspond to positive and negative phases of each orbital, respectively.

We next analyzed the relationship between the electron configuration and bond order expected in these complexes. The single-reference ground state of HfCO (X³Σ⁻) has a 1σ²2σ²1π_x¹1π_y¹ (= 1σ²2σ²1π²) dominant electron configuration. Based on its 1σ²1π² electron population, a double bond (i.e., Hf=C) can be expected for the ground state (Table 1 and Fig. 2). Indeed, its calculated effective bond order based on the CASSCF weights is ∼1.8, which agrees well with our qualitative expectation. The first excited state of HfCO (a¹Δ) has the same electronic configuration as the ground state but with a different spin configuration. Indeed, the coupling of two π electrons in symmetrically equivalent 1π_x and 1π_y orbitals produces this multi-reference singlet electronic state. The destabilization of a¹Δ compared to X³Σ⁻ is rationalized by Hund's rule. The highly excited c¹Σ⁺ state carries the same multi-reference electronic configuration as a¹Δ but is likely higher in energy due to potential mixing with other excited states. The heavily multi-reference A³Φ, B³Π, and d¹Π states are generated by promoting one electron from a 1π orbital of X³Σ⁻ to 1δ orbitals (1d_xy and 1d_x²–y²) with various compositions (Table 1). The only single-reference excited state of HfCO is the b⁵Δ state that has the unique 1σ²2σ¹1π²1δ¹ configuration. Similar to the ground state, we can expect a bond order of two for b⁵Δ, which in fact is almost identical to its effective bond order, 1.9, as could be expected based on its 98% single-reference nature.

Because the dominant single-reference wavefunctions of X³Σ⁻ and b⁵Δ can be constructed at the Hartree–Fock level, we performed highly accurate CCSD(T) calculations on them with triple-ζ, quadruple-ζ, and quintuple-ζ correlation consistent basis sets to quantify basis set sensitivity. For single-reference states, we expect CCSD(T) to be more accurate than MRCI. These CCSD(T) results are reported in Table 2 with the MRCI(+Q) values. MRCI underestimates the D_e of HfCO by 6.6 kcal mol⁻¹ compared to the CCSD(T) result (23.00 vs. 29.65 kcal mol⁻¹ with the AQZ basis set). However, the discrepency between MRCI+Q vs. CCSD(T) is minor at only 2.7 kcal mol⁻¹. In general, MRCI+Q T_e values are slightly higher than MRCI except in two cases (i.e., a¹Δ and c¹Σ⁺), but all T_e differences are 4 kcal mol⁻¹ or less. Similar to D_e, the MRCI+Q T_e values agree well with CCSD(T) where the comparison can be made. After inclusion of spin–orbit coupling at the MRCI level for the MRCI optimized ground state using the Breit-Pauli Hamiltonian, the dissociation energy of HfCO is lowered significantly. Specifically, the D_e of the 0⁺ ground state of the ³Σ⁻ is 13.83 kcal mol⁻¹. This D_e is ∼9 kcal mol⁻¹ lower compared to the spin–orbit untreated ground state (Table 2). The zero-point energy (ZPE) of the ground state of HfCO is 4.21 kcal mol⁻¹ and thus the ZPE-corrected D_e at the MRCI level is 9.62 kcal mol⁻¹. This somewhat smaller D_e might be the reason for its lack of extensive experimental analysis. Note that there is a 108 cm⁻¹ second order spin–orbit splitting for the ground ³Σ⁻ state . The ¹Δ₂ spin–orbit sublevel of the first excited state, a¹Δ, lies 8.77 kcal mol⁻¹ above the , which is only 1.19 kcal mol⁻¹ different from the spin–orbit untreated first excitation energy of HfCO (9.96 kcal mol⁻¹, Table 2).

Table 2 Dissociation energy with respect to ground state fragments D_e (kcal mol⁻¹), equilibrium bond lengths r_e (Å), excitation energy T_e (kcal mol⁻¹), and harmonic vibrational frequencies ω_e (cm⁻¹), for the lowest seven electronic states of ¹⁷⁸Hf¹²C¹⁶O

State	Method^a	D e	r e		T e	ω e
			Hf−C	C−O
a Davidson-corrected MRCI is given as MRCI+Q. All multi-reference calculations are performed with the AQZ basis set (cc-pVQZ-PP/Hf, aug-cc-pVQZ/C,O). The CCSD(T) with TZ (cc-pVTZ-PP/Hf, cc-pVTZ/C,O), ATZ (cc-pVTZ-PP/Hf, aug-cc-pVTZ/C,O), AQZ (cc-pVQZ-PP/Hf, aug-cc-pVQZ/C,O), and A5Z (cc-pV5Z-PP/Hf, aug-cc-pV5Z/C,O) basis sets are performed for only single-reference X^3Σ^− and b^5Δ electronic states. b CCSD(T)-A5Z calculations carried out using the CCSD(T)-AQZ geometry.
X^3Σ^−	CCSD(T)-TZ	28.04	2.107	1.171	0	322, 322, 404, 1895
	CCSD(T)-ATZ	29.11	2.098	1.170	0
	CCSD(T)-AQZ	29.65	2.094	1.166	0
	CCSD(T)-A5Z ^b	29.96			0
	MRCI+Q	26.97	2.085	1.169	0
	MRCI	23.00	2.092	1.147	0
a^1Δ	MRCI+Q				8.86
	MRCI				9.96
A^3Φ	MRCI+Q				14.77
	MRCI				13.20
B^3Π	MRCI+Q				18.57
	MRCI				17.58
b^5Δ	CCSD(T)-TZ		2.094	1.166	22.36	319, 319, 406, 1877
	CCSD(T)-ATZ		2.111	1.175	22.09
	CCSD(T)-AQZ		2.107	1.171	22.10
	CCSD(T)-A5Z ^b				22.13
	MRCI+Q				21.84
	MRCI		2.102	1.152	17.84
c^1Σ^+	MRCI+Q				16.69
	MRCI				19.10
d^1Π	MRCI+Q				22.70
	MRCI				21.92

---

We carried out CCSD(T) geometry optimizations with the ATZ and AQZ basis sets, but to overcome the great computational cost of larger basis set calculations, we performed only a single-point quintuple-ζ basis set calculation at the geometry obtained with our quadruple-ζ basis set. We do not expect a significant structure variation from AQZ to A5Z basis sets, as the difference in bond length observed is less than 0.005 Å by improving from the ATZ to AQZ basis sets. Generally, it is expected that larger basis sets will predict shorter bond lengths,^49,50 and our results are consistent with that expectation (Table 2). Our MRCI-optimized bond distances and CCSD(T) values agree within 0.02 Å. In terms of experimental references values, the experimental C–O stretching frequency has been reported with co-deposition of laser-ablated hafnium with 0.1% CO in neon at 4–10 K.¹² This value¹² is 1869 cm⁻¹, which differs by 26 cm⁻¹ from our CCSD(T) value for X³Σ⁻ (Table 2).

We also evaluated the ionization energy (IE) of the single-reference ground state at the CCSD(T) level of theory. The removal of an electron from the doubly occupied 2σ molecular orbital of HfCO (X³Σ⁻) produces the HfCO⁺ (X²Σ⁻) cation. At the CCSD(T) level this IE is 6.974 eV, which is slightly higher compared to the first IE of the Hf atom (i.e., 6.825 eV) and could be due to the slightly ionic Hf^+0.53–[CO]^−0.53 charge distribution. No experimental reference value is available for comparison to this quantity, so we will use the CCSD(T) reference value for subsequent evaluation of the accuracy of various DFAs in DFT.

Recently, we observed that the D_e of the HfX (X = O, S, Se, Te, Po) series decreases moving from lighter HfO to heavier HfPo, which correlates with the binding elements’ electronegativity (i.e., O > S > Se > Te > Po).⁴⁴ Building upon our prior work, we also analyze the HfCX (X = O, S, Se, Te, Po) series. Because the metal-binding atom C is the same for all HfCX, it is somewhat difficult to make a prediction on the D_e trend for the HfCX series based on our knowledge of the HfX series. Indeed, we observe a reversed trend for HfCX, where D_e increases moving from lighter HfCO to heavier HfCPo. This opposite trend does follow expectations based on the dipole moment (μ) of the ligand (i.e., CO < CS < CSe < CTe < CPo). The relationship between D_e(Hf−CX) vs. μ(CX) is linear with R² = 0.991 (Fig. 3). To confirm the generality of this relationship, at the same level of theory we computed the D_e of the isovalent ZrCX and TiCX series and related it to the dipole moment of the CX species. Note that the ground states of ZrCO and TiCO are both ⁵Δ.^51,52 For both metals (M = Zr and Ti), the D_e increases moving from MCO to MCPo and importantly the near-linear D_evs. μ relationship is preserved (R² values 0.992 and 0.995, Fig. 3). Notably, the slope of the relationship decreases moving from HfCX to TiCX by about 25%. This can be rationalized by the fact that a decrease of D_e was observed moving from HfCX to TiCX, while the dipole moments are necessarily unchanged.

Fig. 3 Dissociation energy (D_e, in kcal mol⁻¹) of M–CX (M = Hf, Zr, Ti; X = O, S, Se, Te, Po) vs. dipole moment (μ, in D) of the free CX computed with CCSD(T) with the cc-pVTZ(-PP) basis for all atoms. The MCX species marked in the plot are MCO, MCS, MCSe, MCTe, and MCPo (left to right).

III B. DFT analysis of HfCO

Having identified single-reference electronic states of HfCO, we next used their CCSD(T) D_e, IE, and the ΔE_T-Q (i.e., the energy between X³Σ⁻ and b⁵Δ) values to evaluate DFT functional errors. Specifically, we assessed 23 functionals that span multiple rungs (GGAs, meta-GGAs, global GGA hybrids, meta-GGA hybrids, range-separated hybrids, and double hybrids) of “Jacob's ladder” using a previously developed workflow (see Computational details).³¹ We generally expect higher accuracy, albeit with higher computational cost, from the functionals that are in higher rungs of the ladder. However, our most expensive double hybrids have large D_e errors compared to less expensive functionals (Fig. 4). Specifically, in comparison to reference values from CCSD(T)-A5Z, the errors of double hybrids were large (24–52%). The global GGA hybrid B3LYP and the range-separated hybrid ωB97X are the best among the studied functionals with ∼4% and 1% errors in D_e, respectively. Compared to CCSD(T), all functionals predicted higher D_es except for the meta-GGA hybrid M06-2X, which underbinds HfCO by 2.3 kcal mol⁻¹.

Fig. 4 Dissociation energy (D_e, in kcal mol⁻¹) of HfCO calculated with different DFT functionals using the def2-QZVP basis set (blue dots). Each class of density functionals is separated with vertical purple lines and ordered by the rung on Jacob's ladder (left to right: GGA, meta-GGA, global GGA hybrid, meta-GGA hybrid, range-separated hybrid, and double hybrid). The horizontal green dashed line represents the CCSD(T)/A5Z dissociation energy of HfCO. Red dots are the dissociation energy between a CO molecule and a Hf surface.

We next extended our comparison to the evaluation of errors in HfCS and HfCSe (ESI† Fig. S3 and S4). The DFA errors are 7–24% higher for HfCO compared to the errors of HfCS and HfCSe for all functionals except for the B3LYP, SCAN0, and ωB97X functionals that have <∼7% errors (ESI† Tables S1–S3). On the whole, the dissociation of all three systems is predicted well by the meta-GGA hybrid M06-2X with less than 8% error. This observation is in line with our past DFT analysis of HfB, which had almost identical D_e predictions for M06-2X and CCSD(T).⁴⁴ Although it should be noted in contrast to HfCX (X³Σ⁻) and HfB (X⁴Σ⁻), which have unpaired valence electrons, M06-2X was previously observed to perform the worst at predicting D_e for the closed-shell ground state of HfO (X¹Σ⁺).⁴⁴ For both HfCS and HfCSe, the two best-performing functionals of HfCO (B3LYP and ωB97X) also provided small errors (∼5–8%), reinforcing our suggestion of the use of these functionals. Notably, all the studied functionals preserved the linear D_evs. μ trend that we observed earlier (Fig. 3 and ESI† Fig. S5). At first the agreement of CCSD(T) and B3LYP on HfCX dissociation is surprising. For example, according to the work by Wilson et al., B3LYP underestimated dissociation of CO from coordination complexes with first row transition metals, including Fe(CO)₅ and Cr(CO)₆, and this trend was reversed in FeCO, where B3LYP overestimated dissociation.⁵³ Nevertheless, it was suggested B3LYP predicted a higher error for the dissociation with Fe due to stronger multireference character of this system, whereas the HfCO molecule studied here is single reference.⁵³ In our previous work we observed B3LYP errors of 7 and 35.4% for dissociations of single-reference ionic HfO (closed-shell X¹Σ⁺) and HfB (X⁴Σ⁻) systems, respectively.⁴⁴ Note that it is rather difficult to come to a firm conclusion of which family or functional is best for all Hf-containing complexes especially if the degree of metal–organic bonding is quite different in each case.

We also evaluated the DFT errors associated with the ΔE_T-Q and IEs of HfCO (Fig. 5 and 6 and ESI† Tables S4, S5). Among all the functionals considered, for these quantities, the double hybrids (DHs), with the exception of PBE0-DH, yield the closest results to CCSD(T). This result differs from the D_e analysis where DHs performed poorly. Specifically, the DHs display values within 0–3 kcal mol⁻¹ of the CCSD(T) reference for ΔE_T-Q and 0–0.11 eV for IE. Surprisingly, the second-best performing family is GGA. Despite the well-known underestimation of ΔE_T-Q by GGA,⁵⁴ the average deviation from the CCSD(T) value is smaller than the other functional families. Overall, several functionals performed well for predicting the ΔE_T-Q and IE of HfCO (Fig. 5 and 6). Importantly, the best-performing functionals for HfCO D_e, the GGA hybrid B3LYP and the range-separated hybrid ωB97X, also had small errors for the IE (1.7 and 5.1%, respectively). The same cannot be said for ΔE_T-Q, where these functionals had 24.4 and 42.3% errors, respectively. Thus, our overall recommendation remains for the use of B3LYP and ωB97X based on D_e and IE despite some caveats regarding the spin state ordering. Nevertheless, it is noteworthy that all functionals correctly identify the ground state, despite in some cases underestimating (e.g., B3LYP) or overestimating (e.g., ωB97X) the quantitative value of the ΔE_T-Q gap (Fig. 5).

Fig. 5 Triplet–quintet adiabatic energy gap (ΔE_T-Q, in kcal mol⁻¹) of HfCO calculated with different DFT functionals using the def2-QZVP basis set (blue dots). Each class of density functionals is separated with vertical purple lines and ordered by the rung on Jacob's ladder (left to right: GGA, meta-GGA, global GGA hybrid, meta-GGA hybrid, range-separated hybrid, and double hybrid). The horizontal green dashed line represents the CCSD(T)/A5Z ΔE_T-Q of HfCO.

Fig. 6 Adiabatic ionization energy (IE, in eV) of HfCO calculated with different DFT functionals using the def2-QZVP basis set (blue dots). Each class of density functionals is separated with vertical purple lines and ordered by the rung on Jacob's ladder (left to right: GGA, meta-GGA, global GGA hybrid, meta-GGA hybrid, range-separated hybrid, and double hybrid). The horizontal green dashed line represents the CCSD(T)/TZ IE of HfCO.

Finally, although the complexes studied in this work are relevant as models for single-atom catalysis, it is also useful to note the extent to which our observations of DFT performance on molecules apply to the solid state, such as in small-molecule adsorption on surfaces relevant to catalysis. In the comparison of the interaction between a CO molecule and the Hf atom versus the Hf surface, functionals within the GGA family yield a D_e value approximately 8 kcal mol⁻¹ smaller for binding the Hf surface compared to the isolated atom. On the other hand, the hybrid functionals B3LYP and PBE0 exhibit very similar D_e values for the molecule and surface, deviating by only 0.67 and 1.11 kcal mol⁻¹, respectively (Fig. 4). If the Hf–C interaction is as short-range as we expect based on the prior success of CO-focused corrections for CO adsorption on Pt(111),^55–57 then the molecular and solid-state D_e values should be comparable. This comparison would lead us to suggest the global hybrid B3LYP as preferable to GGAs to reproduce higher-level CCSD(T) values in molecules and provide a consistent result between the surface and the molecule. Nevertheless, a more comprehensive analysis is necessary in future work, including through presently cost-prohibitive CCSD(T) calculations on the surface, to understand the deviation between molecules and solids with GGAs.

IV. Conclusions

In conclusion, we studied potential energy curves, dissociation energies, excitation energies, harmonic vibrational frequencies, and chemical bonding patterns of several low-lying electronic states HfCO. We found seven bound electronic states of HfCO with respect to the Hf(³F) + CO(X¹Σ⁺) dissociation. Among these, all but the X³Σ⁻ and b⁵Δ are multi-reference in nature. We adopted larger augmented quadruple- and quintuple-ζ basis sets to obtain highly accurate CCSD(T) results for the two single-reference electronic states. Our analysis on MCX (M = Hf, Zr, Ti and X = O, S, Se, Te, Po) demonstrated a linear relationship between the dissociation energies of these complexes and the dipole moments of the CX species, where D_e increases moving from lighter MCO to heavier MCPo.

We used CCSD(T) data to benchmark DFT functionals. Specifically, we used 23 functionals that span multiple rungs of the “Jacob's ladder”. The global GGA hybrid B3LYP and range-separated hybrid ωB97X predicted the HfCO D_e and IE with small errors (less than 8%). Importantly, these two functionals displayed higher errors for predicting ΔE_T-Q hence they might not be ideal for predicting excited-state properties. We further studied the interaction of a CO molecule with a Hf surface with BP86, BLYP, PBE, B3LYP, and PBE0 functionals and observed that the GGAs consistently predicted smaller chemisorption energies (by ∼8 kcal mol⁻¹) compared to their gas-phase D_es but the global hybrids yielded more consistent results. Taken together with our analysis of benchmarks on the molecular systems, this encourages us to suggest B3LYP as a promising functional for both molecular and solid-state models of Hf–C interactions.

Data availability

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

Conflicts of interest

The authors declare no competing financial interest.

Acknowledgements

This work was supported by the United States Department of Energy under grant number DE-NA0003965. C. D. was partially supported by a seed fellowship from the Molecular Sciences Software Institute under NSF grant OAC-1547580. H. J. K. holds a Sloan Fellowship in Chemistry, which supported this work. We would like to thank Shahriar N. Khan for his help with the spin–orbit calculations. The authors thank Adam H. Steeves for providing a critical reading of the manuscript.

References

1. B. M. Gardner, C. C. J. Seechurn and T. J. Colacot, Organomet. Chem. Ind., 2020, 1–22.
2. A. Monney and M. Albrecht, Coord. Chem. Rev., 2013, 257, 2420–2433.
3. B. Cornils, W. A. Herrmann, M. Beller and R. Paciello, Applied homogeneous catalysis with organometallic compounds: a comprehensive handbook in four volumes, John Wiley & Sons, 2017.
4. D. Steinborn, Fundamentals of organometallic catalysis, John Wiley & Sons, 2011.
5. R. Ugo and R. Psaro, J. Mol. Catal., 1983, 20, 53–79.
6. R. R. Ford, in Advance Catalyst, ed. D. D. Eley, H. Pines and P. B. Weisz, Academic Press, 1970, vol. 21, pp. 51–150.
7. J. J. Turner, M. W. George, M. Poliakoff and R. N. Perutz, Chem. Soc. Rev., 2022, 51, 5300–5329.
8. C. Cesari, J.-H. Shon, S. Zacchini and L. A. Berben, Chem. Soc. Rev., 2021, 50, 9503–9539.
9. I. Ciabatti, C. Femoni, M. C. Iapalucci, S. Ruggieri and S. Zacchini, Coord. Chem. Rev., 2018, 355, 27–38.
10. F. A. Cotton, G. Wilkinson, C. A. Murillo and M. Bochmann, Advanced inorganic chemistry, John Wiley and Sons, Inc., 1999.
11. C. R. Groom, I. J. Bruno, M. P. Lightfoot and S. C. Ward, Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater., 2016, 72, 171–179.
12. M. Zhou and L. Andrews, J. Am. Chem. Soc., 2000, 122, 1531–1539.
13. A. Brathwaite and M. Duncan, J. Phys. Chem. A, 2013, 117, 11695–11703.
14. G. Wang, J. Zhao, H. S. Hu, J. Li and M. Zhou, Angew. Chem., 2021, 133, 9420–9424.
15. C. Chi, J. Q. Wang, H. Qu, W. L. Li, L. Meng, M. Luo, J. Li and M. Zhou, Angew. Chem., Int. Ed., 2017, 56, 6932–6936.
16. I. R. Ariyarathna and E. Miliordos, Phys. Chem. Chem. Phys., 2019, 21, 24469–24477.
17. The National Institute of Occupational Safety and Health, https://www.cdc.gov/niosh/idlh/7440586.html, (accessed 08/23/ 2023).
18. U.S. Department of Labor Occupational Safety & Health Administration, https://web.archive.org/web/20080313003040/https://www.osha.gov/SLTC/healthguidelines/hafnium/index.html, (accessed 08/23/ 2023).
19. X. Jin, Y. Zhou, G. Wang and M. Zhou, J. Phys. Chem. A, 2023, 127, 4483–4491.
20. X. Jin, G. Wang and M. Zhou, Phys. Chem. Chem. Phys., 2023, 25, 7697–7703.
21. M. Gao, Q. Zhao, H. Yu, M. Fu and Q. Li, Molecules, 2022, 27, 2885.
22. X. Jin, Y. Bai, Y. Zhou, G. Wang, L. Zhao, M. Zhou and G. Frenking, Angew. Chem., Int. Ed., 2021, 60, 13865–13870.
23. G. Frenking, I. Fernández, N. Holzmann, S. Pan, I. Krossing and M. Zhou, JACS Au, 2021, 1, 623–645.
24. H. Lischka, D. Nachtigallova, A. J. Aquino, P. G. Szalay, F. Plasser, F. B. Machado and M. Barbatti, Chem. Rev., 2018, 118, 7293–7361.
25. K. D. Vogiatzis, M. V. Polynski, J. K. Kirkland, J. Townsend, A. Hashemi, C. Liu and E. A. Pidko, Chem. Rev., 2019, 119, 2453–2523.
26. K. Raghavachari, G. W. Trucks, J. A. Pople and M. Head-Gordon, Chem. Phys. Lett., 1989, 157, 479–483.
27. J. K. Nørskov and T. Bligaard, Angew. Chem., Int. Ed., 2013, 52, 776–777.
28. W. Kohn, Rev. Mod. Phys., 1999, 71, 1253–1266.
29. P. Geerlings, F. De Proft and W. Langenaeker, Chem. Rev., 2003, 103, 1793–1874.
30. L. R. Maurer, M. Bursch, S. Grimme and A. Hansen, J. Chem. Theory Comput., 2021, 17, 6134–6151.
31. C. Duan, S. Chen, M. G. Taylor, F. Liu and H. J. Kulik, Chem. Sci., 2021, 12, 13021–13036.
32. P. J. Knowles and H.-J. Werner, Chem. Phys. Lett., 1988, 145, 514–522.
33. K. Shamasundar, G. Knizia and H.-J. Werner, J. Chem. Phys., 2011, 135.
34. H. J. Werner and P. J. Knowles, J. Chem. Phys., 1988, 89, 5803–5814.
35. H. Werner, P. Knowles, G. Knizia, F. Manby, M. Schütz, P. Celani, W. Györffy, D. Kats, T. Korona and R. Lindh, MOLPRO, version 2015.1, a package of ab initio programs, University of Cardiff Chemistry Consultants (UC3), Cardiff, Wales, UK, 2015.
36. K. Huber and G. Herzberg, Molecular Spectra and Molecular Structure, Van Nostrand Reinhold Company, New York, USA, 1979.
37. S. R. Langhoff and E. R. Davidson, Int. J. Quantum Chem., 1974, 8, 61–72.
38. D. Figgen, K. A. Peterson, M. Dolg and H. Stoll, J. Chem. Phys., 2009, 130.
39. R. A. Kendall, T. H. Dunning Jr and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796–6806.
40. D. E. Woon and T. H. Dunning Jr, J. Chem. Phys., 1993, 98, 1358–1371.
41. A. K. Wilson, D. E. Woon, K. A. Peterson and T. H. Dunning Jr, J. Chem. Phys., 1999, 110, 7667–7676.
42. K. A. Peterson, D. Figgen, E. Goll, H. Stoll and M. Dolg, J. Chem. Phys., 2003, 119, 11113–11123.
43. J. M. Turney, A. C. Simmonett, R. M. Parrish, E. G. Hohenstein, F. A. Evangelista, J. T. Fermann, B. J. Mintz, L. A. Burns, J. J. Wilke and M. L. Abrams, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 556–565.
44. I. R. Ariyarathna, C. Duan and H. J. Kulik, J. Chem. Phys., 2022, 156, 184113.
45. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
46. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
47. D. R. Hamann, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 085117.
48. A. Kramida, Y. Ralchenko and J. Reader, https://physics.nist.gov/asd, 2016.
49. I. R. Ariyarathna and E. Miliordos, Phys. Chem. Chem. Phys., 2018, 20, 12278–12287.
50. E. E. Claveau and E. Miliordos, Phys. Chem. Chem. Phys., 2021, 23, 21172–21182.
51. H. Tan, M. Liao and K. Balasubramanian, J. Phys. Chem. A, 1998, 102, 1602–1607.
52. C. Koukounas, S. Kardahakis and A. Mavridis, J. Chem. Phys., 2005, 123, 074327.
53. N. J. DeYonker, T. G. Williams, A. E. Imel, T. R. Cundari and A. K. Wilson, J. Chem. Phys., 2009, 131.
54. A. Fouqueau, S. Mer, M. E. Casidaa, L. M. L. Daku, A. Hauser, T. Mineva and F. Neese, J. Chem. Phys., 2004, 120.
55. S. E. Mason, I. Grinberg and A. M. Rappe, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 161401.
56. I. Grinberg, Y. Yourdshahyan and A. M. Rappe, J. Chem. Phys., 2002, 117, 2264–2270.
57. G. Kresse, A. Gil and P. Sautet, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 68, 073401.

---

Footnote

† Electronic supplementary information (ESI) available: vbL diagrams for HfCO; dissociation energy of HfCO with DFT; dissociation energy of HfCS with DFT; dissociation energy of HfCSe with DFT; dissociation energy of HfCX (X = O, S, Se) vs. dipole moment with DFT. See DOI: https://doi.org/10.1039/d3cp03550f

---