Errors in the Calculation of 27Al Nuclear Magnetic Resonance Chemical Shifts

Computational chemistry is an important tool for signal assignment of 27Al nuclear magnetic resonance spectra in order to elucidate the species of aluminum(III) in aqueous solutions. The accuracy of the popular theoretical models for computing the 27Al chemical shifts was evaluated by comparing the calculated and experimental chemical shifts in more than one hundred aluminum(III) complexes. In order to differentiate the error due to the chemical shielding tensor calculation from that due to the inadequacy of the molecular geometry prediction, single-crystal X-ray diffraction determined structures were used to build the isolated molecule models for calculating the chemical shifts. The results were compared with those obtained using the calculated geometries at the B3LYP/6-31G(d) level. The isotropic chemical shielding constants computed at different levels have strong linear correlations even though the absolute values differ in tens of ppm. The root-mean-square difference between the experimental chemical shifts and the calculated values is approximately 5 ppm for the calculations based on the X-ray structures, but more than 10 ppm for the calculations based on the computed geometries. The result indicates that the popular theoretical models are adequate in calculating the chemical shifts while an accurate molecular geometry is more critical.


Introduction
Studies with aluminum(III) (Al(III)) complexes are important in bioinorganic chemistry, geochemistry, environmental science, and material science due to the toxic effects of many aluminum compounds and their industrial usage as catalysts and coagulation agents. Because various kinds of species of Al(III) might differ significantly in their toxic effects, speciation of aqueous Al(III) complexes is of great value to study their biological and environmental roles [1,2]. Nuclear magnetic resonance spectroscopy (NMR) has been routinely used in studying the coordination chemistry and speciation of Al(III) in aqueous solutions due to its noninvasive character [3]. 1 H, 13 C, 17 O and 27 Al are commonly used nuclei for observing NMR spectra of Al(III) complexes.
Among them, 27 Al possesses advantages over the others, including 100% abundance, high signal sensitivity, a wide range of chemical shifts, and sensitivity to coordination environment. For example, the 27 Al chemical shift (δ Al ) of the tetrahedral [Al(OH) 4 ] − species is 80 ppm relative to the octahedral [Al(OH 2 ) 6 ] 3+ species. Therefore, solution-state 27 Al NMR spectroscopy has been a very powerful tool to characterize the structures of Al(III) complexes, to monitor the hydrolysis of Al(III), to identify Al(III) species in biological and environmental samples, and to quantify the concentrations of Al(III) species [3]. It is particularly useful for the characterization of those species which are difficult to isolate and/or to crystallize or inert for other spectroscopic methods. 27 Al nucleus has a spin of 5/2 and with a moderately large electric quadrupole moment. Both the chemical shift and line width of 27 Al signals are sensitive to aluminum coordination geometry. For a highly asymmetric coordination geometry, the resonance peaks of 27 Al may be too broad to observe. In aqueous solutions, the coordination number of Al 3+ monomers fluctuates between 4 and 6. Many aqueous Al(III) complexes are octahedral-coordinated and some may be distorted. As a result, δ Al of these complexes are not far away from 0 ppm relative to the aqueous [Al(OH 2 ) 6 ] 3+ species and the chemical shift spectra peaks are not sharp. Furthermore, due to the tendency of Al(III) to hydroxylize, to oligomerize and to form mixed complexes, the solutions are often a complicated system with many coexisting species. Often, many 27 Al chemical shift signals are crowded together in the vicinity of 0 ppm. Although tentative structural information can be made by comparing the spectroscopic features of the unknown species with model complexes, it is usually difficult to make unique assignments. Therefore, computational modeling, along other experimental methods, is often used to complement 27 Al NMR experiments.
With the rapid development of computational chemistry, its application has become almost ubiquitous in all branches of chemistry. In the study of Al(III) complexes, computational methods have been used to predict NMR spectra. For example, Tossell [4] calculated the chemical shieldings of 27 Al, along with other nuclei, by applying the gauge-including atomic orbital (GIAO) method [5,6] with the Hartree-Fock (HF) theory and to elucidate the species formed during the hydrolysis and the oligomerization process of Al(III) in aqueous solution. Amini et al. [7,8] applied the HF theory and density functional theory (DFT) to calculate the solid-state 27 Al NMR spectra of aluminum acetylacetonate and other similar complexes for the clarification of the crystal structure of the compounds. Mirzaei et al. [9,10] used DFT calculations to calculate 27 Al and 31 P chemical shieldings in the investigation of aluminum phosphide and nitride nanotubes. The gauge including projected augmented wave (GIPAW) method [11] developed by Mauri and Pickard in 2001 [12] which enabled the calculation of all-electron NMR parameters in solids was also used to predict the solid-state 27 Al NMR spectra [13,14].
Although the calculations of chemical shieldings have played an important role in the determination of structures, there have been a limited number of systematic studies concerned with the accuracy of computational models that predict the chemical shifts of 27 Al. Mulder et al. [15] reviewed several theoretical methods for the calculation of NMR chemical shifts and they evaluated the applications of these methods in protein structure determination. Jensen [16] evaluated the basis set convergence problem in calculating NMR shielding constants calculated by two DFT methods and proposed to use pcS-n basis sets for calculating shielding constants. In 1999, Kubicki et al. [17] performed a thorough computational study on Al(III) hydrolysis products and Al(III)-carboxylate complexes to determine how accurately the 27 Al chemical shifts can be calculated using ab initio methods. The study was based on the self-consistent isodensity polarized continuum model (SCIPCM) of the aqueous species using HF and second-order Møller-Plesset perturbation (MP2) theories. By comparing the calculated chemical shifts with the experimental values, they suggested that the 27 Al NMR spectra should be reinterpreted in some previous publications. Ten years later, Bi et al. [18] assessed the accuracy of 27 Al chemical shift calculations that used 29 DFT functionals, HF and MP2 models. They concluded that of all the models, the HF and MP2 models gave the most accurate results. Their conclusion was made on the basis of the calculations of only a few Al(III) hydrolysis products. It remains unanswered if the conclusion holds for more general Al(III) complexes.
Previous studies found that the error involved in the ab initio calculation of 27 Al chemical shifts is typically approximately 10 ppm and for some DFT models it could be as large as 30 ppm [18]. Errors of such magnitude pose serious limitations for the application in assigning experimental spectra, since the 27 Al chemical shifts of many octahedral-coordinated complexes are typically distributed in the small range from −10 to 30 ppm. To distinguish between the small differences found in similar complexation species, a computational method must be able to reproduce the chemical shifts with an error within a few ppm and at least be able to put the resonance peaks in the correct order.
There are two major parts to the errors in calculated 27 Al chemical shifts: the error in the prediction of the geometry and the error in the chemical shielding determination. Previous studies presenting the calculations of 27 Al chemical shifts did not differentiate between these two kinds of errors. There have been many publications concerning the accuracy of ab initio models including those that use DFT in predicting complex geometries [19]. The purpose of this work is to evaluate the performance of several commonly used theoretical models, including HF, MP2 and DFT, in predicting the chemical shifts of Al(III) complexes. Our approach is to use the X-ray crystallographic structures to calculate the chemical shifts and then to compare them with the experimentally measured values. It has been well established that X-ray diffraction structures represent well the molecular geometries in solutions except for the rotation-flexible functional groups.

Results and Discussion
More than 200 single-crystal X-ray diffraction determined geometries  of 114 Al(III) complex species were obtained from the Cambridge Crystal Database. (There is more than one crystal structure available for some complexes.) In most species, Al(III) is found in an octahedral-coordination environment binding with six ligands or functional groups. Some Al(III) are found in tetrahedral coordination environments and a limited number of species are in penta-coordination environments. Usually Al(III) binds with ligands via O atoms. Al-N ranks second in the binding modes. A few species of Al(III) coordinated with fluorides and phosphates were also examined. This covers most coordination forms commonly encountered in the study of solution aluminum chemistry. Other coordination forms, such as planar tri-coordination and complexation with other ligands, e.g., H and Cl, were not studied in this work.

Chemical Shielding Constants
The chemical shielding tensors of each structure were calculated with the GIAO formalism [5,6] at the following four popular model levels, B3LYP/6-31G(d), HF/6-31G(d), B3LYP/6-311+G(d,p) and HF/6-311+G(d,p). In addition, MP2/6-311+G(d,p) was also applied to species of [Al(OH 2 ) 6 ] 3+ , [AlF 6 ] 3− and [AlF 5 (OH 2 )] 2− . Figure 1 shows how the calculated isotropic shielding constants at different levels are correlated. The x-axis represents the values calculated at the B3LYP/6-311+G(d,p) level. On the y-axis direction, the shielding constants obtained at the B3LYP/6-31G(d), and HF/6-311+G(d,p) levels are represented as black squares and red circles, respectively. There is a strong linear correlation in a wide range from 350 to 600 ppm for both data sets, though the chemical shieldings calculated at different levels of theory could differ as large as 60 ppm. In general, the hybrid B3LYP functional gives consistently smaller shielding constants than the HF theory and so does the larger basis set (6-311+G(d,p)) than the smaller (6-31G(d)) basis set.
Comparing the two differently sized basis sets, 6-31G(d) versus 6-311+G(d,p), both using the B3LYP model, the coefficient of determination (R 2 ) of the linear regression is 0.95. The best-fit straight line is Equation 1 (black dashed line in Figure 1). The majority of the absolute fitting residuals (74%) are less than 5 ppm. There is one outlier point with a fitting residual of 24 ppm, while all the rest are less than 13 ppm. The outlier point is an Al-phosphate species. The exact reason for the outlier is unclear.
The shielding constants obtained at the HF/6-311+G(d,p) level are slightly higher than those at the B3LYP/6-31G(d) level. They also have strong linear correlations with the values obtained at the B3LYP/6-311+G(d,p) level (R 2 = 0.98). The best-fit straight line is Equation 2 (red dashed line in Figure 1). Most of the absolute fitting residuals (85%) are less than 5 ppm, while the rest are within 11 ppm except the Al-phosphate species which has a large residual of 24 ppm. For the calculations using the 6-31G(d) basis set, the correlation between the B3LYP and HF models is also strongly linear (R 2 = 0.99). Most of the fitting residuals are also within 5 ppm.
A more comprehensive study including the calculations with Dunning's correlation consistent basis sets [146,147], aug-cc-pVXZ (X = D, T and Q), was carried out for a limit number of complexes. The calculated shielding constants were shown in Figure 2. It is clear that the values did not converge well even at the aug-cc-pVTZ or aug-cc-pVQZ level for many geometries. However, the shielding constants calculated with different basis sets are highly correlated, which is consistent with the above observation. Jensen [16] proposed to use pcS-n basis sets with DFT methods to obtain the converged shielding constants. Here, we showed that even though there are large errors involved in the absolute shielding constants calculated with the commonly used basis sets, such as 6-31G(d) and 6-311+G(d,p), the errors are not random and a linear relationship exists between the values obtained with these small size basis sets and basis set limits. Therefore, these basis sets are still valuable in predicting 27 Al chemical shifts.
The MP2 model has been considered a more accurate model in electronic structure calculations. We also calculated the chemical shielding constants using this model with the 6-311+G(d,p) basis set. Because the large memory cost required and CPU time needed for the model, the calculations were done only for three species, [Al(OH 2 ) 6 ] 3+ , [AlF 6 ] 3− and [AlF 5 (OH 2 )] 2− . The calculated chemical shielding constants are also strongly correlated with the values obtained at the B3LYP/6-311+G(d,p) level, R 2 = 0.97. The slope of the best-fit equation is 0.98. The calculations indicate a trend, but more than three data sets (species) will have to be investigated before more firm quantitative relationships can be established.  4 ] − (B) and other complexes (C). The basis sets were labeled on the graphs. In Panels A and B, there are three data points for each basis set which represent the values calculated in vacuum, in the PCM environment of water and methanol. Twelve geometries in A were obtained at the B3LYP/aug-cc-pVQZ, B3LYP/aug-cc-pVTZ, B3LYP/aug-cc-pVDZ and B3LYP/6-311++G(d,p) levels in vacuum, in the PCM environment of methanol and water, and the data sets are colored in black, red, blue, green, cyan, magenta, yellow, brown, orange, pink, purple and gray, respectively. The last geometry in A, which is in light orange color, was obtained at the B3LYP/6-31G(d) level in the PCM environment of water. Ten geometries in B were obtained at the B3LYP/aug-cc-pVQZ, B3LYP/aug-cc-pVTZ, B3LYP/aug-cc-pVDZ and B3LYP/6-311++G(d,p) levels in vacuum, at the B3LYP/ aug-cc-pVQZ and B3LYP/aug-cc-pVTZ levels in the PCM environment of methanol, and at the B3LYP/aug-cc-pVQZ, B3LYP/aug-cc-pVTZ, B3LYP/aug-cc-pVDZ and B3LYP/ 6-31G(d) levels in the PCM environment of water, and the data sets are colored in black, red, blue, green, cyan, magenta, yellow, brown, orange and pink, respectively. Six 3+ , respectively. All the values obtained with a same geometry are connected with straight lines.

Comparison between Experimental and Calculated Chemical Shifts
From the previous subsection, we see that the shielding constants have a strong linear correlation among five levels of theory, though the absolute numbers differ from each other significantly. In practice, the chemical shifts are our primary concern, since the absolute chemical shielding constants are difficult to measure experimentally. Chemical shifts are related to shielding constants through the approximate equation, δ species = σ ref -σ species . 27 Al chemical shifts are usually referred to the aqueous [Al(OH 2 ) 6 ] 3+ species in acidic solutions. We use the averaged chemical shielding constant of the [Al(OH 2 ) 6 ] 3+ species at the corresponding level to calculate the chemical shifts of other species. The results are given in Table 1.     How do these models perform when compared with the experimental measurements? We searched the literature for experimental 27 Al NMR studies of Al(III) species and picked the species whose 27 Al chemical shifts were assigned. The measurements include both solution-state and solid-state magic-angle-spin (MAS) 27 Al NMR spectra. For the latter, the isotropic chemical shifts were calculated by fitting the simulated spectra with the experimental spectra. The solution-state NMR measurements were for the Al(III) species in both aqueous and organic solvents. Therefore, the experimental chemical shifts are measured at a wide range of different chemical environments. But the calculations were done for the isolated molecules. How strongly are the 27 Al chemical shifts dependent on the long-range chemical environment? The answer might be different for different species, but the environment effect might be small for most Al(III) species because the first coordination layer dominates the shielding effect to the 27 Al nucleus. For example, the chemical shifts of the tetrahedral-coordinated [Al(OC(CF 3 ) 3 ) 4 ] − with various kinds of counter ions in different solvents, range from 38.8 to 33.8 ppm [140][141][142][143]. On the other hand, the chemical environment has a significant effect on the small [AlF 6 ] 3− species and fast ligand exchange occurs in aqueous solutions [148,149]. The solid-state MAS NMR studies found the isotropic chemical shifts of [AlF 6 ] 3− range from −0.1 to −17.9 ppm for which the counter ions include H + , K + , NH 4 + and Rb + [150]. The values for [AlF 5 (OH 2 )] 2− obtained from the solid-state MAS NMR spectra range from −6.3 to −14.4 ppm [150]. These results reflect the sensitivity of chemical shifts of small size species to chemical environment. The environment could affect the 27 Al chemical shifts through two aspects. The first is that the geometry parameters, such as bond lengths and bond angles, of the first coordination layer could be perturbed by the packing forces, by the electrical field and/or by solvation. The other factor is that the chemical environment other than the first coordination layer could contribute to the shielding or to the deshielding effect to the metal center. For the isolated molecule models used to compute the chemical shieldings, counter ions and solvents were not included. Therefore, only the first effect was accounted for in the calculations. The calculated chemical shifts vs. the experimental chemical shifts are shown in Figure 3. Black squares are the calculated shifts at the GIAO-B3LYP/6-31G(d) level, and red circles and blue triangles are the calculated values at GIAO-B3LYP/6-311+G(d,p) and GIAO-HF/6-311+G(d,p), respectively.

Tetrahedral-coordinated Al(III) with C's and N's
Comparing the calculated chemical shifts at the GIAO-B3LYP/6-31G(d) level with the experimental values, the absolute errors for most species are less than 12 ppm except for three outliers. About one third (31%) of the absolute errors are less than 2 ppm, while more than 79% in total are less than 6 ppm. The overall root-mean-square difference (RMSD) is 6.5 ppm. If the three outliers whose absolute errors are larger than 12 ppm are removed, the RMSD decreases to 4.4 ppm. From this result, we see that the popular GIAO-B3LYP model with the 6-31G(d) basis set performs quite well for predicting the 27 Al chemical shifts using the X-ray measured geometries.
The three outliers are [AlF 6 ] 3− , [Al(NCCH 3 ) 6 ] 3+ and [Al 4 (μ 2 -OH) 2 (μ 2 -H -1 malato) 4 ] 2− . In calculating δ Al for [AlF 6 ] 3 the large error is probably due to its small size [150]. Calculations including the ion's environment beyond the binding ligands would probably significantly improve the result. For the acetonitrile complex, the ligand is disordered in the crystals and not accounting for the distribution in geometries might lead to a significant error [118]. For the polymeric malate complex, the 27 Al NMR spectra are not well defined, so there might be a significant error associated with the experimental determination of δ Al [73].  27 Al chemical shifts versus the calculated values using the X-ray crystallographic geometries at the GIAO-B3LYP/6-31G(d) (black squares), the B3LYP/6-311+G(d,p) (red circles) and the HF/6-311+G(d,p) (blue triangles) levels. The black solid line, the red dashed line and the dotted blue lines are the best-fit linear regression of the three data sets, respectively, and the best-fit equations are given in Equations 3-5.
When the larger basis set, 6-311+G(d,p), is used, the accuracy improves slightly. The RMSD between the calculated chemical shifts at the GIAO-B3LYP/6-311+G(d,p) level and the experimental values for all the species is 6.6 ppm and it decreases to 3.9 ppm after removing the three outliers whose absolute errors are larger than 12 ppm. About 40% of the absolute errors are less than 2 ppm, while 80% in total are less than 6 ppm.
The accuracy of the chemical shifts calculated at the GIAO-HF/6-311+G(d,p) level is slightly worse than that of the GIAO-B3LYP/6-31G(d) level. Although the RMSD for all the data points is 5.7 ppm which is less than the other two models, it decreases to 4.8 ppm after removing two outliers whose absolute errors are larger than 12 ppm. One third (33%) of the absolute errors are less than 2 ppm and 72% in total are less than 6 ppm.
By comparing the results of the three levels of theory, we see that they are similar in terms of the capability in predicting the 27 Al chemical shifts, though the GIAO-B3LYP/6-311+G(d,p) model is slightly better than the other two.
As seen from Figure 3, the calculated chemical shifts at three levels were well correlated with the experimental results, except for a few outliers. The best-fit equations are given in Equations 3~5 for the correlation between the calculated chemical shifts at the GIAO-B3LYP/6-31G(d), GIAO-B3LYP/6-311+G(d,p) and GIAO-HF/6-311+G(d,p) levels, respectively, versus the experimental values (black solid line, red dashed line and blue dotted line, respectively in Figure 3).

Geometry Optimization
In a practical application scenario, the structure of the proposed complex is unknown. Geometry optimization is used to obtain the structure in the first place. Does the error of the calculated geometries contribute significantly to an error in predicting 27 Al chemical shifts? To investigate this question, the experimental geometries of 33 complexes were optimized at the B3LYP/6-31G(d) level, followed by a shielding tensor calculation using the GIAO-B3LYP/6-31G(d) model.
The calculated δ Al are plotted against the experimental values in Figure 4. On one hand, the calculated δ Al are systematically larger than the experimental values; on the other hand, the slope of the best-fit equation which is given in Equation 6 and plotted as the black line in Figure 4 is 1.10. The 95% confidence interval for the slope is 1.04 to 1.17. The systematic larger calculated chemical shifts may be because the shielding constant of the reference [Al(OH 2 ) 6 ] 3+ species is overestimated in the optimized geometry. The calculated σ Al in the optimized geometry is 16.7 ppm larger than that in the X-ray geometry, which reflects the elongation of Al-O bonds by 0.06 Å in the optimized geometry compared with the mean value in the X-ray geometry. The latter deviation reflects that the single-molecule calculated geometry at the B3LYP/6-31G(d) level causes significant errors in predicting δ Al . Until such time as experimental techniques can unequivocally determine structures of molecules in liquid solvents, this issue cannot be unambiguously settled.
The RMSD for all the data points is 10.9 ppm. Even after removing the points with errors larger than 12 ppm, the RMSD is about 8.4 ppm. If we offset all the calculated δ Al by 15.4 ppm, i.e. to let the best-fit equation cross the origin point, the calculated RMSD is still as large as 7.7 ppm. Therefore, an accurate geometry is critical in calculating 27 Al chemical shifts.

Solvent Effect
Solvent effect plays an important role in solution-state 27 Al NMR. Without the constraint of the periodic boundary condition in the crystals, the geometries will be slightly different, the long-range interactions will be dramatically different and dynamical effect will be also different. In this study, we only examined how the commonly used implicit solvation model affects the prediction of 27 Al chemical shifts. The self-consistent reaction field (SCRF) method [151] with the polarizable continuum model (PCM) [152] of water and/or methanol was applied to the calculation of the shielding constants and the geometry optimization. In Figure 2, we have seen that the calculated shielding constants of [Al(OH 2 ) 6 ] 3+ and [Al(OH) 4 ] − in the PCM environment of water or methanol are nearly the same as those obtained in vacuum using the same geometry. However, the calculated shielding constants using the geometries optimized in the PCM environment of water or methanol differ by 5~10 ppm from those obtained with the optimized geometries in vacuum. The differences in [Al(OH 2 ) 6 ] 3+ are larger than those in [Al(OH) 4  Two more groups of calculations were carried out on 10 representative complexes. In the first group of calculations, the influence of the PCM model to the calculation of the shielding constants was examined using the X-ray geometries at the GIAO-B3LYP/6-31G(d) level. In the second group of calculations, both the geometry optimization and the shielding constants calculation were carried out in the PCM environment of water at the B3LYP/6-31G(d) level. The calculated σ Al values were compared with the values obtained in vacuum using the X-ray geometries. The differences are plotted in Figure 5, where the blue and red bars are for the first and second groups of calculations, respectively. Again, we see that the PCM model only causes a minor change to the calculated σ Al if the same geometry was used. However, if the geometry optimization was applied, significant changes were observed. This result again says that the geometry parameters of the complexes have a big influence on the calculated σ Al .

Computational Details
The crystal structures of 114 Al(III) complexes, determined by single-crystal X-ray diffraction, were obtained from the Cambridge Crystal Structural Database (Cambridge Crystallographic Data Centre, Cambridge, UK). They were used to build the molecular structures of the species. If there was more than one Al(III) species in the asymmetric unit, a molecular structure was built for each structure. For the disordered crystal structures, only the structure with major population was used to build the molecular model. In total, 213 molecular models were created from 206 crystal structures for 114 species. 33 complexes of small to intermediate sizes were optimized using the X-ray structures as the starting point. Figure 5. Difference between the calculated σ Al for 10 Al(III) complexes. The blue bars are the differences between the values in the PCM environment of water and in vacuum using the X-ray geometries. The red bars are the differences between the values calculated using the geometries optimized at the B3LYP/6-31G(d) level in the PCM environment of water and those calculated using the X-ray geometries in vacuum. All the σ Al 3+ and Al(maltol) 3 for 1-10, respectively. Label of the complexes Diff. between the calcd. tensors, ppm Only a single molecule was included in each molecular model. Hydrogen positions were put in their ideal positions with the aid of GaussView 3.0 (Gaussian, Inc., Wallingford, CT, USA) if they were missing from the crystal structures. All the counter ions and solvent molecules were excluded from the models. Chemical shielding tensors were calculated for these X-ray geometry models and the models with further geometry optimization at the B3LYP/6-31G(d) level.

Conclusions
The main purpose of this work is to give a quantitative evaluation of the errors involved in the calculation of 27 Al chemical shifts. We differentiated the source of the errors due to the limitations of the computational models in calculating the chemical shielding tensors and the inaccuracy in geometry prediction by using the X-ray measured geometries and the calculated geometries to compute chemical shifts. The RMSD between the calculated shifts and the experimental values is approximately 5 ppm if the X-ray crystallographic structures were used, while the RMSD is more than 10 ppm for the optimized geometries at the B3LYP/6-31G(d) level. Although we could not isolate the error due to the structural uncertainties in the X-ray crystallography, the popular GIAO-B3LYP/6-31G(d) model produces quite accurate chemical shifts.
There are other error sources not considered in this study, such as the quality of DFT functionals and dynamical effects [161]. In addition, the structure of a complex in solution would be different from its crystal structure; the solvent effect is also more complicated than what has been studied here. However, the results presented in this study shows that the error associated with the shielding constants calculation is small compared with the error in the prediction of geometry. In particular, Jensen showed that the error in the shielding constants calculation can be controlled to less than 1 ppm.
In the future, we should focus on predicting a more accurate geometry in solution or an ensemble of structural configurations generated using molecular dynamics.