Revisiting Gauge-Independent Kinetic Energy Densities in Meta-GGAs and Local Hybrid Calculations of Magnetizabilities

In a recent study [J. Chem. Theory Comput. 2021, 17, 1457–1468], some of us examined the accuracy of magnetizabilities calculated with density functionals representing the local density approximation (LDA), generalized gradient approximation (GGA), meta-GGA (mGGA), as well as global hybrid (GH) and range-separated (RS) hybrid functionals by assessment against accurate reference values obtained with coupled-cluster theory with singles, doubles, and perturbative triples [CCSD(T)]. Our study was later extended to local hybrid (LH) functionals by Holzer et al. [J. Chem. Theory Comput. 2021, 17, 2928–2947]; in this work, we examine a larger selection of LH functionals, also including range-separated LH (RSLH) functionals and strong-correlation LH (scLH) functionals. Holzer et al. also studied the importance of the physically correct handling of the magnetic gauge dependence of the kinetic energy density (τ) in mGGA calculations by comparing the Maximoff–Scuseria formulation of τ used in our aforementioned study to the more physical current-density extension derived by Dobson. In this work, we also revisit this comparison with a larger selection of mGGA functionals. We find that the newly tested LH, RSLH, and scLH functionals outperform all of the functionals considered in the previous studies. The various LH functionals afford the seven lowest mean absolute errors while also showing remarkably small standard deviations and mean errors. Most strikingly, the best two functionals are scLHs that also perform remarkably well in cases with significant multiconfigurational character, such as the ozone molecule, which is traditionally excluded from statistical error evaluations due to its large errors with common density functionals.


Introduction
Molecular magnetic properties such as nuclear magnetic resonance (NMR) shielding constants, NMR spin-spin coupling constants, and magnetizabilities are useful probes of molecular and electronic structure, and they can be studied to pinpoint the location of atoms in a molecule, for instance.However, interpreting the measured spectra often requires computational modeling.][5][6][7][8][9][10] However, the accuracy of DFT for some of these properties has not been thoroughly established in the literature.Some of us recently employed the benchmark set of ref. 7 to assess the accuracy of density functionals 11,12 with an emphasis on newer functionals from rungs 2-4 of the usual Jacob's ladder hierarchy; 13 this benchmark set consists of magnetizabilities for 28 small main-group molecules computed with coupled-cluster theory with singles, doubles and perturbative triples [CCSD(T)], which is a highly accurate wave function theory (WFT).
It was found in ref. 11 that the BHandHLYP 14 global hybrid (GH) functional and some rangeseparated hybrid (RSH) functionals provided the closest agreement with the CCSD(T) reference data, with the smallest mean absolute deviations (MADs) being slightly above 3 × 10 −30 J/T 2 , while Hartree-Fock gave 7.22 × 10 −30 J/T 2 and was ranked 29 th best out of 52 methods evaluated. 11,12Some functionals, in particular of the highly parameterized Minnesota functionals, were found to reach MADs with large errors above 10 × 10 −30 J/T 2 .
Our study in refs.11 and 12 employed Turbomole, 15,16 which until recently relied on the Maximoff-Scuseria (MS) approach 17 to turn the kinetic energy density τ used in metageneralized gradient approximation (mGGA) functionals into a gauge-independent quantity in the presence of a magnetic field.However, recent work by some of us has shown that the MS approach that is also used in the Gaussian 18 program may lead to artefacts for nuclear magnetic shielding constants, such as artificial paramagnetic contributions to the shielding constant of spherical atoms, which can be avoided by employing Dobson's current-density extension of τ , 19 instead; [20][21][22][23] this approach is nowadays available in Turbomole for timedependent density functional theory (TDDFT) calculations [24][25][26] and other properties. 27,28The approach is also used in other packages as well for TDDFT and NMR calculations. 29,30][33][34][35] Holzer et al. 27  Several local hybrid (LH) functionals that employ a position-dependent exact-exchange (EXX) admixture 36,37 and that have been shown to exhibit promising accuracy for nuclear magnetic shielding constants [20][21][22][23] performed excellently for the theory-based test set, with relative mean absolute deviations in magnetizabilities below 1 %.However, Holzer et al. 27 put more emphasis on the comparison to experimental data, in which deviations for most functionals including LH functionals were generally above 5 % and were less systematic than in the theory-based benchmark.
We argue that the comparison to these experimental data does not afford a reliable assessment of the accuracy of density functionals, as many of these data exhibit large error bars.As LH functionals are of interest for many kinds of properties, we revisit their performance for magnetizabilities in this work, based on the theoretical dataset of Lutnaes et al. 7 which was also used in previous studies. 11,12,27Although the limitations of this kind of static DFT benchmarks need to be recognized, 38 such studies often give helpful guidance in terms of the suitable physical contents of density functional approximations.
Going beyond the LH functionals studied in ref. 27, we also include the first modern rangeseparated local hybrid (RSLH) ωLH22t 39 in our assessment, as it has been shown to provide remarkable accuracy for quasiparticle energies of a wide variety of organic chromophores of interest in molecular electronics and organic photovoltaics, 40 while also performing well for many other ground-and excited-state properties. 39We also investigate the performance of recent strong-correlation-corrected LH functionals (scLH) such as scLH22t and scLH22ta, 41 as well as the most recent models with simplified constructions of the sc-correction terms. 42e will show that these functionals can more reliably reproduce magnetizabilities of systems with large static correlation effects such as O 3 , whose magnetizabilities predicted by regular LH functionals significantly deviate from the CCSD(T) reference value. 27n addition, we will also study the Dobson formulation for τ on a wider set of τ -dependent functionals than those studied in ref. 27, including various older and newer mGGA functionals and mGGA-based global and range-separated hybrid functionals.
The layout of this work is as follows.Next, in section 2, we will discuss the employed methodology for computing the magnetizability (section 2.1), gauge origin problems (section 2.2), and LH functionals (section 2.3), and then the employed computational methodology in section 3. We present the results of this study in section 4, and conclude in section 5. Atomic units are used throughout, unless specified otherwise.

Methods for Calculating Magnetizabilities
Magnetizabilities are commonly calculated as the second derivative of the electronic energy with respect to the external magnetic field 9,43-46 The magnetic interaction can also be expressed as an integral over the magnetic interaction energy density ρ B (r), which is the scalar product of the magnetically induced current density J B (r) with the vector potential A B (r) of the external magnetic field B 11,12,47-53 (2) The second derivatives of the magnetic interaction energy with respect to the components of the external magnetic field, together forming the elements of the magnetizability tensor ξ, can be obtained from eq. ( 2) as an integral over the scalar product of the first derivatives of the vector potential of the external magnetic field with the magnetically induced current-density susceptibility (CDT), ∂J B γ (r)/∂B β B=0 , in the limit of a vanishing magnetic field [52][53][54] where the vector potential A B (r) of a homogeneous external magnetic field is and R O is an arbitrary gauge origin.][57][58] The isotropic magnetizability (ξ) is obtained as one third of the trace of the magnetizability tensor where the magnetizability density tensor is defined in terms of the CDT as where ϵ αδγ is the Levi-Civita symbol, and α, β, γ, δ and r δ ∈ {x, y, z} are Cartesian directions.The ξ αα , α ∈ {x, y, z} elements of the magnetizability tensor can be obtained by quadrature as where ρ ξ i;αα is a diagonal element of the magnetizability density tensor at quadrature point i and w i is the corresponding quadrature weight.

Gauge-Origin Problems
The use of a finite one-particle basis set introduces issues with gauge dependence into quantum chemical calculations of magnetic properties.The CDT can be made gauge-origin independent 55,57,58,69 by using gauge-including atomic orbitals (GIAOs), 4,70-73 also known as London atomic orbitals (LAOs), 43,46,74 where i is the imaginary unit, χ µ (r) is a basis function centered at R µ , and c is the speed of light which has the value c = α −1 ≈ 137.036 in atomic units, where α is the fine-structure constant.
The mGGA approximation for the exchangecorrelation energy density contains a dependence on the kinetic energy density τ , which reads in the field-free case as where D µν is the AO density matrix.However, this dependence requires additional care, as the (uncorrected) τ (r) of eq. ( 9) is clearly not a priori gauge invariant even when using GIAOs.
A widely used model to render τ (r) gaugeinvariant was introduced by Maximoff and Scuseria 17 (MS) as where j p (r) is the paramagnetic current density defined as Advantages for coupled-perturbed KS (CPKS) calculations with mGGAs arise from the diagonality of the Hessian for τ M S (r), which allows solving the CPKS equations in a single step.However, the semi-local exchangecorrelation (XC) contribution does not produce a linear response, even though a genuine current-density functional is expected to provide such a response. 75,76τ M S (r) also does not constitute a proper iso-orbital indicator, 77 is non-universal, 24,33,78 and introduces paramagnetic artefacts in shielding calculations. 20 formulation that avoids the disadvantages of the MS model was proposed by Dobson 19,79 τ D (r) in eq. ( 12) leads to gauge independence, while also introducing a current-response of the semi-local XC contribution and thereby a correct physical behavior for τ -dependent density functionals.As the electronic Hessian corresponding to eq. ( 12) is non-diagonal, an iterative solution to the CPKS equations is now necessary, rendering the calculations slightly more expensive than when the physically incorrect eq. ( 10) is used.

Local Hybrid Functionals
In addition to studying the importance of the Dobson formulation (eq.( 12)) of the kineticenergy density in reproducing accurate magnetizabilities, we will also evaluate the accuracy of LH functionals.A form for an LH functional that emphasizes the inclusion of nonlocal correlation terms (often considered to cover nondynamical correlation, NDC) together with full exact exchange and a (semi-local) dynamical correlation (DC) functional is 13,37 where g(r) is the local mixing function (LMF) controlling the fraction of exact exchange included at r.In most of the LH functionals con-sidered here we use a so-called t-LMF defined as the scaled ratio between the von Weizsäcker 80 and KS kinetic energy densities 2][83][84] In the LH functionals with a CF from the Berlin group, 39,41,85 the semi-local CFs are currently derived within the partial integration gauge (pig) approach. 83e also evaluate two recent extensions of LH functionals: the so-called strong-correlationcorrected LH functionals (scLH) 41,42,86 and a range-separated local hybrid (RSLH) functional (ωLH22t). 39In the scLH functionals, a strongcorrelation factor q AC (r) is introduced into the LMF of the LH: This approach is based on the local adiabatic connection approach, and is adapted from the KP16/B13 87 and B13 88 sc-models.The most recent models employ simplified real-space measures to detect strong correlations, as well as modified damping functions to avoid doublecounting of NDC contributions in more weakly correlated situations. 42n the absence of strong correlations, q AC −→ 0.5 as a lower bound, and the underlying LH functional is restored.Whenever the quantities underlying q AC detect locally the presence of strong correlations, q AC is increased, maximally up to 1.0.In the exchange picture, this means that the EXX admixture is locally diminished; in some cases it may even become negative. 41,42his enhances the simulation of nonlocal correlation contributions and is crucial for reducing fractional spin errors and for improving spinrestricted bond dissociation curves. 41,42SLH functionals like ωLH22t 39 may be writ-ten as where e sl X,sr and e ex X,sr are short-range exchangeenergy densities, controlled by the rangeseparation parameter ω.In consequence, the ωLH22t functional has full long-range EXX admixture, 39 like the RSH functionals evaluated here as well, while the short-range EXX admixture is determined by the LMF.

Computational Methods
As in refs.11 and 12, the unperturbed and magnetically perturbed density matrices are generated with the nuclear magnetic shielding module of the Turbomole program (mpshift). 16 locally modified version of Turbomole 7.6 has been used for the present calculations.Turbomole employs Libxc 89 to evaluate many of the presently considered density functionals.A detailed description of the implementations of the LH functionals, 20 the Dobson model, 23 and higher derivatives of the density used in the pig2 CF 21 can be found in the respective publications.The necessary equations for the magnetic-field derivatives of scLHs and RSLHs are outlined in Section S1 of the Supporting Information.
LAOs. 43,46 The selfconsistent field convergence criteria were set as 10 −9 for the energy and 10 −7 for the density.Large numerical integration grids per the approach of Becke 97 were used with the Turbomole 15,16 setting gridsize 7, following the original work of Treutler and Ahlrichs 98 with later extensions documented in the Turbomole manual. 99The non-standard exactexchange integrals occurring in LH, scLH and RSLH functionals are calculated by efficient semi-numerical integration techniques, [100][101][102][103] using standard DFT grids.
The resolution-of-the-identity (RI) approximation was used to evaluate the Coulomb contribution (J), using Turbomole's "universal" auxiliary basis set by Weigend. 104Although one needs to be careful about mixing auxiliary basis sets for different families, we tested the accuracy of this RI-J approximation with a few of the functionals considered in this work, and the resulting magnetizabilities with the RI-J approximation coincided with values obtained without it to the reported number of digits, which is consistent with results recently computed for NMR shieldings. 105,106We note that accurate auxiliary basis sets with controllable accuracy for RI calculations can nowadays be easily generated automatically, 107,108 and recommend such autogenerated auxiliary basis sets to be used when tailored basis sets are not available.We refer the reader to ref.
109 for a recent review of further automatical auxiliary basis generation techniques.
6][57][58] The integral in eq. ( 6) is calculated in Gimic using Becke's 97 multicenter quadrature scheme, 11,12 employing the Numgrid library 110 to generate the atomic quadrature grids with a hardness parameter of 3 for the atomic weight partitioning, radial grids from Lindh et al. 111 and Lebedev's angular grids. 112Both Gimic and Numgrid are free and open-source software. 113,114n this work, we study the accuracy of magnetizabilities reproduced by the 31 functionals listed in table 1 for a dataset of 28 molecules: HCN, HCP, HF, HFCO, HOF, LiF, LiH, N 2 , N 2 O, NH 3 , O 3 , OCS, OF 2 , PN, and SO 2 , which have also been used as benchmark molecules in other studies. 7,11,12The obtained DFT magnetizabilities are compared to the CCSD(T) reference values of ref. 7. Since the magnetizability for O 3 introduces significant uncertainties due to large static correlation effects, it was omitted from the statistical analysis, as was also done in refs.11 and 12.

Results
We will begin the discussion of the results in section 4.1 by examining the accuracy of all the considered functionals with the Dobson and MS models for τ in the subset of data without O 3 , which exhibits strong correlation effects.We then discuss the accuracy of the various functionals O 3 in section 4.2.The magnetizabilities for all studied molecules and functionals with the MS and Dobson formulations of τ can be found in the Supporting Information (SI), accompanied with violin plots of the corresponding error distributions.

Accuracy of Density Functionals with Various Models for τ
Table 2 summarizes the overall statistical evaluations and the ranking of the various functionals within the Dobson and the MS models for τ .We begin by noting that our results agree with those of Holzer et al. 27 for the subset of functionals also studied in ref. 27.
In the mean absolute error (MAE) analysis, the largest effects of using τ D instead of τ MS are obtained for some Minnesota functionals, e.g., −5.47 × 10 −30 J/T 2 for MN15, −2.38 × 10 −30 J/T 2 for M06-2X, and +1.69×10 −30 J/T 2 for M06.Other functionals that exhibit large (more than 0.5 × 10 −30 J/T 2 ) effects on the MAEs include PW6B95 (−0.72 × 10 −30 J/T 2 ) and MN15-L (+0.52 × 10 −30 J/T 2 ).We note here that many Minnesota functionals, including M06 and M06-2X, have been recently found to be numerically ill-behaved, while MN15 and MN15-L appear to behave better. 137,138he MAE values suggest that ensuring proper gauge invariance may either improve or worsen the agreement with the CCSD(T) reference data.Some of us have found a similar behavior for NMR chemical shifts, 22,23 and attributed it to massive error compensation in some of the cases.
However, in most cases the differences between the MS and Dobson formulations are small, and the rankings of the best-performing functionals are not affected very much by    switching from using the MS expression to the Dobson expression.Some exceptions do occur; for instance the change of −0.40 × 10 −30 J/T 2 going from τ MS to τ D leads to an improved ranking by three positions for ωLH22t.
Interesting differences can be seen between the behavior of some first-generation LH functionals based on LSDA exchange-energy densities (e.g., LH12ct-SsirPW92, LH12ct-SsifPW92, and the related scLH21ct-SVWNm) and that of the more advanced functionals (LH, scLH, and the ωLH22t RSLH functional).The accuracy of functionals of the former type deteriorates somewhat after switching from the computationally convenient τ MS to the physically more correct τ D , while the accuracy of the latter type of functionals improves when τ D is used.
A closer analysis of the origin of the differences between the MS and Dobson formulations of τ is beyond the scope of this work.However, a recent study of the Dobson-based gaugeinvariance contributions to TDDFT excitation energies has been able to link the magnitude and even the sign of the effect to the way in which τ enters the enhancement factor of the mGGA functionals and other τ -dependent functionals. 26H20t and its sc-corrected extensions scLH22t and scLH23t-mBR occupy ranks 3, 2 and 1, respectively, with the scLHs improving slightly over their parent LH functional.The reduced MAE of scLH23t-mBR and scLH22t as compared to LH20t arises from some of the molecules that are expected to exhibit larger static correlation effects, as evidenced by larger deviations between CCSD(T) and CCSD results in Ref. 7. Apart from the true static correlation case O 3 (see below), scLH22t gives notable improvements of more than 2×10 −30 J/T 2 for PN and SO 2 .scLH23t-mBR gives large improvements for PN and HCP.The performance of both functionals would be even more impressive if not for the somewhat larger deviations of −6.4 × 10 −30 J/T 2 (scLH22t) and −5.3 × 10 −30 J/T 2 for the LiH molecule as compared to the LH20t value of +3.2 × 10 −30 J/T 2 .scLH23t-mBR-P, with its Padé-based q AC also improves particularly on PN and HCP but deteriorates on LiH (−10.6 × 10 −30 J/T 2 ), hampering the functionals's overall statistical performance compared to the other two scLHs and leaves it slightly behind LH20t in the overall ranking.Effects of the sc-corrections are much smaller for the other molecules, which is consistent with an efficient damping of the sc-factor for weakly correlated systems.It is presently unclear why the effects of the sc-corrections are below 1 × 10 −30 J/T 2 for several other systems with larger differences between CCSD and CCSD(T) (H 2 CO, OF 2 , HOF).
The two best-performing functionals in refs.11 and 12 (BHandHLYP and ωB97X-V) are ranked 8 th and 9 th when considering the Dobson formulation of τ .Strikingly, several LH and scLH functionals, as well as the ωLH22t RSLH functional occupy the seven top positions in the ranking, and several further LH functionals follow in the top 15 of the ranking (table 2).A similar trend was observed earlier when considering a smaller selection of LH functionals. 27However, somewhat different conclusions were then reached because the experimental data used as reference values in that work exhibit very large error bars and are therefore ill-suited for benchmarking purposes.
Individual changes from LH20t to scLH22ta which lacks the damping factor are overall somewhat more pronounced and lead to a slightly larger MAE for the scLH functionals, indicating some deterioration of the accuracy of the magnetizability for more weakly correlated systems, which places scLH22ta behind LH20t, but scLH22ta is still ranked 5 th best.
Further top-performing functionals include LH14t-calPBE (rank 6) and the ωLH22t RSLH (rank 7), but their MAE in the 3.0-3.3×10 −30 J/T 2 range is already comparable to those of the best-performing functionals from refs.11 and 12 (BHandHLYP, ωB97X-V, and CAM-QTP functionals).In agreement with refs.11 and 12, B97M-V remains the highestranked non-hybrid functional in the evaluations, its MAE is reduced from 5.19×10 −30 J/T 2 to 4.78 × 10 −30 J/T 2 by using the Dobson formulation of τ .The normal distributions of the magnetizabilities of the studied functionals with the Dobson formula for τ are shown in fig.1; analogous plots for the MS formula for τ are given in the SI, along with violin plots of the errors for both descriptions, both with and without the inclusion of O 3 in the data set.

Accuracy on O 3
As was already discussed above, due to its large static correlation effects, O 3 has been excluded from statistical evaluations in previous density functional assessments. 7,11,12Its large deviations would otherwise dominate the statistics and it is unclear whether the CCSD(T) reference value is sufficiently accurate for this molecule, given that the inclusion of the perturbational triples (T) contributions reduces the magnetizability by more than 45×10 −30 J/T 2 as compared to the CCSD value. 7In the absence of experimental data, we may estimate a reasonable range of values by also considering earlier complete active space self-consistent field (CASSCF) GIAO (97.8 × 10 −30 J/T 2 ) 139 and multiconfigurational (MC) individual gauge for localized orbitals (IGLO) (89.7×10 −30 J/T 2 ) 140 results.These values are somewhat smaller than the CCSD(T) reference value of 121.5 × 10 −30 J/T 2 and suggest that magnetizabilities around 100 × 10 −30 J/T 2 with an error margin of about ±20 × 10 −30 J/T 2 seem to define the most likely range.We review the various DFT results in light of this range in table 3.
Most functionals overestimate the magnetizability of O 3 dramatically.Among the standard functionals studied here and in refs.11 and 12, the specialized GGA functionals like KT1, KT2, and KT3 get closest to the reference value, with KT1 providing the smallest value of 131.9 × 10 −30 J/T 2 .These DFAs are not good performers in general according to the previous magnetizability benchmark, 11,12 and would enter at ranks 19, 21, and 31 with MAE of 5.87 × 10 −30 J/T 2 , 6.42 × 10 −30 J/T 2 and 9.19 × 10 −30 J/T 2 , respectively, in the statistical evaluation of table 2. Other standard GGA functionals like BLYP or BP86 give magnetizabilities of O 3 around 180 × 10 −30 J/T 2 . 11Similar results for simple GGA functionals were also obtained in the original evaluation of Lutnaes et al.. 7 However, these functionals perform poorly for the entire test set, which is why we will not discuss them further here.Some mGGA functionals like M06-L and B97M-V also attain closer agreement when the τ MS prescription is used, but using the more appropriate τ D results in larger magnetizability and leads to a poorer agreement.When the Dobson formulation is used, most mGGA functionals perform comparably to the simpler GGA functionals, with TPSS and τ -HCTH giving the lowest values of 173.5 × 10 −30 J/T 2 and 178.5 × 10 −30 J/T 2 , respectively.Including exact exchange in GH and RSH functionals significantly increases the calcu-lated magnetizability of ozone, and thereby deteriorates the agreement with the reference value, in some cases they are more than 300 × 10 −30 J/T 2 (BHandHLYP, M06, M06-2X).We find a similar trend for all LH functionals without sc-corrections as well as for the ωLH22t RSLH functional that yield values larger than 200 × 10 −30 J/T 2 .The strikingly good performance of various scLH functionals is thus particularly notable: scLH22ta without damping factor in the sc-corrections gives the lowest value, results for the other sc-functionals are in the range of 112 − 135 × 10 −30 J/T 2 .
Given that four of the five scLHs are also among the five best-performing functionals for the entire test set, these results are significant in suggesting that the scLH functionals to some extent indeed escape the usual zerosum game between delocalization errors and strong-correlation errors, where larger EXX admixtures improve on the former but deteriorate the latter. 141Such an "escape" has been found recently for fractional spin errors and the related spin-restricted dissociation curves of diatomics. 41,42It is gratifying to see this here for a very different property.

Conclusions
This work extends in two directions previous studies of DFT functionals for the computation of molecular magnetizabilities.We considered using a gauge-independent local kinetic energy τ ingredient in a wide variety of mGGA functionals within Dobson's current-DFT formalism (τ D ), which had not been considered in as much detail so far.We examined the effects of the Dobson formalism by comparing the obtained magnetizabilities to values calculated with the Maximoff-Scuseria (MS) formulation (τ MS ) used in previous works.We also extended the assessment to local hybrid functionals i.e., with position-dependent exactexchange admixtures, in particular to their recent strong-correlation corrected and rangeseparated variants.
Regarding gauge invariant formulations of τ , we find that going from the previously used, computationally convenient τ MS to the more physically correct τ D leads in some cases to dramatic changes in the magnetizability, while in other cases the differences between τ D and τ MS are small.The largest effects are seen for some of the highly parameterized mGGA and mGGA hybrid functionals from the Minnesota group, whose numerical behavior has also been recently investigated and found wanting for many functionals. 137,138hile τ D leads to improved agreement with the CCSD(T) reference data for some functionals, it also deteriorates the agreement for other functionals.Notably, the effects of making τ gauge-invariant by the Dobson procedure tend to be smaller for the overall betterperforming functionals, which include many local and range-separated hybrid functionals.
The overall statistical evaluation of a wide variety of different functionals provided evidence that local hybrid functionals can yield particularly accurate magnetizabilities.Indeed, the seven best-performing functionals in the present evaluation are newer local hybrid functionals, their strong-correlation corrected variants, and the recently reported range-separated local hybrid functional ωLH22t.The overall statistical improvement compared to the so far best-performing functionals is moderate but notable, e.g., MAEs of 2.25 × 10 −30 J/T 2 , 2.35 × 10 −30 J/T 2 and 2.48 × 10 −30 J/T 2 for scLH23t-mBR, scLH22t and LH20t, respectively, compared to 3.11×10 −30 J/T 2 for BHandHLYP and 3.23 × 10 −30 J/T 2 for ωB97X-V.
The most striking result is the dramatic improvement obtained with several of the strongcorrelation corrected local hybrid functionals for the static-correlation case O 3 .Importantly, this improvement is achieved while retaining the overall highest accuracy for the weakly correlated systems relevant for the statistical evaluations.This observation for a totally different property than evaluated so far for such functionals is a further indication that strongcorrelation corrected local hybrid functionals offer an escape from the usual zero-sum game between achieving low fractional charge errors and low fractional spin errors.

Figure 1 :
Figure 1: Normal distributions of the magnetizability data calculated with the Dobson formulation of τ .
recently employed the method-ology of refs.11 and 12 to compare the effect of the MS and Dobson formulations of τ on the magnetizabilities of two test sets: the aforementioned theory-based test set of ref. 7, and a test set based on experimental data.Their results appeared at first sight somewhat inconclusive, as the two test sets resulted in different rankings of various functionals.

Table 1 :
The employed local-hybrid functionals (LH), range-separated LH functionals (RSLH), strong correlation LH functionals (scLH), functionals at the meta-generalized gradient approximation (mGGA), global hybrid (GH) functionals as well as range-separated (RS) GGA and mGGA functionals, and one GGA functional AC based on error function (scLH23t-mBR) or Padé functions (scLH23t-mBR-P); simplified construction of underlying function to identify regions of strong-correlation; see ref. 42 for details b See ref. 27 for further details of the a1 and noa2 models c This functional is called B3LYP in Turbomole, and as discussed by Hertwig and Koch, 133 it is slightly different from the B3LYP functional used in ref. 11.
a Use of simplified q

Table 2 :
Mean absolute errors (MAEs), mean errors (MEs), and standard deviations (STDs) of the errors in the magnetizabilities of the 27 studied molecules in units of 10 −30 J/T 2 from the CCSD(T) reference values with the studied functionals