Benchmarking Magnetizabilities with Recent Density Functionals

We have assessed the accuracy of the magnetic properties of a set of 51 density functional approximations, including both recently published and already established functionals. The accuracy assessment considers a series of 27 small molecules and is based on comparing the predicted magnetizabilities to literature reference values calculated using coupled-cluster theory with full singles and doubles and perturbative triples [CCSD(T)] employing large basis sets. The most accurate magnetizabilities, defined as the smallest mean absolute error, are obtained with the BHandHLYP functional. Three of the six studied Berkeley functionals and the three range-separated Florida functionals also yield accurate magnetizabilities. Also, some older functionals like CAM-B3LYP, KT1, BHLYP (BHandH), B3LYP, and PBE0 perform rather well. In contrast, unsatisfactory performance is generally obtained with Minnesota functionals, which are therefore not recommended for calculations of magnetically induced current density susceptibilities and related magnetic properties such as magnetizabilities and nuclear magnetic shieldings. We also demonstrate that magnetizabilities can be calculated by numerical integration of magnetizability density; we have implemented this approach as a new feature in the gauge-including magnetically induced current (GIMIC) method. Magnetizabilities can be calculated from magnetically induced current density susceptibilities within this approach even when analytical approaches for magnetizabilities as the second derivative of the energy have not been implemented. The magnetizability density can also be visualized, providing additional information that is not otherwise easily accessible on the spatial origin of magnetizabilities.


INTRODUCTION
Computational methods based on density functional theory (DFT) are commonly used in quantum chemistry because DFT calculations are rather accurate despite their relatively modest computational costs. Older functionals such as the Becke'88− Perdew'86 1,2 (BP86), Becke'88−Lee−Yang−Parr 1,3 (BLYP), and Perdew−Burke−Ernzerhof 4,5 (PBE) functionals at the generalized gradient approximation (GGA) as well as the B3LYP 6 and PBE0 7,8 hybrid functionals are still often employed, even though newer functionals with improved accuracy for energies and electronic properties have been developed.
The accuracy and reliability of various density functional approximations (DFAs) have been assessed in a huge number of applications and benchmark studies. 9−17 It is important to note that functionals that are accurate for energetics may be less suited for calculations of other molecular properties. 16 In specific, the accuracy of magnetic properties calculated within DFAs has been benchmarked by comparing magnetizabilities and nuclear magnetic shieldings to those obtained from coupled-cluster calculations using large basis sets, 18,19 although modern DFAs have been less systematically investigated. 16,20−23 The same also holds for nuclear independent chemical shifts 24−28 and magnetically induced current density susceptibilities, 29 −36 which have been studied for a large number of molecules, but whose accuracy has never been benchmarked properly.
Magnetizabilities are usually calculated as the second derivative of the electronic energy with respect to the external magnetic perturbation 37−41 (1) Such analytic implementations for magnetizabilities exist in several quantum chemistry programs. However, as the magnetic interaction energy can also be written as an integral over the magnetic interaction energy density ρ B (r) given by 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 30,31,42−45 (2) an approach based on quadrature is also possible. As shown in Section 2, the numerical integration approach for the magnetizability provides additional information about its spatial origin that is not available with the analytic approach based on second derivatives: the tensor components of the magnetizability density defined in Section 2 are scalar functions that can be visualized, and the integration approach can be used to provide detailed information about the origin of the corresponding components of the magnetizability tensor. Similar approaches have been used in the literature for studying spatial contributions to nuclear magnetic shielding constants. 46−53 We will describe our methods for numerical integration of magnetizabilities using the current density susceptibility in Sections 2 and 3. Then, in Section 4, we will list the studied set of density functionals and present the results in Section 5: the functional benchmark is discussed in Section 5.1 and magnetizability densities and spatial contributions to magnetizabilities are analyzed in Section 5.2. The conclusions of the study are summarized in Section 6. Atomic units are used throughout the text unless stated otherwise, and summation over repeated indices is assumed.

THEORY
The current density J B (r) in eq 2 is formally defined as the real part ( ) of the mechanical momentum density where p = −i∇ is the momentum operator. Substituting eq 2 into eq 1 straightforwardly leads to The current density susceptibility tensor 29−31 (CDT) is defined as the first derivative of the magnetically induced current density with respect to the components of the external magnetic field in the limit of a vanishing magnetic field 32−35 The vector potential A B (r) of an external static homogeneous magnetic field is expressed as where R O is the chosen gauge origin. The αβ component of the magnetizability tensor can then be obtained from eqs 4−6 as where the magnetizability density is defined as where ϵ αδγ is the Levi−Civita symbol, α, β, γ, and δ are one of the Cartesian directions (x, y, z), and r δ also denotes one of (x, y, z). The components of the magnetizability density tensor ρ αβ ξ (r) are scalar functions that can be visualized to obtain information about the spatial contributions to the corresponding element of the magnetizability tensor ξ αβ .
As the isotropic magnetizability (ξ̅ ) is obtained as the average of the diagonal elements of the magnetizability tensor we introduce the isotropic magnetizability density ρ ξ̅ (r) defined as which yields information about the spatial origin of the isotropic magnetizability, as we demonstrate in Section 5.2.
Although there is freedom with regard to the choice of the gauge origin of A B (r), the magnetic flux density B is uniquely defined via eq 6, because B = ∇ × (A(r) + ∇f(r)) holds for any differentiable scalar function f(r). The exact solution of the Schrodinger equation should also be gauge invariant. However, the use of finite one-particle basis sets introduces gauge dependence in quantum chemical calculations of magnetic properties. The CDT can be made gauge origin independent by using gauge-including atomic orbitals (GIAOs), also known as (a.k.a.) London atomic orbitals (LAOs) 32,54,55 where i is the imaginary unit and χ μ (0) (r) is a standard atomicorbital basis function centered at R μ . GIAOs eliminate the gauge origin from the expression used for calculating the CDT; the expression we use is given in the Supporting Information (SI). Since the expression for the magnetizability density in eqs 7 and 8 can be computed by quadrature, magnetizabilities can be obtained from the CDT even if the corresponding analytical calculation of magnetizabilities as the second derivative of the energy has not been implemented.

IMPLEMENTATION
The present implementation is based on the gauge-including magnetically induced current (GIMIC) program 56 and the NUMGRID library, 57 which are both freely available opensource software. Gauge-independent CDTs can be calculated with GIMIC 32−35 using the density matrix, magnetically perturbed density matrices, and information about the basis set.
To evaluate eq 7, a molecular integration grid is first generated from atom-centered grids with the NUMGRID library, as described by Becke. 58 In NUMGRID, the grid weights are scaled according to the Becke partitioning scheme using a Becke hardness of 3; 58 the atom-centered grids are determined by a radial grid generated as suggested by Lindh et al., 59 and angular grids by Lebedev 60 are used.
Given the quadrature grid, the diagonal elements of the magnetizability tensor are calculated in GIMIC from the Cartesian coordinates of the n grid points multiplied with the CDT calculated in the grid points. For example, the ξ xx element of the magnetizability tensor is obtained from eq 7 as x are the products of the z and y components of the CDT calculated in grid point i with the Cartesian coordinates y and z of the grid point, respectively, and the external magnetic field perturbation is along the x-axis, B x . The ξ yy and ξ zz elements are obtained analogously. was omitted from the analysis since it is an outlier and due to the fact that the reliability of the CCSD(T) level of theory is not guaranteed for this system: the perturbative triples correction to the magnetizability of O 3 is −46.2 × 10 −30 J/T 2 , indicating that the CCSD(T) result might still have large error bars. 18 The results of this work thus only pertain to the 27 other molecules, as in ref 18.

COMPUTATIONAL METHODS
Electronic structure calculations were performed with Hartree−Fock (HF) and the functionals listed in Tables 1 and 2 using TURBOMOLE 7.5. 110 Several rungs of Jacob's ladder were considered when choosing the functionals listed in Tables 1 and 2: local density approximations (LDAs), generalized gradient approximations (GGAs), and meta-GGAs (mGGAs). Several kinds of functionals are also included: (pure) density functional approximations, global hybrid (GH) functionals with a constant amount of HF exchange, and range-separated (RS) hybrids with a given amount of HF exchange in the short range (SR) and the long range (LR). As can be seen in Tables 1 and 2, the evaluated functionals consist of 1 pure LDA, 8 pure GGAs, 8 global hybrid GGAs, 10 range-separated hybrid GGAs, 12 mGGAs, 8 global hybrid mGGAs, and 4 range-separated mGGAs, in addition to HF.
The Dunning aug-cc-pCVQZ basis set 111−115 (with aug-cc-pVQZ on the hydrogen atoms) and benchmark quality integration grids were employed in all calculations. Universal auxiliary basis sets 116 were used with the resolution-of-theidentity approximation for the Coulomb interaction in all TURBOMOLE calculations. All density functionals were evaluated in TURBOMOLE with LIBXC, 117 except the calculations with the recently published CAMh-B3LYP functional for which XCFun was used. 118 Magnetizabilities were subsequently evaluated with GIMIC by numerical integration of eq 7. The data necessary for evaluating the CDT in GIMIC were obtained from TURBOMOLE calculations of nuclear magnetic resonance (NMR) shielding constants employing GIAOs. 54,55,110,119,120 Although response calculations are not possible at the moment in the presence of the non-local correlation kernel used in ωB97X-V, B97M-V, and ωB97M-V, we have estimated the importance of the van der Waals (vdW) effects on the magnetic properties by comparing magnetizabilities obtained with orbitals optimized with and without the vdW term in the case of SO 2 . The magnetizability obtained with the vdWoptimized orbitals differed by only 0.4 × 10 −30 J/T 2 (0.14%) from that obtained from a calculation where the vdW term was omitted in the orbital optimization. Thus, the vdW term appears to have very little influence on magnetizabilities, as is already well-known in the literature for other properties. 121 The vdW term was therefore not included in the calculations using the ωB97X-V, B97M-V, and ωB97M-V functionals in this study.
The accuracy of the numerical integration in GIMIC was assessed by comparing the TURBOMOLE/GIMIC magnetizability data to analytical values from PySCF, 122 in which LIBXC 117 was also used to evaluate the density functionals. Since PySCF does not currently support magnetizability calculations with mGGA functionals or range-separated functionals, further calculations were undertaken with Gaussian 16. 123 The analytical magnetizabilities from PySCF and Gaussian were found to be in perfect agreement for the studied LDA and GGA functionals available in both codes (LDA, BP86, PBE, PBE0, BLYP, B3LYP, and BHLYP). A comparison of the data from PySCF to the GIMIC data revealed the numerically integrated magnetizabilities to be accurate, as the magnetizabilities agreed within 0.5 × 10 −30 J/T 2 for all molecules using the B3LYP, B97-2, B97-3, BLYP, BP86, KT1, KT2, LDA, PBE, and PBE0 functionals; the small discrepancy may arise from the use of the resolution-of-identity approximation 124 in TURBO-MOLE or from the numerical integration of the magnetizability density. A comparison of the raw data for BP86 and B3LYP is given in the SI. The magnetizabilities calculated with Gaussian and TURBO-MOLE using the meta-GGA functionals were found to differ. The discrepancies between the magnetizabilities obtained with the two programs are due to the use of different approaches to handle the gauge invariance of the kinetic energy density in meta-GGAs, which are described in refs 125 and 126 for Gaussian and TURBOMOLE, respectively. We found the TURBOMOLE data to be significantly closer to the CCSD(T) reference values.
Finally, since we found the implementation of the KT3 functional in LIBXC version 5.0.0 used by TURBOMOLE to be flawed, the KT3 results in this study are based on calculations with PySCF with a corrected version of LIBXC.  Figure 1. The visualization shows the idealized distribution of the error in the magnetizability for each functional, based on the computed mean errors (ME) and standard deviation of the error (STD) given in Table 3. The raw data on the magnetizabilities and the differences from the CCSD(T) reference are available in the SI. Although the error distributions in Figure 1 are instructive, we will employ mean absolute errors (MAEs) to rank the functionals studied in this work in a simple, unambiguous fashion. The MAEs are also given in Table 3.
Examination of the data in Table 3 shows that range-separated (RS) functionals generally yield accurate magnetizabilities. Judged by the mean absolute error, the best performance is obtained with the BHandHLYP GH functional. BHandHLYP is followed by 10 RS functionals, which have much sharper distributions than the rest of the studied functionals. The best performing RS functionals are three of the six Berkeley RS functionals (ωB97X-V, ωB97, and ωB97M-V) and the three RS functionals from the University of Florida's Quantum Theory Project (QTP) CAM-QTP-00, CAM-QTP-01, and CAM-QTP-02. Five of these functionals have 100% long-range (LR) HF exchange, while the CAM-QTP-00 functional has 91% LR HF exchange. The two other RS Berkeley functionals with 100% LR exchange are ranked 11th (ωB97X) and 21st (ωB97X-D) among the studied functionals. The NDs of the studied RS GGA functionals are shown in Figure 1a,b, whereas the NDs of the studied RS mGGA functionals are shown in Figure 1c.
The CAM-B3LYP (65% LR HF exchange) and CAMh-B3LYP (50% LR HF exchange) functionals are among the top 10 functionals (ranked 8th and 10th, respectively). CAM-B3LYP was designed for the accurate description of charge transfer excitations in a dipeptide model, 76 while CAMh-B3LYP functional is aimed at excitation energies of biochromophores. 77 The best Minnesota functional, MN12-SX, is ranked 9th. MN12-SX is a highly parameterized functional with 58 parameters that is known to require the use of extremely accurate integration grids. 13 Furthermore, since MN12-SX is an RS functional with HF exchange only in the short range (SR), it may have problems modeling magnetic properties of antiaromatic molecules sustaining strong ring currents in the paratropic (nonclassical) direction. 127−129 We illustrate this with calculations on the strongly antiaromatic tetraoxa isophlorin molecule in the Supporting Information: MN12-SX yields a magnetizability that is 4 times larger than the local second-order Møller−Plesset perturbation theory (LMP2) reference value, while the magnetizabilities from BHandHLYP and CAM-B3LYP are in good agreement with LMP2. The N12-SX functional ranked 32nd is also an RS functional with 0% LR exchange. The RS Minnesota functionals with 100% LR HF exchange (M11 and revM11) have large MAEs of 9.93 × 10 −30 J/T 2 and 8.87 × 10 −30 J/T 2 and are ranked 44th and 35th, respectively.
The best global hybrid (GH) functional is BHandHLYP, which is ranked first among all functionals of this study, as was already mentioned above. Among GHs, BHandHLYP is followed by QTP-17, which is ranked 12th. Old and established GH functionals like BHLYP a.k.a. BHandH, B3LYP, and PBE0 perform almost as well as QTP-17 and are ranked 13th, 16th, and 20th, respectively. The performance of revB3LYP is Two numbers indicate the exchange and correlation functionals, respectively. A single number indicates an exchange−correlation functional. b Revised version. c Regularized version. d The notation is the same as in Table 1. practically the same as that of B3LYP; the same holds for revTPSSh and TPSSh. The other established GH functionals like B97-2, B97-3, and TPSSh and newer ones like revTPSSh and M08-HX are found in the beginning of the second half of the ranking list, whereas M08-SO, M06, revM06, M06-2X, MN15,  Table 3 are Minnesota functionals. Five other Minnesota functionals are also ranked in the lower half, placing 30th (M08-HX), 32th (N12-SX), 35th (revM11), 38th (M11-L), and 39th (revM06).
The KT1 and KT2 functionals are the best GGA functionals, ranking 18th and 23rd, respectively; both KT1 and KT2 have been optimized for NMR shieldings. 67 The older commonly used GGAs i.e., BLYP, BP86, and PBE are ranked 31st, 37th, and 40th, respectively, which is only slightly better than KT3 ranked 41st and LDA ranked 42nd. The CHACHIYO and N12 functionals, which are newer GGAs, are ranked 43rd and 48th,  respectively. The NDs of the GGA functionals and the LDA are shown in Figure 1k,l.
The magnetizabilities calculated at the HF level are significantly less accurate and have a much larger MAE-STD than those obtained at the DFT levels, and we cannot recommend the use of HF for magnetic properties.
5.2. Magnetizability Densities. Spatial contributions to the magnetizability densities, i.e., the integrand in eq 7, are illustrated for H 2 O, NH 3 , and SO 2 in Figure 2, with Figure 3 showing the corresponding CDTs. The magnetizability densities are calculated with the gauge origin of the external magnetic field (R O ) at (x, y, z) = (0, 0, 0). In the calculations on H 2 O and SO 2 , the magnetic field perturbation is perpendicular to the molecular plane, while for NH 3 , the perturbation is parallel to the C 3 symmetry axis. In the case of H 2 O, the current density flux around the whole molecule (Figure 3) leads to the ring-shaped contribution shown in Figure 2. The magnetic field along the symmetry axis of NH 3 also results in a current density flux around the molecule at the hydrogen atoms (Figure 3), giving rise to a similar ring-shaped contribution shown in Figure  2.
The isotropic magnetizability density of SO 2 shown in Figure  2 has positive (green) and negative (pink) values. Calculations of the CDT show that the oxygen atoms sustain a strong diatropic atomic CDT that flows around the atom, whereas the atomic CDT of the sulfur atom is much weaker (Figure 3). The p-orbital shaped contributions to the magnetizability density of SO 2 around the oxygen atoms in Figure 2 originate from the atomic CDTs. The patterns of the CDT of H 2 O and SO 2 lead to the different magnetizability densities shown in Figure 2a,b, respectively. The positive magnetizability densities in H 2 O and NH 3 are extremely localized close to the atomic nuclei, also because of the vortices of the atomic CDT.  The magnetizability density depends on the gauge origin of the vector potential of the external magnetic field, even though the magnetizability is independent of the gauge origin. 43 The magnetizability densities for H 2 O, NH 3 , and SO 2 calculated with the gauge origin at R O = (1, 1, 1) a 0 are shown in the SI. The contribution of the choice of the gauge origin to the magnetizability computed from eq 7 vanishes when the CDT fulfills the charge conservation condition 29 Calculating the magnetizability for NH 3 with a gauge origin set to R O = (100, 100, 100) a 0 yielded a value that differs by 0.32% from the one computed for R O = (0, 0, 0). When the gauge origin is set to R O = (1, 1, 1) a 0 , the deviation is 2 orders of magnitude smaller because the change in the magnetizability depends linearly on the relative position of the gauge origin. The magnetizabilities of H 2 O and SO 2 also change by only 0.46 and 0.03% when moving the gauge origin from (0, 0, 0) a 0 to (100, 100, 100) a 0 , respectively, showing that charge conservation is practically fulfilled in our calculations. All other positions than (0, 0, 0) for the gauge origin lead to a small, spurious CDT contribution to the magnetizability density. The GIAO ansatz modifies the atomic orbitals leading to a magnetic response of an external magnetic field that is correct to the first order for the one-center problem. 30,130 Even though GIAOs do not guarantee that the integral condition for the charge conservation of the CDT is fulfilled, 131 the basis set convergence is faster and the leakage of the CDT is much smaller when GIAOs are used. 32

CONCLUSIONS
We have calculated magnetizabilities for a series of small molecules using both recently published density functionals, as well as older, established density functionals. The accuracy of the magnetizabilities predicted by the various density functional approximations has been assessed by comparison to coupledcluster calculations with singles and doubles and perturbative triples [CCSD(T)] reported by Lutnaes et al. 18 Our results are summarized graphically in Figure 4: the top functionals afford both small mean absolute errors and standard deviations, but the same is not true for all recently suggested functionals.
Numerical methods for calculating magnetizabilities based on the quadrature of the magnetizability density have been implemented. We have shown that this method allows studies of spatial contributions to the magnetizabilities by visualization of the magnetizability density. The method has been employed to calculate magnetizabilities from magnetically induced current density susceptibilities, which were obtained from TURBO-MOLE calculations of nuclear magnetic shielding constants. Thus, magnetizabilities can be calculated in this way with TURBOMOLE even though analytical methods to calculate magnetizabilities as the second derivative of the energy are not yet available in this program. Further information about spatial contributions to the magnetizability could be obtained in the present approach by studying atomic contributions and investigating the positive and negative parts of the integrands separately in analogy to our recent work on nuclear magnetic shieldings in ref 53, which may be studied in the future work.
Our calculations show that the most accurate magnetizabilities (judged by the smallest MAE) for the studied database are obtained with BHandHLYP, which is an old global hybrid with 50% HF exchange and 50% B88 exchange. The calculations also show that the modern range-separated functionals with 100% long-range HF exchange developed by Head-Gordon and co-workers and by Bartlett and co-workers yield accurate magnetizabilities for the database. Calculations with other range-separated functionals like CAM-B3LYP and CAMh-B3LYP as well as with global hybrid functionals like QTP-17, BHLYP a.k.a. BHandH, B3LYP, and PBE0 yield relatively accurate magnetizabilities for the studied molecules. Meta-GGA functionals are found to yield somewhat better magnetizabilities than GGA and LDA functionals.
However, functionals developed by Truhlar and co-workers do not appear to be well-aimed for calculations of magnetizabilities and other magnetic properties that involve magnetically induced current densities. Magnetizabilities calculated using the popular M06-2X functional are found to be unreliable, and we do not recommend the use of the M06-2X functional in calculations of nuclear magnetic shieldings, magnetizabilities, ring-current strengths, and other magnetic properties that depend on magnetically induced current density susceptibilities. Previous studies have also suggested that the M06-2X functional sometimes underestimates magnetizabilities and ring-current strengths. 128,129,132 Revised versions of Minnesota functionals have been studied in this work and found to yield somewhat more accurate magnetizabilities than the original parameterizations. However, the revised versions also still appear on the second half of the ranking list. (2) Perdew, J. P. Density-functional approximation for the correlation energy of the inhomogeneous electron gas. Phys. Rev. B 1986, 33, 8822. (