Calculating the Aqueous pKa of Phenols: Predictions for Antioxidants and Cannabinoids

We aim to develop a theoretical methodology for the accurate aqueous pKa prediction of structurally complex phenolic antioxidants and cannabinoids. In this study, five functionals (M06-2X, B3LYP, BHandHLYP, PBE0, and TPSS) and two solvent models (SMD and PCM) were combined with the 6-311++G(d,p) basis set to predict pKa values for twenty structurally simple phenols. None of the direct calculations produced good results. However, the correlations between the calculated Gibbs energy difference of each acid and its conjugate base, ΔGaq(BA)°=ΔGaqA−°−ΔGaq(HA)°, and the experimental aqueous pKa values had superior predictive accuracy, which was also tested relative to an independent set of ten molecules of which six were structurally complex phenols. New correlations were built with twenty-seven phenols (including the phenols with experimental pKa values from the test set), which were used to make predictions. The best correlation equations used the PCM method and produced mean absolute errors of 0.26–0.27 pKa units and R2 values of 0.957–0.960. The average range of predictions for the potential antioxidants (cannabinoids) was 0.15 (0.25) pKa units, which indicates good agreement between our methodologies. The new correlation equations could be used to make pKa predictions for other phenols in water and potentially in other solvents where they might be more soluble.


Introduction
Acid dissociation constants (K a , pK a = −log K a ) are crucial physico-chemical quantities that impact chemical, environmental, and biochemical research [1][2][3][4][5][6]. Biochemical kinetic and thermodynamic studies involving acids require the calculation of molar fractions or Gibbs free energies of reaction at physiological pH for which aqueous pK a values are necessary [7][8][9][10][11]. Accurate predictions of aqueous pK a values can also be used to predict pK a values in non-aqueous environments [12]. The quest for determining reliable aqueous pK a values for complex phenolic compounds, including cannabinoids, has been motivated by our antioxidant studies on this family of compounds.
Choosing the best methods for obtaining reliable experimental pK a values can be challenging (due to low solubility, difficulty isolating, high reactivity, and variable ionic strength solutions) and time consuming [13][14][15]. Hence, theoretical calculations are a promising alternative. A traditional methodology uses thermodynamic cycles [14][15][16][17][18][19], which combine experimental or calculated ab initio gas phase Gibbs free energies with calculated solution Gibbs free energies. Another approach uses the dissociation equilibrium, HA (aq) A − (aq) + H + (aq) , and requires experimental data for H + , which is quite variable [14,15,20,21]. Other acid-base equilibria can be applied as well, relative to a reference acid whose experimental pK a is required [22,23]. Alternatively, various linear correlations between calculated properties (in the gas phase or in solution) and experimental pK a values have shown to have important predictive value [24][25][26][27][28]. The application of density functional theory methods combined with continuum solvation methods cosmetics, dyes, and photographic developers [50,51], and stilbenes related to resveratrol (41 and 42) [49]. Another group of ten phenols ( (21)(22)(23)(24)(25)(26)(27)(28)(29)(30), displayed in Figure 2) is used as an inde ent test set to compare our predictions to previously reported experimental or theo aqueous pKa values. This test group includes six complex phenols. The best meth gies are later used to predict the aqueous pKa values of complex phenols with po antioxidant properties that are currently under study by our group [45][46][47]. This gr      Furthermore, aqueous pK a predictions will also be made for of a set of nine cannabinoids displayed in Figure 4. Cannabinoids are phytochemicals found in the Cannabis plant [52]; nevertheless, this term is also used for any substance which interacts with the endocannabinoid system, including drugs that bear no resemblance to plant-derived cannabinoids [53]. Given the current pharmacological interest, inherent legalities, and little data available for cannabinoids, we thought that it would be appropriate to investigate these molecules that also have a phenol ring in their basic structure. Increasing evidence indicates that certain cannabinoids are effective antioxidants, in addition to their therapeutic uses [54][55][56][57][58][59]. For this study, we have chosen eight phytocannabinoids (29,30,(43)(44)(45), and 48-50), which are important components in the Cannabis sativa plant, and two synthetic cannabinoids (46 and 47), all of which are being investigated for potential therapeutic uses. The test set includes ∆ 9 -tetrahydrocannabinol (∆ 9 -THC, 29) and cannabidiol (CBD, 30), which have an experimental and a predicted aqueous pK a value reported, respectively. For molecules with stereocenters, labels have been added to identify the stereoisomer used in each case (see Figures 2 and 4), although the calculated pK a is not affected by this. Furthermore, aqueous pKa predictions will also be made for of a set of nine canna binoids displayed in Figure 4. Cannabinoids are phytochemicals found in the Cannabi plant [52]; nevertheless, this term is also used for any substance which interacts with th endocannabinoid system, including drugs that bear no resemblance to plant-derived can nabinoids [53]. Given the current pharmacological interest, inherent legalities, and little data available for cannabinoids, we thought that it would be appropriate to investigate these molecules that also have a phenol ring in their basic structure. Increasing evidence indicates that certain cannabinoids are effective antioxidants, in addition to their thera peutic uses [54][55][56][57][58][59]. For this study, we have chosen eight phytocannabinoids (29,30,(43)(44)(45), and 48-50), which are important components in the Cannabis sativa plant, and two syn thetic cannabinoids (46 and 47), all of which are being investigated for potential therapeu tic uses. The test set includes ∆ 9 -tetrahydrocannabinol (∆ 9 -THC, 29) and cannabidiol (CBD 30), which have an experimental and a predicted aqueous pKa value reported, respec tively. For molecules with stereocenters, labels have been added to identify the stereoiso mer used in each case (see Figures 2 and 4), although the calculated pKa is not affected by this.
Geometry optimizations were followed by a frequency calculation, both including solvent effects, at the same level of theory to determine the nature of the stationary points, and all structures were confirmed to be local minima in which no imaginary frequencies were present. Given that the conformation used for an acid and its conjugate base can greatly affect the calculated pK a value [67], every effort was made to ensure the most stable conformer was optimized in each case. For example, in 2-substituted halogen groups (molecules 6, 9, 21, and 22), SMD typically favoured the OH facing the halogen while PCM favoured the opposite case. Additionally, the pK a of the (amino) protonated form of molecules 15, 17, and 24-27 was computed to assess the possibility of concurrent microequilibria, of which none were considered significant. The optimized structures of all the species considered in this study at the M06-2X(SMD)/6-311++G(d,p) level of theory are provided in the Supporting Information. The absolute aqueous Gibbs free energies of the different species considered in this study at the various levels of theory at 298.15 K are reported in Tables S1-S5 of the Supporting Information.
This study explores three aqueous acid-base dissociation equilibria (Reactions (R1)-(R3), shown below as (R1)-(R3)) for calculating the absolute pK a value of an acid, HA, whose conjugate base is A − . Reactions (R1) and (R3) are standard ones used in theoretical pK a determinations and, together with Reaction (R2), have produced reasonable results for phenolic molecules [23].
RTln (10) (R1)  [20,23,30,31,38], because it has been corrected for the 1 M reference state and it has sometimes led to a good reproducibility of experimental aqueous pK a values.
Reaction (R2) includes explicit water molecules solvating the anions, which attempts to simulate the solute-solvent interactions that are not fully modelled in the implicit solvation models used. Hydrogen bonding between a water molecule and the conjugate base in A − (H 2 O) can significantly stabilize the charged species. The species OH − (3H 2 O) is the OH − ion solvated by three water molecules, while the species indicated as 3H 2 O is a water cluster of three molecules. The concentration of water used is 55.55 mol/L [23].
Reaction (R3) involves a reference acid similar in structure to the acid studied for which a reliable pK a value is available. For all molecules except phenol, phenol was used as the reference acid (experimental pK a = 9.99). For the determination of phenol, 3-methoxyphenol was used as the reference acid (experimental pK a = 9.65). The uneven distribution of charges across the equilibrium of Reaction (R1) makes this approach more prone to errors. Reactions (R2) and (R3) do not encounter this difficulty since both sides of their equilibria are balanced in terms of charges, thus contributing to better results due to the cancellation of errors. Reactions (R1)-(R3) contain the same number of computed reactant and product species; thus, reference state conversions are not needed.
Molecules 30 and 50 both contain two equivalent deprotonation sites through rotation of the sigma bond connecting the benzene ring to its substituent. As a result, the acid equilibrium constant should be doubled [68]. Accordingly, pK a values for molecules 30 and 50 must be corrected by subtracting log(2) from their respective calculated values using Reactions (R1)-(R3) or the correlation equations [38]. The deprotonation of both sites was tested, and the most stable conformer was used in each case.

Exploring Several Methodologies for the Direct Calculation of pK a Values
Eleven phenols of varying structural complexities are used to initially test the aqueous pK a calculations using five functionals, two solvation models, and three acid-base equilibria. The training set is then increased to twenty phenols (see Figure 1) for all functionals other than PBE0 because of significant technical issues. Mean absolute errors (MAEs), used to assess prediction quality, have been compiled in Table 1. The direct calculation results using the SMD solvation method are always much better than the corresponding results using PCM in each case considered, except when using the M06-2X and TPSS functionals with Reaction (R2). While there are differences between the results obtained with the different functionals using each reaction scheme, the reaction scheme used is much more impactful on the overall accuracy of the methodology applied, and that is how we have organized the discussion that follows. Table 1. Summary of mean absolute pK a errors (MAE) for the phenols in the training set (obtained from direct calculations and using the correlations between the experimental pK a values and calculated ∆G • values) at several levels of theory using Reactions (R1)-(R3) at 298.15 K.   Tables S6 and S7 (Tables S8 and S9) of the Supporting Information when using Reaction (R1) (Reaction (R2)). The MAE ranged between 3.09 and 5.24 (4.73 and 7.38) pK a units when using Reaction (R1) (Reaction (R2)) and the SMD solvation method for the set of twenty (eleven) phenols. Larger errors were usually obtained with the PCM solvation method, between 4.76 and 7.05 (4.75 and 6.83) pK a units when using Reaction (R1) (Reaction (R2)). With very few exceptions, the calculated aqueous pK a values were always overestimated (with positive errors).

Solvent
It is worth mentioning that Reaction (R1) was shown to drastically improve the pK a calculations for phenolic molecules (as well as alcohols, hydroperoxides, and thiols) when three explicit water molecules are added to the HA and A − species, while using the value of −270.29 kcal/mol for the aqueous Gibbs energy of H + after corrections [30,31]. This approach produced MAE of 0.45 pK a units for phenols at the best level of theory reported, B3LYP(SMD)/6-311++G(d,p) [31]. The M06-2X errors were still significantly large using this approach with mean signed errors of −1.40 pK a units.
Reaction (R2), using the Gaussian03 implementation of PCM through single-point energy calculations, was reported to produce much smaller MAE values which were generally in the 0.77-0.86 pK a unit range using similar functionals [23]. Our approach, including PCM as implemented in Gaussian16 in both geometry optimizations and frequency calculations, leads to much larger errors for direct pK a predictions. Given the results obtained for eleven phenols using Reaction (R2), which requires the additional calculation of the explicitly solvated conjugate base, we decided not to take it into account any further.

Results Obtained from the Direct Application of Reaction (R3)
The best direct calculation results are obtained using Reaction (R3), and M06-2X was the best-performing functional with both solvation models. The calculated individual pK a values and their errors are displayed in Table 2 and Table S10, respectively. Table 2 also displays the experimental values used. The best predictive methodologies are expected to have the lowest MAEs and mean signed errors (MSEs, taking their sign into account). When considering the training set of twenty phenols, MAEs between 1.34 and 1.61 (between 1.99 and 2.28) pK a units are obtained when using SMD (PCM). In general, direct calculations with Reaction (R3) led to underestimated (overestimated) pK a values for compounds more (less) acidic than phenol. The more acidic the phenol, the larger and more negative the error in these calculations. The calculated pK a values for the three nitrophenols and nitrosophenol, compounds 1-4, have the largest MAEs (e.g., between −3.36 and −5.02 pK a units with M06-2X(SMD)). These results indicate that the direct methodology tested is insufficient to predict nitrophenols.
When the calculated pK a values of compounds 1-4 are excluded, the MAEs become significantly reduced with values between 0. 50

Exploring Various Correlations with Experimental pK a Values and the Training Set
Correlations between calculated properties (descriptors), including ∆G • values related to the acid dissociation equilibrium, and experimental aqueous pK a values have been previously reported [12,[24][25][26][27][28]38,69,70]. When a correlation has a significantly large R 2 value and the MAE is very small, the fitted equation can be directly used to predict new aqueous pK a values of similar compounds.
Correlations between the experimental aqueous pK a values and the calculated ∆G • (aq) for Reactions (R1)-(R3) produced very high R 2 values (>0.90) in all but one case and much smaller and more consistent MAEs when the corresponding correlation equation is used at a given level of theory (see Table 1 Table 3 and Table S11, respectively, for the various levels of theory considered (only eleven phenols were calculated with the PBE0 functional). The MAEs for these correlations are shown in Table 1, and the associated correlation equations and R 2 values are listed in Table 4.    Table 4; d The calculated aqueous pK a values are reported in Table S11; e Using the pK a (exp) vs. ∆G The correlations using the SMD solvation model have lower MAEs and better R 2 values than the equivalent correlations using PCM, which is the same trend seen from the direct calculations previously discussed (see Table 4 It is important to note that when using these correlations (see Table 3), the calculated errors for the nitrophenols (1,3,4) and the nitrosophenol (2), which were very large when considering direct pK a calculations, are very small, in agreement with the calculated errors for the other compounds. This indicates that our correlations correctly adjust for the previously underestimated pK a predictions. Another observation is that the molecules with intramolecular hydrogen bonding affecting the most stable conformation of the acid form (1,4,5) tend to have a slightly higher error, as seen in Figure 5. We suspect that the additional stabilization in the acidic form may lead to an underestimated pK a prediction.

Predicting Aqueous pKa Values of Complex Phenols
A significant number of phenolic molecules that we are interested in studying (see Figures 2-4) are of greater structural complexity than the twenty molecules included in our training set. Hence, it is essential to test the performance of our correlation equations with more complex phenols. This is a largely underexplored area, partly because most

Correlations between the experimental aqueous pK a values and ∆G
• aq(BA) were reported for simple phenols by Galano et al. at several levels of theory using only the SMD solvation method [38]. Their training set of twenty simple phenols covered the experimental pK a range from 6.33 to 10.31. Using their reported correlation equations, we calculated slightly larger MAE (from 0.21 to 0.46) for the compounds in our training set (see Table 1). Apart from the fact that our correlation was built from these data, this difference is possibly due to the larger number of molecules more acidic than phenol in our training set, compared to theirs.

Predicting Aqueous pK a Values of Complex Phenols
A significant number of phenolic molecules that we are interested in studying (see  are of greater structural complexity than the twenty molecules included in our training set. Hence, it is essential to test the performance of our correlation equations with more complex phenols. This is a largely underexplored area, partly because most complex phenols lack experimental aqueous pK a values and no previous theoretical studies have verified the quality of aqueous pK a predictions for complex phenols using reliable experimental data.

Checking the Predictions with a Test Set
To check the accuracy of our correlations from the training set of twenty simple phenols, we collected ten phenolic molecules of varying complexity (displayed in Figure 2). Seven of them have experimental aqueous pK a values reported (21-24 and 27-29) [42,71,72], and three other ones (25,26, and 30) only have predicted values; 25 and 26 were predicted by the ACD/Laboratories Software [39], and 30 has a minimum experimental aqueous pK a value reported [73]. However, the same experimental methodology approximated the pK a of 29 within 0.1 pK a units [73]. While phenols 21-24 are simple, phenols 25-30 are of significant structural complexity, and 27-29 have experimental pK a values reported. Hence, for the first time to our knowledge, we will be assessing aqueous pK a predictions of complex phenols using correlations involving the experimental values of simpler ones.
, new pK a values should be calculated as pK a (calc) = m∆G • aq(BA) + n; b Using the 6-311++G(d,p) basis set; c The 7 phenols of the test set (which includes 3 large phenols) with experimental pK a values (see Table 5) have been added to the initial set of 20 phenols; d Calculated using 11 phenols.
Using the correlation equations (reported in Table 4) for the training set of twenty phenols employing the functionals M06-2X, B3LYP, BHandHLYP, and TPSS (with both the SMD and PCM solvation methods), the aqueous pK a values of the training set are calculated. Experimental pK a values and errors in their prediction, with MAE and MSE values, are shown in Table 5, while the calculated pK a values are reported in Table S12. MAE and MSE values are reported for the seven phenols with experimental data and for the entire set (including predicted values). Predictions using the correlation equations of Galano et al. [38] and employing the method suggested by Thapa and Schlegel [31] are also reported for comparison.
Unlike the trends seen in the MAE and MSE values previously reported (see, for example, Tables 2 and 3), the correlations using SMD usually exhibit larger errors in the predicted pK a values of the test set than when using PCM. The lowest MAE of 0.24 using the whole test set is achieved with M06-2X(PCM); however, when only experimental pK a values are used, the lowest MAE of 0.23 is produced by B3LYP(PCM). The MSE values are very similar between both solvation models and indicate that our correlations slightly underestimate the test set's pK a values.
Comparing our SMD results with Ref. [38], our values almost always produce smaller MAEs and MSEs. To compare with the predictive ability of the method reported in Ref. [31], a few molecules in our test set were selected (24, 28, and 29). In all cases (except for (R)-Trolox using B3LYP(SMD) and TPSS(SMD)), our correlations produce more accurate values. Moreover, the method reported in Ref. [31] is incompatible with 2-chlorophenol (and in general, with 2-substituted phenols) since the water molecules would not equilibrate around the -OH group in the acid species. Given that the pK a values of nitrophenols are difficult to predict directly, we included 2-nitrophenol in Table 5 (experimental pK a = 7.23) [42]. The pK a predictions of 6.13 and 3.92 by Refs. [31,38], respectively, had much larger errors than when using our M06-2X(PCM) correlation (7.34, see Table S12).
When considering only the complex phenols with experimental aqueous pK a values, the MAEs increased with SMD, but they generally decreased with PCM. Moreover, the MSE values increased for all functionals in both solvation models when only considering the complex phenols with experimental values. The MAEs ranged between 0.45 and 0.53 (0.27 and 0.35) with SMD (PCM). B3LYP(PCM), BHandHLYP(PCM), and TPSS(PCM) produce small MAE values for the complex phenols with experimental pK a values of 0.28, 0.29, and 0.27 pK a units, while their MSE values are −0.23, −0.23, and −0.22 pK a units, respectively. Compared to Ref. [38], our correlations produced significantly lower MAE and MSE values, especially with PCM. The performance disparity between the SMD and PCM solvent models can in part be attributed to the PCM exclusion of explicit non-electrostatic energy contributions. The excellent results with the test set, especially when using PCM, indicate that our correlations can be confidently applied to our prediction sets. Table 5. Predicted pK a errors and experimental pK a values for the phenols in the test set (21)(22)(23)(24)(25)(26)(27)(28)(29)(30) using the corresponding pK a (exp) vs. ∆G  To increase the statistical value of our work, we included the seven phenols with experimental aqueous pK a values from our test set into new correlated equations at all levels of theory. The best correlation graph obtained for the twenty-seven phenols is shown in Figure 5. The associated equations that will be used for predictions are listed in the second half of Table 4. While the MAEs of the new correlations with SMD slightly increased between 0.02 and 0.05 pK a units when the seven phenols were added, the MAEs with PCM decreased between 0.02 and 0.08 pK a units. The MAE range with twenty-seven phenols became 0.26-0.27 (0.26-0.32) using SMD (PCM), with the lowest MAE of 0.26 shared by M06-2X(PCM) and BHandHLYP(PCM). These are excellent results.

Predicting Aqueous pK a Values of Phenols with Potential Antioxidant Activity
Our group has studied the antioxidant properties of molecules 31-42, shown in Figure 3, and molecule 24, to repair damaged leucine residues under physiological conditions (pH 7.4). We require accurate pK a values to understand the biological mechanisms of action of these potential antioxidants. Additionally, we provide pK a predictions for molecules 25 and 26 which are opioid analgesics. Table 6 presents the aqueous pK a predictions using PCM, while additional predictions using SMD are displayed in Table S13. The predicted aqueous pK a values could be used as reference values for approximate pK a predictions in other computationally available implicit solvents, as previously reported [12]. This potential solvent transferability is highly useful as many of the species in the prediction sets have poor aqueous solubility.  [39], and while the predicted value for 27 was 10.09, the experimental value reported is 10.45 [71].
To our knowledge, there are no previous pK a predictions that we can compare our results with for molecules 31-42. However, the small average (median) range of 0.15 (0.12) for values predicted with M06-2X, B3LYP, BHandHLYP, and TPSS functionals using PCM is a promising result. Evidently, there is good agreement between these levels of theory.

Predicting Aqueous pK a Values of Cannabinoids
Finally, we used the same methodologies to predict the pK a for nine phenolic cannabinoids, shown in Figure 4. Of these, 43-46 and 49 present stereoisomerism. In all cases, the naturally occurring isomer was used (see Figure 4). In addition, two synthetic cannabinoids (47 and 48) were considered. Compound 46 contains only one stereocenter at C9, which gives a pair of enantiomers. For consistency, the R enantiomer was used in each case. In the case of 47, which is usually commercialized as a racemic mixture of (S,S)-(+) and (R,R)-(−) isomers, the latter was used in our calculations.
We could only find a reliable experimental pK a value for molecule 29 and an estimated minimum experimental value for 30 [71,73]. These cannabinoids were part of the test set previously discussed. Our average PCM prediction for 29 is 10.16 (with a spread of 0.09 pK a units, see Table 6), which is in good agreement with the reported experimental value of 10.60 [71]. The pK a predictions for 29 using the methodologies described in Refs. [31,38] were 9.23 and 9.69, respectively. Similarly, our average PCM pK a prediction for molecule 30 is 9.96, which is also in good agreement with the minimum estimated experimental value of 9.7 [73]. Using the methodology of Ref. [38], the pK a prediction of 9.55 for 30 seems to be slightly underestimated. Further, these results give us confidence in the accuracy of the methodology followed to calculate the macroscopic pK a value for molecule 50.
Again, due to the structural similarities of the cannabinoids, the average PCM (SMD) predicted pK a range for these molecules was between 9.21 and 10.31 (8.98 and 10.22). An interesting structural trend between molecules 46, 48, and 49, when compared to the other cannabinoids, is their lower expected pK a values (with values of 9.87, 9.26, and 9.21, respectively, if using average PCM predictions; see Table 6) because of the increased conjugate base stabilization from substituent conjugation. All of the cannabinoids were predicted to be similar or slightly more acidic than molecule 29, if considering its experimental value of 10.60; this is supported by molecule 43, an isomer of 29, having an average predicted pK a of 10.15 using PCM (just 0.01 pK a units from 29's average prediction).
The spread of the predicted pK a values for each of these systems using PCM (SMD) was 0.51 (0.69) pK a units or less, with the exception of 45 with PCM predictions spreading over 0.86 pK a units. Similar to what was reported in the previous section, the PCM average (median) range between the levels of theory for each molecule was 0.25 (0.17) which gives us confidence in our pK a predictions for these molecules. Furthermore, since the pK a values of cannabinoids were well reproduced in the test set, our prediction methodology could be extended to other molecules of this family.

Conclusions
Working with an initial training set of eleven structurally simple phenols, which was later expanded to twenty molecules, direct aqueous pK a calculations (using three acid dissociation equilibria) were perform with of five DFT functionals (M06-2X, B3LYP, BHandHLYP, PBE0, and TPSS), using the 6-311++G(d,p) basis set and the SMD and PCM solvent models. Much better and more consistent results were produced from the correlations between the calculated Gibbs energy difference between each acid and its conjugate base, ∆G , and the experimental aqueous pK a values, as previously reported [38]. The correlations using SMD (PCM) produced MAEs between 0.22 and 0.27 (0.28 and 0.40) and R 2 s between 0.947 and 0.975 (0.898 and 0.946). In general, the correlations using twenty phenols with SMD produced more accurate results than PCM.
A new set of ten phenols of varying complexities with experimental and/or predicted pK a values (separated accordingly) was used to test the performance of our correlations. In this case, PCM performed significantly better than SMD and the theoretical methodologies previously reported [31,38] for the entire test set and when the complex phenols were isolated. The best performance (for the set with experimental pK a values) was achieved by B3LYP(PCM) with an MAE (MSE) of 0.23 (−0.08) pK a units. The best performance for the complex phenols with experimental values were achieved by B3LYP(PCM), BHandHLYP(PCM), and TPSS(PCM) with MAE values of 0.28, 0.29, and 0.27 pK a units, respectively. These three functionals are expected to produce the most accurate pK a predictions when combined with the PCM solvent model; however, we have included the remaining levels of theory to form a range of predicted values. Furthermore, we developed new correlations, including the seven molecules from the training set (working with twentyseven phenols in total) to increase the statistical value of our work. The best MAE for the new correlations was shared by M06-2X(PCM), B3LYP(SMD), and BHandHLYP(PCM) with an MAE of 0.26 and R 2 s between 0.955 and 0.960 (see Table 4).
Our correlations were used to predict the pK a values of twelve molecules with potential antioxidant activity and of nine phenolic cannabinoids. The average prediction range with the PCM (SMD) solvation model was 0.15 (0.21) and 0.25 (0.34) pK a units, respectively, which indicates a very good agreement between our methodologies. These aqueous pK a predictions could be used as reference values for predictions in other solvents [12]. In the future, when more experimental data are available, it would be ideal to extend these correlations to a larger set of complex phenolic molecules to create an even better pK a predictive tool.
Supplementary Materials: The following information is available online at https://www.mdpi.com/ article/10.3390/antiox12071420/s1, The absolute Gibbs free energies of the different species considered in this study at the various levels of theory (Tables S1-S5); individual pK a values and their errors are displayed in Tables S6 and S7 (Tables S8 and S9) when using Reaction (R1) (Reaction (R2)); individual pK a errors when using Reaction (R3) (Table S10); correlated pK a values (Table S11); individual pK a values in test set predictions (Table S12); predicted pK a values for prediction sets not shown in Table 6 (Table  S13)