1 Introduction
Molybdenum (Mo) is a versatile transition metal discovered by Karl Scheele in 1778 and is essential for organisms.[1] It forms the catalytic center of numerous enzymes, for example, the aldehyde oxidoreductase (AOR)[2] used for formaldehyde oxidation and the DMSO reductase used for dimethyl sulfoxide (DMSO)[3] reduction. As a transition metal, Mo is positioned between the early and late transition metals in the periodic table of elements, typically exhibiting strong Lewis acidity and demonstrating high reactivity toward various organics,[4] such as alkene metathesis,[5] ammonia oxidation,[6] epoxidation,[7] nitrogen reduction[8] and hydrogenation.[9] Given its unique reactivity profile, Mo has recently emerged as a promising catalyst in transition-metal-catalyzed carbene transfer reactions, a powerful strategy for forming chemical bonds in organic synthesis. These in situ-formed transition-metal carbene intermediates are often generated via the metal-catalyzed decomposition of diazo substrates. Recently, Zhuo’s group[10] reported Mo-catalyzed deoxygenative cyclopropanation and carbonyl-carbonyl olefination reaction[11] (Schemes 1a and 1b). This strategy has broadened the scope of carbene precursors to carbonyl compounds. From the mechanistic point of view, two successive carbonyl oxidative addition steps are involved in Mo-carbene formation. The first one starts from the carbonyl σ-coordination of o-quinone[12] species to the Mo(0) center, followed by the concerted oxidative addition of 1,2-dicarbonyl to the Mo(0) through a five-membered cyclic transition state TS-I (Scheme 1c). This process forms an Mo(II) complex with concomitant release of carbon monoxide. Subsequent oxidation occurs through the π-coordination of carbonyl group to Mo(II) center and the concerted oxidative addition of the carbonyl group to the Mo(II) through a three-membered cyclic transition state TS-II. This step cleaves the C=O double bond directly and generates the key Mo(IV) carbene intermediate. Therefore, these two oxidative addition processes were chosen for the following benchmark study of DFT functionals in this work.
Density functional theory (DFT) calculations enable one to gain more insights into the reaction mechanism.[13] However, the complexity of the electronic structure of Mo poses a challenge to the selection of appropriate DFT methods.[14] Therefore, a comprehensive evaluation of different DFT methods is critical to determine which methods can provide reliable and consistent results for Mo-catalyzed reactions. The core distinction among different exchange- correlation functionals[15] lies in their treatment of electron density and its higher-order derivatives, as well as the degree of incorporating exact exchange energy.[16] For example, the Generalized Gradient Approximation (GGA)[17] accounts for the gradient of electron density (i.e., spatial variations in density), addressing the limitation of the Local Density Approximation (LDA),[18] which relies solely on local density to better describe inhomogeneous electron gases (e.g., electron delocalization effects in molecules). The meta-Generalized Gradient Approximation (meta-GGA)[19] builds upon GGA by further incorporating kinetic energy density, enabling a more precise description of strong electron correlation and nonlocal effects. In this work, fourteen functionals, including density functionals of GGA, meta-GGA, hybrid GGA, and hybrid meta-GGA, were selected to evaluate their capabilities in structural optimization and energy calculation of TS-I and TS-II (Scheme 1c). This benchmark study will provide practical guidance for future computational studies of Mo catalysis.
2 Computational details
All DFT calculations in this study were performed using Gaussian 16[20] or ORCA 4.2.1[21] software packages. The bond lengths from the X-ray crystal structures were used as references in the benchmark for structure optimization. Commonly used density functionals, such as pure GGA functional (BP86-D3), meta-GGA functional (M06-L, M11-L, MN15-L), Hybrid GGA functional (B3LYP-D3, PBE0-D3, ωB97X-D, X3LYP), and hybrid meta-GGA functional (TPSSh, M06, M06-2X, M11, MN15, M06-HF) were chosen for structure optimization.
Density functional theory (DFT) functionals exhibit diverse performance depending on their design and parametrization.[22] GGA functionals like BP86-D3[23] offer computational efficiency but suffer from inaccuracies in non-covalent interactions, charge-transfer excitations, and self-interaction errors. Hybrid functionals such as B3LYP- D3[24] and PBE0-D3[25] address these limitations by incorporating exact Hartree-Fock (HF) exchange[26] and Grimme’s D3 dispersion corrections,[27] improving descriptions of dispersion forces and thermochemistry. However, their reliance on fixed HF ratios restricts adaptability across diverse systems. Range-separated hybrids like ωB97X-D[28] further enhance long-range interaction modeling but increase computational cost. Meta-GGA functionals (M06- L,[29] M11-L,[30] MN15-L[31]) and hybrid-meta variants (M06,[32] M06-2X,[32] M11,[33] MN15[34]) balance accuracy and efficiency. M06-L excels for transition metals but underestimates non-covalent forces, while M06-2X, with higher HF exchange, improves main-group thermochemistry but struggles with delocalization errors. MN15, a modern hybrid-meta functional, demonstrates broad applicability but retains sensitivity to specific systems. Recent trends prioritize dispersion-corrected hybrids (e.g., B3LYP-D3) and system-specific parametrization (MN15, ωB97X-D) to balance accuracy and scalability. Functionals like MN15-L and X3LYP[36] aim for transferability across reaction types, though validation remains ongoing.
High-level ab initio calculations are often used as a reference for DFT calibration. For closed-shell systems without severe multi-configuration characteristics, such as many closed-shell systems containing 4d and 5d transition metals, DLPNO-CCSD(T)[37] (Domain-Based Local Pair Natural Orbital CCSD(T)) with the perturbative treatment of triplet excitations is the “gold standard” of modern quantum chemistry and is often used in the benchmark study. Compared with the CCSD(T) method, the local natural orbital method reduces computational complexity and improves efficiency, which is suitable for this larger-scale molecular system. Therefore, this method and the def2-TZVPP[38] basis set were applied in toluene (SMD model,[39] ε=2.3741). The RI approximation for both Coulomb- and exchange-type integrals (RI-JK) and associated auxiliary basis sets (def2/JK and def2-TZVPP/C) are employed.[40] The resulting single-point energies were used as energy references. The Gibbs free energies discussed below are obtained by adding the thermal free energy correction of the most stable conformations optimized in the gas phase using different generalized functionals combined with the def2- SVP basis set plus the electronic energies obtained in solvent using different generalized functional combinations of the def2-TZVP basis set. The 3D structures of the optimized intermediates or transition states were demonstrated using CYLView software.[41]
3 Results and discussion
3.1 Benchmark for geometry optimization of Mo(0) catalysts
The benchmark study of DFT methods for geometry optimization was initially carried out. Ke’s group[42] developed a new type of strong hydride-donating NHC-Mo(0) catalyst for Mo-catalyzed direct N-alkylation of anilines and nitroarenes. The X-ray crystal structures of NHC-Mo(0) complexes 1 and 2 reported in Ke’s work were used as the reference for the benchmark study of geometry optimization (Figure 1a). These two structures both exhibit an octahedral configuration with two carbon monoxide molecules occupying opposite vertex positions. The bond lengths of seven critical Mo—C and Mo—N bonds are shown in Figure 1a. Fourteen different functionals, including BP86-D3(BJ), B3LYP-D3, PBE0-D3, TPSSh, X3LYP, ωB97X-D, M06- L, M11-L, M06, M06-2X, M06-HF, M11, MN15, and MN15-L were employed to optimize the structures of 1 and 2. The comparison of these seven bond lengths between the DFT-optimized structure and the X-ray crystal structure was performed to evaluate the performance of each DFT method in geometry optimization. The mean absolute error (MAE) and the mean squared error (MSE) were calculated using Eq. 1 and Eq. 2, respectively.
$\mathrm{MAE}=\frac{1}{n} \sum_{i=1,2,3, \ldots}^{n}\left|d_{i}-\hat{d}_{i}\right| \quad(n=7)$
$\mathrm{MSE}=\frac{1}{n} \sum_{i=1,2,3, \ldots}^{n}\left|d_{i}-\hat{d}_{i}\right|^{2} \quad(n=7)$
where di denotes the bond length of the Mo—C and Mo—N bonds in DFT-optimized structure, and $\hat{d}_{i}$ denotes the corresponding bond length in X-ray crystal structure.
The benchmark results of geometry optimization are shown in Figure 1b. These results suggest that the differences between the optimized structures and the crystal structures for different functionals are relatively small (MAE<0.0300 Å, MSE<0.0020 Å2). Among all the functionals used for geometry optimization, TPSSh performs the best as it has the lowest MAE and MSE (MAE=1.00×10－2 Å, MSE=1.11×10－4 Å2). In addition, the performance of B3LYP-D3(BJ) and ωB97X-D is also com- mendable. The calculated MAE values are 1.38×10－2 and 1.09×10－2 Å, respectively. The MSE values are 2.61× 10－4 and 4.86×10－4 Å2, respectively. In contrast, M06-2X performs the worst due to the largest MAE and MSE, thus it is not recommended for geometry optimization in molybdenum catalysis. Therefore, the functionals that ranked within the top three, i.e. B3LYP-D3(BJ), ωB97X-D, and TPSSh, were used to optimize the structures of molybdenum-catalyzed carbonyl oxidative addition transition states in the subsequent studies.
3.2 Benchmark for energy calculation of carbonyl oxidative addition transition states
The carbonyl oxidative addition is a critical step for molybdenum catalysis, which not only triggers the activation of the catalyst to start the catalytic cycle, but also initiates the transformation of substrates. Thus, it is necessary to ascertain a suitable density functional for calculating the activation free energy of carbonyl oxidative addition. Based on this idea, three representative carbonyl oxidative addition reactions under molybdenum catalysis are investigated, including the formation of Mo(II) through the oxidative addition of o-quinone to Mo(0) (Figure 2a), the formation of Mo(IV) through the oxidative addition of 1,2-dicarbonyl compound 3 to Mo(II) (Figure 2b), and the formation of oxo-Mo(IV)-carbene through the oxidative addition of 1,2-dicarbonyl compounds 3 to Mo(II) (Figure 2c). The oxidative addition transition states were respectively optimized for each step, including the five-membered cyclic transition states TS-1, TS-2, and the three-membered cyclic transition state TS-3.
Figure 2 DFT study of three representative carbonyl oxidative addition steps in Mo catalysisThe activation free energies were calculated using DLPNO-CCSD(T) based on the structures optimized by B3LYP-D3, ωB97X-D, and TPSSh, respectively. The 3D structures were optimized using B3LYP-D3(BJ). |
According to the benchmark results in Figure 1, three functionals (B3LYP-D3(BJ), ωB97X-D, and TPSSh) were used for geometry optimization in the gas phase with def2- SVP as the basis set. The single-point energies were calculated using DLPNO-CCSD(T) in toluene solvent with the SMD continuum solvation model based on the gas-phase optimized geometries. The basis set of def2-TZVPP was employed for single-point energy calculation. The activation free energies of carbonyl oxidative addition transition states obtained at the DLPNO-CCSD(T)//TPSSh, DLPNO- CCSD(T)//ωB97X-D, and DLPNO-CCSD(T)//B3LYP-D3 (BJ) levels of theory were respectively used as the standards in the subsequent benchmark studies. 14 functionals, including BP86-D3(BJ), B3LYP-D3(BJ), PBE0-D3(BJ), TPSSh, X3LYP, ωB97X-D, M06-L, M11-L, M06, M06- 2X, M06-HF, M11, MN15, and MN15-L, with a basis set of def2-TZVP were tested for single-point energy calculations. Based on structures optimized using three functionals (B3- LYP-D3(BJ), ωB97X-D, and TPSSh) in gas phase, three sets of energy data will be obtained. The corresponding MAE and MSE with respect to DLPNO-CCSD(T) results were calculated using Eqs. 3, 4:
$\mathrm{MAE}=\frac{1}{n} \sum_{i=1,2,3}^{n}\left|\Delta G_{i}^{\neq}-\Delta \hat{G}_{i}^{\neq}\right| \quad(n=3)$
$\mathrm{MAE}=\frac{1}{n} \sum_{i=1,2,3}^{n}\left|\Delta G_{i}^{\neq}-\Delta \hat{G}_{i}^{\neq}\right|^{2} \quad(n=3)$
where $\Delta G_{i}^{\neq}$ denotes the free energy calculated using 14 functionals based on the structures optimized using TPSSh, ωB97X-D and B3LYP-D3(BJ), and $\Delta G_{i}^{\neq}$ denotes the free energy calculated using DLPNO-CCSD(T) based on the structures optimized using TPSSh, ωB97X-D and B3LYP- D3(BJ).
These benchmark results are summarized in Figures 3~5, respectively. As shown in Figure 3, the MAE and MSE values of 14 functionals based on TPSSh optimized structures are presented. The comparison suggests that M11-L, MN15 and MN15-L have relatively lower MAE and MSE values (MAE=3.3, 3.1, and 3.2 kcal/mol, MSE=11.0, 14.2, and 20.3 (kcal/mol)2, respectively). On the contrary, the TPSSh functional, which performs well in geometry optimization, doesn’t afford good results. This contrast indicates that TPSSh functional is not an ideal option for theoretical study of molybdenum catalysis. Moreover, the BP86-D3(BJ) functional is also not recommended as the MAE and MSE reach up to 7.2 kcal/mol and 81.8 (kcal/mol)2, respectively, which may lead to huge deviations in activation free energy.
When ωB97X-D-optimized structures were employed in the energy benchmark (Figure 4), the overall performance of functionals M06, M11, and ωB97X-D ranks in the top three. The MAE values are 3.0, 2.5, and 2.5 kcal/mol, respectively. The MSE values are 12.2, 10.5, and 11.7 (kcal/ mol)2, respectively. In contrast, the performances of M11-L and BP86-D3(BJ) are worse due to higher MAE and MSE. According to the results shown in Figure 5, the combination of geometry optimization using the B3LYP-D3(BJ) functional with single-point energy calculations using M06, M06-L, or MN15-L delivers the best performance. The MAE values are 1.0, 1.6, and 0.8 kcal/mol, respectively. MSE values are 1.0, 2.5, and 0.9 (kcal/mol)2, respectively. These values all fall within a desired range and are lower than those in Figures 3 and 4. Therefore, B3LYP-D3(BJ) functional is highly recommended for geometry optimization in molybdenum catalysis, and M06, M06-L, or MN15- L functionals are recommended for single-point energy calculation. More specifically, our benchmark study suggests that the combination of MN15-L and B3LYP-D3(BJ) is the most reliable option for computational study of molybdenum catalysis.
Figure 4 Benchmark results of DFT functionals for energy calculation of ωB97X-D optimized carbonyl oxidative addition transition statesThe activation free energy calculated by DLPNO-CCSD(T) is used as the standard of energy. |
4 Conclusion
This study systematically evaluates the performances of 14 DFT functionals for geometry optimization and energy calculation in three types of Mo-catalyzed carbonyl oxidative addition reactions. Benchmark study against X-ray crystal structures reveal that B3LYP-D3(BJ), ωB97X-D, and TPSSh exhibit superior accuracy in structure optimization, achieving mean absolute errors (MAE)≤1.38×10－2 Å and mean squared errors (MSE)≤4.86×10－4 Å2. For free energy calculations, the single-point energy calculation using M06, M06-L, or MN15-L in solvent based on the B3LYP-D3(BJ)-optimized structure yields the lowest deviations (MAE≤1.6 kcal/mol, MSE≤2.5 (kcal/mol)2) relative to the DLPNO-CCSD(T) calculated data. Notably, hybrid-meta functionals such as M06 and M06-L, and meta functionals like MN15-L demonstrate robustness in modeling activation barriers, while pure GGA methods like BP86 show significant errors. These findings highlight the importance of functional selection tailored to specific tasks, which are structure optimization versus energy evaluation in Mo-catalyzed systems. The recommended workflow, i.e., geometry optimization using B3LYP-D3(BJ) with single- point energy calculation using MN15-L, balances the accuracy and computational efficiency, offering the practical guidance for mechanistic studies of molybdenum catalysis. This benchmark study underscores the necessity of functional selection and enhances the reliability of DFT in modeling complex organometallic reactions.
Acknowledgments The author would like to express his gratitude to the supercomputing system of the Supercomputing Center of Wuhan University.
Supporting Information Cartesian coordinates and energies. The Supporting Information is available free of charge via the Internet at http://sioc-journal.cn/.
(Lu, Y.)