EPR Spectroscopy of Cu(II) Complexes: Prediction of g-Tensors Using Double-Hybrid Density Functional Theory

Computational electron paramagnetic resonance (EPR) spectroscopy is an important field of applied quantum chemistry that contributes greatly to connecting spectroscopic observations with the fundamental description of electronic structure for open-shell molecules. However, not all EPR parameters can be predicted accurately and reliably for all chemical systems. Among transition metal ions, Cu(II) centers in inorganic chemistry and biology, and their associated EPR properties such as hyperfine coupling and g-tensors, pose exceptional difficulties for all levels of quantum chemistry. In the present work, we approach the problem of Cu(II) g-tensor calculations using double-hybrid density functional theory (DHDFT). Using a reference set of 18 structurally and spectroscopically characterized Cu(II) complexes, we evaluate a wide range of modern double-hybrid density functionals (DHDFs) that have not been applied previously to this problem. Our results suggest that the current generation of DHDFs consistently and systematically outperform other computational approaches. The B2GP-PLYP and PBE0-DH functionals are singled out as the best DHDFs on average for the prediction of Cu(II) g-tensors. The performance of the different functionals is discussed and suggestions are made for practical applications and future methodological developments.


Introduction
Copper is among the most abundant transition metals in bioinorganic systems [1][2][3][4][5][6][7]. Electron paramagnetic resonance (EPR) spectroscopy is a uniquely useful approach for Cu(II) systems because it provides direct insights into the electronic and geometric structure of Cu(II) centers and the nature of their ligand sphere [4,8]. Relationships between structure and EPR parameters have been previously discussed based on ligand-field theory, as well as experimental correlations [8][9][10][11][12][13]. However, such correlations are often hard to extend to complex structural cases such as those frequently encountered in enzymatic active sites. Therefore, in most cases, interpretation of EPR spectroscopy requires input from quantum chemical calculations that correlate expected EPR parameters with molecular electronic and geometric structure [14][15][16][17][18][19][20][21][22][23][24]. A primary target property for Cu(II) systems is the g-matrix (commonly referred to as g-tensor), which describes the interaction of the unpaired electron with the effective magnetic field. The limitations of classical density functional theory (DFT) methods to accurately predict g-tensors of copper systems have been recognized in multiple studies [20,21,[25][26][27][28][29][30] and have been attributed to deficiencies in the description of the covalency of Cu-ligand bonds [29,31] and of the electronic excitation energies [29]. Notably, wavefunction-based methods, although systematically improvable and able to achieve arbitrary levels of accuracy, cannot at this point in time be applied to copper-based systems without the introduction of approximations that are detrimental for the accuracy of the predicted parameters [25,[32][33][34][35], as well as for hyperfine coupling constants [36], thus severely restricting the utility of these methods. Hence, density functional theory (DFT) currently provides the best compromise between accuracy and cost and serves as the most promising and pragmatic platform for establishing practically useful computational protocols with predictive value.
Over the last thirty years, a plethora of density functionals (DFs) has been developed [37,38], with different levels of approximation and empirical parameters optimized for the prediction of different chemical properties. Notably, EPR parameters practically never feature in training sets of new DFs. Experience shows that local density approximation (LDA)-and generalized-gradient approximation (GGA)-type functionals systematically underestimate g-values [20,21]. Hybrid density functionals (HDFs) include a Hartree-Fock (HF) exchange component (also called exact exchange, EEX). Inclusion of this exchange term in common hybrid functionals is beneficial [20,21,28,30], but standard HDFs still cannot provide sufficient accuracy. Increasing the admixture of exact exchange (EEX) to 30-40% has been recognized as a good strategy to improve the accuracy of g-tensor calculations for 3d, 4d and 5d transition metal complexes [39][40][41][42], while the effect of changing the correlation functional is minimal in comparison. However, increasing the Hartree-Fock exchange percentage introduces spin contamination in unrestricted calculations [43], which may be detrimental for SOMOs that have metal-ligand antibonding character. Local hybrid functionals with position-dependent EEX admixture have been reported to alleviate the problem of spin polarization [44].
Double-hybrid density functionals (DHDFs) additionally introduce unoccupied Kohn-Sham orbitals via perturbation theory. Several benchmark studies have confirmed the superior performance on DHDFs on a wide range of properties, such as relative energies, thermochemistry and electronic excitation energies [45][46][47]. Based on the general idea that "pure" DFT functionals allow a better description of static correlation at the expense of dynamic correlation effects [48][49][50], and MP2 contributes a better description of dynamic correlation [45], the judicious combination of the two is expected to achieve a good balance of dynamic and non-dynamic correlation effects. DHDFs have exhibited excellent performance on response properties calculations, such as g-tensors of organic radicals [51], electric field response properties of water nanoclusters [52], hyperfine coupling constants [53,54] and NMR nuclear shielding constants [55]. However, their performance for g-tensors of transition metal systems in general and for the particularly demanding case of Cu(II) in particular remains unexplored. To determine the scope of applicability of these approximations for the title problem, in this study we investigate the performance of sixteen distinct DHDFs for the prediction of g-tensors of a reference set of eighteen Cu(II) complexes with experimentally known structures and EPR parameters. The performance of DHDFs is contrasted with the best performing hybrid functionals for g-tensor prediction for Cu(II) complexes as reported by a recent benchmark study by Sciortino et al. [26]. Based on our results, we establish that most of the recent DHDFs outperform standard global HDFs and we are able to suggest specific DHDFT-based protocols that maximize the accuracy for g-tensor calculations of copper-based enzymes.
It is noted that the g-tensor is expressed as: where each component Δ , Δ , Δ describes the shift of the system's -value from the -value of the free electron ( 2.002319) and I is the identity matrix. In this work, we express the principal components of the -tensors as -shifts in parts per thousand (ppt), computed as Δ • 1000. The largest component is referred to as Δ or as the parallel -tensor component, Δ ∥ . We define the perpendicular -tensor component as the average of the two smaller components: Δ (Δ Δ )/2. The -tensor components of complexes 1-18 have been experimentally determined from EPR measurements and are given in Table 1. Since Cu(II) complexes have a d 9 configuration, square planar and square pyramidal structures usually exhibit strongly axial EPR signals, i.e., Δ < Δ << Δ , consistent with the high d character of the singly occupied molecular orbital (SOMO). The g-shifts arise predominantly due to the interaction of the ligating atoms with the unpaired electron, so it depends on the nature of the ligands, the nature of the Cu-ligand bond, and the corresponding bond lengths. In this benchmark set we include complexes with variable coordination spheres and a wide range of Δ values, from 83 ppt to 283 ppt. The presence of at least one S atom on the copper coordination sphere decreases the magnitude of Δ . It should be noted that It is noted that the g-tensor is expressed as: where each component ∆g xx , ∆g yy , ∆g zz describes the shift of the system's g-value from the g-value of the free electron (g e = 2.002319) and I is the identity matrix. In this work, we express the principal components of the g-tensors as g-shifts in parts per thousand (ppt), computed as ∆g = (g − g e )·1000. The largest component is referred to as ∆g zz or as the parallel g-tensor component, ∆g . We define the perpendicular g-tensor component as the average of the two smaller components: ∆g ⊥ = (∆g xx + ∆g yy )/2. The g-tensor components of complexes 1-18 have been experimentally determined from EPR measurements and are given in Table 1. Since Cu(II) complexes have a d 9 configuration, square planar and square pyramidal structures usually exhibit strongly axial EPR signals, i.e., ∆g xx < ∆g yy << ∆g zz , consistent with the high d x 2 −y 2 character of the singly occupied molecular orbital (SOMO). The g-shifts arise predominantly due to the interaction of the ligating atoms with the unpaired electron, so it depends on the nature of the ligands, the nature of the Cu-ligand bond, and the corresponding bond lengths. In this benchmark set we include complexes with variable coordination spheres and a wide range of ∆g zz values, from 83 ppt to 283 ppt. The presence of at least one S atom on the copper coordination sphere decreases the magnitude of ∆g zz . It should be noted that mixed ligand complexes exhibit magnetic parameters that are intermediate between the respective complexes with four of the same ligands of each type [8]. This pattern can be observed in the experimental data shown in Table 1. Table 1. Experimental ∆g values (ppt) of the Cu(II) complexes 1-18 studied in this work. Molecular structures of the complexes are shown in Figure 1. The complexes are arranged in groups according to the copper-coordinating atoms for each subgroup, which are also shown for convenience.

Complex
Coord. Core ∆g xx ∆g yy ∆g zz Ref.

Overview of Double Hybrid Density Functionals
In double-hybrid density functional theory, a set of Kohn-Sham orbitals is first obtained by a classical HDFT calculation. Subsequently, an MP2 calculation is performed based on this set of orbitals. A general expression for the total exchange-correlation term for standard DHDFs is: where c X is the amount of DFT exchange, and c DFT C and c MP2 C are the amounts of DFT and MP2 correlation, respectively. Usually, but not always, c MP2 The parameters c for all functionals studied in this work are shown in Table 2. The first DHDF was B2PLYP [47], proposed by Grimme and parametrized in order to reproduce heats of formation. The success of B2PLYP inspired further developments. mPW2PLYP [69] has a modified exchange part and higher accuracy towards the calculation of thermodynamic properties and barrier heights. B2GP-PLYP [70] has increased amounts of exact exchange and perturbative correlation, suggested for general purpose (GP) applications. B2K-PLYP and B2T-PLYP [71] were parametrized for reaction kinetics and thermochemistry. Some researchers attempted to minimize or eliminate the use of empirically defined parameters, striving instead to include parameters that have theoretical justification. Such functionals include PBE-QIDH [72], PBE0-DH [73] and the range-separated RSX-QIDH [74] and RSX-0DH [75], optimized by enforcing reproduction of the total energy of the hydrogen atom.  (2)-(4).  (7) 0.52 (9) As a variation of the standard MP2 approach, some functionals may treat the samespin and opposite-spin correlation contributions separately and with unequal weights. This is referred to as a spin-component scaling (SCS). In case the same-spin contribution is entirely omitted, the approach is referred to as spin-opposite scaled (SOS). The general expression for the total exchange-correlation term for the spin-component scaled doublehybrid functionals with dispersion correction (DSD-DFT) is: where c O is the amount of opposite-spin MP2, c S of same-spin MP2, and s 6 of the dispersion correction. DSD-BLYP [76] and DSD-PBEP86 [77] have been parametrized to reproduce thermochemistry, kinetics, and dispersion forces, and were found to be more accurate than the previously most successful B2GP-PLYP functional.
In range-separated DFs [78][79][80], the interelectronic repulsion, r 12 −1 , is separated using an error function, erf, into a short-range and a long-range component: where the first term accounts for the short-range interaction and the second term accounts for the long-range interaction. ωB2PLYP and ωB2GP-PLYP [81] functionals were parametrized to reproduce electronic excitation energies. In addition, ωB88PP86 and ωPBEPP86 [82] were very recently presented, parametrized to reproduce experimental singlet-singlet as well as singlet-triplet excitations. The first range-separated DHDF that also includes spin component scaling, ωB97X-2 [83], was parametrized for thermochemistry, kinetics, and noncovalent interactions. The parameters µ from Equation (4) for the range separated functionals used herein are shown in Table 2.

Evaluation Criteria
Evaluation of the methods used for the prediction of the g-shifts was based on the differences (D) of the calculated from the experimental parallel and perpendicular components: The mean values of the above parameters for each function for the set of eighteen Cu(II) complexes under study are used to describe the performance of the functional. The mean absolute difference (MAD) of the ∆g component is defined as: and straightforwardly describes the ability of the method to reproduce the experimental values, while the mean difference (MD) reflects the possible systematic over-or underestimation of g-shifts by the specific method. The standard deviation of the above parameters was also considered. The standard deviation of the differences (SDD) is estimated as: and the standard deviation of the absolute differences (SDAD): Finally, the mean absolute percentage differences (MAPD) are examined:

Performance of Functionals
The mean values of the evaluation criteria for each functional are shown in Table 3. To begin with, we observe that the MAPDs for ∆g ⊥ are very close to those for ∆g , therefore we focus our analysis on the parallel components for simplicity and the conclusions largely reflect also the performance for the perpendicular components of the g-tensor. Besides, since the parallel component is the largest, it also has the largest D values. We start the evaluation of the relative performances by examining the MADs and MDs. The MAD ∆g and MD ∆g values for the 19 functionals are plotted in Figure 2, top. Complete results are provided in the Supporting Information, Tables S1-S19.
Focusing first on conventional hybrid functionals, we see that B3LYP is the worst method considered here. It shows the largest MAD ∆g value of 67 ppt, while PBE0 and BHandHLYP follow with MAD ∆g values of 54 and 48 ppt, respectively. Thus, among conventional HDFs the best performance is observed for PBE0 and BHandHLYP, which is in agreement with Sciortino et al. [26]. Only two double-hybrid functionals perform similarly to HDFs. These are B2PLYP and RSX-QIDH with MAD ∆g values of 57 and 55 ppt, respectively. They are also described by the SDD value of each method, given in Table 3. ωB2GP-PLYP, ωPBEPP86 and ωB88PP86 have larger SD values than the B2GP-PLYP and PBE0-DH functionals, which means that there is a larger spread of D values with these functionals. Table 3. Mean difference (MD) and standard deviation of the differences (SDD), mean absolute difference (MAD), standard deviation of the absolute differences (SDAD) and mean absolute percent difference (MAPD), from the experimental values of the parallel and perpendicular g-shifts components calculated with each functional. Values are expressed in ppt. Inspection of the individual D ∆g i values for the complexes 1-18 obtained with each functional allows a more in-depth analysis of our results (Figure 2 bottom). This diagram shows clearly that complexes 11-18, which have at least one S coordinated ligand on copper, are described differently than complexes 1-10, which have only N and O ligating atoms. Specifically, most functionals underestimate the value of parallel g-tensor component of structures 1-10, while they overestimate it for 11-18. This behavior is even more pronounced for the HDFs. Figure 2 bottom shows clearly that even though DHDFs also treat differently the Cu-S bond, they not only have smaller MADs, but also the spread of the ∆g values is smaller. This implies that DHDFs achieve a more balanced description of Cu-ligand bond covalency and this behavior is transferable among different systems, at least to a larger extent than HDFs. In Figure 3, the D(∆g ) obtained by each functional is plotted against the calculated copper Löwdin spin population of the MP2 relaxed density, ρ s (Cu), for six complexes chosen as representative of different copper coordination spheres, i.e., N4, N2O2, O4, N2S2, N3OS and S4 (see Figure S1 for additional examples). Several conclusions can be extracted from these diagrams. First, it can be clearly observed that the value of D(∆g ) increases following a near linear trend along with ρ s (Cu). Second, each system has a different ρ s (Cu) value for which the experimental ∆g is reproduced. Even though no single functional predicts this most favorable ρ s (Cu) value for all complexes, structures with similar coordination spheres are optimally described by the same density functionals ( Figure S1). Among the DHDFs included in this study, the B2PLYP and RSX-QIDH predict the smallest and largest ρ s (Cu) values, respectively, for all complexes, which directly correlates with the observed systematic negative and positive differences from the experimental values, reflected in the large MD(∆g ) values shown in Figure 3. In addition, B3LYP and PBE0 in almost all cases predict a small spin population on Cu. Notably, accurate g-tensor prediction does not necessarily imply prediction of the "correct" spin density, since other factors, such as higher order relativistic contributions, could also alter the D(∆g ) values. We note that the possible lack of correspondence between the accuracy of DFs in the prediction of various observable properties, such as relative energies or ionization potentials, and the quality of the computed densities has been extensively debated recently [84][85][86]. In Figure 3, the D(Δ ∥ ) obtained by each functional is plotted against the calculate copper Löwdin spin population of the MP2 relaxed density, ρs(Cu), for six complexes cho sen as representative of different copper coordination spheres, i.e., N4, N2O2, O4, N2S2 N3OS and S4 (see Figure S1 for additional examples). Several conclusions can be extracte from these diagrams. First, it can be clearly observed that the value of D(Δ ∥ ) increase following a near linear trend along with ρs(Cu). Second, each system has a different ρs(Cu value for which the experimental Δ ∥ is reproduced. Even though no single functiona predicts this most favorable ρs(Cu) value for all complexes, structures with similar coo dination spheres are optimally described by the same density functionals ( Figure S1 Among the DHDFs included in this study, the B2PLYP and RSX-QIDH predict the smal est and largest ρs(Cu) values, respectively, for all complexes, which directly correlate with the observed systematic negative and positive differences from the experimental va ues, reflected in the large MD(Δ ∥ ) values shown in Figure 3. In addition, B3LYP and PBE in almost all cases predict a small spin population on Cu. Notably, accurate -tensor pre diction does not necessarily imply prediction of the "correct" spin density, since othe factors, such as higher order relativistic contributions, could also alter the D(Δ ∥ ) value We note that the possible lack of correspondence between the accuracy of DFs in the pre diction of various observable properties, such as relative energies or ionization potential and the quality of the computed densities has been extensively debated recently [84][85][86]  From the above analysis we can see that overall, the best performing functionals for this set of Cu(II) complexes are the B2GP-PLYP and PBE0-DH, because they have both the lowest MAD ∆g values, and at the same time among the lowest MD ∆g and SD values. PBE0-DH represents a definite improvement over the parent PBE0 functional, which is the best performing HDF. Interestingly, these two DHDFs represent different construction philosophies (PBE0-DH is "non empirical") yet converge to similarly accurate behavior. Notably, even though the parameters of B2GP-PLYP (and B2PLYP) were empirically optimized, they were also suggested to be justified from a purely theoretical standpoint [87]. It is also important to draw attention to the observation that the functionals ωB2PLYP, ωB88PP86, ωB2GP-PLYP, ωPBEPP86 are among the best performing DHDFs. What these functionals have in common is that they were parametrized to reproduce experimental excitation energy reference values [81,82]. From the above analysis we can see that overall, the best performing functionals for this set of Cu(II) complexes are the B2GP-PLYP and PBE0-DH, because they have both the lowest MAD Δ ∥ values, and at the same time among the lowest MD Δ ∥ and SD values. PBE0-DH represents a definite improvement over the parent PBE0 functional, which is the best performing HDF. Interestingly, these two DHDFs represent different construction philosophies (PBE0-DH is "non empirical") yet converge to similarly accurate behavior. Notably, even though the parameters of B2GP-PLYP (and B2PLYP) were empirically optimized, they were also suggested to be justified from a purely theoretical standpoint [87]. It is also important to draw attention to the observation that the functionals ωB2PLYP, ωB88PP86, ωB2GP-PLYP, ωPBEPP86 are among the best performing DHDFs. What these functionals have in common is that they were parametrized to reproduce experimental excitation energy reference values [81,82]. Another interesting point concerns the calculated SDD values in Table 3. These show that the B2PLYP functional has the smallest SDD value of 28 ppt. This is visually represented in Figure 2b, where it can be observed that even though B2PLYP performs poorly, systematically underestimating the g-shifts, its behavior is distinct in that the D values of complexes 11-18 with S coordinating ligands are not very different from complexes 1-10 which have N and O ligands only. By contrast, this is not observed for RSX-QIDH, the other DHDF of this set that performs poorly, systematically overestimating g-shifts. Therefore, we can conclude that B2PLYP achieves a more balanced description of the Cu-ligand bond covalency between different ligand donors.

Functional
Having compared the relative accuracies of the examined DHDFs on predicting the experimental g-tensor components, one may consider whether correlations exist between their construction parameters ( Table 2) and their performance. However, at this point it is not possible to discern any obvious correlation. Perhaps a way to approach this problem is to investigate how reference spin densities, computed with high-level ab initio methods, are reproduced by density functionals with respect to the parameters included in their definition [86]. Although clear choices emerge from the present study for practical applications, further systematic studies are required for the development of improved functionals optimized for spin-density dependent properties.

Conclusions
We investigated the performance of sixteen DHDF approximations for the prediction of experimental g-shifts of a set of eighteen Cu(II) complexes. The DHDFs results were compared with the best performing HDFs for Cu(II) complexes [26]. Our results show that most DHDFs perform significantly better than the best performing standard HDFs, establishing a new standard for the prediction of g-tensors of copper systems. With the exception of B2PLYP and RSX-QIDH, all other DHDFs perform significantly better than HDFs. Based on the criteria defined in Section 3.1, we found that the ranking of the DHDFs with respect to g-tensor prediction is: The best-performing functionals represent different approaches toward construction of DHDFs. B2GP-PLYP is a general purpose empirically fitted functional based on B2PLYP, while PBE0-DH is based on PBE0 and obtained by purely theoretical arguments. Nevertheless, they converge to similarly good performance. Moreover, the equally good performance of the range-separated ωB2PLYP, ωB88PP86, ωB2GP-PLYP and ωPBEPP86 directs our attention to functionals optimized for excitation energies. We suggest that such functionals can be used as a starting point for further refinement, taking into account a more extensive set of transition metal reference compounds for g-tensors, and perhaps also for other spin-dependent properties.
Overall, the double-hybrid methodologies significantly outperform conventional hybrid approaches, offering a good balance between accuracy and computational cost. We note that due to the steep scaling of MP2 (O(N 5 ) where N is a measure of the system size), the computational cost of DHDFs on large (>100 atoms) systems increases more steeply than HDFs in terms of both time and required memory. Hence, for large systems we propose the use of multilevel approaches, where the metal and first coordination sphere ligands are treated with a DHDFT method and the surrounding protein matrix can be treated with a cheaper DFT method. Based on our results, we recommend the use of B2GP-PLYP and PBE0-DH functionals for g-tensor calculations on Cu(II) complexes bearing N, O and S ligands, which are usually encountered in bioinorganic systems. Evidently, the uncertainty of even the best functionals may exceed the uncertainty of experimental values. This suggests that the successful use of these approaches will depend on the type of chemical problem under investigation and in the combination with complementary data. Nevertheless, the present study clearly defines the current state-of-the-art in quantum chemical calculations of g-tensors for Cu(II) systems and encourages further developments in refining double-hybrid DFT for spin density dependent properties.

Computational Details
All calculations were performed with Orca 5 [88]. Geometry optimizations of the eighteen Cu(II) complexes of the benchmark set was performed with the B3LYP [89][90][91] functional, starting from the crystallographic structures obtained from the Cambridge Structural database [92]. Scalar relativistic effects were treated using the zeroth-order regular approximation (ZORA) [93,94]. The resolution of identity (RI) approximation [95] and the chain-of-spheres approximation to exact exchange (COSX) [96] were used throughout to reduce computational time. Tight convergence criteria and high-quality grids (DefGrid3 in Orca convention) were used throughout. The ZORA-recontracted [97] version of the def2-TZVPP basis set [98] was used for Cu and the ZORA-def2-TZVP for all other atoms, along with fully decontracted def2/J auxiliary basis sets [99]. Since MP2 convergence with basis set size can be slower than DFT, the basis set dependence of DHDFs may increase with increasing fraction of the PT2 component [46,47]. To investigate the dependence of calculated g-shifts on the basis set size, we carried out systematic studies on complexes 3, 10 and 17, which have N4, O4 and S4 coordination spheres, respectively, using the B2PLYP and DSD-PBEP86 functionals. In line with previous studies on Cu(II) g-tensor calculations [25], the dependence on the basis set is very weak past the polarized triple-zeta level and the use of the largest available ZORA-def2-QZVPP basis set leads to negligible differences (of the order of 1-2 ppt) in the results compared to the triple-zeta basis sets. The maximum difference in computed g-shifts upon use of the quadruple-zeta basis set (5 ppt) was observed for 17 with the DSD-PBEP86 functional. Since this difference is of the same order of magnitude as the experimental uncertainty, the basis sets used here are judged to be essentially converged for all practical purposes. For the calculation of g-tensors, the DFT functionals tested in this work are: the hybrid functionals PBE0 [100][101][102], BHandHLYP [91], B3LYP [89][90][91], the double-hybrid functionals: B2PLYP [47], mPW2PLYP [69], B2GP-PLYP [70], B2K-PLYP [71], B2T-PLYP [71], PBE-QIDH [72], PBE0-DH [73], DSD-BLYP [76], DSD-PBEP86 [77], and the range-separated double-hybrid functionals: ωB2PLYP [81], ωB2GP-PLYP [81], RSX-QIDH [74], RSX-0DH [75], ωB88PP86 [82], ωPBEPP86 [82], ωB97X-2 [83]. Spin contamination values for HDFs were smaller than 0.01 for all complexes, which suggests that the ground state of all complexes can be adequately described by a single determinant. For the calculations with double-hybrid functionals the NoFrozenCore option was used. In terms of computational costs, we mention as an example that for complexes 1 (17 atoms), 18 (43 atoms) and 12 (63 atoms) the PBE0 calculation needs 1, 16 and 48 min, respectively, using 10 processing cores, while the PBE0-DH calculation needs 3, 184 and 927 min, respectively. The part of the calculation that is responsible for most of the additional computational effort is the calculation of the relaxed MP2 response density, which incorporates orbital relaxation and is consistent with first order properties as analytic derivatives. The spin-orbit coupling was treated using the spin-orbit mean-field (SOMF) operator [102] with the 1X-approximation [103,104] (SOCType 3 in ORCA convention). For the construction of the potential one-electron terms were included, the Coulomb term was computed using the RI approximation, exchange terms were incorporated via one-center exact integrals including the spin-other orbit interaction, and without local DFT correlation terms (SOCFlags 1,3,3,0 in ORCA). Picture change effects were also taken into account.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/magnetochemistry8040036/s1, Figure S1: Correlation of the differences, D(∆g ), of the calculated g-tensor from the experimental value of the parallel g-shift with the Cu spin populations computed with the respective functional for pairs of complexes with the same copper-coordinating atoms 1 and 3, 14 and 16, and with similar coordination spheres 12 and 13. Tables S1-S19: Detailed results for each functional.

Conflicts of Interest:
The authors declare no conflict of interest.