Calculating bond dissociation energies of X−H (X=C, N, O, S) bonds of aromatic systems via density functional theory: a detailed comparison of methods

In this study, the performance of 17 different density functional theory functionals was compared for the calculation of the bond dissociation energy (BDE) values of X−H (X=C, N, O, S) bonds of aromatic compounds. The effect of the size of the basis set (expansions of 6-31(G)) was also assessed for the initial geometry and zero-point energy calculations, followed by the single-point BDE calculations with different model chemistries with the 6-311 + (3df,2p) basis set. It was found that the size of the basis set for geometry optimization has a much smaller effect on the accuracy of BDE than the choice of functional for the following single-point calculations. The M06-2X, M05-2X and M08−HX functionals yielded highly accurate BDE values compared to experimental data (with the average mean unsigned error MUE = 1.2–1.5 kcal mol−1), performing better than any of the other functionals. The results suggest that geometry optimization may be performed with B3LYP functional and a small basis set, whereas the M06-2X, M05-2X and M08-HX density functionals with a suitably large basis set offer the best method for calculating BDEs of ArX−H (X=C, N, O, S) bonds.


Introduction
The Kohn-Sham density functional theory (DFT) offers one of the most widely used methods in computational and theoretical chemistry, with significant progress achieved in developing and validating highly accurate exchange and correlation functionals during the last decade [1][2][3][4]. The broad range of model chemistries offers various strengths and weaknesses, where higher accuracy usually means increasingly, often prohibitively high cost in terms of CPU time. As the development effort reached the final frontiers within the limitations of the method, DFT has increasingly turned into a routine tool for practical applications [5][6][7][8]. With a parallel increase in desktop computing power, the DFT method now advanced from the realm of specialists to become a mainstream, routine method for any chemist seeking fast validation for reaction pathways and energetics or predictions of various physico-chemical properties of compounds. That necessitates moving away from the 'higher is better' attitude towards the level of theory and a kind of cost-benefit analysis: for practical purposes, what level of theory is sufficient to address a particular problem [9].
Bond dissociation energy (BDE) is defined as the enthalpy change of the homolytic cleavage of a specific bond, providing one of the essential characteristics of the breakdown reactivities of compounds in the context of e.g. their antioxidant activity [10,11]. Experimental methods to determine BDE are limited to measuring energy flows, such as calorimetry, and rely on breaking a known weak bond. Thus, experimental measurements of BDE are limited to simple and/or small molecules, hence they are unsuitable in many real-world situations where large and complicated molecules are involved, which is often the case with natural products. The recent advancement of first-principlebased computational chemistry techniques has provided an alternate method for obtaining precise BDE for any molecule of interest [12]. However, only at highly connected post Hartree-Fock levels of theory have ab initio approaches proven accurate among the known computational chemistry methods. The G2 and CCSD(T) approaches, for example, have been shown to be accurate but only for tiny compounds [13][14][15][16][17][18][19]. It is not feasible to perform all such calculations on supercomputers, and the above-stated advantages of DFT justify the less accurate but more practical method, especially for predicting the characteristics of very large molecules [5,10,11,19,20].
Many previous studies demonstrated that the radical scavenging activity of aromatic compounds such as phenolics [23][24][25][26][27][28][29], aromatic amines [30][31][32] or thiols [33,34], particularly in lipid media, follows the hydrogen transfer mechanism characterized by BDE(X-H, X=C, N, O, S) values. While several studies addressed the optimization of model chemistries for calculating specific physicochemical parameters of specific compound families [35][36][37][38][39][40][41][42][43][44][45], these works also highlight that the performance of model chemistries cannot be evaluated in the absolute and has to be repeated for each practical problem sufficiently different from the ones addressed before. One such problem is the computation of BDE(X-H, X=C, N, O, S) of compounds containing aromatic rings. Thus, in this study, the 17 DFT functionals identified above have been evaluated for computing BDE values of the X−H (X=C, N, O, S) bond of ArX−H compounds to find the most convenient method for computing this thermochemical parameter for this specific case.

Computational details
The Gaussian 16 package of programs was used for all calculations in our work. BDE calculations were performed according to the protocol of Wright [11] that is widely used in the literature [10, [18][19][20]46]. In royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 220177 this study, the geometry is optimized, and zero-point energy corrections are calculated with the same model chemistry for all compounds and for all of the following by single-point calculations of BDE values with a range of model chemistries. The calculated BDEs were compared with the experimental BDE values obtained from the literature [47].
BDE is the total energy required for homolytic breakage of a specific bond, here the X−H of the ArX−H compound. BDE is thus the enthalpy change of the reaction For the compounds that have multiple conformers, all of these were screened by the Spartan Software with the MMFF/molecular mechanics module [48] and the conformer with the lowest electronic energy was included in the analysis. ArX−H (X=C, N, O, S) substances were chosen as model compounds for performance testing of the 17 density functionals on the condition that experimental values of BDEs exist for them.

Effect of the size of the basis set
Previous studies showed that the B3LYP functional is sufficient and accurate for geometry optimization and frequency calculations [20,21,[49][50][51][52]. Thus, calculations were performed with B3LYP functional and five typical basis sets including 6-31G(d), 6-31 + G(d), 6-31 + G(d,p), 6-311G(d,p) and 6-311 ++ G(d,p)) for the geometry optimization and frequency calculations, followed by single-point calculations of BDEs using 17 different functionals and the 6-311 + G(3df,2p) level [21]. This is a protocol intended for modest computational resources, where the larger basis set is only used where it is necessary; the aim here is to validate this approach against experimental data.
In the first section, we assess the variation affected by the size of the basis sets by cross-comparison of the data. Arguably any improvement in the numerical values should manifest as a difference between the values calculated with different basis sets, and therefore the analysis does not require referencing experimental data to reveal the ideal size of the basis set. The detailed data including differences in numerical values between basis sets in comparison to the smallest 6-31G(d) level are presented in the electronic supplementary material, table S3-S6 for C−H, S8, S10 and S12 for N−H, O−H and S−H bonds, respectively. Given the very small variation observed for C−H bonds, for N−H, O−H and S−H bonds only the smallest and largest studied basis sets were compared.
As per these data, variations of the BDE values based on geometry optimization with different size basis sets are mostly lower than 0.3 kcal mol −1 (corresponding to 0.3%). The only exception was observed for BDE calculated with B2PLYP functional at the N−H bonds (electronic supplementary material, table S8), and the S−H bond of compounds SH1 and SH3 (electronic supplementary material, table S12). The larger error here is likely related to the chemical environment of these bonds that requires more accurate geometry optimization and is beyond the scope of this work to analyse in detail. Overall it was found that the geometry optimization with larger basis sets does not yield significant differences compared to smaller basis sets; hence even without comparison to experimental data, it is clear that they cannot deliver much improvement of accuracy (range from 0.1% to 0.3%). Importantly the data does not suggest, per se, that the size of the basis set in geometry optimization is irrelevant for BDE calculations, rather than the larger basis set is not expected to deliver substantial improvement, conversely for very large molecules reducing the size of the basis set is a viable compromise to obtain sensible predictions of BDE. As the few specific cases mentioned above imply, for higher accuracy, it is still desirable to perform calculations with a larger basis set. Therefore, in the next section, we focus our analysis on the geometry optimized with the smallest basis set, i.e. the DFT/6-311 + G(3df,2p)//B3LYP/6-31G(d) chemistry model, to evaluate the accuracy of different functionals in the calculation of BDE values of the studied compounds against experimental data.
As shown in table 1, the best-performing functionals for calculating BDE values of C−H substances are M08-HX, M06-2X, BMK and M05-2X, in that order. Remarkably, the mean unsigned errors (MUE) of M06-2X and BMK are equal to M08-HX (1.3 kcal mol −1 ). The maximum absolute error (MaxAE) of M08-HX has its minimum at 3.0 kcal mol −1 . The MaxAE values of BMK and M05-2X are higher than those of the M08-HX and M06HX by about 1.0−1.6 kcal mol −1 . A similar trend was observed at absolute values; the relative MUE of the four listed functionals above are not significantly different (range from 1.4% to 1.6%), whereas the relative MaxAE of BMK and M05-2X (4.8% and 5.5%, respectively) are much higher than that of M08-HX and M06-2X (3.9% and 3.8%, respectively). The rest of the tested functionals for the C−H bond BDE calculations also have the MaxAE lower than 10.0 kcal mol −1 except B2PLYP (12.0 kcal mol −1 ).
These results correlate well with reports on comparing the accuracy of various functionals for energetic calculations on specific classes of compounds. Yan Zhao et al. [21] tested nine functionals against three experimental databases of Izgorodina et al. [62], and found that M06-2X and M05-2X showed excellent performance for calculating energetics involving radicals as well as the alkyl oxygen BDEs of nitroxides. Similar results were also reported in other studies [38][39][40][41][42], once again asserting the superior performance of M06-2X and M05-2X over any other density functionals.

Conclusion
In this study, the BDE(X−H, X=C, N, O, S) values of typical aromatic compounds have been computed and compared using 17 different DFT functionals, based on the protocol by Wright [11]. The results showed that M06-2X, M05-2X and M08-HX yield higher accuracy BDE values than any other tested functionals. The 6-31G(d) basis set with the B3LYP functional is sufficient for geometry optimization, with acceptable absolute and relative value deviation, whereas the larger basis sets 6-31 + G(d), 6-31 + G(d,p), 6-311G(d,p) and 6-311 ++ G(d,p) did not deliver a significant difference. Hence, it is possible to obtain reasonable predictions of BDE of very large molecules based on geometry optimization with a small basis set without substantial loss of accuracy. It should be emphasized that this does not apply to the single-point BDE calculation, where a large basis set was used for accurate results. When compared to experimental data, the M06-2X, M05-2X and M08-HX density functionals were found to offer the most affordable performance for calculating the BDEs of ArX−H (X=C, N, O, S) bonds.