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]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 ρ(r) given by the scalar product of the magnetically induced current density J(r) with the vector potential A(r) of the external magnetic field B(30,31,42−45)an approach based on quadrature is also possible. As shown in Section , 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 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 and 3. Then, in Section , we will list the studied set of density functionals and present the results in Section : the functional benchmark is discussed in Section 5.1 and magnetizability densities and spatial contributions to magnetizabilities are analyzed in Section . The conclusions of the study are summarized in Section . Atomic units are used throughout the text unless stated otherwise, and summation over repeated indices is assumed.

Theory

The current density J(r) in eq is formally defined as the real part of the mechanical momentum densitywhere p = −i∇ is the momentum operator. Substituting eq into eq straightforwardly leads toThe 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(r) of an external static homogeneous magnetic field is expressed aswhere RO is the chosen gauge origin. The αβ component of the magnetizability tensor can then be obtained from eqs –6 aswhere the magnetizability density is defined aswhere ϵαδγ 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 tensorwe introduce the isotropic magnetizability density ρξ̅(r) defined aswhich yields information about the spatial origin of the isotropic magnetizability, as we demonstrate in Section .

Although there is freedom with regard to the choice of the gauge origin of A(r), the magnetic flux density B is uniquely defined via eq , because B = ∇ × (A(r) + ∇f(r)) holds for any differentiable scalar function f(r). The exact solution of the Schrödinger 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 atomic-orbital 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 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 open-source 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 , 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 ξ element of the magnetizability tensor is obtained from eq aswhere the xx component of the magnetizability density tensor at grid point i with quadrature weight w iswhere and 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. The ξ and ξ elements are obtained analogously.

Computational Methods

Calculations are performed for the set of 28 molecules studied in ref (18) that also provides our molecular structures and the CCSD(T) reference values: AlF, C2H4, C3H4, CH2O, CH3F, CH4, CO, FCCH, FCN, H2C2O, H2O, H2S, H4C2O, HCN, HCP, HF, HFCO, HOF, LiF, LiH, N2, N2O, NH3, O3, OCS, OF2, PN, and SO2. However, as in ref (18), O3 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 O3 is −46.2 × 10–30 J/T2, 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).

Electronic structure calculations were performed with Hartree–Fock (HF) and the functionals listed in Tables and 2 using TURBOMOLE 7.5.[110] Several rungs of Jacob’s ladder were considered when choosing the functionals listed in Tables 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 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.

Table 1

Functionals at the Local Density Approximation (LDA) and the Generalized Gradient Approximation (GGA) Considered in This Workf

functional	hybrid	type	notes	LIBXC IDa	references
LDA		LDA		1 + 7	(61−63)
BLYP		GGA		106 + 131	(1, 3), and (64)
BP86		GGA		106 + 132	(1) and (2)
CHACHIYO		GGA		298 + 309	(65) and (66)
KT1		GGA		167	(67)
KT2		GGA		146	(67)
KT3		GGA	PySCF data used	587	(68)
N12		GGA		82 + 80	(69)
PBE		GGA		101 + 130	(4) and (5)
B3LYP	GH	GGA	20% HF	402	(6)
revB3LYPb	GH	GGA	20% HF	454	(70)
B97-2	GH	GGA	21% HF	410	(71)
B97-3	GH	GGA	26.9% HF	414	(72)
BHLYPc	GH	GGA	50% HF	435	(61, 62), and (73)
BHandHLYPd	GH	GGA	50% HF	436	(1) and (73)
PBE0	GH	GGA	25% HF	406	(7) and (8)
QTP-17	GH	GGA	62% HF	416	(74)
N12-SX	RS	GGA	25% SR, 0% LR	81 + 79	(75)
CAM-B3LYP	RS	GGA	19% SR, 65% LR	433	(76)
CAMh-B3LYPe	RS	GGA	19% SR, 50% LR		(77)
CAM-QTP-00	RS	GGA	54% SR, 91% LR	490	(78)
CAM-QTP-01	RS	GGA	23% SR, 100% LR	482	(79)
CAM-QTP-02	RS	GGA	28% SR, 100% LR	491	(80)
ωB97	RS	GGA	0% SR, 100% LR	463	(81)
ωB97X	RS	GGA	15.8% SR, 100% LR	464	(81)
ωB97X-D	RS	GGA	22.2% SR, 100% LR	471	(82)
ωB97X-V	RS	GGA	16.7% SR, 100% LR	531	(83)

Two numbers indicate the exchange and correlation functionals, respectively. A single number indicates an exchange–correlation functional.

Revised version.

Following King et al. in refs (84−86), BHLYP is defined as 50% LDA exchange, 50% HF exchange, and 100% LYP correlation. It is sometimes also known as BHandH, which is its keyword in Gaussian.

BHandHLYP is 50% Becke’88 exchange, 50% HF exchange, and 100% LYP correlation.

CAMh-B3LYP is defined using the XCFun library with α = 0.19, β = 0.31, and μ = 0.33.

GH stands for global hybrid and RS for range-separated hybrid. The amount of Hartree–Fock (HF) exchange or exact exchange in the short range (SR) and the long range (LR) is also given.

Table 2

Meta-GGA Functionals (mGGA) Considered in This Workd

functional	hybrid	type	notes	LIBXC IDa	references
B97M-V		mGGA		254	(87)
M06-L		mGGA		449 + 235	(88)
revM06-Lb		mGGA		293 + 294	(89)
M11-L		mGGA		226 + 75	(90)
MN12-L		mGGA		227 + 74	(91)
MN15-L		mGGA		268 + 269	(92)
TASK		mGGA		707 + 13	(93) and (94)
MVS		mGGA		257 + 83	(95) and (96)
SCAN		mGGA		263 + 267	(97)
rSCANc		mGGA		493 + 494	(98)
TPSS		mGGA		457	(99) and (100)
revTPSSb		mGGA		212 + 241	(96) and (101)
TPSSh	GH	mGGA	10% HF	457	(102)
revTPSShb	GH	mGGA	10% HF	458	(96, 101), and (102)
M06	GH	mGGA	27% HF	449 + 235	(103)
revM06b	GH	mGGA	40.4% HF	305 + 306	(104)
M06-2X	GH	mGGA	54% HF	450 + 236	(103)
M08-HX	GH	mGGA	52.2% HF	295 + 78	(105)
M08-SO	GH	mGGA	56.8% HF	296 + 77	(105)
MN15	GH	mGGA	44% HF	268 + 269	(106)
M11	RS	mGGA	42.8% SR, 100% LR	297 + 76	(107)
revM11b	RS	mGGA	22.5% SR, 100% LR	304 + 172	(108)
MN12-SX	RS	mGGA	25% SR, 0% LR	248 + 73	(75)
ωB97M-V	RS	mGGA	15% SR, 100% LR	531	(109)

Two numbers indicate the exchange and correlation functionals, respectively. A single number indicates an exchange–correlation functional.

Revised version.

Regularized version.

The notation is the same as in Table .

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-the-identity 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 . 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 SO2. The magnetizability obtained with the vdW-optimized orbitals differed by only 0.4 × 10–30 J/T2 (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/T2 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 TURBOMOLE 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 TURBOMOLE 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.

Results

Functional Benchmark

The deviations of the DFT magnetizabilities from the CCSD(T) reference values of ref (18) are visualized as ideal normal distributions (NDs) in Figure . 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 . 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 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 .

Figure 1

Normal distributions (ND) representing the errors in the magnetizabilities for the 27 benchmark reproduced by the studied functionals, obtained by plotting the data presented in Table . The curves are ordered in each figure by increasing standard deviation. The NDs of RS functionals are shown in (a)–(c). The NDs of the GH functionals are shown in (d)–(g). The NDs of the mGGA functionals are shown in (h)–(j). The NDs of the LDA and GGA functionals are shown in (k) and (l).

Table 3

Mean Absolute Errors (MAEs), Mean Errors (MEs), and Standard Deviations (STDs) for the Magnetizabilities of the 27 Studied Molecules in Units of 10–30 J/T2 from the CCSD(T) Reference with the Studied Functionalsa

rank	functional	MAE	ME	STD	rank	functional	MAE	ME	STD
1	BHandHLYP	3.11	2.15	4.65	27	revTPSSh	7.14	7.05	5.94
2	CAM-QTP-00	3.22	0.88	4.67	28	TPSSh	7.20	7.07	6.02
3	ωB97X-V	3.22	2.51	4.36	29	B97-2	7.24	7.07	6.40
4	CAM-QTP-01	3.23	0.59	4.49	30	M08-HX	7.34	5.17	10.27
5	CAM-QTP-02	3.28	–0.23	4.36	31	BLYP	7.91	5.69	8.75
6	ωB97	3.54	2.44	4.75	32	N12-SX	8.04	7.89	7.48
7	ωB97M-V	3.61	0.41	4.75	33	revTPSS	8.20	7.86	6.68
8	CAM-B3LYP	3.73	2.38	4.86	34	TPSS	8.22	7.85	6.85
9	MN12-SX	3.80	0.22	5.34	35	revM11	8.23	6.83	10.03
10	CAMh-B3LYP	4.23	3.22	5.17	36	TASK	8.27	7.31	7.43
11	ωB97X	4.25	3.71	5.22	37	BP86	8.59	7.30	8.75
12	QTP-17	4.58	3.77	5.45	38	M11-L	8.92	5.20	9.26
13	BHLYP	4.73	0.10	6.47	39	revM06	8.94	8.67	10.27
14	B97M-V	5.19	4.13	5.58	40	PBE	9.13	7.07	9.42
15	revB3LYP	5.45	4.34	6.13	41	KT3	9.19	8.38	8.08
16	B3LYP	5.47	4.72	5.97	42	LDA	9.55	5.37	11.36
17	MN12-L	5.79	–2.03	8.02	43	CHACHIYO	9.76	9.17	8.88
18	KT1	5.87	1.15	7.11	44	M11	9.93	7.61	13.77
19	rSCAN	5.91	5.00	6.06	45	M06-2X	10.15	9.01	13.12
20	PBE0	5.96	5.56	6.81	46	MVS	10.35	9.92	9.20
21	ωB97X-D	6.22	5.89	6.35	47	M08-SO	10.40	8.09	14.34
22	SCAN	6.30	5.89	5.96	48	N12	10.89	10.01	9.58
23	KT2	6.42	5.58	7.21	49	MN15	11.45	10.45	12.82
24	MN15-L	6.57	–5.27	6.94	50	M06-L	12.49	12.45	9.42
25	B97-3	6.61	6.61	6.26	51	M06	13.34	13.11	13.16
26	revM06-L	7.00	6.23	5.98	52	HF	18.40	7.48	61.81

The functionals are ordered in increasing MAE.

Examination of the data in Table 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 a,b, whereas the NDs of the studied RS mGGA functionals are shown in Figure c.

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/T2 and 8.87 × 10–30 J/T2 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 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, and M06 are ranked between 39th and 51st. The NDs of the GH functionals are compared in Figure d–g.

B97M-V, at the 14th place, is the best pure mGGA functional. The rSCAN and SCAN functionals are ranked 19th and 22th, respectively, whereas revTPSS and TPSS appear at positions 33 and 34, respectively. The pure mGGA functionals of the Minnesota series are ranked 17th (MN12-L), 24th (MN15-L), 26th (revM06-L), and 50th (M06-L). The performance of the Minnesota pure mGGA functionals, excluding M06-L, is about the same as that of TASK and the other mGGA functionals. The magnetizabilities calculated with the revised M06-L (revM06-L) functional are more accurate than those with M06-L. The MVS mGGA functional is ranked 46th. The NDs for the mGGA functionals are shown in Figure h–j.

The magnetizabilities calculated with several of the Minnesota functionals are inaccurate. Seven of the eight worst performing functionals (M11, M06-2X, MVS, M08-SO, N12, MN15, M06-L, and M06) in Table 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 k,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.

Magnetizability Densities

Spatial contributions to the magnetizability densities, i.e., the integrand in eq , are illustrated for H2O, NH3, and SO2 in Figure , with Figure showing the corresponding CDTs. The magnetizability densities are calculated with the gauge origin of the external magnetic field (RO) at (x, y, z) = (0, 0, 0). In the calculations on H2O and SO2, the magnetic field perturbation is perpendicular to the molecular plane, while for NH3, the perturbation is parallel to the C3 symmetry axis. In the case of H2O, the current density flux around the whole molecule (Figure ) leads to the ring-shaped contribution shown in Figure . The magnetic field along the symmetry axis of NH3 also results in a current density flux around the molecule at the hydrogen atoms (Figure ), giving rise to a similar ring-shaped contribution shown in Figure .

Figure 2

Visualization of the isotropic magnetizability density ρ̅ξ(r) (eq ) shown in the molecular plane of H2O (a) and SO2 (b) as well as in the plane formed by the hydrogen atoms of NH3 (c), positioned 0.06 a0 away from the N atom toward the hydrogen atoms. Negative contributions are shown in pink and positive ones in green. The gauge origin RO is (0, 0, 0) a0.

Figure 3

Streamline representation of the CDT (eq ) of H2O (a), SO2 (b), and NH3 (c). The CDT is calculated with the magnetic field perpendicular to the molecular plane of H2O and SO2 as well as with it along the symmetry axis of NH3. The color scale represents the strength of the CDT in nAT–1a0–2.

The isotropic magnetizability density of SO2 shown in Figure 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 ). The p-orbital shaped contributions to the magnetizability density of SO2 around the oxygen atoms in Figure originate from the atomic CDTs. The patterns of the CDT of H2O and SO2 lead to the different magnetizability densities shown in Figure a,b, respectively. The positive magnetizability densities in H2O and NH3 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 H2O, NH3, and SO2 calculated with the gauge origin at RO = (1, 1, 1) a0 are shown in the SI. The contribution of the choice of the gauge origin to the magnetizability computed from eq vanishes when the CDT fulfills the charge conservation condition[29]Calculating the magnetizability for NH3 with a gauge origin set to RO = (100, 100, 100) a0 yielded a value that differs by 0.32% from the one computed for RO = (0, 0, 0). When the gauge origin is set to RO = (1, 1, 1) a0, 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 H2O and SO2 also change by only 0.46 and 0.03% when moving the gauge origin from (0, 0, 0) a0 to (100, 100, 100) a0, 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 coupled-cluster calculations with singles and doubles and perturbative triples [CCSD(T)] reported by Lutnæs et al.[18] Our results are summarized graphically in Figure : the top functionals afford both small mean absolute errors and standard deviations, but the same is not true for all recently suggested functionals.

Figure 4

Mean absolute errors (the blue solid line) as well as the errors’ standard deviations (red crosses) of the magnetizabilities in 10–30 J/T2 of the 27 studied molecules obtained with the 51 functionals compared to the CCSD(T) reference.

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 TURBOMOLE 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.