# Accurate Calculation of Electron Paramagnetic Resonance Parameters for Molybdenum Compounds

Maria Drosou, Iris Wehrung, Dimitrios Pantazis, Maylis Orio

► To cite this version:

Maria Drosou, Iris Wehrung, Dimitrios Pantazis, Maylis Orio. Accurate Calculation of Electron Paramagnetic Resonance Parameters for Molybdenum Compounds. ChemPhysChem, 2025, 10.1002/cphc.202500317. [https://doi.org/10.1002/cphc.202518994](https://doi.org/10.1002/cphc.202518994)

HAL Id: hal-05318994

https://hal.science/hal-05318994v1

Submitted on 17 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Accurate Calculation of Electron Paramagnetic
Resonance Parameters for Molybdenum
Compounds

Maria Drosou,<sup>1,2*</sup> Iris Wehrung,<sup>3</sup> Dimitrios A. Pantazis,<sup>1*</sup> Maylis Orio<sup>3*</sup>

E-mail: drosou@kofo.mpg.de; dimitrios.pantazis@kofo.mpg.de; maylis.orio@univ-amu.fr

ORCIDs:

Maria Drosou: 0000-0002-4550-710X

Dimitrios A. Pantazis: 0000-0002-2146-9065

Maylis Orio: 0000-0002-9317-8005

<footnote n="1">¹ Max-Planck-Institut für Kohlenforschung, Kaiser-Wilhelm-Platz 1, 45470 Mülheim an der Ruhr, Germany.</footnote>

<footnote n="2">² Department of Chemistry, Quantum Chemistry, TU Darmstadt, Peter-Grünberg-Str. 4, 64287, Darmstadt, Germany.</footnote>

<footnote n="3">³ Aix-Marseille Univ, CNRS, Centrale Med, iSm2, Marseille, France.</footnote>
**Abstract**

Paramagnetic molybdenum compounds are of great interest in inorganic chemistry and metalloenzyme catalysis. Electron paramagnetic resonance (EPR) spectroscopies that determine hyperfine coupling constants (HFCs) and *g*-tensor values are essential for investigating the electronic structure of these compounds, but require support from accurate quantum chemical approaches. Here we present a database of Mo(V) complexes with well-defined structures and EPR parameters and investigate optimal quantum chemical protocols for $^{95}$Mo HFCs and *g*-values. It is shown that the SARC all-electron basis sets can be used without modification to produce converged results for HFCs and *g*-values with the exact-2-component (X2C) Hamiltonian. The dependence of EPR parameters on the functional is studied in detail. Double-hybrid functionals and global hybrids with high exact exchange are top performers for $^{95}$Mo HFCs, with PBE0-DH achieving the best agreement with experiment. Comparison of DFT-derived HFCs with values obtained by coupled cluster theory with the domain-based local pair natural orbital approach (DLPNO-CCSD) shows that DFT remains the method of choice for the present set of compounds. Smaller differentiation among functionals is observed for *g*-tensors, although PBE0-DH is still a top performer and can be recommended as the most reliable approach overall for describing both valence and core properties of Mo compounds.

**KEYWORDS:** Molybdenum, EPR, quantum chemistry, hyperfine couplings, *g*-tensors
**Introduction**

Molybdenum is the only known second-row transition metal essential for all forms of life. Despite a relatively limited set of coordination possibilities encountered in biology, the chemical versatility of the metal means that Mo enzymes are characterized by a wide range of reactivity, catalyzing diverse reactions of central importance in biogeochemical cycles and having important implications in the fields of environment, energy and human health. [1-4] Specific examples include clusters that take advantage of multimetallic cooperativity, such as the multinuclear Fe<sub>7</sub>Mo cluster of nitrogenase[5] and the binuclear MoCu center of CO dehydrogenase,[6] but Mo is also found independently in several mononuclear active sites often coordinated by a pyranopterin-dithiolene ligand (PDT). Molybdopterin enzymes are conventionally categorized into three families based on the coordination sphere of the Mo center, namely xanthine oxidase, sulfite oxidase, and dimethyl sulfoxide (DMSO) reductase families (Figure 1), named according to their archetypical member. Formate dehydrogenases,[7-9] to single out a prominent type of enzyme, belong to the DMSO reductase family and are of particular interest because they catalyze the reduction of CO<sub>2</sub>, a reaction of fundamental significance in the current context of climate crisis mitigation and solar fuels research.[10]

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Structures of the Mo coordination sphere for three families of molybdenum enzymes.</figcaption></figure>

A main challenge in the field is to understand the structural factors that control the reactivity and substrate specificity of Mo cofactors. Toward this end, it is crucial to reliably map structural features at the atomic level to the experimentally observed spectroscopic properties of these cofactors. This can be facilitated by the combined and complementary approaches of synthetic model chemistry, which provides well-defined biomimetic molecular compounds,[11-12] and by quantum chemistry, which provides the mapping between specific structural features and spectroscopic observables.[13-14] The latter can only succeed if reliable and accurate computational methodologies and protocols are available. This represents a continuing quest for applied quantum chemistry in the field of open-shell transition metal
ions. In the case of molybdenum, of particular interest is the oxidation state Mo(V) with a
formal 4d¹ electron configuration (total spin S = 1/2) that lends itself to probing by electron
paramagnetic resonance (EPR) spectroscopy.[¹², ¹⁵]

The two quantities of interest in this context are the molybdenum hyperfine coupling constants (HFCs, A) and the g tensors of the Mo(V) compounds. Both parameters have been targeted by quantum chemical studies in the past,[¹⁶⁻²⁵] with density functional theory (DFT) mostly used to study the sensitivity of computational predictions with respect to the functional, basis set, and the different frameworks to account for relativistic effects. Important observations have been made concerning the capabilities and limitations of representative approaches in capturing the essential physics, leading to practical suggestions such as the requirement for Hartree–Fock exchange admixture (HFX, often referred to also as exact-exchange EXX) in the range of 30–40% for global hybrid functionals.[²⁶⁻²⁸] Most of these studies have focused on small molecules or only a few representative examples of larger coordination compounds, therefore an extensive reference set with sufficient diversity and direct relevance to biological systems is still lacking.

In the present work our first goal is to define such an extensive set of medium-sized Mo(V)
compounds that cover diverse coordination environments with relatively large ligands. After
extensive investigation and validation of the primary literature, we compiled a set of 22
Mo(V) compounds that are structurally well-defined, pose no potential complications due to
non-innocent ligands, and have well-resolved and numerically reliable experimental EPR
parameters. Subsequently, we evaluate a large collection of modern quantum chemical
approaches for the calculation of ⁹⁵Mo HFCs and g-tensors. Given our aim of establishing a
“universal” approach that will be easily applicable to biological systems, typically in the
context of multiscale quantum-mechanics/molecular-mechanics (QM/MM) calculations, we
conduct our investigation using the modern exact 2-component (X2C) Hamiltonian[²⁹⁻³¹] in its
scalar (one-center) form with a finite nucleus model. After determining an optimal choice of
basis sets for Mo, we evaluate a large selection of functionals, sampling several double-hybrid
variants that have not been investigated before in this context. In addition, we examine the
performance of a very recent implementation of the domain-based local pair natural (DLPNO)
approximation to coupled cluster with singles and doubles (CCSD) for ⁹⁵Mo HFCs, a wave
function approach that is fast enough to be applied to realistic-size molecules. Our results
update the state of the art in terms of practical methodologies for accurate calculations of EPR
parameters in Mo systems, highlighting the PBE0-DH double-hybrid functional as a reliable
higher-rung alternative to modified global hybrid functionals for both core and valence
properties.

Methodology

All calculations were performed with the Orca 6.0 program.[32] Geometry optimizations were
carried out with the TPSSh[33] functional and the x2c-TZVPall basis sets.[34] To reduce
computational time the resolution of identity (RI) approximation was used[35] along with the
chain of spheres (COSX) in combination with automatically generated auxiliary basis sets
with the “AutoAux” procedure in Orca.[36] Tight convergence criteria were imposed using the
TightSCF keyword along with the default (DefGrid2) integration grids. Relativistic effects
were included using the exact two-component (X2C) Hamiltonian.[29-31] For the complexes of
the benchmark set that are crystallographically characterized, crystallographic coordinates
were obtained from the Cambridge Structural Database[37] and were used as starting
coordinates for the geometry optimizations, after being individually edited for chemical
correctness and to remove solvent molecules and non-coordinating counter-ions. At first, only
the positions of the hydrogen atoms were optimized, and subsequently, all atoms were
relaxed.

The DFT methods evaluated in the present work for the computation of g-tensors and ⁹⁵Mo
hyperfine coupling tensors include the GGA functionals BP86[³⁸⁻³⁹] and BLYP,[³⁹⁻⁴⁰] MN15-
L,[⁴¹] the meta-GGA TPSS,[⁴²] M06-L,[⁴³] and r²SCAN,[⁴⁴] the hybrid meta-GGA TPSSh,[³³]
the hybrid B3LYP,[³⁹, ⁴⁵] PBE0,[⁴⁶] and B3PW91[⁴⁵, ⁴⁷⁻⁴⁹] with a range of HFX percentages
from 20 to 50%, M05,[⁵⁰] M06,[⁵¹] M06-2X,[⁵¹] and MN15,[⁵²] the range-separated hybrid
ωB97,[⁵³] CAM-B3LYP,[⁵⁴] LC-BLYP, and LC-PBE,[⁵⁵] and the double-hybrid functionals
B2PLYP,[⁵⁶] PBE0-DH,[⁵⁷] B2GP-PLYP,[⁵⁸] ωB2PLYP,[⁵⁹] ωB88PP86,[⁶⁰] and Pr²SCAN50.[⁶¹]
For calculations with double-hybrid functionals, all electrons were included in the MP2
correlation treatment (NoFrozenCore keyword in Orca), and “relaxed” densities were used for
the EPR parameter calculations. All EPR parameter calculations were carried out using
increased integration grids (DefGrid3 settings in Orca). Functionals M05, MN15, and MN15-
L were provided by LibXC.[⁶²]

For the calculation of EPR parameters, we examined different methods for scalar relativity,
including the ZORA[⁶³⁻⁶⁵] and X2C[²⁹⁻³¹] Hamiltonians, along with suitable basis sets. Note
that we use a scalar relativistic version of the X2C Hamiltonian, which is effectively a one-
component method. The X2C Hamiltonian is different from ZORA mainly in that the
decoupling of the one-electron Dirac Hamiltonian is exact, whereas in ZORA it is
approximate. Nevertheless, ZORA and X2C Hamiltonians predict very similar Mo HFCs,
with differences smaller than 2 MHz (Tables S2-S4).

The Mo basis set used for the EPR parameter calculations was carefully selected based on extensive testing as discussed in the next section. Calculations with ZORA were performed using the def2-ZORA-TZVP basis sets[66] on the ligands and the SARC-ZORA-TZVP basis set on Mo.[67] Calculations with X2C were performed using the x2c-TZVPall basis sets[34] on the ligands and we examined the x2c-TZVPall and SARC-DKH-TZVP basis set on Mo. We note that neither of them were designed in their original form for calculations of EPR properties. A SARC basis set for Mo was also explicitly created for the X2C Hamiltonian. To derive such a basis set, starting from the fully decontracted (universal) SARC-TZVP basis set we calculated contraction coefficients using the X2C Hamiltonian, which resulted in a new basis set that we note SARC-X2C-TZVP in the main text. In all cases, the s-decontracted SARC-DKH-TZVP basis set with three additional tight functions[68] is used as a reference. The type of nuclear model used can be important for the calculation of HFCs.[69] Therefore, we compared the point nucleus and finite nucleus[70] approximations for the X2C Hamiltonian, opting for the latter in all method comparisons. In calculations with double-hybrid functionals and in coupled cluster calculations, the SARC-DKH-TZVPP basis set was used for Mo.

The spin- orbit coupling (SOC) operator was treated by a mean- field approximation to the Breit–Pauli operator (SOType 3 in ORCA).[71] Construction of the effective potential included one-electron terms, the Coulomb term was computed exactly, exchange was incorporated via one-center exact integrals including spin-other orbit interactions, and local DFT correlation was included with VWN5 (all these choices are set with the keyword “SOCFlags 1,4,3,1” in ORCA).

The ⁹⁵Mo hyperfine coupling tensors were also computed with the DLPNO-CCSD method using the “unrelaxed” densities.[72-73] For these calculations, the Kohn–Sham orbitals were used to construct the reference determinant. The default DLPNO thresholds[73] for hyperfine coupling calculations were used (HFC1 settings in Orca). These settings involve a tighter PNO threshold for the core electrons, *i.e* $T_{\text{CutPNOCore}} = 1.0×10^{-3}$ instead of $1.0×10^{-2}$ in the NormalPNO and TightPNO settings. Other notable threshold values are $T_{\text{CutMKN}} = 1.0×10^{-4}$, $T_{\text{CutPairs}} = 1.0×10^{-4}$, $T_{\text{CutPNOSingles}} = 0.00$, and $T_{\text{CutPNO}} = 1.0×10^{-7}$. In these calculations, all electrons were treated as active (NoFrozenCore).
The HFC arises from three distinct contributions with different physical origins, which are the Fermi contact term ($A^{\text{FC}}$), the spin-dipolar term ($A^{\text{SD}}$), and the spin-orbit coupling term ($A^{\text{SO}}$). The $A^{\text{FC}}$ and $A^{\text{SD}}$ terms are first-order properties. The $A^{\text{FC}}$ is an isotropic term arising from spin polarization at the *s* orbitals of Mo and is proportional to the spin density at the nucleus. The $A^{\text{SD}}$ is purely anisotropic and arises from dipole-dipole interactions between the nuclear spin and the spin density in other orbitals (Mo *p*, *d*, *f* and orbitals of other atoms). The $A^{\text{SO}}$ term arises predominantly from second-order contributions, while higher-order SOC effects contribute less than 2% of HFC components and are therefore less important for Mo than for heavier elements,[74] such as W.[20] Here, the $A^{\text{SO}}$ term has been computed up to second-order. The $A^{\text{SO}}$ term contains both an isotropic and an anisotropic contribution. For better comparison with experimental values, the $A^{\text{SO}}$ term obtained from Orca ($A_{ii}^{\text{SO}}$, $i = 1, 2, 3$) has been separated into an isotropic term, known as pseudocontact term $A^{\text{PC}}$, and an anisotropic term, $A_{ii}^{\text{SO,an}}$, such that: $A_{ii}^{\text{SO}} = A^{\text{PC}} + A_{ii}^{\text{SO,an}}$. Using these definitions, the components $A_{ii}$ of the complete HFC tensor – up to second-order perturbation theory – can be written as: $A_{ii} = A_{\text{iso}} + A_{\text{an}}$, where the isotropic HFC contribution $A_{\text{iso}}$ corresponds to: $A_{\text{iso}} = A^{\text{FC}} + A^{\text{PC}}$, and the anisotropic contribution $A_{\text{an}}$ is: $A_{\text{an}} = A_{ii}^{\text{SD}} + A_{ii}^{\text{SO,an}}$.
# Results and Discussion

## Definition of a Reference Set of Compounds

Although the subject has been approached in past literature,[16-25] there seems to be no standard collection of Mo(V) complexes to serve as a comprehensive reference set. Therefore, our first goal was to define a sufficiently diverse and representative set of Mo-based mononuclear complexes with reliable experimentally resolved g-tensors and hyperfine coupling constants. After careful examination of the primary sources, we selected twenty-two compounds to serve as our reference set. The optimized structures of all compounds are shown in Figure 2.

The complexes encompass various Mo coordination environments with S, O, N and Cl ligand atom types that render our conclusions broadly applicable to a wide range of systems. All complexes have $S = \frac{1}{2}$ ground states, with molybdenum being in its Mo(V) oxidation state (4d¹). Tris-dithiolene monoanions **1** and **2** feature a distorted trigonal prismatic coordination sphere. Compounds **3**, **4**, **19** and **20** are square pyramidal Mo-oxo species in which Mo is five-coordinate, slightly elevated from the equatorial plane. The remaining compounds are oxo-molybdenum complexes where Mo is six-coordinate; among these, complexes **7**, **16**, and **17** feature sulfido-ligands on the Mo center. Complexes **10** and **12**, along with their protonated forms **11** and **13**, respectively, resemble the sulfite oxidase family of enzymes and their proposed catalytic intermediates. Structures **21** and **22**, on the other hand, represent the *cis* and *trans* geometric isomers of the same compound with an heteroscorpionate ligand. The compounds of the benchmark set that have been crystallographically characterized are: **1**,[75] **3**,[76] **4**,[77] **5** and **6**,[78] **9**,[79] **15-17**,[80] **19**,[81] **20**,[82] and **21-22**[25]. For these complexes, the root mean square deviations (RMSDs) of the geometrically optimized (TPSSh) structures from the crystallographic structures are smaller than 0.08 Å, which confirms that the optimization protocol should be sufficiently reliable for all optimized structures of our reference set.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Structures of the Mo(V) complexes considered in the present study. Mo centers are colored light blue; S atoms yellow; O red; N dark blue; B purple; Cl green; C grey; H white. Definition of ligand abbreviations: bdt = benzenedithiolene; qdt = quinoxaline-2,3-dithiolene; L3S = (2-dimethylethanethiol)bis(3,5-dimethylpyrazolyl)methane; L1 = N,N'-dimethyl-N,N'-bis(2-mercaptophenyl)-1,2-ethylenediamine; L2 = N,N'-dimethyl-N,N'-bis(2-hydroxyphenyl)-1,2-ethylenediamine; Tp* = tris-(3,5-dimethylpyrazolyl)hydroborate; mdt = 1,2-dimethyl-1,2-dithiolene; L1O = bis(3,5-dimethylpyrazolyl)-3-tert-butyl-2-hydroxy-5-methylphenyl)methane.</figcaption></figure>

The *g*-tensors and HFCs for all complexes in the benchmark set, as determined by EPR, are presented in Table 1. The *g* tensor components are expressed as shifts (Δ*g*) from the free electron *g*-value (*g*<sub>e</sub> = 2.0023), in parts per thousand (ppt). It is important to acknowledge the inherent limitations[19] in the reliability of molybdenum HFCs due to the low natural abundance of molybdenum isotopes with magnetic nuclei (⁹⁵Mo: 15.9% and ⁹⁷Mo: 9.6%, both with *I* = 5/2), arising from the use of natural-abundance samples in the EPR experiments.
Nevertheless, the selected compounds have been carefully chosen to ensure the highest possible data quality. Complexes **1-2**, **7-11**, and **14-18** were measured using multifrequency EPR, and for complexes **14** and **19** refined values[19] have been considered. The robustness of the dataset is further evidenced by the absence of significant and persistent outliers in the comparison of different methods, even in light of the structural diversity of the complexes, underscoring the reliability of the data. The *g*-tensors of the selected compounds cover a broad range of values and symmetries, with *g*-shifts spanning from -57 to 22 ppt and from -248 to -8 ppt for $\Delta g_{11}$ and $\Delta g_{33}$, respectively, while *g*-tensor anisotropy, $\Delta g_{11} - \Delta g_{33}$, ranges from 2 to 225 ppt. In all complexes, the Mo HFC tensor is predominantly axial, except for complex **1**, where it exhibits rhombic symmetry. The largest HFC component ($A_{11}$) varies within the range of 117 to 249 MHz.

<table><caption>**Table 1.** List of the Mo(V) complexes considered in the present study and their experimentally determined *g*-shifts ($\Delta g_{ii} = g_e - g_{ii}$ in parts per thousand, ppt) and $^{95}$Mo hyperfine coupling constants ($A_{ii}$ in MHz).</caption><thead><tr><th>Complex</th><th>Δg<sub>11</sub></th><th>Δg<sub>22</sub></th><th>Δg<sub>33</sub></th><th>A<sub>11</sub></th><th>A<sub>22</sub></th><th>A<sub>33</sub></th><th>Ref.</th></tr></thead><tbody><tr><td>1 [Mo(bdt)<sub>3</sub>]<sup>-</sup></td><td>18</td><td>5</td><td>-8</td><td>117</td><td>87</td><td>39</td><td>[83]</td></tr><tr><td>2 [Mo(qdt)<sub>3</sub>]<sup>-</sup></td><td>10</td><td>4</td><td>-9</td><td>120</td><td>117</td><td>21</td><td>[83]</td></tr><tr><td>3 [MoO(SPh)<sub>4</sub>]<sup>-</sup></td><td>18</td><td>-20</td><td>-20</td><td>157</td><td>66</td><td>66</td><td>[84]</td></tr><tr><td>4 [MoO(bdt)<sub>2</sub>]<sup>-</sup></td><td>20</td><td>-18</td><td>-26</td><td>-</td><td>-</td><td>-</td><td>[85]</td></tr><tr><td>5 (L3S)MoO(SPh)<sub>2</sub></td><td>12</td><td>-36</td><td>-56</td><td>165</td><td>87</td><td>81</td><td>[78]</td></tr><tr><td>6 (L3S)MoO(bdt)</td><td>14</td><td>-30</td><td>-53</td><td>156</td><td>78</td><td>72</td><td>[78]</td></tr><tr><td>7 [(L1)MoOS]<sup>-</sup></td><td>14</td><td>-69</td><td>-114</td><td>160</td><td>70</td><td>68</td><td>[86]</td></tr><tr><td>8 (L1)MoO(SH)</td><td>13</td><td>-43</td><td>-50</td><td>153</td><td>69</td><td>70</td><td>[86]</td></tr><tr><td>9 (L1)MoO(SEt)</td><td>22</td><td>-38</td><td>-47</td><td>175</td><td>71</td><td>67</td><td>[79]</td></tr><tr><td>10 [(L1)MoO<sub>2</sub>]<sup>-</sup></td><td>-16</td><td>-87</td><td>-192</td><td>206</td><td>94</td><td>81</td><td>[86]</td></tr><tr><td>11 (L1)MoO(OH)</td><td>-22</td><td>-55</td><td>-59</td><td>192</td><td>87</td><td>73</td><td>[86]</td></tr><tr><td>12 [(L2)MoO<sub>2</sub>]<sup>-</sup></td><td>-23</td><td>-105</td><td>-248</td><td>-</td><td>-</td><td>-</td><td>[87]</td></tr><tr><td>13 (L2)MoO(OH)</td><td>-22</td><td>-55</td><td>-59</td><td>-</td><td>-</td><td>-</td><td>[87]</td></tr><tr><td>14 [(Tp*)MoO(S<sub>2</sub>CNEt<sub>2</sub>)]<sup>+</sup></td><td>-22</td><td>-32</td><td>-48</td><td>193</td><td>90</td><td>85<sup>a</sup></td><td>[88]</td></tr><tr><td>15 (Tp*)MoO(bdt)</td><td>2</td><td>-30</td><td>-68</td><td>180</td><td>78</td><td>72</td><td>[80]</td></tr><tr><td>16 (Tp*)MoS(bdt)</td><td>-5</td><td>-34</td><td>-86</td><td>177</td><td>79</td><td>78</td><td>[80]</td></tr><tr><td>17 (Tp*)MoS(cat)</td><td>-38</td><td>-43</td><td>-105</td><td>202</td><td>90</td><td>87</td><td>[80]</td></tr><tr><td>18 (Tp*)MoO(cat)</td><td>-34</td><td>-36</td><td>-83</td><td>192</td><td>81</td><td>78</td><td>[80]</td></tr><tr><td>19 [MoOCl<sub>2</sub>(mdt)]<sup>-</sup></td><td>2</td><td>-35</td><td>-59</td><td>186<sup>a</sup></td><td>75<sup>a</sup></td><td>75<sup>a</sup></td><td>[81]</td></tr><tr><td>20 [MoOCl<sub>4</sub>]<sup>-</sup></td><td>-37</td><td>-37</td><td>-56</td><td>227</td><td>103</td><td>103</td><td>[89]</td></tr><tr><td>21 trans-(L1O)MoOCl<sub>2</sub></td><td>-37</td><td>-45</td><td>-56</td><td>226</td><td>104</td><td>90</td><td>[25]</td></tr><tr><td>22 cis-(L1O)MoOCl<sub>2</sub></td><td>-57</td><td>-61</td><td>-73</td><td>226</td><td>105</td><td>96</td><td>[25]</td></tr></tbody></table>

<footnote><sup>a</sup>Simulated values were obtained from Fritscher et al.[19]</footnote>
decontracted versions. For the present purposes, we created an adapted version of the SARC basis set for Mo, in which we retained the common primitive GTFs but contracted the CGTFs in each angular momentum subspace in a way that is fully adapted to the X2C Hamiltonian. Subsequently we discuss all variants, but note that the optimal contraction coefficients for the X2C Hamiltonian are extremely similar to those of the DKH2 Hamiltonian, to the point where the two basis sets are almost indistinguishable. This is not the case with the ZORA Hamiltonian, where CGTFs have obviously distinct contraction coefficients.

The impact of this difference is evident when comparing the SARC-DKH-TZVP basis set for Mo with the SARC-ZORA-TZVP, both used with the X2C Hamiltonian: The ZORA version produces a significantly different $A^{FC}$ term, although the anisotropic components are very similar. Decontracting the s-CGTF of the SARC-DKH-TZVP and even supplanting it by three additional tight (high angular momentum) s functions produced by multiplying the tighter exponent by 2.5, 6.25 and 15.625, in the spirit of the “core property” adaptation of older basis sets for lighter elements,[68] does not provide perceptible changes in the computed properties. Substitution of the DKH-adapted with the X2C-adapted higher angular momentum CGTFs has no effect either. We note that the s-decontracted x2c-TZVPall basis sets are closer to the default (and modified) SARC results, confirming that the x2c-TZVPall basis sets in their default contraction cannot be used for calculations of the target properties.

<table><caption>**Table 2.** Calculated HFCs (total values and individual components, in MHz) of compound 1 obtained using the X2C Hamiltonian with the B3PW91 functional with 40% HFX and finite nucleus using different basis sets for Mo.</caption><thead><tr><th rowspan="2">Mo Basis Set</th><th colspan="3">Total</th><th colspan="3">Isotropic</th><th colspan="3">Anisotropic</th></tr><tr><th>A<sub>11</sub></th><th>A<sub>22</sub></th><th>A<sub>33</sub></th><th>A<sup>FC</sup></th><th>A<sup>PC</sup></th><th>A<sub>11</sub><sup>SD</sup></th><th>A<sub>22</sub><sup>SD</sup></th><th>A<sub>33</sub><sup>SD</sup></th><th>A<sub>11</sub><sup>SO,an</sup></th><th>A<sub>22</sub><sup>SO,an</sup></th><th>A<sub>33</sub><sup>SO,an</sup></th></tr></thead><tbody><tr><td>x2c-TZVPall</td><td>113.1</td><td>94.0</td><td>17.8</td><td>63.5</td><td>11.5</td><td>34.4</td><td>16.6</td><td>-51.0</td><td>3.8</td><td>2.5</td><td>-6.3</td></tr><tr><td>x2c-TZVPall s-decontracted</td><td>134.9</td><td>115.6</td><td>39.0</td><td>85.1</td><td>11.5</td><td>34.6</td><td>16.7</td><td>-51.3</td><td>3.9</td><td>2.5</td><td>-6.3</td></tr><tr><td>SARC-DKH-TZVP</td><td>124.4</td><td>107.1</td><td>38.2</td><td>83.2</td><td>6.7</td><td>31.7</td><td>15.4</td><td>-47.1</td><td>2.9</td><td>1.9</td><td>-4.8</td></tr><tr><td>SARC-ZORA-TZVP</td><td>116.4</td><td>99.1</td><td>30.2</td><td>75.2</td><td>6.8</td><td>31.7</td><td>15.4</td><td>-47.1</td><td>2.9</td><td>1.9</td><td>-4.8</td></tr><tr><td>SARC-DKH-TZVP s-decontracted +3s</td><td>123.9</td><td>106.7</td><td>37.7</td><td>82.7</td><td>6.7</td><td>31.7</td><td>15.4</td><td>-47.1</td><td>2.9</td><td>1.9</td><td>-4.8</td></tr><tr><td>SARC-X2C-TZVP s-decontracted +3s</td><td>124.0</td><td>106.7</td><td>37.7</td><td>82.7</td><td>6.7</td><td>31.7</td><td>15.4</td><td>-47.1</td><td>2.9</td><td>1.9</td><td>-4.8</td></tr></tbody></table>

The choice of a point nucleus versus a finite nucleus, the latter typically modeled by a Gaussian charge distribution, is another methodological choice that can have an impact on the computed HFC values, specifically on the $A^{FC}$ term. Here we compared the two options for
the X2C Hamiltonian and observed that the differences are less than 0.5 MHz for the Fermi
contact term, an effect that is rather insignificant in view of other sources of uncertainty in
these calculations. Nevertheless, we decided to use the finite nucleus model in all subsequent
calculations because it does not incur additional computational cost and is physically more
correct. Under normal circumstances the use of a finite nucleus model may require specific
adaptation of the basis set in the form of tighter *s* functions, i.e. primitive GTFs with higher
exponents, but the fact that here we opt for the *s*-decontracted version of an already highly
flexible basis set ensures that the basis set is also converged when a finite nucleus model is
employed. We also confirmed that the point vs finite nucleus properties are converged with
respect to the extension of the *s*-decontracted space through extra-tight *s* functions. The
results presented in Table 2 and Tables S2-S3 indicate that the default SARC-DKH-TZVP
basis set for Mo with its default contraction pattern is fully sufficient for all practical
calculations of Mo EPR properties with the X2C Hamiltonian and a finite nucleus model.

The above discussion focused on the HFC, a property that is very sensitive to the flexibility of
the basis set, particularly in the core region. The *g*-tensor depends more sensitively on the
description of the valence electronic structure and thus can be sensitive to higher angular
momenta. The basis sets we use are already sufficiently polarized, nevertheless we tested the
effect of further extending the high-angular momentum subspace of the basis set. These
calculations confirmed that addition of more and higher angular momentum functions
compared to the default polarization functions of the SARC-DKH-TZVP basis set has
negligible effect on the computed *g*-tensors (Table S4). Therefore, the default version of the
SARC-DKH-TZVP basis set as used for <sup>95</sup>Mo HFCs is also converged in terms of size for *g*-
tensor calculations of Mo compounds. This allows us to use a common basis set for all EPR
properties and this is what will be applied in all calculations reported in the following
sections.
⁹⁵Mo Hyperfine Coupling Constants from DFT

After having established convergence of DFT calculated hyperfine coupling constants and g-tensors with respect to the Mo basis set, we now evaluate the performance of various density functional approaches. To this end, we calculated the ⁹⁵Mo HFCs of 19 complexes in our benchmark set, excluding complexes 4, 12, and 13, as their HFCs have not been experimentally determined (Tables S5-S31). We employed a diverse selection of DFT approaches, encompassing various functional families: GGA (BP86, BLYP, M06-L), meta-GGA (TPSS, r²SCAN), hybrid (TPSSh, B3LYP, B3PW91 with a range of HFX percentages from 20 to 50%), range separated hybrid functionals (ωB97, CAM-B3LYP, LC-BLYP, LC-PBE), and double-hybrid (B2PLYP, PBE0-DH, B2GP-PLYP, ωB2PLYP, ωB88PP86, Pr²SCAN50). For the calculations, we employed the X2C Hamiltonian with the SARC-DKH-TZVP basis set on molybdenum, except for the double-hybrid functionals, where the SARC-DKH-TZVPP basis set was used to maximize electron correlation recovery during the MP2 calculation.

The performance of each method for the molybdenum hyperfine coupling constants is evaluated using the absolute percent deviation (APD) relative to the experimental values. This assessment focuses on three key parameters; the isotropic HFC ($A_{iso}$), which is the average of the three HF tensor components; the anisotropic difference ($A_{11} - A_{33}$), which is the difference between the largest and smallest HF tensor components; and the largest component of the HF tensor ($A_{11}$). The APDs for each of these parameters are calculated as follows:

$$APD = \left| \frac{A^{\text{calc}} - A^{\text{exp}}}{A^{\text{exp}}} \right| \times 100$$

The mean APDs (MAPDs) averaged across the 19 examined complexes for each functional, along with their standard deviations, are given in Table 3, and the MAPDs of $A_{iso}$ and $A_{11} - A_{33}$ are graphically represented in Figure 3. The functionals are arranged in ascending order of MAPD($A_{iso}$), and in Table 3 they are further categorized according to their respective rungs on the DFT ladder, to enable a more straightforward comparison of methods within the same rung. Among the three criteria, MAPD($A_{iso}$) emerges as the most significant metric due to its significant differentiation among the examined methods. In contrast, as illustrated in Figure 3, MAPD($A_{11} - A_{33}$) exhibits less variability, as a consequence of the largely systematic nature of the deviations from the experimental values.
<table><caption>**Table 3.** MAPDs of A<sub>iso</sub>, A<sub>11</sub> – A<sub>33</sub> and A<sub>11</sub> and their standard deviations (SD in MHz) for calculated HFCs using different DFT functionals with respect to experimental HFC values for 19 Mo(V) complexes of the benchmark set.</caption><thead><tr><th rowspan="2">Functional</th><th colspan="2">A<sub>iso</sub></th><th colspan="2">A<sub>11</sub> - A<sub>33</sub></th><th colspan="2">A<sub>11</sub></th></tr><tr><th>MAPD</th><th>SD</th><th>MAPD</th><th>SD</th><th>MAPD</th><th>SD</th></tr></thead><tbody><tr><td><strong>GGA</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>BLYP</td><td>37.6</td><td>3.8</td><td>19.5</td><td>7.9</td><td>31.3</td><td>3.0</td></tr><tr><td>BP86</td><td>39.7</td><td>3.9</td><td>21.7</td><td>7.6</td><td>33.4</td><td>2.9</td></tr><tr><td>MN15-L</td><td>61.2</td><td>8.4</td><td>20.7</td><td>8.9</td><td>46.4</td><td>4.3</td></tr><tr><td><strong>meta-GGA</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>r<sup>2</sup>SCAN</td><td>18.5</td><td>6.2</td><td>22.8</td><td>7.6</td><td>20.2</td><td>4.4</td></tr><tr><td>M06-L</td><td>20.0</td><td>8.3</td><td>21.9</td><td>7.9</td><td>19.4</td><td>8.1</td></tr><tr><td>TPSS</td><td>25.0</td><td>4.4</td><td>21.7</td><td>7.6</td><td>24.1</td><td>3.5</td></tr><tr><td><strong>Hybrid</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>B3PW91 40%</td><td>4.9</td><td>4.3</td><td>11.7</td><td>5.7</td><td>5.2</td><td>4.1</td></tr><tr><td>B3PW91 50%</td><td>10.4</td><td>6.2</td><td>11.1</td><td>5.4</td><td>5.9</td><td>4.5</td></tr><tr><td>B3PW91 30%</td><td>11.2</td><td>5.2</td><td>12.6</td><td>6.6</td><td>11.1</td><td>4.7</td></tr><tr><td>PBE0</td><td>12.8</td><td>5.2</td><td>13.7</td><td>7.9</td><td>13.1</td><td>4.4</td></tr><tr><td>TPSSh</td><td>15.4</td><td>5.0</td><td>17.8</td><td>8.2</td><td>16.5</td><td>4.1</td></tr><tr><td>B3LYP</td><td>19.2</td><td>4.7</td><td>12.7</td><td>7.5</td><td>16.7</td><td>4.2</td></tr><tr><td>M05</td><td>20.2</td><td>12.6</td><td>11.4</td><td>7.5</td><td>15.8</td><td>8.5</td></tr><tr><td>B3PW91 20%</td><td>21.3</td><td>4.6</td><td>13.9</td><td>8.3</td><td>18.7</td><td>3.9</td></tr><tr><td>M06-2X</td><td>66.4</td><td>19.8</td><td>29.5</td><td>58.7</td><td>50.5</td><td>37.3</td></tr><tr><td>MN15</td><td>78.4</td><td>10.0</td><td>16.6</td><td>18.4</td><td>49.9</td><td>9.5</td></tr><tr><td>M06</td><td>80.4</td><td>12.7</td><td>19.0</td><td>21.4</td><td>57.5</td><td>17.8</td></tr><tr><td><strong>Range-separated</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ωB97</td><td>13.8</td><td>7.2</td><td>11.7</td><td>6.3</td><td>8.5</td><td>5.7</td></tr><tr><td>CAM-B3LYP</td><td>14.0</td><td>5.3</td><td>11.2</td><td>5.9</td><td>11.5</td><td>5.3</td></tr><tr><td>LC-PBE</td><td>21.1</td><td>6.0</td><td>12.4</td><td>6.1</td><td>15.9</td><td>5.4</td></tr><tr><td>LC-BLYP</td><td>25.8</td><td>5.0</td><td>11.4</td><td>6.5</td><td>18.9</td><td>5.2</td></tr><tr><td><strong>Double-hybrid</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>PBE0-DH</td><td>4.2</td><td>3.5</td><td>11.6</td><td>6.6</td><td>5.5</td><td>3.7</td></tr><tr><td>ωB2PLYP</td><td>7.6</td><td>7.7</td><td>11.3</td><td>7.2</td><td>6.2</td><td>5.1</td></tr><tr><td>Pr2SCAN50</td><td>9.2</td><td>8.5</td><td>14.2</td><td>7.3</td><td>9.5</td><td>5.8</td></tr><tr><td>ωB88PP86</td><td>10.9</td><td>10.8</td><td>12.2</td><td>7.2</td><td>9.0</td><td>7.2</td></tr><tr><td>B2GP-PLYP</td><td>11.6</td><td>16.4</td><td>14.8</td><td>12.6</td><td>9.3</td><td>8.5</td></tr><tr><td>B2PLYP</td><td>15.0</td><td>11.1</td><td>13.8</td><td>7.7</td><td>12.9</td><td>5.8</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Graphical representation of the MAPDs of A<sub>iso</sub> and A<sub>11</sub> – A<sub>33</sub> for the calculated values of the set of 19 Mo(V) complexes obtained with different functionals with respect to experimental results.</figcaption></figure>

Looking at Figure 3, we observe that double-hybrid functionals demonstrate the best performance. Among these, PBE0-DH stands out with the lowest MAPD(A<sub>iso</sub>) of 4.2%, while ωB2PLYP and Pr<sup>2</sup>SCAN50 also perform well, maintaining MAPD(A<sub>iso</sub>) values below 10%. Among all other methods, only the hybrid functional B3PW91 with 40% HFX achieves a comparable MAPD(A<sub>iso</sub>) of 4.9%. In terms of standard deviations, PBE0-DH again excels, achieving the lowest SD(A<sub>iso</sub>) of 3.5%, outperforming B3PW91 40% which has a SD(A<sub>iso</sub>) of 4.3%. As shown in Table 3, both PBE0-DH and B3PW91 40% demonstrate similarly strong
performance in the MAPD($A_{11}$ – $A_{33}$) and MAPD($A_{11}$) criteria, making them the top-performing methods among those examined.

Regarding hybrid functionals, their MAPDs depend heavily on the HFX percentage. The results confirm previous suggestions that the optimal HFX percentage for HFCs of ⁹⁵Mo complexes is between 30 and 50%, and particularly 40% is strongly favorable in the case of Mo(V) systems. All functionals with HFX percentage smaller than 30% show MAPD($A_{iso}$) values that exceed 12%. These methods tend to underestimate the HF tensor components. In addition, range-separated functionals, perform similarly to hybrid functionals with lower HFX percentages.

GGA and meta-GGA functionals occupy the bottom half of the plot in Figure 3. Meta-GGA functionals, however, show a notable improvement over GGA functionals, achieving MAPDs close to 20%, whereas GGA functionals exhibit MAPDs exceeding 30% for both the MAPD($A_{iso}$) and MAPD($A_{11}$) criteria. GGA functionals systematically underestimate the HF tensor components, a consequence of their well-documented tendency to overestimate the covalency of the M–L bond, leading to an underestimation of the spin density on the Mo nucleus. The systematic nature of the deviations is evident in the relatively small standard deviations of BLYP and BP86, which are comparable to those of PBE0-DH. This consistency suggests that these very low-cost methods could be effectively paired with a scaling factor, making them suitable for handling very large systems or for efficiently screening extensive datasets of compounds.

The hybrid Minnesota functionals M06-2X, MN15, and M06 show by far the worst performances for the HFC calculations, with MAPDs higher than 50% for $A_{iso}$ and $A_{11}$. Interestingly, the anisotropic parameters $A_{11}$ – $A_{33}$ for MN15 and M06 are not notably larger than those of other methods. However, they exhibit very large standard deviations, indicating that the differences from the experimental values are not systematic. This result mirrors previous reports regarding the poor performance of these functionals for transition metal HFCs.[26, 95-96] In view of the good performance of Minnesota functionals in thermochemistry, thermochemical kinetics, and noncovalent interactions, the present limitation likely arises from the large number of empirical parameters used in their construction, which constrains the spectrum of properties that they can accurately target.[97] Notably, the degree of spin contamination in these functionals is similar to that of B3PW91 40%, with the deviation from the ideal value of 0.75 ranging between 0.003 and 0.070. Moreover, the computed total Mo spin populations are not substantially different from those of other hybrid functionals (see
Supplementary Information Tables S33-S60), which suggests that the problem likely lies in the poor description of spin polarization in the (semi)-core region.[24, 26]

The PBE0-DH results for all complexes are shown in Table 4. We observe that the largest contributor of the ⁹⁵Mo HF components is the Fermi contact term, $A^{FC}$. Spin-orbit coupling contributions to the isotropic term ($A^{PC}$) are up to 15 MHz and account on average for 7% of the $A_{iso}$ term. The rhombic character for compound **1** and the axial character of all other complexes are correctly predicted by the calculations. The predicted anisotropic spin-orbit contributions also participate in the axial character, especially for **10** in which the $A_{11}^{SO}$ term is 19.7 MHz.

<table><caption>**Table 4.** Calculated HFCs (total values and individual components, in MHz) obtained using the PBE0-DH functional for the 22 Mo(V) complexes of the benchmark set.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Total</th><th colspan="3">Isotropic</th><th colspan="3">Anisotropic</th></tr><tr><th>A<sub>11</sub></th><th>A<sub>22</sub></th><th>A<sub>33</sub></th><th>A<sup>FC</sup></th><th>A<sup>PC</sup></th><th>A<sub>11</sub><sup>SD</sup></th><th>A<sub>22</sub><sup>SD</sup></th><th>A<sub>33</sub><sup>SD</sup></th><th>A<sub>11</sub><sup>SO,an</sup></th><th>A<sub>22</sub><sup>SO,an</sup></th><th>A<sub>33</sub><sup>SO,an</sup></th></tr></thead><tbody><tr><td>1</td><td>119.3</td><td>107.8</td><td>37.0</td><td>83.2</td><td>4.8</td><td>28.3</td><td>17.5</td><td>-45.8</td><td>3.0</td><td>2.3</td><td>-5.3</td></tr><tr><td>2</td><td>115.5</td><td>113.6</td><td>36.5</td><td>83.6</td><td>4.9</td><td>24.5</td><td>22.6</td><td>-47.1</td><td>2.5</td><td>2.5</td><td>-5.0</td></tr><tr><td>3</td><td>150.0</td><td>64.9</td><td>64.9</td><td>88.2</td><td>5.1</td><td>54.7</td><td>-27.4</td><td>-27.3</td><td>2.2</td><td>-1.1</td><td>-1.1</td></tr><tr><td>5</td><td>159.0</td><td>75.8</td><td>72.1</td><td>96.2</td><td>6.1</td><td>54.2</td><td>-25.6</td><td>-28.5</td><td>2.7</td><td>-1.0</td><td>-1.8</td></tr><tr><td>6</td><td>160.1</td><td>72.4</td><td>69.9</td><td>95.3</td><td>5.5</td><td>56.2</td><td>-27.1</td><td>-29.1</td><td>3.3</td><td>-1.4</td><td>-1.9</td></tr><tr><td>7</td><td>151.2</td><td>74.3</td><td>71.5</td><td>88.2</td><td>10.8</td><td>44.6</td><td>-22.7</td><td>-21.9</td><td>7.8</td><td>-2.1</td><td>-5.6</td></tr><tr><td>8</td><td>162.3</td><td>76.0</td><td>72.7</td><td>97.3</td><td>6.4</td><td>56.3</td><td>-26.7</td><td>-29.6</td><td>2.5</td><td>-1.0</td><td>-1.5</td></tr><tr><td>9</td><td>155.7</td><td>73.8</td><td>68.9</td><td>93.2</td><td>6.2</td><td>54.3</td><td>-25.1</td><td>-29.2</td><td>2.1</td><td>-0.6</td><td>-1.5</td></tr><tr><td>10</td><td>205.7</td><td>99.1</td><td>90.4</td><td>116.3</td><td>15.4</td><td>54.4</td><td>-26.8</td><td>-27.6</td><td>19.7</td><td>-5.9</td><td>-13.8</td></tr><tr><td>11</td><td>185.2</td><td>91.1</td><td>85.8</td><td>111.9</td><td>8.8</td><td>59.8</td><td>-27.7</td><td>-32.1</td><td>4.9</td><td>-2.0</td><td>-2.9</td></tr><tr><td>14</td><td>183.7</td><td>86.1</td><td>83.7</td><td>110.8</td><td>7.0</td><td>61.7</td><td>-30.4</td><td>-31.3</td><td>4.4</td><td>-1.5</td><td>-2.9</td></tr><tr><td>15</td><td>173.7</td><td>81.2</td><td>77.5</td><td>104.2</td><td>6.6</td><td>58.2</td><td>-27.5</td><td>-30.7</td><td>4.9</td><td>-2.2</td><td>-2.7</td></tr><tr><td>16</td><td>166.3</td><td>77.2</td><td>67.0</td><td>96.4</td><td>7.0</td><td>57.4</td><td>-23.9</td><td>-33.5</td><td>5.6</td><td>-2.5</td><td>-3.1</td></tr><tr><td>17</td><td>186.9</td><td>85.0</td><td>82.7</td><td>107.8</td><td>10.4</td><td>62.1</td><td>-27.6</td><td>-34.5</td><td>6.8</td><td>-5.7</td><td>-1.1</td></tr><tr><td>18</td><td>187.0</td><td>81.0</td><td>80.7</td><td>107.0</td><td>9.2</td><td>63.7</td><td>-33.1</td><td>-30.6</td><td>7.2</td><td>-2.3</td><td>-4.9</td></tr><tr><td>19</td><td>179.5</td><td>84.3</td><td>81.7</td><td>108.1</td><td>7.0</td><td>60.6</td><td>-29.0</td><td>-31.6</td><td>4.0</td><td>-2.0</td><td>-2.0</td></tr><tr><td>20</td><td>206.4</td><td>100.5</td><td>100.5</td><td>126.1</td><td>9.7</td><td>65.9</td><td>-33.0</td><td>-33.0</td><td>4.8</td><td>-2.4</td><td>-2.4</td></tr><tr><td>21</td><td>207.2</td><td>105.6</td><td>103.8</td><td>129.2</td><td>9.6</td><td>62.3</td><td>-30.1</td><td>-32.2</td><td>6.2</td><td>-3.2</td><td>-3.0</td></tr><tr><td>22</td><td>209.9</td><td>107.5</td><td>105.0</td><td>129.5</td><td>11.3</td><td>63.8</td><td>-31.4</td><td>-32.4</td><td>5.5</td><td>-2.1</td><td>-3.5</td></tr></tbody></table>

The $A^{FC}$, $A_{33}^{SD}$, and $A^{PC}$ components of complex **1** calculated using all examined DFT functionals are plotted in Figure 4. The plot shows that the differences between the results obtained using different functionals can be traced back to the calculated Fermi contact term.
All functionals predict $A_{11}^{SD}$ and $A^{PC}$ components, in the range 22.9 to 33.1 MHz and 1.4 and 7.9 MHz, respectively, except for some Minnesota functionals that are outliers. Therefore, achieving accurate $A^{FC}$ is the biggest methodological challenge for Mo(V) compounds.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Decomposition of the calculated HFCs for compound 1 obtained using different functionals.</figcaption></figure>

To investigate the origins of the different behaviors of various DFT functionals on the ⁹⁵Mo HFCs, we analyze the molecular orbital contributions to the Fermi contact term, focusing on the core spin polarization. Figure 5 shows the molecular orbital analysis of the Mo $A^{FC}$ term of complex 1 obtained using a selection of different DFT approaches. As ⁹⁵Mo has a negative nuclear *g*-factor, this induces negative spin density at the nucleus from the 2s and 4s orbitals. The poor performance of functionals M06 and BLYP can be correlated with significantly smaller 4s contributions to the Fermi contact term relative to the best-performing functionals PBE0-DH and B3PW91 with 40% HFX. Moreover, M06 and M05 predict relatively larger 2s and 3s contributions. This results in lower ratios between the 4s contributions and those of the 2s and 3s orbitals compared to better-performing methods, resembling previous observations in 3d systems.[24, 26]
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Core orbital contributions for compound 1 using different functionals.</figcaption></figure>

The performance of double-hybrid functionals is likely correlated to the balance between HFX and MP2 contributions. B2PLYP (53% HFX and 27% MP2) has previously demonstrated superior performance for predicting HFCs of transition metal complexes compared to global hybrid functionals,[98] suggesting that the MP2 contribution may compensate for exaggerated spin polarization caused by the high HFX percentage. Of course, the balance between HFX and MP2 components is likely system-dependent, while both components may eventually limit the performance of any given double-hybrid functional in cases of transition metal systems with significant static correlation (relevant examples have been noted in the case of NMR chemical shifts for 3d transition metal systems)[26]. For the present Mo(V) systems it appears that PBE0-DH with 50% HFX and 12.5% MP2 (notably, PBE0-DH was constructed on a purely theoretical basis without empirical fitting[99]) achieves the best compromise, while the other double-hybrid functionals with HFX percentages of 50% or higher and MP2 contributions exceeding 25% lag behind PBE0-DH.
## Coupled Cluster Calculations of $^{95}$Mo HFCs

The coupled-cluster method with singles, doubles, and perturbative triples, CCSD(T), commonly regarded as the “gold standard” for quantum chemistry, provides highly accurate HFCs for transition metal complexes.[27] However, its steep computational scaling significantly restricts its application to small, simplified systems. In contrast, most of the systems analyzed in this study involve more than 1000 basis functions, making CCSD(T) impractical for such cases. The DLPNO-CCSD(T) method leverages the locality of electron correlation to reduce the information content of the coupled-cluster wavefunction. This approach enables the calculation of properties of transition metal complexes that have been notoriously challenging for DFT methods, extending coupled-cluster accuracy to larger systems, including enzyme models. A recent development is the implementation of the DLPNO-based Lagrangian scheme, which enables the computation of spin-densities that approximate the canonical CCSD densities, allowing the calculation of hyperfine coupling constants.[73] While these methods are computationally more demanding than DFT, they have the significant advantage of avoiding system-dependent biases, thereby offering robust and transferable computational protocols. Aiming to establish a universal approach for enzymes and biomimetic complexes, we assessed the DLPNO-CCSD method for the prediction of $^{95}$Mo HFCs. We note that the DLPNO-CCSD method in its current implementation can only provide the Fermi contact ($A^{FC}$) and spin dipolar ($A^{SD}$) terms. In the absence of the spin-orbit contribution, the total HFCs cannot be obtained directly at this level. The DLPNO-CCSD results ($A^{FC} + A^{SD}$) are given in Table 5, excluding complexes 16 and 21, where convergence issues were encountered.

<table><caption>**Table 5.** Calculated HFCs (including the $A^{FC}$ and $A^{SD}$ terms, in MHz) obtained with DLPNO-CCSD for 17 Mo(V) complexes of the benchmark set.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Isotropic</th><th colspan="3">Anisotropic</th></tr><tr><th>A<sub>11</sub></th><th>A<sub>22</sub></th><th>A<sub>33</sub></th><th>A<sup>FC</sup></th><th>A<sub>11</sub><sup>SD</sup></th><th>A<sub>22</sub><sup>SD</sup></th><th>A<sub>33</sub><sup>SD</sup></th></tr></thead><tbody><tr><td>1</td><td>127.6</td><td>115.2</td><td>42.7</td><td>95.2</td><td>32.5</td><td>20.0</td><td>-52.5</td></tr><tr><td>2</td><td>140.6</td><td>71.0</td><td>50.1</td><td>87.3</td><td>53.4</td><td>-16.2</td><td>-37.2</td></tr><tr><td>3</td><td>168.2</td><td>77.5</td><td>77.2</td><td>107.7</td><td>60.6</td><td>-30.1</td><td>-30.4</td></tr><tr><td>5</td><td>171.8</td><td>83.1</td><td>79.7</td><td>111.5</td><td>60.3</td><td>-28.4</td><td>-31.9</td></tr><tr><td>6</td><td>175.6</td><td>84.2</td><td>82.2</td><td>114.0</td><td>61.6</td><td>-29.8</td><td>-31.8</td></tr><tr><td>7</td><td>141.4</td><td>72.1</td><td>71.7</td><td>95.0</td><td>46.3</td><td>-23.0</td><td>-23.3</td></tr><tr><td>8</td><td>176.3</td><td>85.3</td><td>82.1</td><td>114.5</td><td>61.7</td><td>-29.3</td><td>-32.5</td></tr><tr><td>9</td><td>170.6</td><td>82.4</td><td>77.9</td><td>110.3</td><td>60.3</td><td>-27.9</td><td>-32.4</td></tr></tbody></table>
<table><tr><td><b>10</b></td><td>140.3</td><td>112.1</td><td>41.8</td><td>98.1</td><td>42.2</td><td>14.1</td><td>-56.3</td></tr><tr><td><b>11</b></td><td>192.0</td><td>96.3</td><td>91.3</td><td>126.5</td><td>65.5</td><td>-30.3</td><td>-35.2</td></tr><tr><td><b>14</b></td><td>192.6</td><td>93.4</td><td>93.1</td><td>126.3</td><td>66.2</td><td>-33.0</td><td>-33.3</td></tr><tr><td><b>15</b></td><td>186.9</td><td>93.3</td><td>90.2</td><td>123.4</td><td>63.4</td><td>-30.2</td><td>-33.3</td></tr><tr><td><b>17</b></td><td>192.3</td><td>96.4</td><td>94.8</td><td>127.9</td><td>64.4</td><td>-31.4</td><td>-33.0</td></tr><tr><td><b>18</b></td><td>191.1</td><td>89.2</td><td>86.3</td><td>122.2</td><td>68.9</td><td>-33.0</td><td>-35.9</td></tr><tr><td><b>19</b></td><td>194.7</td><td>97.3</td><td>94.1</td><td>128.7</td><td>66.0</td><td>-31.4</td><td>-34.6</td></tr><tr><td><b>20</b></td><td>218.5</td><td>113.4</td><td>112.9</td><td>148.3</td><td>70.2</td><td>-34.9</td><td>-35.3</td></tr><tr><td><b>22</b></td><td>216.9</td><td>113.2</td><td>112.3</td><td>147.5</td><td>69.4</td><td>-34.3</td><td>-35.1</td></tr></table>

To assess the performance of the DLPNO-CCSD against experiment, we can combine the $A^{FC}$ and $A^{SD}$ terms from DLPNO-CCSD with the spin-orbit coupling term ($A^{SO}$) derived from another method. Notably, the spin-orbit contributions to the total HFCs vary less than the Fermi contact and spin dipolar terms among the different functionals for all complexes of the benchmark set (see Figure 4 and Tables S6–S31). We define a composite approach denoted DLPNO-CCSD*, where the total HFCs are computed as a sum of the $A^{FC}$ and the $A^{SD}$ terms calculated with the DLPNO-CCSD method, and the $A^{SO}$ term provided by the best-performing DFT functional, PBE0-DH, so that $A_{ii}(\text{DLPNO-CCSD}^*) = A^{FC}(\text{DLPNO-CCSD}) + A_{ii}^{SD}(\text{DLPNO-CCSD}) + A_{ii}^{SO}(\text{PBE0-DH})$. Comparison of experimental and calculated HFC values using DLPNO-CCSD* and PBE0-DH are shown in Figures 6a and 6b, respectively, and numerical values for DLPNO-CCSD* are given in Table S32.

Figure 6a shows that the DLPNO-CCSD* values overestimate the HFC components in most cases. The MAPD($A_{iso}$) for the 17 complexes is 13.5%, with a SD of 6.4%, and the MAPD($A_{11}$) is 7.7%, with a SD of 4.8%. This places the approach in the mid-range of DFT methods plotted in Figure 3, specifically between PBE0 and ωB97. However, the DLPNO-CCSD* values have the lowest MAPD($A_{11} - A_{33}$) of 9.3% with a SD of 6.2%. Therefore, despite systematically overestimating the HFC values, DLPNO-CCSD* most accurately reproduces the anisotropy of the HF tensors.

In Figure 6c, the calculated isotropic HFC contributions ($A_{iso}$), derived as a sum of $A^{FC}$ and $A^{PC}$ terms obtained from DLPNO-CCSD and from PBE0-DH, respectively, are plotted against the experimental isotropic HFC values. The calculated values are larger than the experimental values for all compounds except **10** and **20**. The $A^{FC}$ contributions obtained from DLPNO-CCSD exceed those from PBE0-DH by ca. 3–22 MHz for all compounds except **10**, for which DLPNO-CCSD predicts a lower $A^{FC}$ contribution by 18.2 MHz. Hence, the overestimation of the $A_{iso}$ contribution by the composite DLPNO-CCSD* method results either from the
overestimation of the Fermi contact term by DLPNO-CCSD or from the overestimation of the DFT-derived pseudocontact term. The second case would imply that the success of PBE0-DH may be due to compensation between an overestimated pseudocontact term and an underestimated Fermi contact term, but this cannot be determined at present. Overall, at this stage we can conclude that $^{95}$Mo HFCs obtained with PBE0-DH are more reliable than those utilizing DLPNO-CCSD components.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Correlation of experimental and calculated $^{95}$Mo total HFCs from a composite DLPNO-CCSD* approach (a) and from PBE0-DH (b), and correlation of experimental and calculated isotropic HFCs from the composite DLPNO-CCSD* method (c).</figcaption></figure>

These results indicate that DLPNO-CCSD holds promise as a method for $^{95}$Mo HFCs only if future developments address missing spin-orbit components and current sources of error. Apart from the absence of second-order contributions, other factors that influence the accuracy of DLPNO-CCSD include the use of relaxed versus unrelaxed densities and the incorporation of triple excitation corrections (T), which has been shown to significantly impact the isotropic HFCs of some small transition metal complexes.[27] Addressing these aspects could further enhance the method's accuracy and reliability. Regarding the thresholds used for the DLPNO calculations to select the electron pairs and pair natural orbitals (PNOs) that will be treated at the canonical CCSD level, it is worth noting that the HFC1 settings employed here (see Computational Details) have been optimized specifically for HFC calculations and employ tight thresholds, especially in the core region.[73] At this point, and considering that other wave function-based approaches yield poor estimates for transition metal HFCs,[96] DFT remains the method of choice for accurately predicting molybdenum HFCs. Finally, it is important to note that the use of Kohn-Sham orbitals is imperative as initial input in these calculations; Hartree-Fock produces severe artifacts for the present
systems, and hence is inappropriate as reference, thus also precluding a meaningful
application of pure MP2 methods for the present quantities.

Comparison of DFT Methods for g-Tensors

In this section, we compare the performance of various DFT methods for g-tensors of Mo-
based compounds using the 22 complexes of our benchmark set (Tables S33-S60). The
experimentally measured g-shifts of our benchmark set span a wide range of 2–248 ppt in
absolute values. Consequently, we employ the mean absolute deviations (MADs) as the
evaluation criterion for comparing different methods, rather than MAPDs. The parameters
analyzed include the g tensor asymmetry (Δg₁₁ – Δg₃₃), the isotropic g value (Δgᵢₐ), and the
g-shift of the largest g-tensor component (Δg₁₁). The results are shown in Table 9 and
graphically represented in Figure 7. Here, the tensor asymmetry is considered as the primary
metric, since it is more important for practical applications and is affected to a lesser extent by
higher-order SOC effects, which tend to be systematic.

<table><caption>**Table 9.** MAD and standard deviations (SD in ppt) of Δg₁₁ – Δg₃₃, Δgᵢₐ and Δg₁₁ for calculated g-tensors using different DFT functionals with respect to experimental values for the 22 Mo(V) complexes of the benchmark set.</caption><thead><tr><th rowspan="2">Functional</th><th colspan="2">Δg<sub>11</sub> - Δg<sub>33</sub></th><th colspan="2">Δg<sub>iso</sub></th><th colspan="2">Δg<sub>11</sub></th></tr><tr><th>MAD</th><th>SD</th><th>MAD</th><th>SD</th><th>MAD</th><th>SD</th></tr></thead><tbody><tr><td><strong>GGA</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>BLYP</td><td>12</td><td>15</td><td>19</td><td>11</td><td>17</td><td>10</td></tr><tr><td>BP86</td><td>13</td><td>15</td><td>20</td><td>11</td><td>18</td><td>10</td></tr><tr><td>TPSS</td><td>16</td><td>18</td><td>21</td><td>13</td><td>17</td><td>11</td></tr><tr><td><strong>meta-GGA</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>M06-L</td><td>18</td><td>20</td><td>19</td><td>15</td><td>12</td><td>11</td></tr><tr><td>MN15-L</td><td>19</td><td>20</td><td>17</td><td>15</td><td>11</td><td>10</td></tr><tr><td>r<sup>2</sup>SCAN</td><td>20</td><td>21</td><td>21</td><td>15</td><td>13</td><td>12</td></tr><tr><td><strong>Hybrid</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>B3PW91 30%</td><td>7</td><td>7</td><td>9</td><td>8</td><td>11</td><td>9</td></tr><tr><td>PBE0</td><td>8</td><td>7</td><td>11</td><td>9</td><td>12</td><td>9</td></tr><tr><td>B3PW91 20%</td><td>8</td><td>8</td><td>12</td><td>9</td><td>13</td><td>10</td></tr><tr><td>B3LYP</td><td>8</td><td>8</td><td>12</td><td>9</td><td>13</td><td>10</td></tr><tr><td>B3PW91 40%</td><td>9</td><td>8</td><td>8</td><td>6</td><td>10</td><td>9</td></tr><tr><td>B3PW91 50%</td><td>11</td><td>11</td><td>7</td><td>5</td><td>11</td><td>9</td></tr><tr><td>TPSSh</td><td>13</td><td>15</td><td>18</td><td>12</td><td>14</td><td>11</td></tr><tr><td>M06</td><td>20</td><td>19</td><td>4</td><td>3</td><td>10</td><td>8</td></tr></tbody></table>
<table><tbody><tr><td>M05</td><td>26</td><td>29</td><td>6</td><td>5</td><td>8</td><td>7</td></tr><tr><td>M06-2X</td><td>40</td><td>61</td><td>17</td><td>14</td><td>11</td><td>11</td></tr><tr><td>MN15</td><td>69</td><td>102</td><td>25</td><td>31</td><td>11</td><td>10</td></tr><tr><td colspan="7"><strong>Range-separated</strong></td></tr><tr><td>CAM-B3LYP</td><td>8</td><td>8</td><td>9</td><td>7</td><td>12</td><td>9</td></tr><tr><td>LC-BLYP</td><td>8</td><td>8</td><td>10</td><td>7</td><td>14</td><td>9</td></tr><tr><td>ωB97</td><td>9</td><td>9</td><td>8</td><td>7</td><td>11</td><td>10</td></tr><tr><td>LC-PBE</td><td>10</td><td>10</td><td>8</td><td>6</td><td>11</td><td>11</td></tr><tr><td colspan="7"><strong>Double-hybrid</strong></td></tr><tr><td>PBE0-DH</td><td>7</td><td>6</td><td>9</td><td>7</td><td>10</td><td>8</td></tr><tr><td>B2PLYP</td><td>11</td><td>12</td><td>17</td><td>8</td><td>17</td><td>12</td></tr><tr><td>ωB2PLYP</td><td>13</td><td>12</td><td>8</td><td>5</td><td>11</td><td>10</td></tr><tr><td>B2GP-PLYP</td><td>15</td><td>20</td><td>16</td><td>8</td><td>19</td><td>16</td></tr><tr><td>ωB88PP86</td><td>16</td><td>17</td><td>15</td><td>6</td><td>18</td><td>13</td></tr><tr><td>Pr<sup>2</sup>SCAN50</td><td>23</td><td>51</td><td>18</td><td>9</td><td>14</td><td>10</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Graphical representation of the MADs of $\Delta g_{11} - \Delta g_{33}$ and $\Delta g_{iso}$ for the calculated values of the set of 22 Mo(V) complexes obtained with different functionals with respect to experimental values.</figcaption></figure>

The deviations from the experimental values are in general very small, unlike systems such as Cu systems for example where deviations tend to be more significant.[95, 100-101] However, accuracy is here more important to distinguish between different electronic structures. Overall, the results of the nine best-performing functionals lie in a narrow range with MAD($\Delta g_{11} - \Delta g_{33}$) of 7–9 ppt with SDs between 6 and 9 ppt. These functionals are PBE0-DH, B3PW91 with 30% HFX, PBE0, CAM-B3LYP, B3PW91 with 20% HFX, B2LYP, LC-BLYP, B3PW91 with 20% HFX, and $\omega$B97. Among them, PBE0-DH is the only double-hybrid, whereas all the rest are hybrid functionals. Therefore, the best-performing functional
for the prediction of <sup>95</sup>Mo HFCs (Figure 3) is also the best-performing one for g tensor prediction (Figure 7). PBE0-DH has the smallest MAD(Δg<sub>11</sub> − Δg<sub>33</sub>) of 7 ppt with the smallest SD of 6 ppt. Most of the remaining functionals also show relatively small MADs with values less than 20 ppt, but their reliability is further limited due to larger standard deviations of 10–20 ppt. Only the worst-performing functionals Pr<sup>2</sup>SCAN50, M06, M06-2X, and MN15 show MAD values that exceed 20 ppt and standard deviations above 50 ppt.

Significant differences are observed in the ranking of methods based on their performance for g-tensors and for molybdenum HFC calculations. Among the double-hybrid functionals, the standout performer is PBE0-DH, which excels in both properties. In contrast, B2PLYP ranks in the middle of the plot, while all other double-hybrid functionals fall in the bottom half. Additionally, all range-separated hybrid functionals perform very well for g-tensors, whereas for HFC calculations, LC-BLYP and LC-PBE exhibit larger deviations compared to most other methods, which highlights the different demands of these properties. Here, the optimal HFX percentage for the hybrid functional B3PW91 is 30%, consistent with the fact that g-tensors require lower HFX percentages than metal HFCs, due to the errors from spin contamination.[26]

The GGA and the Minnesota functionals collectively exhibit the poorest performances for both g-tensors and HFCs. The GGA functionals tend to generally and systematically overestimate g-tensor components. In contrast, the nature of the deviations observed with the hybrid Minnesota functionals (M05, M06, M06-2X, and MN15) does not appear to be systematic, which is also reflected on the large standard deviations, which exceed 29 ppt particularly for M05, M06-2X, and MN15. Notably, M05 and M06 exemplify a remarkable error compensation, achieving the smallest MAD(Δg<sub>iso</sub>) of only 4 and 6 ppt, respectively. This metric does not align with their overall performance, as indicated by their significantly larger MAD(Δg<sub>11</sub> − Δg<sub>33</sub>) of 26 and 20 ppt, respectively, and SDs, which place them among the worst-performing functionals.

Overall, PBE0-DH stands out as the only double-hybrid functional that performs optimally for both HFCs and g-tensor calculations of Mo complexes. Notably, this functional is also among the top-performers for the particularly challenging problem of copper EPR parameter calculations.[101] It remains difficult to identify the physical origin of its success right now, or to directly relate the performance to any specific construction feature of these functionals. A possible reason might be the relatively smaller HFX of PBE0-DH, 50%, and at the same time
the small MP2 (12.5%) contribution, as discussed earlier. This directs toward an optimal
construction of double-hybrid functionals for spin-dependent properties.
Conclusions

In this work, we compiled an extensive and varied reference set of Mo compounds with reliable EPR data, and we investigated theoretical methodologies for the calculation of metal (⁹⁵Mo) hyperfine coupling constants and *g*-values. Our investigation began by determining a suitable basis set for Mo, and we found that the SARC-DKH-TZVP basis set can provide converged HFCs and *g*-values without modifications such as decontraction of the core orbitals or augmentation with additional tight *s* functions. Using this basis set, we conducted a balanced evaluation of quantum chemical methodologies, including DFT and wavefunction-based approaches. The assessment encompassed a diverse and representative group of DFT families of functionals, including GGA, meta-GGA, hybrid, range-separated, as well as double-hybrid functionals. Among the approaches examined herein, double-hybrid and range-separated functionals were tested for ⁹⁵Mo HFCs and *g*-tensors for the first time.

Our results indicated that the double-hybrid functional PBE0-DH is the best-performing DFT method for both HFCs and *g*-tensor calculations, as it provides the lowest mean absolute deviations from the experimental data for both properties. It was also shown that other double-hybrid functionals such as ωB2PLYP, Pr²SCAN, ωB88PP86, and B2GP-PLYP also provide small deviations from the experimental HFC values. However, these functionals provide poor estimations of *g*-tensors. Therefore, only PBE0-DH succeeds in both properties simultaneously. Regarding hybrid functionals, consistently with previous observations, the HFX of 40% is optimal for HFCs calculations, whereas for *g*-tensors, a lower HFX of 30% is preferable. Nevertheless, the B3PW91 with HFX 30-50% provides highly accurate and reliable *g*-tensors. Thus, the B3PW91 functional with 40% HFX can be proposed as the second-best method for both properties. The range-separated functionals achieved accuracy comparable to hybrid functionals with low HFX percentage. They yield poor estimations for HFCs, but perform remarkably well for *g*-tensors. In line with the well-known deficiencies of GGA and meta-GGA functionals that lead them to overestimate the covalency of M-L bonds, they were shown to systematically overestimate *g*-tensor components and underestimate HFCs with respect to experimental values. The functionals of the Minnesota group cannot deliver reliable predictions for any of the EPR properties examined, and the deviations from the experimental values are not systematic, especially for the hybrid functionals M05, M06, M06-2X, and MN15. These results for the present *d*¹ systems are largely consistent with previously reported benchmark studies on EPR parameters of *d*⁹ Cu(II) complexes.[⁹⁶, ¹⁰¹] However, they should not be considered generalizable to all transition metal cases and other
electron configurations, which may be affected by issues such as high spin contamination.[26]

Ultimately, local hybrid approaches that incorporate position-dependent HFX admixture may emerge as the most general practical solution to the problem of balancing the description of core and valence spin polarization,[26, 102] while other promising approaches based on the post-Kohn-Sham random phase approximation (RPA) methods[103] also require further investigation in large transition metal complexes.

The performance of a wavefunction-based method, the local variant of coupled-cluster, DLPNO-CCSD, was also assessed. Currently, this approach provides only the Fermi contact and spin-dipolar terms of the HFCs. To compare the HFCs obtained with DLPNO-CCSD with the experimental values, we used the spin-orbit coupling contribution obtained from the PBE0-DH functional that achieves optimal performance. The accuracy of this composite method is mediocre and, hence, it does not compete with the performance of double-hybrid and hybrid functionals with HFX percentages higher than 30%.

Overall, the PBE0-DH functional is suggested as the best approach for robust and reliable calculations of EPR parameters of Mo complexes that can be used in bioinorganic chemistry applications on enzymes. Further investigation of the construction parameters of PBE0-DH that makes it so accurate would be interesting for future developments.
**Acknowledgment.** The authors acknowledge support by the Max Planck Society. M.D. acknowledges support from the Alexander von Humboldt Stiftung (Humboldt Research Fellowship) and from the Fonds der Chemischen Industrie (Liebig Fellowship). The authors gratefully acknowledge financial support of this work by the French National Research Agency (CODEC, grant no. ANR-19-CE05_0030_01) and the Deutsche Forschungsgemeinschaft (CUBISM, grant no. ANR-18 CE092_0040_01 / DFG project no. 406697875), and from the France-Germany Hubert Curien Program – German Academic Exchange Service (DAAD) (Procope project 42525PB / DAAD project 57445526 / Procope Plus project no. 0185-DEU-22-001 LG 4).

**Additional information**

Supplementary information is available in the online version of the paper.

Complete original data are provided in Edmond, the Open Research Data Repository of the
Max Planck Society, at: https://doi.org/10.17617/3.CVNSRT.

**Competing financial interests**

The authors declare no competing financial interests.
References

[1] Y. Zhang, S. Rump, V. N. Gladyshev, Comparative genomics and evolution of molybdenum utilization, *Coord. Chem. Rev.* **2011**, *255*, 1206-1217.

[2] Y. Zhang, V. N. Gladyshev, Molybdoproteomes and Evolution of Molybdenum Utilization, *J. Mol. Biol.* **2008**, *379*, 881-899.

[3] R. Hille, J. Hall, P. Basu, The Mononuclear Molybdenum Enzymes, *Chem. Rev.* **2014**, *114*, 3963-4038.

[4] M. C. Metzger, P. Basu, in *Microbial Metabolism of Metals and Metalloids* (Ed.: C. J. Hurst), Springer International Publishing, Cham, **2022**, pp. 359-415.

[5] O. Einsle, D. C. Rees, Structural Enzymology of Nitrogenase Enzymes, *Chem. Rev.* **2020**, *120*, 4969-5004.

[6] H. Dobbek, L. Gremer, R. Kiefersauer, R. Huber, O. Meyer, Catalysis at a dinuclear [CuSMo(=O)OH] cluster in a CO dehydrogenase resolved at 1.1-Å resolution, *Proc. Natl. Acad. Sci. U. S. A.* **2002**, *99*, 15971-15976.

[7] S. T. Stripp, B. R. Duffus, V. Fourmond, C. Léger, S. Leimkühler, S. Hirota, Y. Hu, A. Jasniewski, H. Ogata, M. W. Ribbe, Second and Outer Coordination Sphere Effects in Nitrogenase, Hydrogenase, Formate Dehydrogenase, and CO Dehydrogenase, *Chem. Rev.* **2022**, *122*, 11900-11973.

[8] T. Hartmann, N. Schwanhold, S. Leimkühler, Assembly and catalysis of molybdenum or tungsten-containing formate dehydrogenases from bacteria, *BBA - Proteins and Proteomics* **2015**, *1854*, 1090-1100.

[9] L. B. Maia, I. Moura, J. J. G. Moura, in *Enzymes for Solving Humankind's Problems: Natural and Artificial Systems in Health, Agriculture, Environment and Energy* (Eds.: J. J. G. Moura, I. Moura, L. B. Maia), Springer International Publishing, Cham, **2021**, pp. 29-81.

[10] A. M. Appel, J. E. Bercaw, A. B. Bocarsly, H. Dobbek, D. L. DuBois, M. Dupuis, J. G. Ferry, E. Fujita, R. Hille, P. J. A. Kenis, C. A. Kerfeld, R. H. Morris, C. H. F. Peden, A. R. Portis, S. W. Ragsdale, T. B. Rauchfuss, J. N. H. Reek, L. C. Seefeldt, R. K. Thauer, G. L. Waldrop, Frontiers, Opportunities, and Challenges in Biochemical and Chemical Catalysis of CO₂ Fixation, *Chem. Rev.* **2013**, *113*, 6621-6658.

[11] S. J. N. Burgmayer, M. L. Kirk, Advancing Our Understanding of Pyranopterin-Dithiolene Contributions to Moco Enzyme Catalysis, *Molecules* **2023**, *28*, 7456.

[12] C. G. Young, Chemical systems modeling the d¹ Mo(V) states of molybdenum enzymes, *J. Inorg. Biochem.* **2016**, *162*, 238-252.

[13] M. Orio, D. A. Pantazis, Successes, challenges, and opportunities for quantum chemistry in understanding metalloenzymes for solar fuels research, *Chem. Commun.* **2021**, *57*, 3952-3974.

[14] V. Barone, S. Alessandrini, M. Biczysko, J. R. Cheeseman, D. C. Clary, A. B. McCoy, R. J. DiRisio, F. Neese, M. Melosso, C. Puzzarini, Computational molecular spectroscopy, *Nat. Rev. Methods Primers* **2021**, *1*, 38.

[15] L. B. Maia, I. Moura, J. J. G. Moura, in *Future Directions in Metalloprotein and Metalloenzyme Research* (Eds.: G. Hanson, L. Berliner), Springer International Publishing, Cham, **2017**, pp. 55-101.

[16] M. M. Cosper, F. Neese, A. V. Astashkin, M. D. Carducci, A. M. Raitsimring, J. H. Enemark, Determination of the g-Tensors and Their Orientations for cis,trans-(L-N₂S₂)MoVOX (X = Cl, SCH₂Ph) by Single-Crystal EPR Spectroscopy and Molecular Orbital Calculations, *Inorg. Chem.* **2005**, *44*, 1290-1301.

[17] J. Fritscher, P. Hrobárik, M. Kaupp, Computational Studies of Electron Paramagnetic Resonance Parameters for Paramagnetic Molybdenum Complexes. 1. Method Validation on Small and Medium-Sized Systems, *J. Phys. Chem. B* **2007**, *111*, 4616-4629.

[18] S. C. Drew, C. G. Young, G. R. Hanson, A Density Functional Study of the Electronic Structure and Spin Hamiltonian Parameters of Mononuclear Thiomolybdenyl Complexes, *Inorg. Chem.* **2007**, *46*, 2388-2397.

[19] J. Fritscher, P. Hrobárik, M. Kaupp, Computational Studies of EPR Parameters for Paramagnetic Molybdenum Complexes. II. Larger MoV Systems Relevant to Molybdenum Enzymes, *Inorg. Chem.* **2007**, *46*, 8146-8161.

[20] P. Hrobárik, O. L. Malkina, V. G. Malkin, M. Kaupp, Relativistic two-component calculations of electronic g-tensor for oxo-molybdenum(V) and oxo-tungsten(V) complexes: The important role of higher-order spin-orbit contributions, *Chem. Phys. 2009*, **356**, 229-235.

[21] R. G. Hadt, V. N. Nemykin, J. G. Olsen, P. Basu, Comparative calculation of EPR spectral parameters in $[MoVOX_4]$, $[MoVOX_5]^{2-}$, and $[MoVOX_4(H_2O)]^-$ complexes, *Phys. Chem. Chem. Phys.* **2009**, *11*, 10377-10384.

[22] S. Gohr, P. Hrobárik, M. Repiský, S. Komorovský, K. Ruud, M. Kaupp, Four-Component Relativistic Density Functional Theory Calculations of EPR g- and Hyperfine-Coupling Tensors Using Hybrid Functionals: Validation on Transition-Metal Complexes with Large Tensor Anisotropies and Higher-Order Spin-Orbit Effects, *J. Phys. Chem. A* **2015**, *119*, 12892-12905.
[23] A. Wodyński, M. Kaupp, Density Functional Calculations of Electron Paramagnetic Resonance g- and Hyperfine-Coupling Tensors Using the Exact Two-Component (X2C) Transformation and Efficient Approximations to the Two-Electron Spin-Orbit Terms, *J. Phys. Chem. A* **2019**, *123*, 5660-5672.

[24] A. Wodyński, M. Kaupp, Noncollinear Relativistic Two-Component X2C Calculations of Hyperfine Couplings Using Local Hybrid Functionals. Importance of the High-Density Coordinate Scaling Limit, *J. Chem. Theory Comput.* **2020**, *16*, 314-325.

[25] V. N. Nemykin, J. R. Sabin, B. W. Kail, A. Upadhyay, M. P. Hendrich, P. Basu, Influence of the ligand-field on EPR parameters of cis- and trans-isomers in MoV systems relevant to molybdenum enzymes: Experimental and density functional theory study, *J. Inorg. Biochem.* **2023**, *245*, 112228.

[26] C. J. Schattenberg, T. M. Maier, M. Kaupp, Lessons from the Spin-Polarization/Spin-Contamination Dilemma of Transition-Metal Hyperfine Couplings for the Construction of Exchange-Correlation Functionals, *J. Chem. Theory Comput.* **2018**, *14*, 5653-5672.

[27] M. Munzarová, M. Kaupp, A Critical Validation of Density Functional and Coupled-Cluster Approaches for the Calculation of EPR Hyperfine Coupling Constants in Transition Metal Complexes, *J. Phys. Chem. A* **1999**, *103*, 9966-9983.

[28] M. L. Munzarová, P. Kubáček, M. Kaupp, Mechanisms of EPR Hyperfine Coupling in Transition Metal Complexes, *J. Am. Chem. Soc.* **2000**, *122*, 11900-11913.

[29] D. Peng, N. Middendorf, F. Weigend, M. Reiher, An efficient implementation of two-component relativistic exact-decoupling methods for large molecules, *J. Chem. Phys.* **2013**, *138*, 184105.

[30] Y. J. Franzke, N. Middendorf, F. Weigend, Efficient implementation of one- and two-component analytical energy gradients in exact two-component theory, *J. Chem. Phys.* **2018**, *148*, 104110.

[31] Y. J. Franzke, J. M. Yu, Hyperfine Coupling Constants in Local Exact Two-Component Theory, *J. Chem. Theory Comput.* **2022**, *18*, 323-343.

[32] F. Neese, Software update: The ORCA program system—Version 5.0, *WIREs Comput. Mol. Sci.* **2022**, *12*, e1606.

[33] V. N. Staroverov, G. E. Scuseria, J. Tao, J. P. Perdew, Comparative assessment of a new nonempirical density functional: Molecules and hydrogen-bonded complexes, *J. Chem. Phys.* **2003**, *119*, 12129-12137.

[34] P. Pollak, F. Weigend, Segmented Contracted Error-Consistent Basis Sets of Double- and Triple-ζ Valence Quality for One- and Two-Component Relativistic All-Electron Calculations, *J. Chem. Theory Comput.* **2017**, *13*, 3696-3705.

[35] F. Neese, An improvement of the resolution of the identity approximation for the formation of the Coulomb matrix, *J. Comput. Chem.* **2003**, *24*, 1740-1747.

[36] G. L. Stoychev, A. A. Auer, F. Neese, Automatic Generation of Auxiliary Basis Sets, *J. Chem. Theory Comput.* **2017**, *13*, 554-562.

[37] C. R. Groom, I. J. Bruno, M. P. Lightfoot, S. C. Ward, The Cambridge Structural Database, *Acta. Crystallogr. B* **2016**, *72*, 171-179.

[38] J. P. Perdew, Density-functional approximation for the correlation energy of the inhomogeneous electron gas, *Phys. Rev. B Condens. Matter.* **1986**, *33*, 8822-8824.

[39] A. D. Becke, Density-functional exchange-energy approximation with correct asymptotic behavior, *Phys. Rev. A* **1988**, *38*, 3098-3100.

[40] C. Lee, W. Yang, R. G. Parr, Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density, *Phys. Rev. B* **1988**, *37*, 785-789.

[41] H. S. Yu, X. He, D. G. Truhlar, MN15-L: A New Local Exchange-Correlation Functional for Kohn-Sham Density Functional Theory with Broad Accuracy for Atoms, Molecules, and Solids, *J. Chem. Theory Comput.* **2016**, *12*, 1280-1293.

[42] J. Tao, J. P. Perdew, V. N. Staroverov, G. E. Scuseria, Climbing the Density Functional Ladder: Nonempirical Meta-Generalized Gradient Approximation Designed for Molecules and Solids, *Phys. Rev. Lett.* **2003**, *91*, 146401.

[43] Y. Zhao, D. G. Truhlar, A new local density functional for main-group thermochemistry, transition metal bonding, thermochemical kinetics, and noncovalent interactions, *J. Chem. Phys.* **2006**, *125*, 194101.

[44] J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, J. Sun, Accurate and Numerically Efficient r²SCAN Meta-Generalized Gradient Approximation, *J. Phys. Chem. Lett.* **2020**, *11*, 8208-8215.

[45] A. D. Becke, Density- functional thermochemistry. III. The role of exact exchange, *J. Chem. Phys.* **1993**, *98*, 5648-5652.

[46] C. Adamo, V. Barone, Toward reliable density functional methods without adjustable parameters: The PBE0 model, *J. Chem. Phys.* **1999**, *110*, 6158-6170.

[47] A. D. Becke, A new mixing of Hartree-Fock and local density- functional theories, *J. Chem. Phys.* **1993**, *98*, 1372-1377.

[48] J. P. Perdew, Y. Wang, Accurate and simple analytic representation of the electron-gas correlation energy, *Phys. Rev. B* **1992**, *45*, 13244-13249.
[49] J. P. Perdew, Generalized gradient approximations for exchange and correlation: A look backward and forward, *Physica B: Condensed Matter* **1991**, *172*, 1-6.

[50] Y. Zhao, N. E. Schultz, D. G. Truhlar, Exchange-correlation functional with broad accuracy for metallic and nonmetallic compounds, kinetics, and noncovalent interactions, *J. Chem. Phys.* **2005**, *123*, 161103.

[51] Y. Zhao, D. G. Truhlar, The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals, *Theor. Chem. Acc.* **2008**, *120*, 215-241.

[52] H. S. Yu, X. He, S. L. Li, D. G. Truhlar, MN15: A Kohn–Sham global-hybrid exchange–correlation density functional with broad accuracy for multi-reference and single-reference systems and noncovalent interactions, *Chem. Sci.* **2016**, *7*, 5032-5051.

[53] J.-D. Chai, M. Head-Gordon, Long-range corrected hybrid density functionals with damped atom–atom dispersion corrections, *Phys. Chem. Chem. Phys.* **2008**, *10*, 6615.

[54] T. Yanai, D. P. Tew, N. C. Handy, A new hybrid exchange-correlation functional using the Coulomb-attenuating method (CAM-B3LYP), *Chem. Phys. Lett.* **2004**, *393*, 51-57.

[55] H. Iikura, T. Tsuneda, T. Yanai, K. Hirao, A long-range correction scheme for generalized-gradient-approximation exchange functionals, *J. Chem. Phys.* **2001**, *115*, 3540-3544.

[56] S. Grimme, Semiempirical hybrid density functional with perturbative second-order correlation, *J. Chem. Phys.* **2006**, *124*, 034108.

[57] E. Brémond, C. Adamo, Seeking for parameter-free double-hybrid functionals: The PBE0-DH model, *J. Chem. Phys.* **2011**, *135*, 024106.

[58] A. Karton, A. Tarnopolsky, J.-F. Lamère, G. C. Schatz, J. M. L. Martin, Highly Accurate First-Principles Benchmark Data Sets for the Parametrization and Validation of Density Functional and Other Approximate Methods. Derivation of a Robust, Generally Applicable, Double-Hybrid Functional for Thermochemistry and Thermochemical Kinetics, *J. Phys. Chem. A* **2008**, *112*, 12868-12886.

[59] M. Casanova-Páez, M. B. Dardis, L. Goerigk, $\omega$B2PLYP and $\omega$B2GPPLYP: The First Two Double-Hybrid Density Functionals with Long-Range Correction Optimized for Excitation Energies, *J. Chem. Theory Comput.* **2019**, *15*, 4735-4744.

[60] M. Casanova-Páez, L. Goerigk, Time-Dependent Long-Range-Corrected Double-Hybrid Density Functionals with Spin-Component and Spin-Opposite Scaling: A Comprehensive Analysis of Singlet–Singlet and Singlet–Triplet Excitation Energies, *J. Chem. Theory Comput.* **2021**, *17*, 5165-5186.

[61] L. Wittmann, H. Neugebauer, S. Grimme, M. Bursch, Dispersion-corrected $r^2$SCAN based double-hybrid functionals, *J. Chem. Phys.* **2023**, *159*, 224103.

[62] S. Lehtola, C. Steigemann, M. J. T. Oliveira, M. A. L. Marques, Recent developments in LIBXC — A comprehensive library of functionals for density functional theory, *SoftwareX* **2018**, *7*, 1-5.

[63] E. van Lenthe, E. J. Baerends, J. G. Snijders, Relativistic regular two-component Hamiltonians, *J. Chem. Phys.* **1993**, *99*, 4597-4610.

[64] E. van Lenthe, E. J. Baerends, J. G. Snijders, Relativistic total energy using regular approximations, *J. Chem. Phys.* **1994**, *101*, 9783-9792.

[65] E. van Lenthe, J. G. Snijders, E. J. Baerends, The zero- order regular approximation for relativistic effects: The effect of spin–orbit coupling in closed shell molecules, *J. Chem. Phys.* **1996**, *105*, 6505-6516.

[66] F. Weigend, R. Ahlrichs, Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy, *Phys. Chem. Chem. Phys.* **2005**, *7*, 3297-3305.

[67] J. D. Rolfes, F. Neese, D. A. Pantazis, All-electron scalar relativistic basis sets for the elements Rb–Xe, *J. Comput. Chem.* **2020**, *41*, 1842-1849.

[68] F. Neese, Prediction and interpretation of the $^{57}$Fe isomer shift in Mössbauer spectra by density functional theory, *Inorg. Chim. Acta* **2002**, *337*, 181-192.

[69] E. Malkin, I. Malkin, O. L. Malkina, V. G. Malkin, M. Kaupp, Scalar relativistic calculations of hyperfine coupling tensors using the Douglas–Kroll–Hess method with a finite-size nucleus model, *Phys. Chem. Chem. Phys.* **2006**, *8*, 4079-4085.

[70] L. Visscher, K. G. Dyall, Dirac–Fock atomic electronic structure calculations using different nuclear charge distributions, *Atom. Data Nucl. Data Tabl.* **1997**, *67*, 207-224.

[71] F. Neese, Efficient and accurate approximations to the molecular spin-orbit coupling operator and their use in molecular g-tensor calculations, *J. Chem. Phys.* **2005**, *122*, 034107.

[72] D. Datta, S. Kossmann, F. Neese, Analytic energy derivatives for the calculation of the first-order molecular properties using the domain-based local pair-natural orbital coupled-cluster theory, *J. Chem. Phys.* **2016**, *145*, 114101-114118.

[73] M. Saitow, F. Neese, Accurate spin-densities based on the domain-based local pair-natural orbital coupled-cluster theory, *J. Chem. Phys.* **2018**, *149*, 034104.
[74] M. Witwicki, Overcoming Challenges in Density Functional Theory-Based Calculations of Hyperfine Coupling Constants for Heavy Heteroatom Radicals, *ChemPhysChem* **2025**, e202400978.

[75] A. Cervilla, E. Llopis, D. Marco, F. Pérez, X-ray Structure of (BuₕN)[Mo(1,2-Benzenedithiolate)₃]. Trigonal-Prismatic versus Octahedral Coordination in Tris(1,2-Benzenedithiolate) Complexes, *Inorg. Chem.* **2001**, *40*, 6525-6528.

[76] J. R. Bradbury, M. F. Mackay, A. G. Wedd, The crystal and molecular structure of tetraphenylarsonium Tetrakis(benzenethiolato)oxomolybdate(V), *Aust. J. Chem.* **1978**, *31*, 2423-2430.

[77] S. Boyde, S. R. Ellis, C. D. Garner, W. Clegg, Structural comparison of oxobis(benzene-1,2-dithiolato)molybdenum-(V) and -(IV) complexes, *J. Chem. Soc., Chem. Commun.* **1986**, 1541-1543.

[78] K. Peariso, B. S. Chohan, C. J. Carrano, M. L. Kirk, Synthesis and EPR Characterization of New Models for the One-Electron Reduced Molybdenum Site of Sulfite Oxidase, *Inorg. Chem.* **2003**, *42*, 6194-6203.

[79] M. L. Mader, M. D. Carducci, J. H. Enemark, Analogues for the Molybdenum Center of Sulfite Oxidase: Oxomolybdenum(V) Complexes with Three Thiolate Sulfur Donor Atoms, *Inorg. Chem.* **2000**, *39*, 525-531.

[80] S. C. Drew, J. P. Hill, I. Lane, G. R. Hanson, R. W. Gable, C. G. Young, Synthesis, Structural Characterization, and Multifrequency Electron Paramagnetic Resonance Studies of Mononuclear Thiomolybdenyl Complexes, *Inorg. Chem.* **2007**, *46*, 2373-2387.

[81] B. S. Lim, M. W. Willer, M. Miao, R. H. Holm, Monodithiolene Molybdenum(V,VI) Complexes: A Structural Analogue of the Oxidized Active Site of the Sulfite Oxidase Enzyme Family, *J. Am. Chem. Soc.* **2001**, *123*, 8343-8349.

[82] C. D. Garner, L. H. Hill, F. E. Mabbs, D. L. McFadden, A. T. McPhail, Crystal and molecular structure, electron spin resonance, and electronic spectrum of tetraphenylarsonium tetrachloro-oxomolybdenum(V), *J. Chem. Soc., Dalton Trans.* **1977**, 853-858.

[83] S. Sproules, P. Banerjee, T. Weyhermüller, Y. Yan, J. P. Donahue, K. Wieghardt, Monoanionic Molybdenum and Tungsten Tris(dithiolene) Complexes: A Multifrequency EPR Study, *Inorg. Chem.* **2011**, *50*, 7106-7122.

[84] A. V. Astashkin, F. Neese, A. M. Raitsimring, J. J. A. Cooney, E. Bultman, J. H. Enemark, Pulsed EPR Investigations of Systems Modeling Molybdenum Enzymes: Hyperfine and Quadrupole Parameters of Oxo-¹⁷O in [Mo¹⁷O(SPh₄)], *J. Am. Chem. Soc.* **2005**, *127*, 16713-16722.

[85] N. Ueyama, H. Oku, M. Kondo, T.-a. Okamura, N. Yoshinaga, A. Nakamura, Trans Influence of Oxo and Dithiolene Coordination in Oxidized Models of Molybdenum Oxidoreductase: Synthesis, Structures, and Properties of Q₂[Mo VI O₂(1,2-benzenedithiolato)₂] (Q = NEt₄, PPh₄) and Related Complexes, *Inorg. Chem.* **1996**, *35*, 643-650.

[86] G. L. Wilson, R. J. Greenwood, J. R. Pilbrow, J. T. Spence, A. G. Wedd, Molybdenum(V) sites in xanthine oxidase and relevant analog complexes: comparison of molybdenum-95 and sulfur-33 hyperfine coupling, *J. Am. Chem. Soc.* **1991**, *113*, 6803-6812.

[87] G. Peng, J. Nichols, E. A. McCullough, Jr., J. T. Spence, Models for the Molybdenum(VI/V) Centers of the Molybdenum Hydroxylases and Related Enzymes. Geometry, Electronic Structure, and EPR g-Tensor Predictions from ab Initio and Semiempirical Molecular Orbital Studies, *Inorg. Chem.* **1994**, *33*, 2857-2864.

[88] I. K. Dhawan, J. H. Enemark, EPR Studies of Oxo-Molybdenum(V) Complexes with Sulfur Donor Ligands: Implications for the Molybdenum Center of Sulfite Oxidase, *Inorg. Chem.* **1996**, *35*, 4873-4882.

[89] K. K. Sunil, M. T. Rogers, ESR studies of some oxotetrahalo complexes of vanadium(IV) and molybdenum(V), *Inorg. Chem.* **1981**, *20*, 3283-3287.

[90] M. Bühl, C. Reimann, D. A. Pantazis, T. Bredow, F. Neese, Geometries of Third-Row Transition-Metal Complexes from Density-Functional Theory, *J. Chem. Theory Comput.* **2008**, *4*, 1449-1459.

[91] D. A. Pantazis, X.-Y. Chen, C. R. Landis, F. Neese, All-Electron Scalar Relativistic Basis Sets for Third-Row Transition Metal Atoms, *J. Chem. Theory Comput.* **2008**, *4*, 908-919.

[92] D. A. Pantazis, F. Neese, All-Electron Scalar Relativistic Basis Sets for the Lanthanides, *J. Chem. Theory Comput.* **2009**, *5*, 2229-2238.

[93] D. A. Pantazis, F. Neese, All-Electron Scalar Relativistic Basis Sets for the Actinides, *J. Chem. Theory Comput.* **2011**, *7*, 677-684.

[94] D. A. Pantazis, F. Neese, All-electron scalar relativistic basis sets for the 6p elements, *Theor. Chem. Acc.* **2012**, *131*, 1292.

[95] G. Sciortino, G. Lubinu, J.-D. Maréchal, E. Garribba, DFT Protocol for EPR Prediction of Paramagnetic Cu(II) Complexes and Application to Protein Binding Sites, *Magnetochemistry* **2018**, *4*, 55.

[96] R. J. Gómez-Piñeiro, D. A. Pantazis, M. Orio, Comparison of Density Functional and Correlated Wave Function Methods for the Prediction of Cu(II) Hyperfine Coupling Constants, *ChemPhysChem* **2020**, *21*, 2667-2679.
[97] N. Mardirossian, M. Head-Gordon, How Accurate Are the Minnesota Density Functionals for Noncovalent Interactions, Isomerization Energies, Thermochemistry, and Barrier Heights Involving Molecules Composed of Main-Group Elements?, *J. Chem. Theory Comput.* **2016**, *12*, 4303-4325.

[98] S. Kossmann, B. Kirchner, F. Neese, Performance of modern density functional theory for the prediction of hyperfine structure: meta-GGA and double hybrid functionals, *Mol. Phys.* **2007**, *105*, 2049-2071.

[99] M. G. Medvedev, I. S. Bushmarinov, J. Sun, J. P. Perdew, K. A. Lyssenko, Density functional theory is straying from the path toward the exact functional, *Science* **2017**, *355*, 49-52.

[100] R. J. Gómez-Piñeiro, M. Drosou, S. Bertaina, C. Decroos, A. J. Simaan, D. A. Pantazis, M. Orio, Decoding the Ambiguous Electron Paramagnetic Resonance Signals in the Lytic Polysaccharide Monooxygenase from *Photorhabdus luminescens*, *Inorg. Chem.* **2022**, *61*, 8022-8035.

[101] M. Drosou, C. A. Mitsopoulou, M. Orio, D. A. Pantazis, EPR Spectroscopy of Cu(II) Complexes: Prediction of *g*-Tensors Using Double-Hybrid Density Functional Theory, *Magnetochemistry* **2022**, *8*, 36.

[102] A. Wodyński, B. Lauw, M. Reimann, M. Kaupp, Spin-Symmetry Breaking and Hyperfine Couplings in Transition-Metal Complexes Revisited Using Density Functionals Based on the Exact-Exchange Energy Density, *J. Chem. Theory Comput.* **2024**, *20*, 2033-2048.

[103] F. Bruder, F. Weigend, Y. J. Franzke, Application of the Adiabatic Connection Random Phase Approximation to Electron-Nucleus Hyperfine Coupling Constants, *J. Phys. Chem. A* **2024**, *128*, 7298-7310.
TOC GRAPHIC

<figure><img src="image_9.png" /></figure>

Using a curated database of experimentally characterized, biologically relevant Mo(V) compounds, we evaluate the optimal choice of basis sets, density functionals, and wavefunction methods for accurately computing *g*-tensor values and ⁹⁵Mo hyperfine coupling constants, with the goal of establishing reliable protocols for predicting EPR of Molybdenum complexes and enzyme cofactors.