Nuclear magnetic resonance chemical shifts with the statistical average of orbital-dependent model potentials in Kohn–Sham density functional theory

In this paper, an orbital-dependent Kohn–Sham exchange-correlation potential, the so-called statistical average of (model) orbital potentials, is applied to the calculation of nuclear magnetic resonance chemical shifts of a series of simple molecules containing H, C, N, O, and F. It is shown that the use of this model potential leads to isotropic chemical shifts which are substantially improved over both local and gradient-corrected functionals, especially for nitrogen and oxygen atoms. This improvement in the chemical shift calculations can be attributed to the increase in the gap between highest occupied and lowest unoccupied orbitals, thus correcting the excessively large paramagnetic contributions, which have been identified to give deficient chemical shifts with both the local-density approximation and with gradient-corrected functionals. This is in keeping with the improvement by the statitical average of orbital model potentials for response properties in general and for excitation energies in part...


I. INTRODUCTION
In the last decade, density functional theory ͑DFT͒ ͑Refs. 1-3͒ has been increasingly employed for the prediction of molecular properties and for the study of chemical reactivity in many areas of chemistry. 4,5 The success of DFT relies on its good accuracy achieved at a relatively low computational cost, due to the use of various explicit approximations to the exact exchange-correlation ͑XC͒ functional. [6][7][8][9] The main goal of the present paper is to further improve the accuracy which can be obtained from such calculations by improving the approximation to the exchange-correlation potential which is usually the main source of error, 10,11 thus moving closer to the benchmark accuracy obtained in, e.g., coupled-cluster calculations.
It has been demonstrated that the quality of the KS potential is a crucial factor in obtaining accurate response properties. 10 A particular successful modeling of xc is the so-called statistical average of different model orbital potentials ͑SAOP͒, yielding an orbital dependent expression for the potential. 12,13 The distinguishing feature of this new potential is that it incorporates physically well motivated features in both the asymptotic region and the ͑sub͒-valence and core regions. In previous papers, [13][14][15][16][17] the SAOP potential has been shown to provide high quality results for a whole series of static and frequency-dependent molecular response properties and to provide a substantial improvement upon the LDA and GGA potentials.
In view of the improvement provided by the SAOP potential for a range of response properties, 13,[15][16][17] we now want to investigate it for NMR chemical shifts. The main reason is that, although traditional wave-function-based correlated methods achieve an acceptable accuracy, 18 these methods are not attractive for application to large molecules. An alternative route to high accuracy NMR shift calculations for possibly very large molecules would require an improved exchange-correlation potential, which should still be computationally expedient, i.e., not essentially more expensive than the GGA potentials.
In the present work, NMR calculations by means of the SAOP potential have been carried out for a series of 44 simple molecules, only containing H, C, N, O, and F. For comparison purposes this series has been chosen to be the same as the one already studied by Patchovskii et al. 19 These authors implemented the Perdew-Zunger 20 self-interaction correction ͑SIC͒ within a molecular DFT program ͑ADF͒, using the Krieger-Li-Iafrate 21,22 approximation to the optimized effective potential method, and the Vosko-Wilk-Nusair ͑VWN͒ functional ͑SIC-VWN͒. This SIC-VWN potential was shown to achieve a nice improvement with respect to VWN and BP86 potentials. More recently it has been shown 23 that using a SIC-GGA functional gave a slight further improvement, although not universally so. The correction to the Kohn-Sham potential by the introduction of the SIC is much more important than changing the LDA ͑VWN͒ to a GGA ͑BP86 or revPBE͒ functional. We will also compare the SAOP potential to results obtained recently by Wilson and Tozer, 24 using the same series of 21 simple molecules that these authors have used to test various ab initio and DFT levels of theory. Of special interest are the results that they achieved by means of the application of Kohn-Sham potentials, which are determined from theoretical electron densities. They call the local KS potentials multiplicative ͑MKS͒; associated orbitals and eigenvalues are used to calculate the NMR shielding constants. Specifically, Wilson and Tozer 24 achieved the best results with the MKS͑BD͒ potential, which is obtained from Brueckner Doubles calculations as the correlated method to generate the electronic density. The constructed local ͑multiplicative͒ KS potential is able to regenerate this correlated density as the sum of orbital densities of the noninteracting KS system. They obtained results largely improved with respect to the HCTH ͑Ref. 25͒ functional or the B97-1 ͑Ref. 26͒ hybrid functional, which again indicates that the quality of the KS potential is a key factor in the NMR shielding calculations. In addition to the SAOP potential, we have also investigated the potential obtained with the gradient-regulated connection to the asymptotically correct van Leeuwen-Baerends ͑LB͒ ͑Ref. 27͒ potential ͑BP-GRAC-LB͒, see Ref. 28. The BP-GRAC-LB potential shares with the SAOP potential the correct asymptotic behavior, but is in the bulk region an unmodified GGA ͑in the present case BP86͒ potential.
In both series of molecules to be studied here, the SAOP potential presents a nice performance improvement over VWN or BP86 results. The SAOP results are in fact quite close to the SIC ͑Ref. 19͒ results, but they are still inferior with respect to the MKS͑BD͒ ͑Ref. 24͒ results. However, it is cumbersome, and for very large systems computationally very expensive, to generate the MKS potential from highquality correlated calculations. The low computational cost and the ease of use of the SAOP potential compared to the MKS͑BD͒ potential, and to a lesser extent the SIC calculations, recommend it for routine use. The BP-GRAC-LB results do not offer a similar improvement, they appear to be very close to the BP86 ones.

II. COMPUTATIONAL DETAILS
The SAOP ͑Refs. 12-14͒ xc potential xc SAOP is constructed as the statistical average over the occupied KS or- of the model orbital potentials xci mod . The latter are obtained with the interpolation, between the modified potential xc LB␣ of van Leeuwen and Baerends ͑LB͒, 27 which has the proper Coulombic asymptotics Ϫ1/r, and the potential xc GLLB of Gritsenko, van Leeuwen, van Lenthe, and Baerends ͑GLLB͒, 29 which correctly reproduces the atomic shell structure in the inner regions. With Eqs. ͑1͒ and ͑2͒ xc SAOP provides a balanced approximation to the KS potential xc in all regions.
All NMR calculations have been performed by means of the Amsterdam Density Functional program, [30][31][32][33] which already has the SAOP potential implemented. The Slater basis set used has been the standard set VII implemented in the package, and all-electron calculations were performed. This basis set was chosen after having tested a series of eight molecules with five different Slater basis sets, for the VWN potential. The absolute shieldings obtained are tabulated in Table I  •(n f ) basis functions for C, N, O, and F; and 10•(ns), 4 •(np), and 3•(nd) basis functions for H. The results of Table I show that all these basis sets are rather good, the differences between the basis sets are in almost all cases ͑except in the few cases where theory and experiment are very close͒ small compared to the difference between theoretical and experimental results. The basis sets QZ, VII, and VIIIB are slightly better converged than V and PAZ, but the differences between these three bases are insignificant. Basis VIIIB was discarded as being computationally very expensive and not improving on basis VII. The VII basis set was chosen as being not unduly large but still very accurate. The PAZ basis is sufficiently close in quality to allow us to compare the SAOP results of this work with the SIC-VWN ones.
The NMR calculations for the series of 44 molecules were carried out at the same geometries as in Ref. 19, i.e., optimized molecular geometries using the BP86 GGA functional, and the ADF standard basis set IV of TZP quality. Frozen ͓1s͔ cores were used on C, N, O, and F. On the other hand, for the series of Wilson and Tozer, 24 all calculations were performed at near-experimental geometries.
In NMR experiments, shielding constants are measured relative to a reference substance, and are given as the chemical shift. Calculations, on the other hand, produce the absolute tensors , and their traces iso . To a good approximation, the relationship between chemical shifts and shieldings is given by

͑3͒
A straightforward application of Eq. ͑3͒ to the conversion of theoretical absolute shieldings to chemical shifts tends to bias the comparisons between different theoretical techniques by assigning an excessive significance to the error in the calculated absolute shielding of the reference nucleus. Therefore we have followed the same procedure as in Ref. 19, i.e., we have adjusted ref for each theoretical method such that the average signed error in the calculated chemical shifts will be zero. The reference shieldings ( ref ), used to calculate the chemical shifts by Eq. ͑3͒, are given in Table III along with the average absolute errors.

III. RESULTS AND DISCUSSION
NMR calculations have been carried out for two series of simple molecules. The discussion will be focused first on the series of 44 molecules treated in Ref. 19 and next on the series of 21 molecules treated in Ref. 24. Thus, the chemical shifts of a selected set of 44 small molecules, calculated with the methodology detailed above using VWN, BP86, SAOP, and BP-GRAC-LB potentials, are summarized in Table II. Whenever possible, available experimental results are also given for comparison, either obtained in the gas phase, or in nonpolar noncoordinating solvents. Table III gives, apart from the reference shieldings ( ref ), the average absolute errors relative to experiment. Meanwhile, Figs. 1-4 provide a visual comparison of the residual errors in the calculated NMR chemical shifts for C, N, O, and F.
First, we will comment on the chemical shifts obtained for each atom, paying special attention to the average errors relative to experiment ͑see Table III͒. We analyze together the results obtained by the VWN, BP86, and SAOP potentials ͑Table IIIa͒. The BP-GRAC-LB results will be studied separately ͑Table IIIb͒, since the BP-GRAC-LB calculations require as input the ionization energy, for which we have used the experimental ionization potentials which are known for a subset of 29 out of the 44 molecules.
Concerning the average errors relative to experiment for the 44 molecules ͑Table IIIa͒ the nice performance of the SAOP potential can be easily recognized. It considerably improves in general the chemical shifts for carbon, nitrogen, oxygen, and fluorine, with respect to the VWN and BP86 potentials. However, for hydrogen, we note that the SAOP potential ͑0.37 ppm͒ gives an average error equal to the one obtained by the VWN ͑0.36 ppm͒, and slightly larger than the one given by BP86 ͑0.25 ppm͒. Nevertheless, as also noted by Patchkovskii et al., 19 this value of Ϸ0.3 ppm corresponds to a typical vibrational correction to the 1 H NMR chemical shifts. 34,35 In special cases the vibrational correction may even become twice as large. 36 Therefore, evaluating the performance of the SAOP potential for 1 H would require further work on the vibrational corrections.
For 13 C chemical shifts, the average error achieved by the SAOP potential ͑4.7 ppm͒ shows the improvement given by this method with respect to the VWN ͑8.9 ppm͒ and BP86 ͑6.2 ppm͒ potentials. There are only three cases (CO 2 , CF 3 CN, and CH 3 F) that are not improved by the SAOP potential, cf. Table II. Moreover, it must be mentioned that there are some cases ͓H 2 CO, (CH 3 ) 2 CO, C 3 O 2 , CH 3 NC, and CF 3 CN] that, even though they are improved with respect to BP86, present errors relative to experiment larger than 10 ppm. It has been suggested 19 that these deviations may be due to environmental effects. The 13 C correlation plot ͑Fig. 1͒ exhibits the generally nice behavior of the SAOP potential in the cases where it improves the results, while highlighting the still unsatisfactory cases with errors of 10 ppm or larger.
Nitrogen and oxygen chemical shifts present the largest improvement when applying the SAOP potential. For 14 N and 15 N chemical shifts, the VWN and BP86 errors are 41.0 and 33.6 ppm, while the SAOP average error is 20.2 ppm, a significant improvement. It is worth mentioning the case of the N 2 O 3 system ͑see Fig. 2͒, where the error in the calculated nitroso nitrogen's chemical shift is strongly reduced ͑VWN: 735, BP86: 668, SAOP: 544, expt: 366͒, but still quite large. Therefore, if we take this outlier case out of the statistical study, the average error would be reduced to 21.7, 17.8, and 10.9 ppm, for VWN, BP86, and SAOP, respectively. Pyridine is the only case where the chemical shift    The behavior of the SAOP potential when calculating oxygen chemical shifts is displayed in the oxygen correlation plot of Fig. 3, where special attention should be payed to the systems that have been singled out above.
The last atom to be discussed is fluorine. Only twelve systems have been calculated, so the statistics for the average errors is relatively poor. Nonetheless, the improvement given by the SAOP potential is worth noting, the average errors decreasing from 25 Fig. 4͒.
The BP-GRAC-LB chemical shifts have been calculated for a selection of 29 of the previous 44 systems, i.e., those for which the experimental ionization potentials are known. The average errors relative to experiment are tabulated in Table IIIb. We can easily recognize that the BP-GRAC-LB results are very similar to the BP86 ones. And qualitatively the errors found for each atom follow the same trend as for the 44 systems. It has been shown 17 that GGA potentials like BP86 have, compared to asymptotically correct potentials like BP-GRAC-LB, the Rydberg excitations way too low and in addition induce many spurious low-lying Rydberg-type excitations. They also produce a distorted spectral structure of response properties like the polarizability ␣ and the Cauchy coefficient S Ϫ4 . However, the valence excitations are rather similar between the BP86 and BP-GRAC-LB potentials. A similar observation can be made for the energies of the virtual orbitals of valence and Rydberg type which directly enter the expressions for the chemical shielding in the coupled perturbed Kohn-Sham formulation. The observed similarity between BP86 and BP-GRAC-LB chemical shifts can then be understood if the valence excitations have a dominating effect on the chemical shifts, while the Rydberg excitations do not have a large influence.
The deficiencies in the GGA-DFT NMR chemical shifts for nitrogen and oxygen [37][38][39][40] have been attributed to an insufficient separation between the occupied and virtual MOs, leading to excessively large paramagnetic contributions. We can take the HOMO-LUMO gap as a probe of the behavior of an approximate potential in the molecular region ͑as opposed to the asymptotic region͒, which will indicate whether the valence excitations in general will contribute more ͑smaller HOMO-LUMO gap͒ or less ͑larger gap͒ to the paramagnetic term. Table IV contains the HOMO-LUMO difference for the 44 systems studied, at the four levels of theory, VWN, BP86, SAOP and ͑when available͒ BP-GRAC-LB. There is indeed a correlation between the observed behavior of the HOMO-LUMO gap for the various potentials and the results for the shieldings. The SAOP potential gives the largest HOMO-LUMO gap in each system, thus avoiding the overestimation of the paramagnetic contribution, resulting in the most accurate chemical shifts, as previously observed. The HOMO-LUMO gaps for the BP-GRAC-LB potential are very close to the BP86 ones, which is consistent with the similarity in chemical shifts for these two potentials. It is to be noted that this is a rather global argument, which may mask considerable differences in the details for different molecules. For instance, we have shown earlier 17 that in the case of N 2 the lowest excitation shifts down by as much as 2.5 eV when one uses LDA or BP86 instead of SAOP or BP-GRAC-LB, which does not appear to be in agreement with the relatively small HOMO-LUMO gap differences in Table IV. However, this lowest excitation is not a HOMO to LUMO transition, so this large effect is not visible in Table  IV. It happens to be a Rydberg (3 g →3 u ) excitation. When this Rydberg excitation has little contribution to the chemical shift, as presumed earlier, the large shift of the excitation in the BP-GRAC-LB potential as compared to BP86 will have little effect. Discussion of individual cases at this level of detail is outside the scope of the present paper.
From the results above, we can conclude that the SAOP potential improves the chemical shift calculations, at least for the carbon, nitrogen, oxygen and fluorine atoms. For the same series of molecules, Patchovskii et al. 19 achieved a comparable, even slightly better, improvement over the LDA and GGA chemical shifts by means of the SIC-VWN ͑Ref. 19͒ and SIC-GGA ͑Ref. 23͒ potentials. The present SAOP potential has the advantage that it is much easier to implement, and can be used without further precautions in the SCF calculations, with a lower computational cost. We note that the HOMO-LUMO energy difference in the case of the SAOP potential is very similar to that obtained by SIC-VWN. Since this is an important factor in the paramagnetic term in particular, this provides some insight in the origin of the similar accuracy of both potentials. Of course, the asymptotic corrections applied in the SAOP and other asymptotically corrected potentials ͑HCTH-AC, BP-GRAC-LB͒ are a means of incorporating the proper asymptotic Ϫ1/r behavior that the exact xc potential owes to its selfinteraction correction part ͑the xc hole contains Ϫ1 elec-tron͒. However, the difference between BP-GRAC-LB on the one hand, and SAOP and SIC on the other hand, show that the asymptotic behavior is in this case not the important feature of the potential, but rather its shape in the molecular region that determines the HOMO-LUMO gap.
The approximate potentials like SAOP are supposed to approximate the exact Kohn-Sham potential, and their performance for a property should be judged against the performance obtained when the exact Kohn-Sham orbitals and orbital energies would be used. Wilson and Tozer 24 have recently generated such presumably rather accurate Kohn-Sham potentials from correlated densities obtained with ''Brueckner Doubles'' CI calculations. They denote these as multiplicative Kohn-Sham potentials, MKS͑BD͒. Table V contains the isotropic NMR shieldings with the MKS͑BD͒ potentials, together with the VWN, BP86, SAOP, and BP-GRAC-LB values, for the series of molecules studied by Wilson and Tozer. 24 Experimental values are included for comparison. We do not discuss the NMR shieldings obtained for individual molecules, but we focus on the mean ͑calculated minus experiment͒ error, and the mean absolute error, denoted as d and ͉d͉, respectively, which are both included in Table V. These errors are generally dominated by the errors of the O 3 molecule, that is why we also include the errors obtained when the O 3 results are omitted.
Since dϭϪ͉d͉ for all but the MKS͑BD͒ potential, we see that all these other potentials give systematically too low shielding values for this set of molecules. Qualitatively the errors when including the O 3 system are quite similar to those when omitting it, although quantitatively the errors are significantly lower in the latter case. We focus on the results when O 3 is omitted. Again, from the absolute mean errors, it is seen that the SAOP (͉d͉ϭ29.8) potential improves the VWN (͉d͉ϭ55.0) and BP86 (͉d͉ϭ46.1) ones. The BP-GRAC-LB results are again very close to those achieved by the BP86 potential, reinforcing the idea that Rydberg excitations do not affect the NMR calculations. The SAOP gives errors quite similar to those achieved by the hybrid B97-1 This MKS potential also proves to yield improved NMR shifts compared to conventional ͑hybrid͒ GGA functionals, and its density is easily obtained. However, the generation of the KS potential that reproduces a given density, either the correlated ͑BD͒ one or the B97-1 one, is not a routine procedure and in fact prone to numerical instabilities. It has not yet been performed for truly large molecules. It is therefore mandatory to develop potentials that are preferably simple functionals of the density, or or-bital densities, but at the same time approach the true KS potential closely enough so that the potential accuracy of the DFT calculations of chemical shieldings, which is apparent from the results of Wilson and Tozer, 24 can be achieved in a routine manner. The SAOP potential is obviously a step in the right direction.

IV. CONCLUSIONS
In this paper, we evaluate the performance of the SAOP potential for the calculation of NMR chemical shifts. This SAOP potential has been previously successfully tested for other response properties, like excitation energies, ionization potentials, and ͑hyper͒polarizabilities. [13][14][15][16][17] SAOP results show considerable improvement with respect to previous potentials, like VWN or BP86, at least for the carbon, nitrogen, oxygen, and fluorine chemical shifts. Also, a few NMR calculations carried out on third period atoms ͑S, P, and Cl͒ seem to improve when using the SAOP potential. Nevertheless, it is not an universal improvement, as some tests on Se systems did not result in the large improvement found for the first period atoms. Thus, future work will be focused on studying chemical shifts of heavier nuclei. The SAOP potential is not competitive with the nearexact KS potential, as can be seen from the results of Wilson and Tozer with the MKS͑BD͒ potentials. 24 It is more closely comparable to the SIC-LDA and SIC-GGA potentials studied by Patchovskii et al. 19,23 The SAOP potential offers the compromise of being computationally expedient while considerably improving over LDA and GGA potentials.

ACKNOWLEDGMENTS
One of the authors ͑J.P.͒ acknowledges the Departament d'Universitats, Recerca i Societat de la Informació ͑DURSI͒ of the Generalitat de Catalunya for benefiting from the doctoral fellowship No. 2000FI-00582.