Testing of Exchange-Correlation Functionals of DFT for a Reliable Description of the Electron Density Distribution in Organic Molecules

Researchers carrying out calculations using the DFT method face the problem of the correct choice of the exchange-correlation functional to describe the quantities they are interested in. This article deals with benchmark calculations aimed at testing various exchange-correlation functionals in terms of a reliable description of the electron density distribution in molecules. For this purpose, 30 functionals representing all rungs of Jacob’s Ladder are selected and then the values of some QTAIM-based parameters are compared with their reference equivalents obtained at the CCSD/aug-cc-pVTZ level of theory. The presented results show that the DFT method undoubtedly has the greatest problems with a reliable description of the electron density distribution in multiple strongly polar bonds, such as C=O, and bonds associated with large electron charge delocalization. The performance of the tested functionals turned out to be unsystematic. Nevertheless, in terms of a reliable general description of QTAIM-based parameters, the M11, SVWN, BHHLYP, M06-HF, and, to a slightly lesser extent, also BLYP, B3LYP, and X3LYP functionals turned out to be the worst. It is alarming to find the most popular B3LYP functional in this group. On the other hand, in the case of the electron density at the bond critical point, being the most important QTAIM-based parameter, the M06-HF functional is especially discouraged due to the very poor description of the C=O bond. On the contrary, the VSXC, M06-L, SOGGA11-X, M06-2X, MN12-SX, and, to a slightly lesser extent, also TPSS, TPSSh, and B1B95 perform well in this respect. Particularly noteworthy is the overwhelming performance of double hybrids in terms of reliable values of bond delocalization indices. The results show that there is no clear improvement in the reliability of describing the electron density distribution with climbing Jacob’s Ladder, as top-ranked double hybrids are also, in some cases, able to produce poor values compared to CCSD.


Introduction
Due to the relatively low computational cost and generally good accuracy of the results obtained, density functional theory (DFT) [1][2][3] is currently the most frequently used method of describing electronic structure of molecules in computational chemistry [4]. The known problem of this method is the lack of the exact form of the so-called exchangecorrelation functional, as a result of which it is necessary to use its worse or better approximations [5,6]. Unfortunately, one could say that at present there are as many exchangecorrelation functionals as ants in an anthill, which leads to the situation that we currently have hundreds of them, and choosing the most appropriate one for a given problem is one of the biggest concerns when using DFT. Therefore, benchmark calculations in which the reliability of DFT functionals is tested are extremely important . However, the arises here is the correct choice of this reference structure. A good solution would perhaps be the use of the highest possible method to obtain the reference geometry. However, a limitation connected with computational possibilities will always force the compromise between the quality of calculations and the size of the system. Our goal is to test EDD in a real molecular system consisting of various atoms, connected via different types of bonds. Therefore, the mentioned limitation does not allow to obtain a computationally satisfactory reference geometry. Therefore, we naturally directed our attention towards experimental data. For this purpose, we use the good quality X-ray crystal structure of 2,2-dichloro-1-(1H-pyrrol-2-yl)ethan-1-one [57] (see the "Crystal Structure and Methodology" section for details) that we solved (see Figure 1). The 2,2-dichloro-1-(1H-pyrrol-2-yl)ethan-1-one molecule consists of various types of atoms connected by different bonds. These atoms are relatively small and thus give the opportunity to utilize computationally more advanced post-SCF calculations in order to obtain a high-quality reference EDD. Additionally, this molecule is planar due to bonding character and partial intramolecular π-electron delocalization. Therefore, it can be expected that its solid phase and gas phase geometries will not be significantly different due to possible packing effects. The use of geometry from the X-ray experiment has yet another justification. Currently the state-of-the-art X-ray analysis allows to experimentally estimate details of EDD in interatomic space; however, the standard procedure is that experimental results are confronted with theoretical EDD produced for geometry from the crystal (see, for instance, ref. [58] and citations therein). In that context, our analysis gains potentially additional important application, namely, an indication of the DFT functional which reproduces EDD most reliably.

Results and Discussion
As already mentioned in the Methodology, the electron density distribution (EDD) is probed by the most commonly used QTAIM parameters computed either at the bond critical point (BCP) or ring critical point (RCP): the electron density (ρ), Laplacian of the electron density (∇ 2 ρ), and total electronic energy density (H). Additionally, in the case of bonds, the delocalization index (δ) was also determined. The analysis of the impact of the theoretical method on EDD in BCP will be presented in the first subsection, and in RCP in the second. The impact of the method on the value of the bond delocalization index will be presented in the third subsection. Figure 2 shows the difference between the electron density value determined for the BCP (ρ BCP ) of a given bond within the given method and the reference CCSD method, i.e., ∆ρ BCP = ρ BCP (method) − ρ BCP (CCSD). Thus, a positive value for ∆ρ BCP indicates an overestimation, and a negative value for ∆ρ BCP indicates an underestimation relative to the reference value for CCSD. First, it should be noted that the obtained ∆ρ BCP values vary considerably depending on the bond under consideration. Importantly, the greatest discrepancies concern the double and highly polarized C1=O3 bond (see also Figure 1), for which the obtained ∆ρ BCP values are generally positive, indicating an overestimation of the ρ BCP value in relation to CCSD. By far, the largest error in relation to the CCSD value was obtained for the M06-HF functional (0.020 a.u.). Then, functionals BHHLYP (0.013 a.u.), B3LYP, X3LYP, and M11 (0.011 a.u.) followed. Quite the opposite, of all the considered functionals, MN12-SX performed the best, giving a negligible error value of only −0.001 a.u. A clearly greater error in the value of ρ C=O BCP , but also a negative one (−0.004 a.u.), was obtained by M06-L. It is worth mentioning that both of these functionals, i.e., MN12-SX and M06-L, were found to be among the best for general purposes [5]. Nevertheless, for all other functionals, the obtained ∆ρ C=O BCP values were greater than 0.005 a.u. showing that DFT has a particular problem with accurately describing EDD in this bond. The high sensitivity of the EDD of the double and strongly polarized C=O bond to the applied methodology was already noticed earlier [56]. Interestingly, the GGA (0.005-0.009 a.u.) and M-GGA (ca. ±0.005 a.u.) functionals generally performed better than the H-GGA (0.006-0.013 a.u) functionals (HM-GGA gave both the greatest and the lowest ∆ρ C=O BCP values) and double hybrids (0.009 ca.).
At the other extreme, i.e., with ∆ρ BCP < 0, the greatest error was obtained for the nonpolar C1-C4 bond, which as a result of electron density delocalization is intermediate between single and double ( Figure 1). In this case, the largest deviation was obtained for the local SVWN functional (−0.013 a.u.) and all GGA functionals (ca. −0.010 a.u.) except BLYP and B97-D, which gave slightly smaller errors (−0.007 and −0.006 a.u., respectively). For the other functionals, the description of EDD in this bond is clearly better. In particular, VSXC (−0.0028 a.u.), SOGGA11-X (−0.0027 a.u.), and mPW2-PLYP (−0.0025 a.u.) performed the best. Interestingly, BHHLYP as the only functional gave a positive error, which was even slightly smaller (0.0018 a.u.). Our results show well that density functionals may have a particular problem with the correct description of EDD not only in highly polarized double bonds but also in bonds that are formally nonpolar but associated with large electron delocalization. This result should be seen as a warning when describing organic compounds. However, the situation is somewhat similar in the case of the MP2 method, because, although in general, clearly smaller than for DFT, the highest values of ∆ρ BCP were also obtained for the bonds C1=O3 and C1-C4 (0.0038 and −0.0050 a.u., respectively). On the contrary, the difference to CCSD for MP3 and MP4SDQ is negligible.
Tognetti and Joubert reported [26] that the GGA and M-GGA functionals give too small ρ BCP values and that adding the exact exchange, that is the use of the H-GGA and HM-GGA functionals, remedies this defect somewhat. Our results confirmed this, but differences between the values obtained using the GGA or M-GGA functionals and their hybrid counterparts are rather small. Referring to a very recent study by Brémond et al. [50] it is worth paying special attention to the performance of M06-2X and M06-HF, because these HM-GGA functionals gave the smallest and the largest errors for the electron density at critical points, respectively. As can be seen in Figure 2, also, our results show that M06-2X, regardless of the bond, performs fairly well compared to the other functionals. Nevertheless, as our results show, other functionals, and in particular VSXC, M06-L, SOGGA11-X, and MN12-SX, appear to give better or comparable ρ BCP values, though perhaps not for all bond types. Both M06-L and MN12-SX were among the best general purpose functionals in Peverati and Truhlar's earlier studies [5], while the SOGGA11-X functional gave the best Fukui functions [38]. Additionally, as shown in Figure 2, the M-GGA functional VSXC also deserves attention. Moreover, TPSS, TPSSh, and B1B95 also give quite good ρ BCP values, although again perhaps not necessarily for all bond types ( Figure 2). All of these functionals were also previously mentioned as some of the best in the EDD description [26,38,39]. Therefore, our results are in line with the previous results and complement them nicely. On the other hand, when it comes to the poor description of ρ BCP by M06-HF [50], the extremely large error obtained for the strongly polarized C=O bond actually disqualifies this functional, although for some types of bonds (e.g., C7-C8 and C4-C7), it performed very well ( Figure 2). The differences in ∇ 2 ρ BCP values, i.e., Laplacian of the electron density, are shown in Figure 3. As the second derivative of the electron density, ∇ 2 ρ BCP is a quantity quite sensitive to the adopted methodology, e.g., the basis set used [56]. The current research also confirmed this. Namely, for many bonds, the differences in the obtained values of ∆∇ 2 ρ BCP are very large. Again, the C=O bond is clearly distinguished for which the values of ∆∇ 2 ρ BCP can be either significantly negative (e.g., ca. −0. The first three form the imidazole ring, while the last one is a linker to the -(CO)-CHCl 2 group. Thus, they all participate in a conjugated system with a high degree of electron charge delocalization. Together with the previously discussed large errors for the C=O and C-N bonds, this result again demonstrates the difficulty of describing EDD in bonds with large charge delocalization.
In contrast, the ∆∇ 2 ρ BCP values are significantly lower for C1-C2 and both C-Cl bonds, and therefore for the single bonds. While the local functional SVWN and all GGA functionals gave a significantly smaller error in the ∆∇ 2 ρ BCP value for C-Cl (ca. 0.000 a.u.), the trend is generally reversed in the case of hybrid and double hybrid functionals (and, in addition, ∆∇ 2 ρ C−Cl BCP is generally negative), but the differences for both bond types are small. Nevertheless, it can be seen that BHHLYP, M06-HF, and MN12-SX gave the greatest errors. Interestingly, along with B1B95, the BHHLYP functional gave the best ∆∇ 2 ρ BCP values in the earlier calculations by Tognetti and Joubert (however, only 10 functionals were tested then) [26]. Our results presented in Figure 4 clearly show, however, that while this functional indeed produced reasonable ∆∇ 2 ρ BCP values for many types of chemical bonds (e.g., C1-C4, C7-C8, C4-C7), the overall performance of BHHLYP is rather unimpressive.
The obtained ∆H BCP values are shown in Figure 4. Again, it can be seen that the obtained ∆H BCP values can be both positive and negative, not only depending on the bond but also for the selected bond (see the bright blue line for the C4-N5 bond or the red line for C1-O3 (i.e., C=O)). In the former case, the ∆H C4−N5 BCP values are in a wide range, from −0.033 a.u. for VSXC, BHHLYP, and SOGGA11-X, up to 0.038 a.u. for SVWN. However, by far, the largest error value was obtained by the M06-HF functional for the C=O bond (−0.124 a.u.). Much smaller, but also significant, ∆H C=O On the contrary, the largest positive error (up to 0.060 a.u.) was obtained within the SVWN functional for the C8-C9, C7-C8, C4-C7, and C1-C4 bonds. For these bonds, the VSXC (ca. 0.028 a.u.), BHHLYP (ca. 0.020 a.u.), and M06-HF (ca. 0.025 a.u.) functionals performed the best. The smallest error in the H BCP value was obtained for both C-Cl bonds, but in this respect, the VSXC and BHHLYP functionals performed the worst (both −0.011 a.u.).

Testing of DFT Functionals for the Electron Density Distribution in a Ring Critical Point
Another source of information on the reliability of the description of EDD by means of exchange-correlation functionals of the DFT method may be the values of ρ, ∇ 2 ρ, and H obtained in the RCP, which is possible due to the presence of the imidazole ring in the molecule under consideration. The resulting ∆ρ RCP values are shown in Figure 5. The result of a similar analysis but relating to ∆∇ 2 ρ RCP is shown in Figure 6. It is worth noting that in the case of ∆∇ 2 ρ RCP , the H-GGA and DH-GGA functionals behaved similarly, giving a fairly constant value, around 0.008 a.u., although deviations for SOGGA11-X and especially for MN15 are visible. Quite the contrary, in the case of GGA, M-GGA, and especially HM-GGA functionals, the differences obtained are significant.
The results for ∆H RCP are shown in Figure 7. As for ∆ρ RCP and ∆∇ 2 ρ RCP , the scatter of the obtained ∆H RCP values is large. The largest errors in relation to CCSD were obtained by M06-L (0.0071 a.u.) and VSXC (0.0069 a.u.). Slightly smaller ∆H RCP errors (0.005-0.006 a.u.) were obtained for BLYP, HCTH, B97-D, B3LYP, X3LYP, and M11. On the contrary, the smallest deviations were obtained by functionals B1B95 (0.0004 a.u.) and MN12-SX (−0.0006 a.u.), which are the only ones that clearly outperformed the MP2 method (0.0014 a.u.). Thus, these two functionals should be considered worth using for QTAIM-based π-electron delocalization analysis [59]. It is worth noting that in the case of H-GGA functionals, as in the case of ∆∇ 2 ρ RCP (Figure 6), the MN15 functional clearly stands out, giving a negative, not a positive, deviation from the CCSD value. Moreover, even without this functional, the values of ∆H RCP are clearly more different than in the case of ∆∇ 2 ρ RCP .

Testing of DFT Functionals for the Bond Delocalization Index
Bond delocalization index (δ) is one of the most important QTAIM-based parameters because, as a quantity describing the number of electrons shared by two atomic basins, it is directly related to the bond order [60][61][62]. Moreover, it is strongly related to the exchange energy and therefore describes the covalent component of a bond [63,64]. The delocalization index differences determined for each method are presented in Figure 8. The smallest errors ∆δ with a fairly constant value of ca. 0.120 a.u., regardless of the functional used (but double hybrids, of course, with ∆δ amounting to ca. 0.045 a.u. for B2-PLYP and mPW2-PLYP and 0.067 a.u. for PBE0-DH), were obtained for the single C1-C2 bond (brown line). This result shows that in terms of the reliability of the δ value obtained by DFT, a single nonpolar bond such as C-C is the least problematic.
From the results presented in Figure 8, it is clear that DH-GGA functionals significantly outperform the functionals from other rungs of the Jacob's Ladder in terms of the reliability of the δ value. Of these, it is rather difficult to clearly find the best, although perhaps VSXC, BHHLYP, SOGGA11-X, and M06-2X can be recommended. Quite the opposite, in order to determine reliable δ values, the use of SVWN, PBE, MN15, and M11 can certainly be discouraged.

Methodology
The aim of this article is to test the exchange-correlation functionals of DFT in terms of their reliability of the description of the electron density distribution (EDD) in a real molecule. As already mentioned in the Introduction section, the reference structure is the high-quality X-ray crystal structure of 2,2-dichloro-1-(1H-pyrrol-2-yl)ethan-1-one (see Figure 1). This molecule has various types of bonds: nonpolar (e.g., C1-C2 and C1-C4), highly polarized (e.g., C1=O3 and C2-Cl14), single (e.g., C2-Cl14), double (C1=O3), and intermediate (e.g., N5-C4). This makes it possible to carry out the analysis also in terms of the type of bonding. Moreover, the presence of a ring makes it possible to study the influence of the functional selection on the EDD in it. Using this experimentally determined geometry, wave functions used to describe EDDs were determined utilizing the Gaussian 09 [66] and Gaussian 16 [67] programs.
The EDD of the 2,2-dichloro-1-(1H-pyrrol-2-yl)ethan-1-one molecule was probed by the most popular and useful QTAIM parameters characterizing bonds and rings, namely, the electron density (ρ), its Laplacian (∇ 2 ρ), and the total electronic energy density (H)-all computed at either the bond or ring critical point (marked as BCP or RCP, respectively) [54]. Importantly, as noted by Brémond et al. [50], the electron density at critical points well represents the quality of EDD of the whole molecular system. Additionally, the bond delocalization index (δ) was computed as well. All these parameters were obtained using the AIMAll program [68].

Conclusions
The problematic issue when performing computations using the DFT method is the selection of the appropriate exchange-correlation functional. Earlier theoretical studies show that there is no universal functional. Some functionals are better for some purposes and other functionals are better for others. To make matters worse, it follows that a given functional possibly produces acceptably good values for one parameter, while the values obtained for other parameters may have large errors. The worst, however, is when these parameters belong to the same group of physicochemical data; for example, they are the lengths of different bonds or the values of different angles in the same molecule.
The aim of this article was to investigate the reliability of various exchange-correlation functionals of the DFT method in terms of describing the electron density distribution. For this purpose, calculations of several fundamental parameters derived from QTAIM were performed and then compared to the reference values obtained by means of the CCSD method. These calculations were made using 2,2-dichloro-1-(1H-pyrrol-2-yl)ethan-1-one taken from the X-ray crystal structure determined for the purposes of this project.
The presented results show that the DFT method has a particular problem in reliably describing the electron density distribution in multiple and highly polarized bonds such as C=O, as well as conjugated bonds associated with large charge delocalization. In the case of 2,2-dichloro-1-(1H-pyrrol-2-yl)ethan-1-one, these are the bonds belonging to the imidazole ring and formally a single C-C bond connecting this ring to the -(CO)-CHCl 2 group. On the contrary, description of the electron density distribution in single bonds proved to be much more reliable.
The obtained conclusions fit well with the current knowledge on the general behavior of exchange-correlation functionals. Namely, even for the considered QTAIM-based parameters, the functionals generally behave in an unsystematic way, and even if one gives good values for one parameter, it generally gives bad, or at least less reliable, values for another. This is quite a pessimistic side effect of using the DFT method. Moreover, this observation leads to the conclusion that nowadays it is much easier to discourage the use of certain functionals than to recommend some. The presented results show that the use of the M11 functional, in particular, should be discouraged for the purpose of a reliable description of the electron density distribution in molecules. In this respect, SVWN, BHHLYP, M06-HF, and, to a slightly lesser extent, BLYP, B3LYP, and X3LYP, also perform rather poorly. It should be emphasized here that in this group of nonrecommended exchange-correlation functionals there is B3LYP, which, however, is the most used. Therefore, its followers should seriously consider whether this functional really gives reliable results for the quantity they are concerned with, especially when these quantities are closely related to the electron density distribution in chemical bonds.
It is much more difficult to recommend a functional in terms of a reliable description of all the QTAIM-based parameters considered. The results show that none of the tested functionals stands out clearly in this respect. However, in the case of the electron density itself, being the most important QTAIM-based parameter, the VSXC, M06-L, SOGGA11-X, M06-2X, MN12-SX, and, to a slightly lesser extent, TPSS, TPSSh, and B1B95 can be recommended. It should be noted, however, that VSXC gave a very poor description of the electron density distribution in the ring critical point. On the other hand, the worst here was the M06-HF functional, which owes its worst position to the large error for the C=O bond. It was shown that functionals B1B95 and MN12-SX best describe the value of the total electronic energy density at the ring critical point and are therefore best suited for describing π-electron delocalization in an aromatic ring. Particularly noteworthy is the fact that double hybrid functionals clearly outperform other functionals in terms of a reliable value of the bond delocalization index, and, thus, the covalent component of a bond. It is noteworthy that all the functionals mentioned here turned out to be among the best in previous research [5,26,38,39,50].
It is also worth emphasizing that there is no clear overall improvement in the reliability of the description of the electron density distribution when climbing Jacob's Ladder, i.e., from local functionals via gradient and meta-gradient to hybrid and double hybrid ones. The lowest-ranked local functionals in this hierarchy do, in many cases, yield unreliable results, but they may also produce fairly good or good results in other cases. On the contrary, theoretically, the best double hybrid functionals, although giving some of the most reliable results for many parameters (especially the bond delocalization index), fail in the case of others (e.g., the total electronic energy density at the bond critical point of the C=O bond). This result confirms the older conclusion of Medvedev et al. from 2017 [36] that the performance of the latest exchange-correlation functionals in terms of the proper description of the electron density distribution is still not satisfactory.

Conflicts of Interest:
The authors declare no conflict of interest.