On the convergence of zero-point vibrational corrections to nuclear shieldings and shielding anisotropies towards the complete basis set limit in water

ABSTRACT The method and basis set dependence of zero-point vibrational corrections (ZPVCs) to nuclear magnetic resonance shielding constants and anisotropies has been investigated using water as a test system. A systematic comparison has been made using the Hartree–Fock, second-order Møller–Plesset perturbation theory (MP2), coupled cluster singles and doubles (CCSD), coupled cluster singles and doubles with perturbative triples corrections (CCSD(T)) and Kohn–Sham density functional theory with the B3LYP exchange-correlation functional methods in combination with the second-order vibrational perturbation theory (VPT2) approach for the vibrational corrections. As basis sets, the correlation consistent basis sets cc-pVXZ, aug-cc-pVXZ, cc-pCVXZ and aug-cc-pCVXZ with X = D, T, Q, 5, 6 and the polarisation consistent basis sets aug-pc-n and aug-pcS-n with n = 1, 2, 3, 4 were employed. Our results show that basis set convergence of the vibrational corrections is not monotonic and that very large basis sets are needed before a reasonable extrapolation to the basis set limit can be performed. Furthermore, our results suggest that coupled cluster methods and a decent basis set are required before the error of the electronic structure approach is lower than the inherent error of the VPT2 approximation.


Introduction
Accurate prediction of nuclear magnetic parameters, including shieldings, shielding anisotropies and indirect spin-spin coupling constants is nowadays possible with correlated ab initio methods combined with sufficiently large and flexible basis sets [1][2][3]. However, for a realistic comparison of high-level theoretical results, typically obtained for isolated molecules, with experimental results rovibrational averaging of the calculated properties is needed [4,5]. Furthermore, since nuclear magnetic resonance (NMR) experiments are mostly performed in solution, solvent effects should also be considered [6,7], but this is outside the scope of the current article. CONTACT  Early studies of zero point vibrational corrections (ZPVCs) to magnetic parameters of small molecules were published about 30 years ago [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26]. Unfortunately, calculations of vibrational corrections, even in the simplest approximation -second-order vibrational perturbation theory (VPT2) [8], require the calculation of parameters defined as the third (cubic force-constants and shielding gradients) and fourth (shielding Hessians) derivatives of the electronic energy, and the number of such parameters, that must be calculated, increases rapidly with the number of atoms in the molecule. Thus, NMR calculations on medium size or larger molecules are usually performed on either an experimental or a theoretically optimised equilibrium geometry and without taking vibrational effects into account. It has even been suggested to compensate for the errors in B3LYP calculations of chemical shifts by performing these calculations at geometries optimised at the Hartree-Fock level. [27]. The studies that include vibrational effects for larger molecules often limit their investigation to the Hartree-Fock (HF) or density functional theory (DFT) level of theory and employ smaller basis sets. Ruud,Åstrand and Taylor reported e.g. calculations of ZPV corrections to C and H shieldings in 38 molecules at HF level of theory. They found for the H shieldings a degree of transferability of these corrections between similar molecular fragments in contrast to the C shieldings [28].
Since H 2 O is a small molecule, it has previously been studied in relation to the size of the vibrational effects on nuclear shieldings [8,9,19,21,24,26,[29][30][31] and many times for the equilibrium geometry, e.g. in [32]. A pioneering study by Fowler and Raynes [9] used an empirical force-field with HF shielding surfaces to obtain a ZPVC to σ ( 17 O) of −13.1 ppm, whereas a later but similar study by Wigglesworth et al. [24] using instead MCSCF for the shielding surface obtained −9.9 ppm. Other studies mostly obtained values within this range. HF calculations by Fukui et al. [19] and MCSCF calculations by Vaara et al. [21] expanded their force field in moleculedependent internal coordinates. Later studies, such as those of Ruud et al. [26], Auer [29] and Kupka et al. [30], tend to favour more general approaches based on normal coordinate expansions of the surfaces. The most thorough calculation of the nuclear shielding of water was presented by Puzzarini et al. [31] as part of an attempt to redefine the absolute shielding scale of 17 O. They used CCSD(T) with a large basis set and an accurate numerical description of the vibrational problem to compute a ZPVC of −11.7 ppm, which is 0.8 ppm lower than their VPT2 value. A conservative estimate of the error associated with the VPT2 approach is, therefore, around 10 %. The error associated with the method and basis set used for the calculation of the ZPVC should consequently not be larger than that. Whereas originally specialised programs were often used for the calculation of ZPV corrections to magnetic properties, automatic routines are nowadays available in both 'free for academic use' licenced programs such as Dalton [33] or CFOUR [34,35] and commercial programs such as Gaussian [36].
The magnitude of absolute magnetic parameters depends heavily on the quality of the one-electron basis set. In case of the electronic energy and atomic/molecular parameters directly related to the energy, their values change irregularly when changing the basis set size in the traditional Pople type basis sets. However, other hierarchies of basis sets have been constructed so that convergence of the electronic energy towards the complete basis set (CBS) limit is observed when performing calculations with increasingly large basis sets. The most popular such series is the correlation-consistent basis sets of Dunning [37][38][39] which was designed with correlated calculations in mind, whereas a more recent addition is the polarisation-consistent basis sets of Jensen [40,41] optimised for DFT calculations. Later it was shown that other magnetic properties often display similar patterns of convergence [42][43][44][45][46]. Specialised basis sets optimised for the calculation of magnetic properties were also developed based on both the correlation and polarisation consistent basis set [43,44,[47][48][49][50] as e.g. the (aug)-pcS-n series [44], which is the shielding constant version of the polarisation consistent basis sets. On the other hand, due to the high intrinsic cost of calculating vibrational corrections, these are usually calculated using fairly small basis sets and there are only few studies in the available literature thoroughly investigating the basis set dependence of ZPVCs to magnetic properties [31,51,52]. The question remains, therefore, whether also the ZPVCs display a regular convergence and if it is possible to extrapolate these to the CBS limiting values at HF, DFT, secondorder Møller-Plesset perturbation theory (MP2) or coupled cluster (CC) levels of theory. Furthermore, is it possible to approximate large-basis set, high-level CC ZPVC calculations by lower level calculations? In the present study, we have investigated these question for the simple system, water, carrying out calculations of ZPVCs at the HF, MP2, CCSD, CCSD(T) and DFT using the B3LYP exchange-correlation functional levels with the cc-pVXZ, aug-cc-pVXZ, cc-pCVXZ, aug-cc-pCVXZ (X = D, T, Q, 5, 6) and aug-pc-n and aug-pcS-n (n = 1, 2, 3, 4) series of basis sets.

Theory
In NMR experiments, molecules are exposed to an external magnetic field B, which will interact with the nuclear magnetic moments of the atomic nuclei m N . However, due to the surrounding electrons, the magnetic field experienced by each nuclei is slightly changed. The leading term in this effect is the nuclear shielding tensor σ N with elements defined as [1,53,54] where E e is the electronic energy. When dealing with perturbations due to an external magnetic field, like in the case of the nuclear shielding tensor in Equation (1), it is necessary to specify a gauge origin for the external magnetic field. All practical quantum chemical methods project the one-electron Schrödinger equation into a finite basis of 'atomic' orbitals, which creates a spurious dependence on this gauge origin. Today, the standard solution is to use gauge-including atomic orbitals (GIAOs) also known as London orbitals [55][56][57], which has been implemented for most common analytical derivative techniques [58][59][60][61][62], although alternative ways for dealing with the gauge-origin problem are available [47,[63][64][65]. The entire shielding tensor cannot be determined experimentally. In typical liquid or gas phase experiments, only the isotropic nuclear shielding can be obtained An often reported quantity is also the shielding anisotropy, which is defined from the principal elements of the symmetrised shielding tensor (σ N 11 ≤ σ N 22 ≤ σ N 33 ) as [66] The theory of perturbation theory-based approximations to the vibrational corrections is well established [8], so we will present only a quick summary, more thorough accounts can be found in the reviews of Ruden and Ruud [4] and Faber et. al. [5]. The normal mode description of molecular vibrations is obtained by diagonalising the second-order term in the Taylor expansion of the electronic energy as a function of mass-weighted atomic positions. Perturbation theory can then be applied by treating the higher order terms in this expansion as corrections to the harmonic vibrational Hamiltonian typically by taking the nth term of the expansion as the (n − 2)th order correction. This allows for the addition of anharmonic vibrational effects to the vibrational energy levels in the VPT2 approach and for the calculation of vibrational expectation values of a molecular property provided that the functional dependence of the property on the nuclear coordinates is known. This dependence can again be expressed in form of a Taylor expansion in normal coordinates. In order to ensure a balance between the potential energy and property surface, one can employ a double perturbation approach and include orders also in the expansion of the property surface, or one can simply truncate the Taylor expansion of the property surface at second order. Either way, the lowest order approximation to the vibrational correction in the lowest vibrational state is where {q i } are the dimensionsless reduced normal coordinates and k ijj is the third derivative of the potential energy surface w.r.t. these reduced normal coordinates. Basis set extrapolations are based on the assumption, that the given property converges monotonously towards the correct result, if a hierarchy of basis sets has been constructed in a systematic manner. It should then be possible to extrapolate to the result obtained in the CBS limit. Various functional forms have been suggested for this extrapolation. Studies on helium, for example, suggest that an inverse-third-power expansion is appropriate for the correlation energy [67,68], whereas an exponential convergence is usually assumed to hold for HF and DFT results [69,70]. It has also been suggested to fit the HF and correlation contribution to a property separately using an exponential and inverse-third-power expansion, respectively [70]. However, it is not clear how such a combined approach could be employed in the context of vibrational averaging using correlated wave function methods, as this involves combined basis set effects from the geometry optimisations, force constant and property calculations. Since properties such as magnetic shieldings can be described as derivatives of the electronic energy, it is often assumed that the same procedures as for energies work also for properties. Indeed, it has been demonstrated previously that this is often the case [46,71]. Still one might run into problems, when different contributions to a particular property exhibit distinct basis set dependencies preventing an overall monotonic convergence. Vibrational corrections, which have contributions from force constants and property derivatives calculated at optimised geometries (OptGs) might be such a property.

Computational details
Gaussian 09, ver. D.01 [72] has been used to perform the RHF and B3LYP calculations, whereas the MP2, CCSD and CCSD(T) calculations have been carried out with the CFOUR program [34]. GIAOs have been employed in all calculations of the nuclear shieldings. Calculations at the experimental equilibrium geometry use the parameters given in [73]: R OH = 95.72pm, θ HOH = 104.52°. The zero-point vibrational averaging in Gaussian has been performed using Barone's recent implementation of VPT2 methodology [36]. Very tight optimisation criteria as well as a very large density grid for DFT were used in the calculations. The core-valence correlation-consistent basis sets cc-pCVXZ and aug-cc-pCVXZ, which should Table . The isotropic nuclear shielding of  O, σ (  O), in ppm calculated using the aug-cc-pVZ, aug-cc-pCVZ (aug-cc-pCVZ for CCSD) and aug-pcS- basis sets and various methods. ExpG and OptG are results at the experimental and optimised geometry, while ZPVA and ZPVC are the zero-point vibrational averaged values and corrections. ExpC is the difference between the ZPVA and ExpG value. CBS-D or CBS-J are extrapolated values using the inverse third power scheme, (), and the results for the two largest basis sets in the aug-cc-p(C)VXZ and aug-pcS-n series.
give a better description of core-valence correlation, were therefore only employed in the correlated wave function calculations. The CFOUR calculations using the polarisation consistent and valence Dunning-type basis sets were performed using the frozen-core approximation, whereas calculations using the core-valence versions include the core in the correlation treatment. Since the calculation of vibrational corrections inevitably involves numerical differentiation, tight convergence criteria were selected with the exception of the largest basis sets for which the threshold in the SCF calculation was 10 −8 and for the CC amplitudes 10 −7 . For the numerical differentiation, the default step length of 0.05 (dimensionless) reduced normal coordinates were used.
The CBS values were estimated from two-parameter fits to Equation (5) or three-parameter fits to Equation (6) using the results, which were obtained with the largest basis sets, if not specified otherwise.

Method dependence
In Table 1, we compare the 17 O isotropic absolute shielding constant, σ ( 17 O), obtained with various methods and Table . The isotropic nuclear shielding of  H, σ (  H), in ppm calculated using the aug-cc-pVZ aug-cc-pCVZ (aug-cc-pCVZ for CCSD) and aug-pcS- basis sets and various methods. ExpG and OptG are results at the experimental and optimised geometry, while ZPVA and ZPVC are the zero-point vibrational averaged values and corrections. ExpC is the difference between the ZPVA and ExpG value. CBS-D or CBS-J are extrapolated values using the inverse third power scheme, (), and the results for the two largest basis sets in the aug-cc-p(C)VXZ and aug-pcS-n series.

Basis
ExpG OptG ZPVA ZPVC ExpC the largest basis sets employed in this study, the augcc-pCV6Z (aug-cc-pCV5Z for CCSD), aug-cc-pV6Z and aug-pcS-4 basis sets. The calculations at the experimental geometry (ExpG) show that the HF, B3LYP and MP2 results differ by 8-11 ppm or about 3 % from the CC results. Both HF and B3LYP underestimate σ ( 17 O), while MP2 overestimates it implying that MP2 overestimates the correlation correction by almost a factor of 2. When using instead a geometry, which was calculated with the same method as the shieldings, the HF results change by as much as 9 ppm, whereas for the other methods, which are generally considered to be better geometries, the change is less dramatic. Even more remarkable is that the HF result for the oxygen shielding is coincidentally within 1 ppm of the CCSD(T) value both calculated at their respective equilibrium geometries and thus in much better agreement with the CC results than MP2 or B3LYP.
The largest basis results of the vibrational corrections relative to the optimised geometries (ZPVC) calculated with B3LYP and MP2 are off by about 2 ppm as compared to the CCSD(T) value, while the HF value differs by only 0.5 ppm. Furthermore, since the B3LYP ZPVC is 2 ppm more negative than the CCSD(T) value and the B3LYP optimised equilibrium geometry (OptG) value is already too low, the difference between the two methods is further Table . The nuclear shielding anisotropy of  O, σ (  O), in ppm calculated using the aug-cc-pVZ aug-cc-pCVZ (aug-cc-pCVZ for CCSD) and aug-pcS- basis sets and various methods. ExpG and OptG are results at the experimental and optimised geometry, while ZPVA and ZPVC are the zero-point vibrational averaged values and corrections. ExpC is the difference between the ZPVA and ExpG value. CBS-D or CBS-J are extrapolated values using the inverse third power scheme, (), and the results for the two largest basis sets in the aug-cc-p(C)VXZ and aug-pcS-n series.  Table 1, the ZPVC relative to the OptGs are essentially the same for all methods, −0.5 ppm. Also, the 'ordering' of the methods differs from the one observed for oxygen. For hydrogen at both the experimental and OptGs and in the vibrationally averaged values, B3LYP overestimates the CCSD(T) shielding, while MP2 underestimates it. The HF results again depend very much on the geometry such that, at the ExpG they underestimate the CCSD(T) value, while at the OptG they overestimate it. . The nuclear shielding anisotropy of  H, σ (  H), in ppm calculated using the aug-cc-pVZ aug-cc-pCVZ (aug-cc-pCVZ for CCSD) and aug-pcS- basis sets and various methods. ExpG and OptG are results at the experimental and optimised geometry, while ZPVA and ZPVC are the zero-point vibrational averaged values and corrections. ExpC is the difference between the ZPVA and ExpG value. CBS-D or CBS-J are extrapolated values using the inverse third power scheme, (), and the results for the two largest basis sets in the aug-cc-p(C)VXZ and aug-pcS-n series.  Tables 3 and 4 show the corresponding results for the shielding anisotropies σ ( 17 O) and σ ( 1 H), respectively. For the 17 O equilibrium geometry anisotropies, the trends of the methods are mostly the same as for the isotropic shielding, although the relative discrepancies between the methods are different. The differences in the oxygen shielding anisotropy at the experimental geometry (ExpG) of the cheaper methods from the CCSD(T) results are with 4-8 ppm slightly smaller than for the isotropic shielding but percentage wise larger with 8%-15 %. However, the 'ordering' of the methods is different. HF gives now the largest results and MP2 the smallest. This implies that MP2 overestimates the correlation correction, like for the isotropic shieldings, and that B3LYP underestimates it. Optimising the geometry does not change this significantly. Neither does vibrational averaging, as the correlated methods agree approximately on a ZPVC of −2 ppm with the B3LYP value (−2.27 ppm) being closer to the −2.17 ppm of CCSD(T) than MP2. Again like for σ ( 17 O), adding the ZPV correction to the MP2 OptG value increases the discrepancy to the corresponding CCSD(T) value. But contrary to the isotropic shielding, the HF ZPVC to the anisotropy is less than half of the CCSD(T) value and thus also increases the difference to CCSD(T). The differences between the methods in σ ( 1 H) are rather small. At both the experimental and optimzed geometries, B3LYP underestimates the CCSD(T) value by 1 ppm, while MP2 overestimates the it by 0.5 ppm.
The HF values, on the other hand, depend again strongly on the geometry. While at experimental geometry HF only slightly overestimates the CCSD(T) value, the difference is a factor of 4 larger at OptG. Similar to σ ( 1 H), all methods agree on a ZPV correction to σ ( 1 H) of −1 ppm, which is, however, twice as large as for the isotropic shielding.

Trends in basis set convergence
A meaningful basis set extrapolation requires a monotonic convergence of the results to the basis set limit. In order to asses whether the calculated nuclear shieldings actually follow such a trend, plots of the basis set progression of the isotropic shieldings of each of the methods are given for a selection of basis sets in Figures 1-7. All results can also be found in Tables S1-S10 of the supplementary material. Figures 1-3, we see that one cannot expect a monotonic convergence in all cases. If one follows the trend for basis sets without diffuse functions, cc-pCVXZ ( in Figures 1-3 Figures 4-7), one observes a monotonic convergence from above for the two shielding constants and from below for the two anisotropies. However, the convergence is very slow and the progression appears to be almost linear between the triple and quintuple zeta basis sets, only showing signs of convergence beyond the quintuple zeta basis sets. Adding diffuse functions in the augcc-p(C)VXZ basis sets ( or • in Figures 1-3 Figures 4-7), there appears a difference between the results for σ ( 17 O) at the experimental equilibrium geometry on one side and the OptG and vibrationally averaged results on the other. Overall, the results converge quickly, but on closer inspection one can see that only at the experimental equilibrium geometry the results converge monotonically from above. Nevertheless, for the aug-cc-pVXZ series, the difference between the QZ and 5Z ExpG results (0.54 ppm) is still   larger than difference between the TZ and QZ results (0.13 ppm), while for the core-valence basis sets aug-cc-pCVXZ, the differences are 0.38 and 0.27 ppm, respectively, suggesting that core-valence results converge more evenly. On the other hand, at the OptG and when adding the ZPV corrections, a dip shows clearly up in the curve at the TZ level before converging now from below in the case of the aug-cc-pCVXZ series while the aug-cc-pVXZ results continue to oscillate up to the sextuple zeta level. For σ ( 17 O), the results with the augmented correlation consistent basis sets converge also much faster and monotonically but there is no significant difference in the convergence between the two geometries or the vibrational averaged results apart from that the core-valence basis set gives consistently larger values of σ ( 17 O). In case of σ ( 1 H), the diffuse functions are of much less importance. A lot of the improvement on including diffuse functions is due to an improved description of the geometry, since the change on adding them is much smaller for the calculations at the experimental geometry, Figure 1, than The hydrogen results obtained with the polarisation consistent basis sets aug-pc-n ( in Figures 1-3) do not differ significantly from the corresponding aug-ccp(C)VXZ results. However, the aug-pc-n results exhibit a clear minimum for n = 2 before converging from below for σ ( 17 O), while for σ ( 17 O) they go through a maximum for n = 3 and are not converged yet at n = 4.

or filled symbols, dashed lines in
Finally, the specialised shielding basis sets aug-pcS-n (× in Figures 1-3 Figures 4-7) exhibit a completely different behaviour for the oxygen isotropic and anisotropic and the hydrogen isotropic shieldings than the other basis sets. For the oxygen and hydrogen shieldings, they exhibit an oscillatory behaviour, while the shielding anisotropies converge monotonically. The aug-pcS-n results for the hydrogen shielding and for both anisotropies converge faster towards the largest basis set result than the corresponding aug-cc-pCVXZ results. However, this is again not the case for σ ( 17 O). At the experimental equilibrium geometry, aug-pcS-1 is still closer to the largest basis set result than aug-cc-pCVDZ, but this does not hold for the larger basis sets or at the OptG and the vibrationally averaged values. Actually, the largest aug-pcS-n result is for σ ( 17 O) closer to the largest result from the aug-cc-pVXZ series than from the corevalence series. Furthermore, for σ ( 17 O) aug-pcS-n does not converge faster than aug-cc-p(C)VXZ.

or solid lines in
To understand the reasons for the sometimes uneven convergence at the optimised geometry and in the vibrational averaged shieldings, one has to remember first of all that not only shielding constants enter in the ZPVA values, but also force constants and first and second derivatives of the shielding constants, which might exhibit quite a different basis set dependence. Second, the basis sets were also changed in the optimisation of the geometries, which implies that we are looking at the combined effect of basis set convergence of shieldings, shielding derivatives and force constants and their variation with a changing geometry. Actually, it is, therefore, surprising that one can observe a monotonic convergence in some cases. On the other hand, in the case of σ ( 17 O) and the aug-pcSn basis sets, it is primarily a basis set effect, as the calculations at the experimental equilibrium geometry lack already a monotonic convergence.
The trends are very similar for the other methods (Figures 4-7), but one should note two points.
Whereas for the correlated methods, the calculations at the experimental (symbol •) and optimised (symbol ) geometries converge towards almost the same value of the shielding, that is not really the case for B3LYP, and the limiting values for HF at the two geometries differ vastly as already discussed in the previous section. Furthermore in the HF case, the convergence towards the basis set limit is not really faster than in the correlated calculations, as one might have expected. For our present purpose, it is also interesting to look at the convergence of the ZPVCs separately. In Figure 8, the basis set convergence of the ZPVC to σ ( 17 O), σ ( 1 H), σ ( 17 O) and σ ( 1 H) at the CCSD(T) level is plotted. None of the basis sets exhibits a monotonic convergence for the ZPVC to all four properties. Four out of the five basis set series agree on the final value of the ZPV correction to σ ( 17 O), i.e. the value obtained with the largest basis set in the series, whereas aug-pc-n predicts a slightly smaller absolute value. However, only the cc-pCVXZ basis sets exhibit a monotonic but very slow convergent behaviour, while the larger basis sets give results which either oscillate or do not show a convergent behaviour for the three largest basis sets in the series. For σ ( 17 O), the cc-pCVXZ, aug-cc-pVXZ and aug-pcS-n series agree on the final values, while both the aug-pc-n and aug-cc-pCVXZ series predict numerically smaller values. The cc-pCVXZ and aug-pcS-n basis sets exhibit also a montonic convergent behaviour in contrast to the other basis sets, whose results show some oscillations. For hydrogen, the variation in the ZPVC of both the isotropic shielding and the shielding anisotropy due to the basis set is an order of magnitude smaller than for oxygen, although the shielding anisotropy is only a factor 2-3 smaller. The changes are thus rather irrelevant. Nevertheless, only some of the basis set series exhibit a monotonic convergent behaviour. For σ ( 1 H), the aug-pcS-n and aug-cc-pVXZ series of basis sets converge monotonically, the first from below and the second from above in terms of the absolute values, while for σ ( 1 H) none of the basis set series shows a trend which could easily be extrapolated to the CBS limit.

Basis set extrapolation
Having sketched out the characteristics of the basis set progression, we can now attempt to extrapolate towards the basis set limit using both the inverse third power formula in Equation (5) and the exponential formula in Equation (6). The question, we are mostly interested in, is whether one can approximate the results obtained with the largest basis sets in the aug-cc-p(C)VXZ and aug-pcS-n series by extrapolating results from smaller basis sets. We will discuss this explicitly for the isotropic shieldings at the OptG and the ZPVA value both at the HF and CCSD(T) level as these are the two extremes in our series of methods. The results of the CBS extrapolation for all methods and properties using Equation (5) are given in Tables S11-S20 of the supplementary material and for the two largest basis sets, they are shown in Tables 1-4 where they have been discussed previously. Results for the exponential extrapolation, Equation (6), at the HF and CCSD(T) level are given in Tables S21 and S22 in the supplementary material.
In Figures 9 and 10, we have, therefore, compared for σ ( 17 O) and σ ( 1 H) the values obtained with a basis set with cardinal number X/n (× for aug-pcS-n and + for aug-cc-p(C)VXZ) with estimates of the CBS value obtained by extrapolating the results of this basis and the next smaller one(s) (open symbols). For σ ( 17 O), there is no significant difference between the results obtained at the HF and CCSD(T) level or with the aug-cc-p(C)VXZ and aug-pc-n series of basis sets. We observe that the estimates of the CBS values from Equation (5) using two data points exhibit not surprisingly oscillations like the original data. The estimated CBS values are often further away from the results of the largest basis set than the original smaller basis set values. Essentially only fits including the largest basis sets predict CBS values that are close to the largest basis set results. The situation is somewhat different for the exponential fits, Equation (6), where the estimated CBS values are close to the results of the largest basis set employed in this particular fit. Nevertheless, they are still not closer to the result of the largest basis in the series. The conclusion is thus that it is unfortunately not possible for σ ( 17 O) to estimate the CBS value from calculations with the smaller basis sets DZ-QZ or n= 1-3. In the case of hydrogen and the aug-pcS-n series of basis sets, the situation is essentially the same as for oxygen. The original data values oscillate and thus also the estimated CBS values. The variation is, however, much smaller than for oxygen. On the other hand, for the augcc-p(C)VXZ basis set series the situation is quite different. The calculated shieldings exhibit a monotonic convergence, which holds then also for the estimated CBS values. The latter are, therefore, also often closer to the results of the largest basis set than the original smaller basis set values. The differences between the two point inverse-third-power fits and the three points exponential fits are also smaller than for oxygen. For σ ( 17 H) in combination with the aug-cc-p(C)VXZ basis set series, it is thus possible to improve the results of smaller basis sets by extrapolation. For all the other properties, however, one still can extrapolate the results obtained with the largest basis sets towards a CBS value. This has been done for both the aug-cc-p(C)VXZ and aug-pcS-n series and the results are shown in Tables 1-4.

Conclusion
We have investigated the method and basis set dependence of the shielding constants and anisotropies in water and their ZPVCs. For the total shieldings and anisotropies calculated with the largest basis sets either at equilibrium geometries or the vibrationally averaged values, we confirm the well known trend, that MP2 overestimates the correlation effects for these properties and DFT with the B3LYP exchange correlation functional underestimates them. The HF results behave rather erratically compared to the CCSD(T) results leading e.g. to close agreement in case of the ZPV averaged isotropic shielding of oxygen.
With respect to the basis set dependence of both the equilibrium geometry and the vibrationally averaged values independent of the method employed, we observe rather slow convergence for the non-augmented basis sets so that sufficient convergence is not achieved before the sextuple zeta level. Adding the diffuse functions significantly accelerates convergence apart for the anisotropy of the hydrogen shielding, and sufficient convergence is already achieved at the quintuple zeta level or for n = 3 in the case of the aug-pc(S)-n basis sets. While the augcc-p(C)VXZ results for the hydrogen shielding and all series of basis sets for the anisotropies exhibit a reasonable monotonic behavior which allows for extrapolations, the oxygen shielding results approach the value of the largest basis in a way, which prevents a meaningful extrapolation to the CBS limit with the exception of the results of the two or three largest basis sets.
Also for the ZPVCs, one observes significant differences on adding diffuse functions. Furthermore, although the values converge with increasing basis set, the convergence is not monotonic in most cases, which means that one cannot accurately extrapolate to the CBS limit from calculations with smaller basis sets. The variation in the ZPVC is, on the other hand, smaller than the variation in the equilibrium geometry values and smaller than the inherent error in the VPT2 approach. An augmented basis set of TZ or n = 2 quality will, therefore, suffice.
With respect to the method dependence of the ZPVCs, we observe that for the hydrogen shielding and shielding anisotropy basically any method is good enough, while for oxygen there are significant differences. While HF predicts a ZPV correction to the oxygen shielding close to the CCSD(T) result, the B3LYP and MP2 ZPVCs are disappointing. They differ by about 20 % from the CCSD(T) result, which we believe to be larger than the intrinsic error of using the VPT2 approach for calculation of vibrational corrections. Still this might be acceptable in many cases, as ZPVCs are rarely more than 10 % of the equilibrium value and this error may, therefore, not exceed the correlation error of the equilibrium calculation, even if that is carried out with a method such as CCSD(T).
The main conclusion of our study is that one should not expect a priori that extrapolations to the CBS values of magnetic properties and in particular of their vibrational corrections are meaningful. Oscillatory convergence can occur, which inhibits extrapolations, since the commonly used forms of fit functions do not allow for such behaviour. Also the composite nature of VPT2 calculations, where both geometry, force-fields and shielding-derivatives may not exhibit the same converges behaviour, adds to the challenge. When going all the way to the largest basis sets in both the aug-cc-p(C)VXZ and aug-pcs-n series of basis sets, one can obtain CBS values of the ZPVC that one can be confident in. However, we were also able to show, at least for oxygen, that attempts to use extrapolation techniques with basis sets of intermediate sizes did generally not improve the results compared to the value of the largest basis set. While functional forms for the extrapolation functions, able to handle these challenges, could in principle be developed, these will include more than the two or three fitting parameters common in present day schemes, which makes them unsuitable for actual applications due to the rapidly increasing computational cost of ZPVC calculations as the basis set size is increased.
Although we have illustrated our points only for a simple single molecule, water, we believe our main conclusions to be transferable not only to other simple hydrides but also to larger molecules. The rapidly increasing number of vibrational modes and therefore increasing number of contributions to the vibrational corrections should make it even more unlikely that the vibrational corrections exhibit a monotonic convergence towards the CBS limit in larger molecules.