Magnetic properties with multiwavelets and DFT: the complete basis set limit achieved

Multiwavelets are emerging as an attractive alternative to traditional basis sets such as Gaussian-type orbitals and plane waves. One of their distinctive properties is the ability to reach the basis set limit (often a chimera for traditional approaches) reliably and consistently by fixing the desired precision e. We present our multiwavelet implementation of the linear response formalism, applied to static magnetic properties, at the self-consistent field level of theory (both for Hartree–Fock and density functional theories). We demonstrate that the multiwavelets consistently improve the accuracy of the results when increasing the desired precision, yielding results that have four to five digits precision, thus providing a very useful benchmark which could otherwise only be estimated by extrapolation methods. Our results show that magnetizabilities obtained with the augmented quadruple-z basis (aug-cc-pCVQZ) are practically at the basis set limit, whereas absolute nuclear magnetic resonance shielding tensors are more challenging: even by making use of a standard extrapolation method, the accuracy is not substantially improved. In contrast, our results provide a benchmark that: (1) confirms the validity of the extrapolation ansatz; (2) can be used as a reference to achieve a property-specific extrapolation scheme, thus providing a means to obtain much better extrapolated results; (3) allows us to separate functional-specific errors from basis-set ones and thus to assess the level of cancellation between basis set and functional errors often exploited in density functional theory.


Introduction
Density functional theory (DFT) is nowadays the de facto standard for quantum chemistry applications. 1,2There are several reasons for the success of DFT: it is conceptually simple, focussing directly on the observable three-dimensional electron density; its Kohn-Sham formulation 3 allows the problem to be recast as the optimization of a single-determinant wavefunction for the fictitious Kohn-Sham system of independent particles, allowing the tools quantum chemists have employed for decades in connection with Hartree-Fock (HF) theory to be straightforwardly applied in the optimization of the Kohn-Sham state.However, the apparent simplicity also comes with a significant challenge: the definition of the (unknown) exchange-correlation (XC) functional. 4To address this shortcoming, a large library of functionals has over the years been developed, 5,6 allowing an XC functional to be chosen which is best suited for the problem at hand.This last argument is nevertheless unsatisfactory.Despite attempts at creating hierarchies of XC functionals which could be expected to perform better as the complexity or information content of the functional is increased, 7,8 such a universal hierarchy remains elusive: on the one hand improving functionals to yield better results is legitimate and necessary; on the other hand, such a wide choice of XC functionals makes it often possible to achieve the desired result for a specific problem or substrate, though with little to no predictive power when experimental reference data are not available.Avoiding pitfalls arising from spurious error cancellations is not easy, also because scientific literature is often biased towards positive outcomes.† An important way of assessing and improving the quality of currently available XC functionals is to benchmark their performance for different applications. 9,101][12][13] The molecular energy is by far the most important quantity, superseding molecular structure, vibrational properties, reactivity, and dynamics.Nevertheless, the interest This journal is © the Owner Societies 2016 in molecular properties other than the energy has been steadily growing, both in connection with spectroscopic investigations [14][15][16] and with the development of materials with specific properties. 17,18FT is for instance often used in connection with response theory to obtain a wide range of molecular properties. 19This asymmetry (functionals optimized/benchmarked with respect to their energetic performance and employed for property calculations) mandates a thorough assessment of the quality of such functionals for molecular properties, and several such studies have already been presented in the literature, see for instance the review by Laurent and Jacquemin of benchmark studies of XC functionals used in time-dependent DFT studies and references therein. 20nother important aspect to be considered is the choice of basis set.For molecular calculations with DFT, a linear combination of atomic orbitals approach have almost exclusively been used.When a given basis is employed, the error on the final result is invariably a combination of a direct basis set error and a functional error, that is, difference between the value obtained with the given functional and the ideal exact one in the limit of a complete basis set.To complicate matters further, most functionals are parametrized against thermochemical data by making use of a given basis set, 9,12,13 thus introducing an ''indirect'' basis set error.As a consequence, the definition and quality of a basis set limit for molecular properties within DFT is still an open question.In particular, in two recent studies, Lutnaes et al. 21and Teale et al. 22 have extensively addressed such issues for magnetizabilities and nuclear magnetic resonance (NMR) shielding constants of a set of 28 small, closed-shell molecules.
Low-scaling DFT methods 23 combined with modern highperformance computer clusters allow substrates with several thousand atoms to be modelled: ideally, a large basis set should be employed to guarantee convergence in the property value; however, using large basis sets for large systems is challenging due to numerical problems caused by near linear dependencies that will arise in the basis set. 24nother potential problem is given by the integration grids employed in DFT: they are optimized to achieve the best compromise between accuracy and efficiency for energy calculations, but there is no guarantee that such grids will work equally well for the functional derivatives required for the evaluation of molecular properties.
3][34][35][36][37][38][39][40][41] The first three groups have been successfully applied in materials modeling, especially in combination with periodic boundary conditions and pseudopotentials.However, for molecular all-electron calculations, the most promising approach so far seems to be the one using multiwavelets, pioneered originally by Harrison and coworkers. 32everal properties make multiwavelets attractive compared to atom-centered bases: they are by construction orthonormal, avoiding linear dependencies; completeness in the basis is achieved by refinement, with a rigorous, predefined error control; function representations can be refined adaptively, limiting the memory footprint; a separated tensorial representation of integral convolution operators is employed, 42 coupled with the non-standard form of operators, 43 achieving narrow-banded, diagonally dominant matrices that preserve the adaptive refinement.
The multiwavelet formalism provides therefore a reliable route to compute molecular properties with guaranteed precision with respect to the complete basis set (CBS) limit.This has already been illustrated for excitation energies 35 and for electric polarizabilities. 36,37In this work, we consider magnetic properties, 44 and in particular magnetizabilities and NMR shielding constants.
The rest of this paper is organized as follows: in Section 2 we briefly summarize the theoretical foundation for the calculation of molecular energies and magnetic response properties in a multiwavelet basis.In Section 3 we present the evaluation of second-order molecular magnetic response properties, with special attention to the molecular magnetizability and nuclear magnetic shielding tensor.Section 4 details the computational protocol of our calculations.In Section 5 our results are presented and discussed.We end the paper with a brief summary and outlook.

Static linear response equations
The MW formalism has been applied to the ground-state energy 32,33 and geometry, 34 as well as dynamic polarizabilities 36,37 and excitation energies 35,38,39 in a linear response formalism.In the present work, we further extend the MW toolbox to include two important second-order magnetic properties: the magnetizability and NMR shielding tensors. 44,45In the following we present the static linear response equations, as our derivation differs slightly from the dynamic equations presented previously in the literature.
In order to solve the SCF problem, the action of the MW representation is required for only two operators, which are of the same form: where k 2 = 0 corresponds to the Poisson operator (used for electrostatic potentials), and k 2 4 0 is the bound-state Helmholtz (BSH) operator (used in the iterations of the SCF equations).
Our MW implementation for functions and operators has been described previously, 40,41 where the parallel performance as well as the inherent linear scaling of the algorithms have been demonstrated in the case of electrostatic Coulomb potentials.A derivative operator is also required for the kinetic energy operator, density gradients for GGA functionals, and angular momentum operators.We have implemented a derivative operator following the work of Alpert et al. 46

Unperturbed system
The ground-state SCF problem can be written in a general, non-canonical form to allow for localized molecular orbitals (MOs) 47,48 where is the Fock matrix, and F ˆ= T ˆ+ V ˆthe Fock operator.The potential operator in Kohn-Sham DFT and HF theory can be written in a general form as where 0 r l r 1 gives the amount of exact exchange.The electron interaction operators Ĵ j p E ¼ ð rðr 0 ; r 0 Þj p ðrÞ r À r 0 j j dr 0 ; (5) are expressed in terms of the one-particle density matrix r(r,r 0 ) which is constructed from the occupied MO The solution of the unperturbed self-consistent field (SCF) problem follows closely the work of Harrison, Yanai and coworkers, 32,33 where the SCF equations are rewritten in integral form 49 and iterated until convergence, including a Krylov subspace accelerated inexact Newton (KAIN) iterative accelerator. 50The BSH operator G ˆ(eqn (1)) is the inverse of the kinetic energy operator, shifted by the diagonal element of the Fock matrix 2G ˆi = (T ˆÀ F ii ) À1 .The second term in the operator argument in eqn ( 9) corrects for the use of non-canonical (localized) orbitals, and vanishes if the Fock matrix is diagonalized.

Perturbed system
Adding a small static perturbation h ˆ(1) to the unperturbed Hamiltonian h ˆ(0) will lead to small changes in the orbitals, and the perturbation in the density can be expressed in terms of the unperturbed {j (0) } and first-order perturbed {j (1) } orbitals.Introducing the density operator (the projector onto the occupied orbital space) and retaining terms up to first order, we get where r(0) is the ground-state density operator and r(1) the first-order perturbed density.This change in electron density in turn changes the potential operators in eqn (5-7), and to first order we get Ĵð1Þ j p E ¼ ð r ð1Þ ðr 0 ; r 0 Þj p ðrÞ r À r 0 j j dr 0 ; (13) Setting up the SCF problem for the perturbed system to first order leads to the modified Sternheimer equation 51Fð0Þ j where the perturbed Fock operator contains both the explicit perturbation and the induced perturbations in the potential operator (F ˆ(1) = h ˆ(1) + V ˆ(1) ).The perturbed Fock matrix elements are obtained by expanding eqn (3) and collecting all the firstorder terms: The orbitals satisfy the following weak orthogonality condition hj (0) i |j (1) j i + hj (1) i |j (0) j i = 0, (18)   which is equivalent to the idempotency condition on the density operator r2 = r to first order.However, in the diagonal terms the orbital phase factors can be chosen arbitrarily, and the off-diagonal terms do not contribute to the density to firstorder.Therefore one can in both cases impose the stronger orthogonality condition: hj (0) i |j (1) which has been proposed as a way to speed up convergence. 52n the other hand, for larger systems, the orthogonalization procedure becomes prohibitively expensive, and it would be more efficient to ignore these redundant projections. 53Using the strong orthogonality condition, the first and third terms in eqn (17) will vanish, and we get Fð1Þ j The Sternheimer equation (eqn ( 16)) can now be written in integral form, in the same way as the unperturbed SCF This journal is © the Owner Societies 2016 equation (eqn (9)): The diagonal element of the unperturbed Fock matrix appears in the BSH operator as 2G ˆi = (T ˆÀ F (0) ii ) À1 .We remark that eqn (21) defines a set of coupled equations which corresponds to the static limit of the response equations of Sekino et al. 36,37 and closely resembles the working equations of Yanai et al. 35,38 for excitation energies in time-dependent HF/DFT.Kottmann et al. 39 recently used the same equations to compute excitation energies in the configuration interaction singles (CIS) approximation.

Magnetic properties
A general second-order magnetic property M can be expressed as an energy derivative with respect to two parameters a and b 44,45 The property will have two contributions, one from a secondorder interaction operator h ˆ(a,b) , and one from a pair of firstorder interaction operators h ˆ(a) and h ˆ(b) .The former is known as the diamagnetic contribution, and is computed as an expectation value of the second-order operator and the zeroth-order (unperturbed) density (r = r 0 denotes the trace of the density matrix, and will always be assumed in the following) rð0Þ ĥða;bÞ dr ¼ X i j ð0Þ i D ĥða;bÞ j The latter is known as the paramagnetic contribution and is computed by perturbing the system (i.e.solving the response equations) using one of the operators h ˆ(a) to get the corresponding perturbed orbitals j (a) and density r(a) , and computing the expectation value (tracing the density matrix) of the second operator h ˆ(b)

Magnetizability and NMR shielding
The magnetizability tensor n and NMR shielding tensor s are second-order magnetic properties that can be identified as energy derivatives with respect to the external magnetic field B and the nuclear magnetic moment M K associated with nucleus K and are thus computed from the following expressions where m, n = x, y, z are the components of the perturbing fields and r(B) m denotes the density perturbed by the m component of the h ˆ(B) operator.The interaction operators are obtained by differentiating the Hamiltonian with respect to B and/or M K and evaluating at zero perturbational strength.For closed-shell systems we get ‡ 44,45 ĥðBÞ ¼ 1 2 where l ˆjO = Àir jO Â r j is the angular momentum operator, r jO is the position of electron j relative to the gauge origin O, r K is the position of nucleus K, and a E 1/137 is the fine-structure constant.
The perturbed orbitals j (B) are obtained by solving the Sternheimer equations (eqn ( 21)) with the following perturbed Fock operator (there is one for each component m of the perturbing field) Notice that the perturbed density (in real space) vanishes for pure imaginary perturbations, so J ˆ(1) and V ˆ(1) XC (if there is no explicit current dependence in the functional) can be omitted from this Fock operator.5][56] In our case, however, we only have the occupied orbitals, and the equations must still be solved iteratively to sample the (complete) virtual space by application of the BSH Green's function.Notice also that even if the perturbed density vanishes, the paramagnetic expectation value (e.g. with h ˆ(MK) ) does not r ðBÞ m ðr; rÞ 0; The above expressions can be directly related to others, such as Ramsey's 57 original sum-over-states expression for the shielding tensor where h (B) , h (MK) , h (B,B) and h (B,MK) now denote the corresponding matrix representations in the chosen AO basis.As a side note, we note that in the present formulation, due to the interchange theorem, 58,59 the order of the perturbations can be swapped, so that in the case of NMR shielding constants, the response equations are solved using h ˆ(MK) as perturbation operator instead of h ˆ(B) , i.e.
ð rðBÞ A situation where such a swapping of operators is useful is when computing shielding tensors of selected nuclei in large molecules.Whereas h ˆ(B) is a global operator affecting the entire molecule (plus any explicit solvent molecules), the operator h ˆ(MK) is localized around nucleus K (decays as r À2 following Biot-Savart's law for induced magnetic fields), leading to localized perturbations that can be treated much more efficiently than the global perturbations arising from h ˆ(B) .This approach can also be used in AO-based formulations in combination with gauge-including atomic orbitals (GIAOs), however, when the swap is employed to achieve linear scaling, 60 it requires the implementation of a number of new contributions.1][62] However, one has to keep in mind that the number of response equations to be solved increases with the number of selected nuclei, and it eventually becomes beneficial to instead perturb the system with the global h ˆ(B) operator, which requires the solution of only three response equations for all nuclei.

Computational details
All Gaussian-type orbital (GTO) calculations were performed with the Dalton program, 63 using Dunning's correlation-consistent (cc-pVXZ 64 ) and Jensen's polarization-consistent (pcS-n 65 ) basis sets.The latter was specifically optimized for computing NMR shielding constants using DFT.The calculations were performed using GIAOs unless otherwise specified.
The MW calculations were performed with a development version of the MRChem program package. 66The exchangecorrelation functionals and their derivatives to arbitrary order are provided by the XCFun library. 67In the results presented, k denotes the polynomial order of the basis, e is the overall numerical accuracy used in the calculation (the internal threshold used for truncating the MW representations of all functions and operators) and Df is the convergence threshold in the orbital residuals (both ground-state and response).As starting guess for the ground-state calculations we used converged wave functions in small GTO basis sets (e.g.cc-pVDZ).A zero initial guess was used for the response functions, and the convergence was a bit slower than the corresponding ground-state calculation (for the ground state, typically 2-4 iterations are needed to gain one order of magnitude in accuracy, for the response equations one or two iterations more, but with a larger KAIN iterative history).

Results
This section is divided in two parts.In the first we calibrate and benchmark MRChem both for the convergence towards the CBS limit and the origin dependence, through some test calculations of magnetizabilities and NMR shielding constants of the second-row hydrides (geometries given in Table 1).Then we show the MRChem performance on the challenging case of magnesium oxide, which has been shown to ''display pathological behaviors with respect to basis set convergence''. 65 In the second part, we perform a statistical analysis of a larger set of molecules, originally considered by Lutnaes et al. 21for magnetizabilities and later by Teale et al. 22 for shieldings.We will here assess the quality of the density functionals and basis sets that are typically being used for such calculations.2 shows the Hartree-Fock magnetizability of the water molecule computed at two different gauge origins using MW and GTO basis sets of different quality.The MW calculations are grouped in three different overall numerical accuracies with a factor of 100 between them (e = 10 À3 , 10 À5 , 10 À7 ), and for each of these we look at the convergence of the property with respect to the norm of the orbital residuals (Df).The overall precision e should give the maximum relative accuracy that we are able to obtain as the orbitals converge, which is reflected in the total energies when we compare with the reference value of Yanai et al. 33 Even though the convergence of the energy is quadratic in the convergence of the orbitals, the final error in the converged energy will be given by the overall accuracy of the MW calculation.
For magnetizability calculations with a gauge origin within the molecular geometry (r O = (0, 0, 0), close to the center of mass), the error in the diamagnetic contribution is expected to be linear in the error of the ground-state orbitals (Df) and limited by the overall accuracy (e) of the calculation.This can be seen from the numbers in Table 2, and the error in the last (most accurate) MW value is expected to be around 10 À6 atomic units.The paramagnetic contribution depends on both the ground-state and perturbed orbitals, and the absolute error is similar to the error in the corresponding diamagnetic contribution, which means that we get a consistent accuracy in the total magnetizability of the molecule (that is, around 10 À6 a.u.).For comparison we see that the cc-pV6Z basis set (including GIAOs) only yields 10 À2 a.u.precision for the water magnetizability, while augmenting the basis set with extra diffuse functions increases the accuracy by two orders of magnitude.
5.1.2Gauge origin dependence.It is well known that the (arbitrary) placement of the gauge origin will affect the quality of the result whenever an incomplete basis set is used in the calculation of magnetic properties.Several solutions to the problem have been proposed.Nowadays, the most common approach is to use GIAOs, first introduced by London in 1937. 68The method was further developed by Ditchfield in the 1970s, 69 but was not made efficiently applicable until 1990, with the work of Wolinski, Hinton and Pulay. 70In the London orbital approach, each AO is made explicitly dependent on the external magnetic field, and a local gauge origin is defined at the center of each AO.In this way, the magnetic properties become gauge-origin independent by construction, though not gauge invariant and current conserving. 71,72he MW basis is in principle complete up to the predefined truncation threshold.Therefore, we expect a much less severe gauge dependence in the computed magnetic properties (effective gauge origin independence as well as gauge invariance to within the accuracy threshold), and very high accuracy is attainable for small molecules by choosing a common gauge origin within the molecular geometry (e.g.center of mass).However, as the accuracy in the response property is relative to the absolute value of the paramagnetic contribution (the diamagnetic contribution is more accurate because it only depends on the unperturbed system), we also expect that it becomes progressively more difficult to attain the same accuracy when the origin is moved out of the molecular framework, as both the diamagnetic and paramagnetic contributions to the magnetizability will increase in magnitude.
Table 2 shows also the same magnetizability computed with a displaced gauge origin (r O = (5, 5, 5)).For the MW calculations, we see that the diamagnetic part is computed with the same absolute accuracy as before (comparing e = 10 À5 , Df = 10 À5 results with e = 10 À7 , Df = 10 À7 , the variations are on the sixth decimal place in both cases), but for the paramagnetic part, only the relative accuracy is maintained (variations in the fourth digit after the comma).This means that two digits are lost in the total magnetizability in this case, as the magnitude of the paramagnetic part is increased from B1 to B100.The MW basis still provides reliable and systematically improvable results, but this shows that the origin should not be chosen Table 2 Hartree-Fock magnetizability (a.u.) of water computed using two different gauge origins r O , relative to the molecular geometry given in Table 1.GTO: (0, 0, 0) computed using GIAOs, (5, 5, 5)  arbitrarily (and preferably somewhere within the framework of the molecular structure).For comparison, traditional GTOs are dependent on the GIAO parametrization to yield reasonable results.Whereas the GIAO results would be independent of the choice of gauge origin up to the same numerical issues as observed for the multiwavelet basis, 73 without GIAOs, if the gauge origin is moved away from the center of the Gaussians, even the largest basis sets are unable to yield a correct result.This is clearly shown in Table 2, where the magnetizability value for water is correct to B10 À4 a.u.for the aug-cc-pV6Z basis when London orbitals are employed, but becomes 30% larger with a displaced gauge origin (r O = (5, 5, 5)) without London orbitals.
5.1.3NMR shielding constants.Table 3 shows the NMR shielding constants of the second-row hydrides computed using the B3LYP 74 functional and MW and GTO basis sets of different quality.As for the magnetizability, the MW calculations are grouped in three different overall accuracies (e) and we present the convergence of the total NMR shielding constant with respect to the orbital residuals (Df).The gauge origin is chosen close to the center of mass of each molecule (r O = (0, 0, 0) relative to the geometries given in Table 1), whereas the GTO calculations are performed using GIAOs.The final MW shielding constants are expected to be accurate to at least 10 À3 ppm for the second row elements, and 10 À4 ppm for the hydrogen atoms.For comparison, shielding constants computed with the largest non-augmented Gaussian basis set (pcS-4) are correct to B0.1 and B0.01 ppm, respectively, and augmenting with extra diffuse functions does not significantly improve the results.The performance of the pcS-n basis sets are thus, as expected, quite good, because they are specifically optimized for NMR shielding calculations.

Magnesium oxide.
A molecule which has proven difficult to handle with traditional basis sets is magnesium oxide.In the calibration of the pcS-n basis sets, Jensen decided to remove a handful of molecules from the original test set, as their errors were so large that they would have ruined the statistics. 65The worst of these systems was MgO, and Table 4 shows the shielding constants of this molecule computed with three different functionals (B3LYP, Becke's half-and-half 75 functional and Hartree-Fock), using both MW and GTO bases.Whereas B3LYP starts out with a massive overestimation of the shielding constant for the smallest pcS-0 basis set, we do observe a systematic improvement and monotonic convergence for this functional when we increase the cardinal number of the basis.This was also observed and analyzed by Jensen. 65For the BHandH functional, however, the convergence is less systematic and we have to go to the pcS-2 basis and beyond in order to get qualitative agreement between the numbers, and the situation gets even worse for Hartree-Fock, where anything less than pcS-3 gives completely erratic results, and even the biggest basis is only converged to one digit.Again, the augmented basis sets show only slight improvements over the standard pcS-n.
In the MW basis, on the other hand, we see a systematic convergence for all functionals, although a bit slower than for the second-row hydrides presented above.The final MW shielding constants should be accurate to within 1 ppm, except for the Hartree-Fock oxygen value, where the error is on the order of 10 ppm.

Benchmarking HF and DFT with CBS results
In two recent works, Lutnaes et al. 21and Teale et al. 22 investigated the performance of a number of density functionals to assess their ability to reproduce magnetizabilities and shieldings, respectively.In these studies, 28 molecules were considered, comparing the DFT results with coupled cluster results and with experimental data.In order to minimize basis set errors, they used large basis sets (the largest basis set used was the correlation-consistent, core-valence aug-cc-pCVQZ basis set of Woon and Dunning. 76) They also made use of an extrapolation method with a two-point exponential parametrization for HF and a polynomial extrapolation for the correlation part.We report here the extrapolation formula for the HF part: e aY e aX À e aY : In particular, for the exponential extrapolation of the HF values, they used the same parametrization (a = 1.63) that is employed for the molecular energy. 77This choice was justified by the fact that second-order molecular properties are energy derivatives.With a MW approach, it is possible to provide very accurate benchmark results that can be employed to check the quality of the GTO results, both for HF and for the available density functionals.Additionally, we can determine whether the exponential extrapolation procedure leads to an improvement in the results, as well as investigate to what extent it predicts the correct basis-set limit.
We have therefore considered the same set of molecules and performed calculations at increasing precision (e = 10 ÀZ , Z = 3, 4, 5, 6), to compare our MW results with the results obtained using GTOs.We will here summarize our main findings.
Based on the initial test calculations on the first-row hydrides, we will in all the following MW calculations converge the orbitals (both ground-state and perturbed) to 10e.The reason for this is twofold: firstly, the property is usually converged within the expected error bars at this point, and secondly the convergence might be affected by numerical noise when we approach the limit of the guaranteed accuracy of the computation.
5.2.1 Magnetizabilities.The HF magnetizabilities computed with MRChem are reported in Table 5, together with the reference values of Lutnaes et al. 21The progression of the MW results clearly show how the results gain consistently in accuracy when Z is increased, and our most accurate results are converged to the fifth digit.In more detail, the MW3 results have 1-2 correct digits, MW4 have about three and MW5 four correct digits.The only notable exception is ozone for MW3, which is also qualitatively wrong.
The comparison with the best GTO results shows a very good agreement with differences of about 0.1-0.2Â 10 À30 J T À2 .Ozone proves to be a challenging system for GTOs with a deviation of 1.4 Â 10 À30 J T À2 .We can therefore conclude that the aug-cc-pCVQZ basis is able to attain very good accuracy for magnetizabilities, yielding results that are comparable to our MW5 values.
A very similar picture is obtained for the density functionals examined, with aug-cc-pCVQZ performing on par with MW5.Our MW6 results are reported in Table 6 together with the augcc-pCVQZ results from Lutnaes et al. 21 detailed error analysis, and the accuracy of the extrapolation formula are not discussed because of the agreement down to the last available digit in the GTO results, and rounding effects would therefore play a major role.
5.2.2 NMR shielding constants.The NMR shielding constants computed with different basis sets for Hartree-Fock and DFT are presented in Tables 7 and 8, respectively, and the statistical basis set errors of the different methods are reported in Table 9, taking the MW6 values as reference.We begin by considering the progression along the series of MW calculations for Hartree-Fock.The MW3 results yield errors of a few tens of ppm, whereas MW4 are accurate to within 1-2 ppm and MW5 results are on average 0.1-0.2ppm away from the MW6 values.We can therefore conclude that our best MW results are accurate to within 0.1 ppm or better.It is interesting that ozone, which is completely wrong at MW3, becomes in line with the other cases when MW5 and MW6 are compared.This is a clear indication of the systematic nature of the MW basis, which is able to provide the flexibility required to achieve the CBS result throughout the series.
We can now compare MW6 results to the GTO values.In particular, the best aug-cc-pCVQZ results are accurate to within 1-2 ppm (MAE = 1.8/1.2with/without O 3 ), i.e. comparable to MW4.The extrapolated values (a = 1.63) show an improvement over the aug-cc-pCVQZ ones, and the error is reduced (MAE = 1.1/0.7),although not substantially.However, if we change the parametrization to a = 1.05, a much closer agreement is achieved (MAE = 0.38/0.25 with/without O 3 ).Other measures for the error yield a similar picture.This result suggests that the extrapolation procedure, which is originally justified because the shielding is a second derivative of the energy, is a reasonable way to estimate CBS results.However, making use of the same exponent as the molecular energy, is improving over the aug-cc-pCVQZ values only slightly.If the exponent is set to 1.05 (larger corrections than a = 1.63) the improvement is substantial.Considering that a is the only adjustable parameter for 72 shielding values, this result is a confirmation of the validity of the two-point exponential parametrization.
We now turn our attention to DFT, where for clarity we present the raw data only for quadruple-z and MW6 for the chosen functionals in Table 8, while the statistical errors are given for all basis sets in Table 9.The sequence of MW calculations at increasing precision paints a similar picture to that seen for HF: MW3 yields results that are several ppm away from the MW6 values and also worse than triple-z quality numbers, MW4 is a huge leap forward, outperforming the quadruple-z results.MW5 yields results that are in most cases 0.1-0.2ppm away from the MW6 values.This progression is  again an indication that the error in the MW6 values can be considered well below 0.1 ppm from the CBS results, and MW6 values can in practice be as a reference.
The MW6 reference values enable us to apply to DFT the same extrapolation procedure used for HF.With a = 1.63, we observe a reduction of the errors (e.g. for B3LYP the MAE goes from 1.59 to 0.9).With a = 1.05 the reduction is more marked (for B3LYP MAE = 0.28/0.24with/without O 3 ).Slightly better results could be achieved if the parameter a was optimized for each functional separately, e.g. using a = 0.9 leads to a further substantial reduction in the standard deviation in the Hartree-Fock results, but we observe the opposite effect for all DFT functionals.The main point is, however, that neither for HF nor for DFT is a = 1.63 the optimal value, and our results indicate that a lower value, yielding a larger correction, performs much better for HF and DFT alike.

Computation times.
In Table 9, we also report the average timings for the GTO as well as the MW calculations.All timings include the ground-state SCF optimization (to obtain r(0) ) followed by three response calculations (to obtain r(B) m for m = x, y, z).The magnetizability and NMR shielding tensors of all nuclei in the molecule are then computed simultaneously according to eqn (26) and (27).The timings are given in minutes wall-time on a single compute node with 16 CPUs (2.6 GHz Intel E5-2670) using a shared memory parallelization strategy.For the GTO calculations, the timings are obtained from Dalton, on the same computer architecture and also using 16 processors using the message passing interface (MPI).
Whereas the MW code is far from optimized for production calculations, we can draw some interesting conclusions from the numbers.First of all, there is an order of magnitude difference in computation time between pure and hybrid density functionals: for HF, PBE0 and B3LYP, the exact-exchange is a clear bottleneck, especially in the magnetic response solver because the pure functionals do not contribute to the perturbed Fock operator (eqn (31)), but also in the ground-state calculation as the exact exchange scales quadratically with the number of orbitals.Although the current implementation fully supports orbital localization, which should have an effect on the exact exchange, no particular attempt has yet been made in exploiting localization to achieve linear scaling algorithms.The systems treated in the current work are anyway too small for this to have any appreciable effect.
Along the series of MW calculations, we observe that the computational costs roughly doubles at each ten-fold increase of the requested precision.
Concerning the comparison between MW and GTO, and limiting the comparison to pure functionals (hybrid functionals and HF are affected by the exchange bottleneck), we notice how TZ-quality results are quite cheap to achieve, being ten times faster than MW3 results.However MW4 results, which are comparable to the extrapolated [TQ]Z values, are only three times more expensive than QZ values: along the GTO series, the computational cost increases ten times at each step (aug-cc-pCV5Z confirmed the trend, for the 24 molecules where this basis set is available).
In conclusion, GTOs come from decades of developments, where a large effort has been poured into improving the underlying algorithms and fine-tuning the basis set compositions.This explains why a moderate basis set, such as TZ, is able to give reasonable results with good performance.But when the precision requirements are increased, a MW basis becomes eventually superior.
5.2.4Comparison with coupled cluster and experiment.With the proper CBS limits established for the different DFT functionals, we can assess their ability to reproduce the theoretical limit: full configuration interaction (FCI) in a CBS.In practice coupled cluster singles and doubles with perturbative triples corrections (CCSD(T)) including basis set extrapolation has been taken as a reference. 22Additionally, experimental data including zero-point vibrational corrections (ZPVCs) have been considered.Several investigations of this type have been conducted in the past, 22,[78][79][80][81] and our observations agree with the current consensus regarding the CBS limit.It is well known that the employed density functionals systematically underestimate the shielding constant: this has been attributed to a too small HOMO-LUMO gap. 78Some attempts have also been made to fix the problem by simple level shifting of the virtual orbital energies, 78,82,83 although the theoretical justification is questionable.
We present in Table 10 the statistical errors with respect to coupled cluster (CC) the empirical equilibrium (experimental values including ZPVC).In the comparison to CCSD(T), O 3 is omitted, because it is an extreme outlier with errors on the order of 1000 ppm.In the comparison to experiments, around one fourth of the data set is instead omitted because accurate experimental data is lacking.The same data set as used by Teale et al. 22 has been used, and details about the excluded systems can be found in their paper.All GTO calculations, as well as all experimental numbers, are taken from the work of Teale et al. 22 Focusing first on the CBS limit, represented here by the MW6 numbers, we observe an underestimation of 20-30 ppm, without significant variations throughout the set of functionals employed.The magnitude of the mean error (ME) and the mean absolute error (MAE) is about the same, indicating that the errors are systematic in sign.It is interesting to note that the cancellation between functional error and basis set error is rather systematic: a more or less uniform improvement of the results, compared to both CCSD(T) and experiment is obtained when a poorer basis set is employed.In particular, the aug-cc-pCVDZ basis set displays an improvement of around 10-15 ppm in the MAE, compared to the CBS limit, accompanied by a significant reduction in the standard deviation (SD).This trend has been already reported by Kupka and coworkers, [79][80][81] and has also been recommended to NMR spectroscopists. 84inally, taking experimental values without ZPVC as reference (not shown in the tables), the errors for generalized gradient approximation (GGA) and hybrid functionals with the double-z basis are further reduced to around 8 and 20 ppm for the MAE and SD, respectively.The combination of several errors eventually gives the right answer in the case of NMR shielding calculations with DFT.This demonstrates the challenges of optimizing functionals explicitly for a specific property, as for instance done for the Keal-Tozer functionals for magnetic properties. 85

Conclusions and outlook
We have presented a new real-space implementation of a static linear response solver in a multiwavelet framework.We have applied the formalism to compute magnetizabilities and NMR shielding constants of small closed-shell molecules using HF and Kohn-Sham DFT.We have shown that the MW basis provides reliable numerical results for a wide range of molecular systems, including challenging cases such as MgO and O 3 .The accuracy of the MW basis depends only on a single input parameter: the overall relative precision e.By tightening e, the CBS values can be attained reliably and consistently.While high accuracy is attainable also with traditional GTO bases, they rely on a careful preoptimization of the parameters, and a wide variety of customized basis set families are available for different molecular properties, and it is the users responsibility to choose a basis set that is suited for the problem at hand.
We have shown that the MW basis provides magnetic response properties that are gauge invariant within the chosen relative precision, which means that the method does not rely on a GIAO parametrization using field-dependent phase factors for each AO (which the implementation enormously) in order to provide high-accuracy results for small molecules.For larger systems, however, we expect a reduction of the overall accuracy, as the size of both the paramagnetic and diamagnetic contribution will grow faster than the total property, affecting the relative precision that can be achieved.In order to treat larger systems in the future, we plan to implement a local-origin method, such as the individual gauge for localized orbitals (IGLO) method of Kutzelnigg 86 or the localized orbital/local origins (LORG) method of Hansen and Bouman. 87As in the London orbital approach, these methods also use different origins for different orbitals, but at the MO rather than AO level, placing the origin at the centroid of charge of each (localized) orbital.
While these methods also lead to gauge-origin independent results, they rely on orbital localization, and in a GTO calculation also on the use of resolution-of-the-identity to simplify certain complicated integrals.The GIAO approach has therefore dominated the calculation of magnetic properties during the past two decades.However, some of these problems should not be present in the MW framework because the exact placement of the origin does not affect the quality of the results, as long as it keeps the size of the paramagnetic contribution as low as possible, and the MW basis satisfies the resolution-of-the-identity condition to within the selected precision.By employing local origin methods, the paramagnetic part can be decomposed into orbital contributions of much smaller magnitude (as the origin is chosen individually for each localized orbital), leading to good accuracy also for larger molecules.
Table 10 Statistical errors of shielding constants with respect to extrapolated CCSD(T) and empirical equilibrium (experiment minus zero-point vibrational corrections) results for the different methods employed in the paper.We report the basis set convergence of aug-cc-pCVXZ (X = D, T, Q) as well as MW6, which can be considered converged to the CBS limit.For the CCSD(T) comparison O 3 is omitted, while only the molecules where good experimental values are available were included in the comparison with empirical equilibrium (see text for details).All GTO calculations, as well as all experimental numbers, are taken from ref.We have employed our method to probe the CBS limit of HF and DFT for a large test set of small molecules, initially proposed by Lutnaes et 21 for magnetizabilities and later Teale et al. 22 for NMR shielding constants.We found that the aug-cc-pCVXZ basis sets performed very well for magnetizabilities, yielding typical errors of B0.1% for quadruple-z compared to the MW basis-set limit results.As expected, the NMR shielding constants were more challenging for GTOs, with an average error of 1-2 ppm (around 2-3%) for the largest aug-cc-pCVQZ basis.In order to improve on these numbers, a two-point exponential basis set extrapolation formula was employed.We found that the parametrization that is commonly used for total energies does not transfer directly to properties, and a significantly lower exponent (close to 1.0 rather than 1.63) should be used for the estimation of the basis-set limit for the shielding constants.Although the optimal value might differ somewhat between different functionals, a common exponent of 1.05 was able to reduce the errors by almost an order of magnitude for all functionals that were considered, compared to the quadruple-z value.
Finally, we performed a comparison between the CBS limit for HF and DFT and accurate CC and experimental values for the shielding constants.We found that the original study by Teale et al. 22 using GTO basis sets were close enough to the CBS limit for their conclusions with respect to CC and experiment to hold.We found, however, that there is a significant systematic -yet fortuitouscancellation of errors between the functionals and the basis set: with respect to raw experimental values (not zero-point vibrational energy corrected), surprisingly accurate shielding constants can be achieved with a moderate aug-cc-pCVDZ basis set.
The current implementation is limited to SCF levels of theory, and since correlated wave function methods rely on at least sixdimensional numerical representations, they are extremely expensive in a MW framework.][90] The performance of our MW approach, compared to GTO calculations, gives a mixed picture.At moderate accuracy, GTOs deliver reasonable values with limited effort (aug-cc-pCVTZ), whereas our MW3 values are not as accurate and computationally more expensive.The picture changes dramatically if higher accuracy is demanded: increasing the accuracy with GTOs implies a 10-fold increase in the computational cost, whereas MW calculations are only twice as expensive at each step along the MWX (X = 3, 4, 5, 6) series.As a result, aug-cc-pCVQZ calculations are only 2-3 times faster but less accurate than MW4.Some preliminary tests show that MW5 is indeed less expensive than the aug-cc-pCV5Z basis.The main performance bottleneck in the current MW implementation is clearly the exact HF exchange, which yields a 10-fold increase in the timings compared to pure functionals.

Table 1
Geometries (Bohr) of the second-row hydrides used in the benchmark calculations This journal is © the Owner Societies 2016 computed without GIAOs, number of contracted functions in parenthesis.See Section 4 for computational details 33Reference energy from Yanai et al.33.This journal is © the Owner Societies 2016 Phys.Chem. Chm.Phys., 2016, 18, 21145--21161 | 21151

Table 3
B3LYP nuclear shielding constants (ppm) of second-row hydrides.Geometries given in Table1.See Section 4 for computational details

Table 6
Magnetizabilities (10 À30 J T À2 ) computed using a range of density functionals with the aug-cc-pCVQZ basis set and the MW basis at precision e = 10 À6 .All GTO calculations are taken from the ESI of ref.21 This journal is © the Owner Societies 2016

Table 8
Shielding constants (ppm) computed using a range of density functionals with the aug-cc-pCVQZ basis set and the MW basis at precision e = 10 À6 .All GTO calculations are taken from the ESI of ref.22

Table 9
Statistical basis set errors of all computed shielding constants (including O 3 ) for the different methods employed in the paper, MW calculations at e = 10 À6 are taken as reference.For each combination of basis/ functional the following errors are reported: medium error, medium absolute error, medium relative error, medium absolute relative error, standard deviation and maximum error.Average timings in minutes on 16 CPUs aug-cc-pCV[TQ]Z a = 1.05 ME 0.18 À0.13 À0.04 À0.03 À0.01 22.For each combination of basis/functional the following errors are reported: medium error, medium absolute error, medium relative error, medium absolute relative error, standard deviation and maximum error This journal is © the Owner Societies 2016 Phys.Chem.Chem.Phys., 2016, 18, 21145--21161 | 21159