DFT calculations of 1H- and 13C-NMR chemical shifts of 3-methyl-1-phenyl-4-(phenyldiazenyl)-1H-pyrazol-5-amine in solution

Geometries of the 3-methyl-1-phenyl-4-(phenyldiazenyl)-1H-pyrazol-5-amine azo-dye compound and its tautomer were optimized using B3LYP and M06-2X functionals in coupling with TZVP and 6–311 + G(d,p) basis sets. The 1H- and 13C-NMR chemical shifts of all species were predicted using 13 density functional theory (DFT) approaches in coupling with TZVP and 6–311 + G(d,p) basis sets at the different optimized geometries by applying the using GIAO method using the eight geometries. The selected functionals are characterized by having different amount of Hartree–Fock exchange. The selected DFT methods were B3LYP, M06-2X, BP86, B97XD, TPSSTPSS, PBE1PBE, CAM-B3LYP, wB97XD, LSDA, HSEH1PBE, PW91PW91, LC-WPBE, and B3PW91. The results obtained were compared with the available experimental data using different statistical descriptors such as root mean square error (RMSE) and maximum absolute error (MAE). Results revealed that the prediction of the 1H-NMR chemical shifts has more significant dependence on the applied geometry than that of the prediction of the 13C-NMR chemical shifts. Among all the examined functionals, B97D and TPSSTPSS functionals were found to be the most accurate ones, while the M06-2X functional is the least accurate one. Results also revealed that the prediction of NMR chemical shifts using TZVP basis sets results is more accurate results than 6–311 + G(2d,p) basis set.

Azo-Dyes, consisting of 60-70% of dyes, have a great importance in science and technology, such as their applicability in optical materials, molecular data processing, ion sensors, and biological applications 1 . Simultaneously, most of the azo dyes are tautomeric or potentially tautomeric [2][3][4][5] . The tautomerism phenomenon, which plays a significant role in the stability and the color of the dyes, has been extensive studied during the last decades [6][7][8][9][10][11] . The intramolecular 1,3-proton transfer process, which is very fast at the NMR time scale, occurs by the migration of a proton from the amino nitrogen atom to the carbonyl oxygen atom. Thus, to overcome this problem, an averaged signal is typically measured 7,12 . Therefore, reference compounds are needed to evaluate the proportion of the azo and hydrazo tautomers. The selection of these references that represent the pure azo and hydrazo tautomers is necessary. Thus, there is a need for such reference compounds to be identified and used in the NMR investigations of tautomeric azo dyes.
Up to now, nuclear magnetic resonance (NMR) spectroscopy has been one of the most powerful and reliable methods for structure elucidation. 1 H-and 13 C-NMR are the most widely studied nuclei. However, other nuclei such as 15 N, 31 P, 18 O.etc., are also under consideration, which can be attributed to the corresponding routinely available techniques readily yielding an overabundance of structural and kinetical information for the proper identification of novel compounds. Unfortunately, complex compounds suffer from the main problem in the full signal assignment, which is sometimes impossible, even with the utilization of multidimensional NMR experiments.
Calculating the nuclear magnetic shielding tensor (NMST) is required for the theoretical NMR investigations. For this reason, several computational methods were developed such as the GIAO (gauge-independent atomic orbitals) [13][14][15][16][17] , the IGLO (an individual gauge for localized orbital) 18 , and the CGST (continuous set of gauge transformation) 19,20 methods. In general, the GIAO method has been shown to achieve faster convergence of the shielding value with respect to the size of the basis set and, consequently, has been the most commonly employed procedure 21 .
Last past years have witnessed a great development in the growth of quantum mechanical calculations quickly, which aids the determination of compound structures [22][23][24] . Different levels of theory can be effectively used to compute the NMR shielding constants. These levels include coupled-cluster theories (e.g., coupled-cluster singles-doubles (CCSD) and coupled-cluster singles-doubles-perturbative-triples, CCSD(T)) 23 and Møller-Plesset (MP) second-order perturbation theory (MP2) 22,25 . Exactly, the most accurate sets and highly comparable with the experimental data of benchmark NMR shielding constants were performed using MP2, CCSD, and CCSD(T) methods in various basis sets. Unfortunately, these methods are highly expensive and computationally time consumable, and they are only exclusively applicable to small model molecules and not for larger ones [22][23][24][25] .
In spite of being less accurate, DFT methods, on the opposing, can handle structures composed of several hundreds of atoms [26][27][28][29][30][31] . DFT methods allow the calculation of useful chemical properties, structural and electronic analysis of molecules and the theoretical values of the 13 C-and 1 H-NMR chemical shifts with a high degree of accuracy compared to the experimental data 32,33 . Bango et al. 27 determined the 13 C-NMR chemical shifts experimentally and theoretically by DFT calculations B3LYP with the 6-31G(d,p) and cc-pVTZ basis sets. Bally and Rablen 28 proved that B3LYP/6-31G(d,p)u + 1 s (one of the most economical methods) and B3LYP/ccp-pVTZ (ca.8 times more expensive in terms of CPU time) predicted the 1 H-NMR with root mean square of 0.5 Hz, compared to the experimental results. Broadbelt et al. 34 accurately computed the 1 H-and 13 C-NMR chemical shifts using the BMK method with different basis sets in toluene-d 8 solvent using IEFPCM solvation model. The agreement between the experimental data and the DFT computed results was improved by inclusion of vibrational corrections for DFT methods, which contrasts the behavior of the coupled-cluster methods 23,24,29,[34][35][36][37][38][39][40] . This methodology could be applied as a promising alternative in support of traditional NMR experimental techniques. Furthermore, previous studies 40,41 showed that DFT theory brings a respected tool for the prediction of 13 C-NMR properties, even with a small basis set. Additionally, it was recommended to use a larger basis set, including polarization functions, for a more accurate prediction of trends in the experimental 1 H-NMR chemical shifts 41 . Furthermore, most of studies that related to the computational of 1 H-and 13 C-NMR were mainly focused on DFT and ab initio performance for structures in CDCl3, which expected due to the fact that a large number of experimental spectra are taken in deuterated chloroform and there is an abundance of data for numerous molecules and functional groups.
Test set. The 3-methyl-1-phenyl-4-(phenyldiazenyl)-1H-pyrazol-5-amine (T1) may have another conformer (T2) that can rapidly interconvert at the NMR time scale with relative energy of about 9.3 kcal/mol 10 ( Fig. 1). The 1 H-and 13 C-NMR experimental chemical shifts of the main species (T1) in CDCl 3 were taken exclusively from Deneva et al. 10 . All protons that are bonded to the nitrogen atoms were not included in the validation test because they are strongly affected by the ability to form hydrogen bonds as well as the solvent environment 34 . This study is undertaken to present and discuss the following issues: (i) computation of 1 H-and 13 C-NMR chemical shifts for the species under probe at different geometries that optimized using B3LYP and M06-2X approaches with TZVP and 6-311 + G(2d,p) basis sets within the DFT framework using the GIAO 42 method and (ii) statistical validation of the theoretical NMR chemical shifts based on the available experimental data.
Computational methodology. Geometries of the investigated species were optimized using two-hybrid DFT approaches, B: B3LYP 43,44 and the M: M06-2x (meta-GGA) 35,45 in acetonitrile as a solvent using the Polarizable Continuum Model (the integral equation formalism variant, IEFPCM) 48 . Two basis sets were included in the optimization process, which are: (i) the triple-ζ Pople basis sets augmented with additional polarization functions 6-311 + G(2d,p) 46 and (ii) triple zeta for valence electrons plus polarization function, TZVP 47 . Therefore, To ensure that all geometries are minima to the electronic potential energy surface, the vibrational frequencies of each species were calculated at the same levels of theories, and the absence of any imaginary frequencies was noted. The isotropic magnetic shielding tensors were computed using thirteen different XC DFT approaches, applying the gauge-independent atomic orbital (GIAO) 42 approach in coupling with the same basis sets CHCl 3 solvent. The classification of the thirteen XC different functionals are given in supplementary materials. In order to the simplify the following of our study, DFT approaches B3LYP and M06-2X are given the abbreviations T and P, respectively. The two bases sets TZVP and B3LYP/6-311 + G(2d,p) are given the abbreviation T and P, respectively. So that, BT, BP, MT and MP optimization levels were considered. To calculate the NMR chemical shifts, www.nature.com/scientificreports/ XC DFT approach is coupled with the two bases sets at the considered geometries. In total, eight combinations of theoretical approach will be used to compute the NMR chemical shifts (see Table S1 of the supporting information). All quantum chemical calculations were performed using Gaussian 09 computational software package 49 . Because NMR chemical shifts result from a dynamic process, the isotropic magnetic shielding tensors were averaged over all symmetry related carbons and hydrogens where applicable. For each 1 H and 13 C nucleus, the average isotropic magnetic tensor of the two tautomeric forms σ cal was taken. Then, the isotropic chemical shifts δ cal was defined as δ cal = σ TMS − σ cal , where σ TMS is the anisotropic shielding constant of 1 H and 13 C in trimethylsilane (TMS). Notice that σ TMS values of 1 H-and 13 C-NMR of the TMS were computed using the same strategy that used for the azo-dye species and based on the MP2 geometry (all σ TMS are listed in Table S2 of the supporting information). Also, it should be noted that the σ TMS for the 1 H-and 13 C-NMR were calculated as the average for the 12 hydrogen nuclei and 4 carbon atoms, respectively. The δ cal values were collected in Tables S3 and S4 of the supporting information.
Due to the possibility of the tautomerization of the investigated species at the NMR time, the scaled 1 H-and 13 C-NMR chemical shifts of the hydrogen and carbon nuclei of the T1 and T2 species were averaged. The averaged 1 H-and 13 C-NMR chemical shifts were tested for possible outliers using 1.5xIRQ rule, and no outliers were detected. Computation of the different statistical descriptors was performed using Microsoft Excel 2016 using the following equations:

Results and discussion
The optimized structures of the two envisaged tautomers of the azo-dye species as obtained using B3LYP/TZVP level of theory in acetonitrile are shown in Figure S1 of the supporting information. The statistical parameters (MAE, MAX, RMSE and MAPE) are measured for suitable errors. Among all of the statistical descriptors, RMSE is the most powerful descriptor that can be used to compare the different DFT methods. The smaller the RMSE parameter, the more accurate the DFT approach to estimate the 1 H-and 13 C-NMR chemical shifts. Tables S5 and S6 Tables S5 and S6). The poorest regression was obtained for M06-2X functional, while the best one was obtained for TPSSTPSS and B97D functionals. Figure S1 shows graphically by column the statistical descriptors (RMSE, MAX, and MAE) computed for the 1 H-NMR chemical shift deviations. For all combinations, the results show that the MAX descriptor (blue column), which corresponds to the maximum deviation with respect to the experimental data, is much higher than that of RMSE (red column) and MAE (dark red column) descriptors. The corresponding MAX arises from the computation of 1 H-NMR chemical shifts of the hydrogen atom attached to C12 (Fig. 2). The reason may be attributed to the possible interaction (by noncovalent van der Waal) between the hydrogen atom (van der Waal radii = 2.02 Å) and the N(β) (van der Waal radii = 2.25 Å) of the azo-dye group, resulting in a weak hydrogen bond. The geometrical structure of the T1 and T2 tautomers show that the distance between the hydrogen atom and N(β) atom is an average of 2.450 Å, which strongly confirms our suggestion.

Benchmark of 1 H NMR chemical shifts.
The data given in Tables S2 and S5, show that the signed deviations (ME) are typically positive, indicating that the computed 1 H NMR chemical shifts are exclusively higher than that of the experimental ones. It is also found that ~ 93.5% of the ME values are positive. The number of negative ME values are ~ 96 out of 1456, with a percentage of ~ 6.5%. Moreover, most of the negative values belong to the hydrogen atom attached to C10 atom, which may be attributed to the consideration of the average value of the isotropic magnetic tensors of the three hydrogen atoms attached to C10.
Scientific Reports | (2022) 12:17798 | https://doi.org/10.1038/s41598-022-22900-y www.nature.com/scientificreports/ Basis set effect. In order to analyze the effect of the basis set that couples with the XC functional on the 1H-NMR chemical shift, two basis sets TZVP and 6-311 + G(2d,p) were selected. For the sake of comparison, the eight combinations (Table S4) were divided into two parts based on the basis sets, which coupled with NMR XC functional, as follows: • Part (1) in which the NMR functional combined with TZVP basis set, and it includes C1, C3, C5 and C7 combinations, and • Part (2) in which the NMR functional combined with 6-311 + G(d,p) basis set, includes C2, C4, C6, and C8 combinations.
To simplify the discussion, the average absolute differences in RMSE values were considered and plotted as pie charts (Fig. 3). In this case, changing the basis set coupled with the XC functionals to compute the 1H-NMR chemical shifts is considered. Taking into account the same geometry level, the results show that the accuracy of the 1 H-NMR chemical shifts as computed by TZVP basis set is higher than that computed by 6-311 + (2d,p) basis set. Actually, it is found that the results obtained by TZVP basis set are 0.069-0.101 ppm higher than those obtained by 6-311 + G(2d,p) one. Moreover, the highest differences are found in the case of using 6-311 + G(2d,p) basis set in the optimization process (combinations C3, C4, C7 and C8). These results led us to conclude that the  (Table S1). Figure 4 shows the graphical representation by columns for the average differences of the RMSE values obtained by taking into account the same basis sets that are coupled with the different XC functionals. Computation of the 1 H-NMR chemical shifts using DFT/TZVP at B3LYP/6-311 + G(2d,p) geometries (Combination 3) is more accurate than those at B3LYP/TZVP geometries (C1) (Fig. 4a). Whereas the computed 1 H-NMR chemical shifts at M06-2X/6-311 + (2d,p) geometry (combination 5) are more accurate than those at the M06-2X/TZVP geometry (Combination 7). On the other hand, using B3LYP/6-311 + G(2d,p) geometry can compute more accurate 1 H-NMR chemical shifts than the M06-2X/TZVP geometries (combinations C3 and C5). Similar results can also be found for the case of DFT/6-311 + G(2d,p) NMR level at the four geometries (Fig. 4b).
Based on the above results, one may conclude that taking into account either M06-2X or B3LYP functionals, using of the 6-311 + G(2d,p) basis set in the optimization steps computes accurate more accurate 1 H-NMR chemical shifts than the TZVP basis sets.

XC functionals effect.
Besides to the basis set and geometry level effects, one of the main questions of this study is whether the DFT functionals can compute the most accurate 1 H-NMR chemical shifts? The answer of this question can be achieved by a careful looking at the results obtained indicates, in general, that M06-2X, LC-WPBE, and wB97XD functionals compute the poorest 1 H-NMR chemical shifts with very high RMSE values   Tables S5 and S6 and shown in Figure S1). Whereas, the highest degree of accuracy is obtained by B97D, TPSSTPSS, and BP86 functionals, which is in an excellent agreement with the literature 50 . This finding can be attributed to the amount of HF exchange in the functionals used: (M06-2X (54%), LSDA (45%), and wB97XD (25%) 50 . (Table S1) to compute the 1 H-NMR chemical shifts. To shed some light on the effect of the several combinations on the accuracy of the computed 1 H-NMR chemical, Fig. 5a,b shows the RMSE% of the B97D and TPSSTPSS functionals. It is clearly seen that the most accurate results are obtained for the DFT approach with the TBP combination (C3). Indeed, the RMSE% are 9.5% and 10.8% for B97D and TPSSTPSS functionals, respectively. Whereas the least accurate results were computed using the DFT apprroach with PBT combination (C4), with RMSE% of 15.2% and 14.9% for B97D and TPSSTPSS functionals, respectively. These results led us to safely conclude that the 1 H-NMR chemical shifts, with a high degree of accuracy, can be computed by either B97D or TPSSTPSS functionals in coupling with TBP.

Combination effect. Eight combinations (C1-C8) have been considered
To get an exclusive order for the different DFT approaches, the computed descriptors obtained for the TBT, TBP, TMT and TMP combinations were averaged, and the results are listed in Table 1. The most important statistical descriptors (RMSE, MAX, MAE, and MAPE) are graphically shown in Figure S3. It is clearly obvious from Table 1 and Figure S3 that the corresponding RMSE, MAX, MAE, MSE, ME and MAPE descriptors are 0.006 (1.6%), 0.020 (3.5%), 0.006 (1.6%), 0.004 (2.8%), 0.008 (2.7%) and 0.002 (3.4%) ppm, respectively. These results indicate that there are no significant differences in the computed 1 H-NMR chemical shifts based on using either B3LYP or M06-2X geometry. These findings led us to conclude that the main factors that can influence the calculation of the 1 H-NMR chemicals are the DFT approach and the selected basis sets, regardless of the geometry level used.
Benchmark of 13 C-NMR chemical shifts. In the following, we are going through the analysis of the results obtained for the prediction of 13 C-NMR. Tables S7 and S8 list all the statistical descriptors (RMSE, MSE, MAE, MAX, MAPE, R 2 , slope and intercept) computed for 13 C-NMR deviations. An excellent linear regression correlation coefficient (R 2 > 0.986-0.993) was obtained. It is also found that, however, the slope (0.963-1.149 ppm), 84% of these values are in agreement with the ideal values of slope (Fig. 6). The intercept (0.342-3.786 ppm) differed slightly from the ideal values for all combinations (see Tables S5 and S6). The best regressions are obtained for TPSSTPSS and B97D functionals, while the worst ones are obtained for the M06-2X functional. Figure 7 shows only the statistical descriptors (RMSE, MAX, MAE, and MAPE) for the best 13 C-NMR results (B97D, TPSSTPSS, BP86, B3PW91, PBE1PBE and PW91PW91 functionals), which agrees with the previous results 39 . With the exception to the results of wB97XD and CAM-B3LYP functionals, the accuracies of the different DFT approaches that were used to compute the 1 H-and 13 C-NMR chemical shifts are almost the same. Therefore, the most accurate 13 C-NMR chemical shifts are computed using B97D, TPSSTPSS, and BP86 functionals. Whereas the least accurate results are obtained using M06-2X, LC-WPBE, and CAM-B3LYP functionals. C-NMR chemical shifts, we apply the same scenario we used for 1 H-NMR results. In comparison with the experimental results 10 , Fig. 8 shows the 3D color pie charts of the RMSE computed for the 13 C-NMR chemical shifts obtained using B97D and TPSSTPSS, BP86, and PW91PW91. For both functionals, the least deviations of 12.13 and 12.22% (TPSSTPSS)) were computed using the B97D/TZVP//M06-2X/TZVP and TPSSTPSS/TZVP// M06-2X/TZVP levels, respectively, while the highest deviations (least accurate) were computed at 12.83 and 12.76% were computed using the B97D/6-311 + G(2d,p)//M06-2X/6-311 + G(2d,p) and TPSSTPSS/TZVP// B3lYP/6-311 + G(2d,p) levels, respectively. Based on these results, one concludes that the prediction of the 13 C NMR chemical shifts can be accurately computed by using the D97D and/or TPSSTPSS in coupling with the TZVP basis set at the M06-2X/TZVP geometry.
To rank the different DFT approaches that we have used to compute the 13 C-NMR chemical shifts, the data in Tables S7 and S8 were averaged and classified into two main combinations. The first one averages the computed descriptors based on the B3LYP geometry with the TZVP and 6-311 + G(2d,p) basis sets, while the second one averages the computed descriptors based on M06-2X geometry with the TZVP and 6-311 + G(2d,p) basis sets. The results are listed in Table 2, and the data are sorted by RMSE descriptors.
For the B97D, TPSSTPSS, BP86, PW91PW91and B3PW91, the average statistical errors will be as follows: RMSE = 4.74 and 4.32 ppm, MAX = 9.60 and 9.10 ppm, MAE = 3.83 and 3.50 ppm, ME = 2.10 and 1.31 ppm, and MAPE = 4.55 and 4.16% for B3LYP and M06-2X geometries, respectively. These results strongly confirm that the prediction of the 13 C-NMR chemical shifts on the base of the M06-2X geometry, regardless of the NMR functional and basis set, is more accurate than the B3LYP geometry.
Finally, the data given in Tables S7 and S8, show that the signed deviations (ME) are typically positive, indicating that the computed 13 C NMR chemical shifts are exclusively higher than the experimental ones. It is found that > 87% of the ME values are positive. Indeed, the number of negative ME equals approximately 340 out of 2496 values, with a percentage of ~ 13%). The maximum number of the negative ME values was computed for TPSSTPSS functionals with a percentage of ~ 40%.

Conclusions
Eight different combinations of DFT functionals, geometry level, and basis sets have been examined for the prediction of 1Hand 13C-NMR chemical shifts in 3-methyl-1-phenyl-4-(phenyldiazenyl)-1H-pyrazol-5-amine. The computed results have been predicted using experimental shifts as a reference using several statistical descriptors. From the calculated data, the following remarked conclusion can be extracted: • Excellent linear correlations have been obtained between the DFT-calculated and experimental 1 H-and 13C NMR chemical shifts for all functional, basis set, and geometry level tested. • The slopes of the linear correlations are within the ideal value for the 13 C-NMR chemical, and it is slightly deviated in the case of 1 H-NMR chemical shifts. • The accuracy of the 1 H-NMR chemical shifts computed using the TVZP basis set on the four suggested geometries was higher than those obtained by the 6-311 + G(2d,p) basis set. • The B97D, TPSSTPSS, and BP86 functionals can compute the 1 H-and 13C-NMR chemical shifts with a reasonable degree of accuracy. • For 13 C-NMR chemical shifts, the least errors have been obtained using the combination of B97D and TPSST-PSS functionals with TZVP basis sets with taken into account M06-2X/TZVP geometry.

Data availability
All data generated or analyzed during this study are included in this published article [supporting informationupdated]. Additionally, the datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.