Better performance of Hartree–Fock over DFT: a quantum mechanical investigation on pyridinium benzimidazolate types of zwitterions in the light of localization/delocalization issues

Context With the advent of fast computing facilities, combined with rapid emerges of many new and intricate quantum mechanical functionals, computations with pure Hartree–Fock (HF) theory are now-a-days regarded as trivial or obsolete, or even considered as not reliable by many researchers. Consequently, current trends in computational chemistry show extensive use of post-HF theories for smaller molecular systems and various DFT methods for organic and inorganic chemistry related problems (larger molecules/systems). In this contribution, I have tried to show that sometimes, HF might be more suitable over DFT methodologies in addressing structure–property correlations. Molecules studied here were previously synthesized by Boyd in 1966 and important experimental data were produced by Alcalde and co-workers in 1987. Comparison of computed and experimental results clearly shows that HF method was more effective in reproducing the experimental data compared to especially the DFT methodologies. Reliability of HF method was further assured from the very similar results shown by the CCSD, CASSCF, CISD and QCISD methods. Current study also indicates that the localization issue associated with HF proved to be advantageous over delocalization issue of DFT based methodologies, in correctly describing the structure–property correlation for zwitterion systems. Methods All computations were performed with Gaussian 09. A wide-range of quantum mechanical methodologies, HF, B3LYP, CAM-B3LYP, BMK, B3PW91, TPSSh, LC-ωPBE, M06-2X, M06-HF, ωB97xD, MP2, CASSCF, CCSD, QCISD, CISD and semi-empirical methods like, Huckel, CNDO, AM1, PM3MM and PM6, were used for investigations. Graphical abstract


Introduction
History of 'Quantum Mechanics' shows that, Hartree-Fock theory (HF) refers to the 1927's original work of D. R. Hartree (published in January 1928).[1] Hartree's work was presented to scientific community immediately after the 1926's profound work of Erwin Schrödinger [2] (famously known as Schrödinger Equation, iℏ d dt �Ψ(t)⟩ = � H�Ψ(t)⟩ ).[3] In his work Hartree offered a systematic procedure, famously known as Self-Consistent Field (SCF) method (also recognised as Hartree's Method), where he presented the solutions to determine energy and wave functions (approximate), for elementary quantum mechanical systems.[1] In 1928, J. C. Slater [4] and J. A. Gaunt [5] (both simultaneously and independently), demonstrated that the variational principle to a trial wave-function (an ansatz) can be considered as an appropriate theoretical basis for the SCF method.Then, in and around 1930, independently J. C. Slater [6] and V. A. Fock, [7,8] applied the anti-symmetry to the electronic solutions.With the use of Slater determinant (of one-particle wavefunctions), [9] which was known to have the essential property of anti-symmetry, ansatz was made suitable to be used in the variational principle (to be noted, earlier in 1926, Heisenberg and Dirac [10] used the principle of anti-symmetry).Considering all these parallel developments, Hartree reformulated his original theory (with the inclusion of Born-Oppenheimer approximation [11]), and made it more suitable for solving time-independent Schrodinger equation, applicable to real physical systems.Later this revised theory gained the popularity as the "Hartree-Fock Method" (HF Method) [3,12,13].
Then it was realized that the non-inclusion of electron correlations in HF Method is one major setback and was needed to be addressed properly.This lead to the development of many well-ordered approaches, which were able to appropriately address the correlation issue associated with HF method.Approaches like, Moller-Plesset Perturbations theory (with 2 ND order consideration, frequently used one is the MP2 method) [14], Configuration Interaction (CI) [15], Multi-configuration Self-Consistent Field (MC-SCF) [16], Complete Active Space Self-Consistent Field (CASSCF) [17], Coupled Cluster (CC) [18][19][20], are the few notable theoretical developments.All these well-established theories are collectively known as the post-HF theories [21,22].Worth to mention here that, the HF theory was not only the bedrock for these post-HF theories [21,22], but also contributed significantly for the developments of the well-known DFT (density functional theory) theory [23,24].Many DFT based methodologies are now being used extensively by almost all the researchers worldwide, for various areas of scientific research [25][26][27][28][29][30][31].
When computational capabilities were still in their emerging stage (around in the 70's and 80's of the 20 TH century), the long well-established HF theory, was in the forefront of computational chemistry research.At present, looking at the trends in the use of methodologies in computational chemistry research, one can see the general tendencies for using post-HF theories (mainly for relatively smaller molecular systems) and DFT (for medium to larger systems) [21,22,28,[30][31][32][33][34].While the post-HF theories have the advantage of being superior to basic HF approach, in producing accurate results (many times very close to experiment), at the same time, being computationally expensive they are generally not favoured for medium to large systems.About the DFT based methodologies, current trends clearly show that they are the most widely used computational tools by all the computational chemists (being applied to small to medium to large systems, or to say all kinds of systems).With wide varieties of functionals being available, no wonder, application of DFT based methodologies became the main choice for the computational chemists in the present era of science [21,[32][33][34].
About the strong dominance of DFT in the field of organic chemistry, one can easily see that besides many profound theoretical articles, now-a-days even most of the experimental articles in organic chemistry are most often augmented with a substantial portion of DFT computations [21,28,30,33].About the HF method, now a days it is generally not the method of choice for many researchers, even in the field of computational organic chemistry.In this contribution, I investigated the structure-property correlation of some wellknown zwitterionic organic molecules, which were synthesized by Boyd, long-time back in 1966 (Scheme 1) [35].Main motivation of this study is to show that sometimes one may be able to achieve equally good or even better performances by HF method compared to DFT-based methodologies, and thus HF method can still be treated as suitable method of choice for computational organic chemistry.Then in 1987 (after 20 years), Alcalde et al. [36], resynthesized these Boyd's zwitterions (pyridinium benzimidazolates), and investigated their crystal structures as well as one of the fundamental molecular properties, the dipole moments.Based on its relatively large dipole moment value, Molecule 1 (Scheme 1) was subjected to some interesting studies in the later years, mainly for the properties directly/indirectly linked to its dipole moment as well as its interesting charge transfer property.
For the Molecule 1, around 10 years later, Abe et al. [37] investigated it for nonlinear optics, and compared their computational results from HF method, with the experimental dipole moment value (10.33D) reported by Alcalde et al. [36].They found that HF method was able to reproduce (also co-authored with J. Abe) [38], on the same molecule, they mentioned about the better agreement between the HF dipole moment with the experiment compared to DFT, but they were a bit unconvinced on the performance of HF method.Hence, the question arises "Are we undervaluing the results produced by HF method by not considering it seriously?"or "by effectively reproducing the experimental results for zwitterionic systems, is it showing cues for its usefulness?"Hence to find an answer to the posed questions, I carried out a detailed investigation of the some already synthesized zwitterionic molecules using a wide range of quantum mechanical methods, ranging from classical semiempirical methods, to DFT to post-HF methods and along with the HF methods.Not only the dipole moment, but also several other fundamental properties were computed to assess the performances of various methodologies.

Molecular structures
The fully optimized geometries of the Pyridinium Benzimidazolate zwitterion (Molecule 1) computed using all the methods mentioned in Sect."Computational methods", were analysed and important structural parameters are shown in Table 1.These computed structural data were compared with the available experimental crystal structure data (provided in Scheme 1) to assess the effectiveness of various quantum mechanical methods in reproducing the experimental structure of Molecule 1.For Molecule 1, one of the important structural parameters is the twist angle between the two aryl units (donor-acceptor junction: D 2,1,7,8 ).Experimental crystal structure predicted the molecule to be fully planar, with a 0.0° twist angle [36].All the computational methods also predicted it to be fully planar like the experiment, except the MP2 method.MP2 computations predicted it to be 13.7° twisted at the junction, quite unusual behavior compared to all other methodologies.To eliminate any possibility of insufficiency of basis set consideration for MP2, I carried out full optimization with larger basis sets.Computations with aug-cc-pVTZ failed due to the insufficient computing infrastructure available in my lab.Then computation was carried out using aug-cc-pVDZ basis set.
Interestingly, the obtained optimized structure was found to be almost fully planar like other methods (other structural parameters obtained using this large basis set were found to be same or close to the values obtained using 6-31 + + G(d,p) basis set).A reduction of twist angle from 13.7° to 0.5°, with the change of basis set from 6-31 + + G(d,p) to aug-cc-pVDZ clearly indicates that, with cc-pVTZ basis set it will possibly show a fully planar conformation as predicted by all other methods.Thus, one can say that while dealing with biaryl types of zwitterionic molecules which are susceptible for internal rotations at the inter-ring junctions, with MP2 methodology a large basis set is essential to address the structures more accurately [58].As other methods (even the well-known HF theory) were found to be more suitable in addressing the structural aspects, even in the lower basis set domains, hence MP2 method is not advisable for similar zwitterionic systems for researchers, who are limited with computing resources.
Like the inter-ring twist angle, another important structural parameter is the junction bond (R 1,7 : Scheme 1) and the reported experimental value is 1.45 Å. Analysis of the computed values from various methodologies (Table 1) it can be observed, while all the DFT methods predicted underestimated values (except LC-ωPBE) compared to experiment [36], the post-SCF (CASSCF, CCSD, QCISD) methods (except MP2) along with the HF methods (with different basis sets) predicted values very much closer (same cases exact values) to the experiment [36].Interesting to note here that almost exact matching values with the experiment were shown by HF methods with larger basis sets.This is again a clear indication that, even simple HF method is quite suitable in reproducing the experimental parameters.Other structural parameters shown in

Dipole moments
Density functional approaches are well-known to be predicting improved energies and are generally believed that they can inevitably be more accurate in predicting electron densities compared to Hartree-Fock method.This belief has been questioned by researchers in the recent times [59][60][61][62][63][64][65][66][67][68].
It is well-known that one of the simple molecular properties is the dipole moment of any molecule, and often it is used to precisely understand the electron density distributions of any polar molecule.Hence, to assess how various wellknown quantum mechanical methods predict the molecular dipole moment (indirect inferences can be obtained about the performances of various methodologies on predictions of appropriateness of the distributions of electron densities) is the main intention of this analysis [59].To test this aspect for the first time, for a zwitterionic molecule, both the crystal structure data and the experimental dipole moment were available for one such molecule (Molecule 1, Scheme 1).
Another molecule for which both the crystal structure and dipole moment were also reported in the same work but owing to its exceedingly large size (multiple phenyl and alkyl substituted derivative of the reference molecule), such an investigation (with all the methods investigated in this work) will be impossible for us (with our limited computing resources) [36].
Computed dipole moment data for the Molecule 1, obtained from various methodologies are shown in Table 2.As can be seen from Table 2, when compared with the experimental dipole moment of 10.33 D (in the same report, Alcalde et al. also computed the dipole moment with semiempirical method and reported the value as 11.06 D [36]), it can be seen that, while the post-HF and HF methods predicted the dipole moment values closer to the experimental values, at the same time the predicted values by the hybrid DFT and also the semi-empirical methods were found to be deviated strongly from the experimental value.To be noted that most of the long-range corrected DFT methods predicted the dipole moment values closer to experimental value.Worth to mention here that, while HF and post-SCF methods predicted either nearly equivalent or slightly overestimated values compared to the experiment, at the same time the long-range corrected methods predicted slightly underestimated values (even the other DFT and semi-empirical methods predicted more underestimated deviations from the experimental dipole moment).
Dipole moment of 10.31 D predicted by simple HF/6-31G was found to be almost same as reported in earlier experimental value of 10.33 D [36].With other basis sets, the HF method predicted slightly larger values compared to experiment.At the same time a value of 10.34 D predicted by post-HF method CISD with 6-31G basis set can be regarded as exact match with the experiment.Other post-SCF methods predicted the dipole moment values slightly higher than experiment.All the long-range corrected DFT methods, predicted the dipole moments comparatively closer to experimental value compared to all the hybrid DFT methods, but still larger differences compared to experiment were observed.This gives a clear indication that HF method is quite suitable in predicting not only the structure and but also dipole moment of the zwitterionic molecule (Molecule 1) more efficiently than all the DFT based methods.Hence, at this stage one can say that even the simple HF method can be considered as suitable in accounting the electron density distribution around of the zwitterionic Molecule 1 and can also be expected to efficient in the cases of other similar zwitterionic systems (definitely need more investigations to validate this claim).
To test the observed suitability of HF, we looked in the same work of Alcalde et al. [36], we found three more molecules for which they reported the experimental dipole moments, hence they can be ideal candidates for the test.Out of three molecules, two molecules are extremely large (with multiple phenyl substitutions) and hence it will not be possible for us to do the line of investigations carried out with those two molecules.Luckily one molecule was of reasonable size and a substituted derivative of Molecule 1.The methyl substituted derivative of Molecule 1 (pyridinium benzimidazolate), which is Molecule 2 is shown in Table 3 (geometry of molecule is shown in its optimized orientation).Computed dipole moments from various methodologies are also shown in Table 3.
Though experimental dipole moment for Molecule 2 was reported, but Alcalde et al. [36], did not report the crystal structure data for it.But, in view of similarity in structure and moreover as the two methyl substituents are remotely located from the bridge, hence we can expect Molecule 2 be structurally very close to Molecule 1.For all the investigated methodologies, we found the Molecule 2 to be in fully planar conformations (even in MP2 with smaller basis sets).Now about the central junction bond, we observed slightly different values in various methodologies compared to Molecule 1 (it was 1.45 Å for Molecule 1).While HF and post-HF methodologies predicted it to be closer to Now about the dipole moments, unlike Molecule 1, here for Molecule 2 surprisingly we observed that in all the methodologies predicted underestimated values of dipole moments compared to experimental value of 13.52 D. B3LYP predicted a value of 6.85 D which is almost half of the experimental value.But with the long-range corrected methodologies, like CAM-B3LYP and ωB97xD, we observed improved values than the B3LYP.Among all the methodologies, HF predicted the largest dipole moment (around 10.1 D).Slightly lower values of observed in all the post-HF methodologies, compared to the HF predicted dipole.These marginal differences may be due to the lower level of basis sets employed for most of the post-HF methodologies.Despite the underestimated values in all methodologies, once again HF was proved itself to be more suitable by predicting close to the experiment the experimental value (only a difference of 3.4 D was observed).Based on the similarity of trends shown by both Molecules 1 and 2, with respect to the performances, we can say that among various methodologies, simple HF methodology is quite useful in predicting structure-property correlations closer to experiment.Though at this stage exclusively we can't advocate the persistence of such a trend for all kind of zwitterions but based on the good performances shown by HF methodology and unusual worst performance by B3LYP for the current zwitterions, we can say that one need to be careful while computationally dealing with zwitterions.

Frontier molecular orbital energies
Besides the observations of the dipole moment values, one of the key molecular properties which is going to be directly affected by the variances in distributions of electron densities (shown by various methodologies) is the energies related to frontier molecular orbitals.Hence frontier molecular energetics data for the Molecule 1, computed using various methodologies were tabulated and are shown in Table 4.As expected, large differences were observed in the energetics of HOMO (highest occupied molecular orbital: E HOMO ), LUMO (lowest unoccupied molecular orbital: E LUMO ) and the ΔE HLG (energy difference between the two: E LUMO -E LUMO ).
From Table 4, it can be seen that the lowest values of ΔE HLG were obtained for the hybrid DFT based methods and larger HOMO-LUMO gaps were observed for the HF method, with nearly similar values for the post-HF methods.At the same time all the long-range corrected DFT methods showed larger HOMO-LUMO gaps compared to hybrid DFT methods and lower gaps compared to HF as well as post-HF methods (except M06-HF which predicted it to be like the hybrid DFT methods).As orbitals play primary roles in charge transfer and excitation properties of a molecule, hence such a situation can be expected to affect many related properties.Based on the near similar values of HF and post-SCF methods, and at the same time approaching values for larger HLG values for the long-range corrected methods, one can say that the hybrid DFT methods may not be properly accounting the frontier molecular orbital energetics for this molecule.
Analysis of the trends in the stabilities of the frontier molecular orbitals indicate that while the HF and post-SCF methods predicted the HOMO to be stabilized, at the same time the DFT based methods predicted it to be less stabilized (whereas LC-ωPBE and ωB97xD methods stabilized HOMO like that of the HF method).Interesting observations were found for the energetics of LUMO.While all the post-HF and HF methods predicted a destabilized LUMO (positive energies), all the DFT based methods predicted it to be highly stabilized with negative energy values.Such a situation can be attributed as being largely responsible for the observed molecular band gap (of HOMO-LUMO gap) for different methodologies.This gives a clear indication that there might be some intrinsic problems associated with the methodologies, which ultimately affects the electron density distributions around the molecule, and consequently the orbital energies as well as the dipole moment values are also getting strongly influenced.Many recent works proposes that such possibilities may be analysed in the light of the localization-delocalization problem [59-61, 65, 66].

Localization and delocalization problems
As discussed previously, the question arises here is that the observed differences in properties can exclusively attributed only to the difference in geometries or some other factors might be affecting these behaviours shown by HF and DFT methods.To test this, I took the optimized geometry of the B3LYP/6-311 + + G(d,p) and computed the dipole moment at HF/6-311 + + G(d,p) level, without doing any further optimizations.Also, to avoid any basis set related adverse effects, I investigated the above situation with a 6-311 + + G(d,p) basis set.Interestingly, the observed dipole moment was found to be 10.85 D, larger than the values obtained from B3LYP/6-311 + + G(d,p) full optimization, and slightly larger than HF/6-311 + + G(d,p) full optimization.This is clear indication that the geometry may not have significant influence on the observed dipole moment.Also, influence of the geometry only, on the dipole moment can be ruled out based on the MP2 results, where it showed a completely different geometry (with a twisted conformation around the aryl-aryl junction) compared to all other methods, but a dipole moment was found to be close to experimental value.Interestingly according to S. Nam et al., [69], "Problematic calculations are density-sensitive and using HF densities fixes these issues."Not only that, but also all the other post-HF approaches produced very similar geometry and dipole moment values like that the traditional HF approach.
Based on all these observations, I tried to move our focus on the nature of problems associated with the HF and DFT.Computed electrostatic potential (ESP) map, frontier molecule orbitals (HOMO & LUMO) and canonical forms of the reference molecule (Fig. 1) are used in the analysis of the localization and delocalization issues associated with the HF and DFT methodologies.One of the problems described in many earlier works [59][60][61][62][63][64][65][66][67][68][69][70][71][72][73][74][75], and recently highlighted in a perspective article by H. J. Kulik, is the localization/ delocalization problem associated with the HF and DFT methodologies.H. J. Kulik stated, "Degree of localization or delocalization is a problem in both approximate density function theory and in Hartree-Fock."[70].The perspective also highlights the pernicious nature of this problem, like, how sometimes over-delocalization (DFT) can turn an insulator to metal and over-localization (HF) can produce wrong results in barrier height estimations or energetics of a reaction [70].In the present context, the reference molecule being zwitterionic nature, this problem looks more probable, as a slightest over-or under-estimation of the localization/ delocalization can strongly affect the structure-property correlations.A delocalization may impact some stabilization to a particular canonical form (Fig. 1) and may ultimately affect the frontier molecule orbital energies.Figure 1 shows that both HOMO & LUMO are exhibiting localized natures of population density distributions, hence any delocalization (representative canonical form) may affect the density distributions as well as the energetics.Also, the ESP map clearly indicates centralized regions of positive and negative potentials in the molecule (Fig. 1) and any delocalization can directly affect not only the potential distributions of the consolidated regions, but also many other directly related properties of the molecule.
Associated to delocalization problem, self-interaction error associated with DFT is discussed in a recent work of Lee et al. [71].They discussed about the difference in properties predicted by DFT and HF method for anions.They even suggested that HF theory can be a naive solution to address the problem, shown by DFT methods in predicting the frontier molecular orbital energies and other related properties of anionic systems [71].As the systems considered here are zwitterionic natures (with the donor part anionic in nature), hence underperformance of DFT compared to HF may also be addressed as explained by Lee et al. [71].In the present case I have also observed very similar problems shown by the DFT methods in predicting the frontier molecular orbital energies of the reference zwitterionic molecule (Table 4 and related discussions in Sect."Frontier molecular orbital energies").As discussed in the beginning of this section, the orbital energies E HOMO = -5.51eV and E LUMO = -3.00eV shown by the B3LYP/6-311 + + G(d,p), got changed to -6.93 eV and + 0.14 eV respectively when computed at HF/6-311 + + G(d,p) level, without doing any further optimization of the geometry obtained at B3LYP/6-311 + + G(d,p) method (these results are similar and close to the results obtained from the HF/6-311 + + G(d,p) full optimizations).This clearly indicates that the localization/delocalization problem associated with the HF and DFT methods plays very significant role.In the case of DFT methods, it was observed that the associated delocalization issue gave a stabilized LUMO (negative energy) and a slightly destabilized HOMO, in all the cases compared to HF method.Thus, the localization issue associated with the HF method was found to be in an advantageous position to accurately address the orbital energetics.This inference can be drawn based on the near equivalent results produced by the post-HF theories.Also, to mention here that DFT methods incorporating the long-range interaction factors and dispersion corrections were found to be in better positions in producing the orbital energetics (still the LUMO energies were found to be negative) and consequently structure and dipole moments were also found to be close to the HF results (LC-ωPBE was found to be performing better than all other methods).Other inference can be drawn is that localization issue associated with HF method showing a preference towards the more zwitterionic behaviour of the molecule, and at the same time the delocalization issue associated with the DFT based methods showing the preferences for a more quinonoid type behaviour for the reference molecule.

Conclusions
This contribution reports the computational results for Boyd's zwitterions (pyridinium benzimidazolates) and compares them with the experimental structural data and properties already reported in the literature.Computational investigations were carried out using various methodologies like, HF, B3LYP, CAM-B3LYP, BMK, B3PW91, TPSSh, LC-ωPBE, M06-2X, M06-HF, ωB97xD, MP2, CASSCF, CCSD, QCISD, CISD, Huckel, CNDO, AM1, PM3MM and PM6, to evaluate the performances of all these well-known computational methods, against the available experimental structural and dipole moment data.A comparison with the experimental data clearly indicates that HF and post-HF methodologies were able to reproduce the experimental data (both the structure and dipole moments).At the same time Hybrid DFT methodologies showed substantial deviations from the experimental data.But, with the inclusion of long-range interaction corrections and/or dispersion corrections in DFT, some improved performances in reproducing the experimental data were observed.Better performances of HF and post-HF methodologies compared to DFT, were explained in the light of localization/delocalization problem associated with the HF and DFT.It was observed that localization issue associated with HF preferring more zwitterionic, and the delocalization issue associated with the DFT based methods, showing the preferences for a more quinonoid type behaviour for the pyridinium benzimidazolate type of zwitterions.These molecules being zwitterionic nature, a slightest over-or under-estimation of the localization/delocalization can strongly affect the structure-property correlations.

Table 1
[36]e1are the directly associated with the two aryl units present in the molecule, Important structural parameters (experimental crystal structure data: Scheme 1[36]) of the Molecule 1, computed using various methodologies observations indicate the DFT based methods were reproducing some of those values closer to the experimental data, compared to all other methods (even the post-SCF methods).Nevertheless, the central aryl-aryl junction bond and the twist angle which are important in establishing effective communication between the donor and acceptor sides, indicated that the values obtained in HF and post-HF methods are closer to experimental results than the DFT based methods.A direct comparison of CASSCF/6-31 + + G(d,p) and HF/6-31 + + G(d,p) methods can be made as both were performed in same basis set combinations.Such a comparison indicates that both the methodologies predicted near equivalent geometries, and very close to the experimental crystal structure data[36].Based on this analogy and comparison with the experimental data, one can say that even HF/6-31 + + G(d,p) is the most efficient methodology to properly account for the structure of this zwitterion and needed to be tested for other zwitterionic systems to see if it is a general trend or not.

Table 3
[36]le moments in Debye (D) for the Molecule 2, computed using various methodologies.Experimental dipole moment shown is from reference[36].Molecule 2 in its optimized geometric orientation is also shown

Table 4
Orbital energies of HOMO, LUMO and HOMO-LUMO gap in electron Volts (eV) of the Molecule 1 computed using various methodologies a represents the values from aug-cc-pVDZ basis set computation