Is HAM / 3 ( Hydrogenic Atoms in Molecules , Version 3 ) a Semiempirical Version of DFT ( Density Functional Theory ) for Ionization Processes ?

Nós calculamos os potenciais de ionização verticais (VIPs) de nove moléculas pequenas, bem como, os potenciais de uracil e do C 2 F 2 usando diferentes métodos: i) semi-empírico HAM/3; ii) semi-empírico AM1; iii) não empírico Teoria do Funcional de Densidade (TFD) com os modelos uDI(B88-P86)/cc-pVTZ e -ε(SAOP)/TZP; iv) ab initio HF/cc-pVTZ. Os resultados numéricos obtidos com HAM/3 são mais próximos dos resultados obtidos com o TFD do que o Hartree-Fock (HF). Nós também calculamos as energias de ligação de elétrons do caroço (CEBE) da anilina, nitrobenzeno, e p-nitroanilina com o HAM/3 e a TFD empregando o método ∆E. O modelo de DFT designado como ∆E KS (PW86-PW91)/TZP produziu resultados precisos de CEBE, com desvio médio absolutos de 0,14 eV. Enquanto que a magnitude absoluta dos CEBEs calculados pelo método HAM/3 tem um erro de menos de 3 eV, os deslocamentos químicos (∆CEBE) têm erros menores que 0,55 eV. Mesmo que os resultados de CEBE não apresentem uma resposta definitiva à pergunta do título, as tendências nos VIPs indicam que o HAM/3 não se aproxima do TFD com potenciais de troca-correlação precisos, mas indicam uma proximidade com funcionais semelhantes ao B88-P86.


Introduction
The semiempirical molecular orbital (MO) method known as Hydrogenic Atoms in Molecules, Version 3 (HAM/3) was developed by group of Lindholm and coworkers in 1977. 1 The method has been successfully applied to the calculation of mainly vertical ionization potentials (VIPs), excitation energies, and electron affinities of wide variety of molecules.Some selected works are listed in references 2-7.However, extensions of HAM/3 to the calculation of X-ray emission spectra, 8 Auger electron spectra, 9 valence-electron shake-up satellites, 10 core-electron shake-up satellites, 11 and relative intensities for valence region of X-ray photoelectron spectra will not be discussed in this paper.The theoretical foundation of HAM/3 has been described in detail in the book of Lindholm and Åsbrink. 12According to the authors, theoretical foundation of HAM/3 method is based on density functional theory (DFT). 13,14In other words, HAM/3 was considered to be a semiempirical version of DFT.In their book, Parr and Yang 15 called attention to the HAM model being a semiempirical version of the Kohn-Sham method.
The increasing popularity of DFT in the last decade is quite remarkable.DFT is one of the most promising theoretical methods to study electronic structure of medium and relatively large molecular system as well as solid state in high accuracy.Although there are empirical parameters in most choices of the exchange-correlation functional E xc , the parameters are assumed to be universal for all systems (and hence not empirical for each element), and consequently DFT calculations are often considered to be a priori, firstprinciple, or nonempirical.The advancement of DFT is due to development of more and more reliable exchangecorrelation functionals E xc , together with availability of powerful computers with relatively low cost.However, there still remains an area in which a semiempirical method such as HAM/3 can play.Computer times of calculation for a large molecular system by HAM/3 are orders of magnitude smaller than nonempirical DFT calculation.There are occasions when calculation of approximate molecular properties of a large number of molecules is desired in short time.In such cases, a semiempirical method is the natural method of choice.The object of the present paper is to investigate numerically to determine whether or not HAM/3 is really a semiempirical version of DFT for ionization processes.We therefore calculate vertical ionization potentials (VIPs) of outer valence electrons as well as coreelectron binding energies (CEBEs) of some molecules for comparable studies.
Before we begin, let us briefly review the theoretical studies of VIPs and CEBEs.Early studies of VIPs, using semiempirical (such as CNDO, MINDO, AM1, etc.) or ab initio Hartree-Fock (HF) methods, rely on Koopmans' theorem (KT).HAM/3 is an exception, using a simplified version of Slater's transition-state model 16 called diffuse ionization (DI) for VIPs. 1 On the other hand, a number of DFT calculations of VIPs have been performed previously.Initial studies 17 showed that the VIPs for the lowest cationic states of each symmetry can be computed quite well, with average absolute deviations (AADs) of 0.3 eV from experiment, regardless of whether one uses local density approximation or generalized gradient approximations such as B88 18 -P86 19 or PW86 20 -PW91, 21 and regardless of whether one uses Slater's transition-state method or total energy difference ∆E.However, the list of molecules has not been extensive enough.Shapley and Chong used the PW86-PW91 functional and calculated 181 VIPs of 41 molecules by ∆E method. 22The overall AAD from experiment was 0.55 eV.The most important impact of that study is that many VIPs, even for the lowest cationic state of each symmetry, for perfluoro molecules were in error by over 2 eV (even with several other functionals).More recently, Chong and his collaborators 23 showed that the energies ε k of the outer-valence occupied Kohn-Sham orbitals gave approximate relaxed VIPs better than KT.Excellent agreement between -ε k of outer-valence Kohn-Sham MOs for N 2 , CO, HF and H 2 O and experimental VIPs of the molecules were obtained, especially when accurate V xc was used, the average absolute deviation from experiment being less than 0.1 eV.In addition, calculations of 64 molecules were performed with the approximate V xc known as the statistical averaging of orbital potentials (SAOP).Reasonable agreement between -ε k (SAOP) and the outer valence VIPs (including many perhalo molecules) was found, with AAD of 0.4 eV, much better than KT with HF MOs.In this study, we shall designate the use of -ε k of Kohn-Sham orbital k (instead of Hartree-Fock orbital) to approximate VIP k as meta-Koopmans' theorem (mKT).
Let us now turn to CEBEs.Since the early days of development of experimental technique of X-ray photoelectron spectroscopy, the electrostatic potential models or the thermochemical equivalent-core approximation, 24 were used to interpret chemical shifts of CEBEs.Since then, several reviews have been published. 25n a potential model, for instance, the chemical shift is related to the atomic charge of the atom concerned and to those of surrounding atoms.There are empirical parameters that have to be determined by a least-squares fit to the experimental data in the model.Hence, the potential model can be called as an empirical model.There is fairly good correlation between observed and calculated chemical shift for many cases.However, there are cases where correlation is poor, depending on type of molecules treated.The semiempirical HAM/3 gives the CEBE directly as the difference between the ground-state molecule and the corehole cation.Chong 26 compared CEBEs of a variety of molecules calculated by HAM/3 with observed ones.A persistent error associated with C-H bonds was also reported. 26Recently, Chong 27 proposed a method that enabled one to calculate accurate CEBEs by the densityfunctional theory.The method employs an unrestricted generalized transition-state (uGTS) model. 28Pulfer et al. have confirmed the reliability of the method with a total of seventy-six cases. 29More recently, a more reliable method, called ∆E KS (PW86-PW91)/cc-pCVTZ, was found by Cavigliasso and Chong. 30In the shorthand notation, ∆E KS is the difference in the total energies E's of coreionized cation and neutral parent molecule calculated by DFT using correlation consistent polarized core-valence triple zeta (cc-pCVTZ) basis set; PW86 is Perdew-Wang 1986 exchange functional, 20 and PW91 is Perdew-Wang 1991 correlation functional. 21They used a modified version of DeMon computational package 31 for the calculations.The initial high accuracy of molecular CEBEs, with AAD of only 0.15 eV for 32 cases, was confirmed later with a total of 78 cases. 32An extensive review of DFT calculation of CEBE has been given in a chapter contributed by Chong. 33In this chapter, he also proposed the use of Amsterdam Density Functional (ADF) package 34 for calculation of accurate CEBEs, because localizing a core hole is much easier with ADF than with other program such as DeMon.

Method
In this work, we calculate the VIPs of nine small molecules, as well as uracil and C 2 F 4 , by several different methods: HAM/3, two models within nonempirical DFT, ab initio HF, 35 and the semiempirical MO method called Austin Model, Version 1 (AM1). 36In the first DFT method studied here, ionization energies of a system (neutral molecule M) can be calculated employing an approximate version of Slater's transition state (TS) method, in which one half electron is removed from the system.Absolute values of occupied MO energies of the species M +0. 5 give the ionization energies.The HAM/3 calculations were performed with the computer program purchased from QCPE. 37In calculations of VIPs of the molecules by HAM/ 3 using the keyword PES, the restricted diffuse ionization (rDI) model was automatically employed.This model is the first model and it is designated as HAM/3 (rDI).The rDI model is a simplification of Slater's transition state method, by distributing the +0.5 charge over all the valence MOs, with electrons with α and β spins occupying the same spatial orbitals.We also calculated negative of k-th molecular orbital energy (-ε k ) of a neutral form of the nine small molecules by HAM/3 for the sake of comparison.This is the second model and it is designated as HAM/3 (-ε k ).The negative of k-th molecular orbital energy (-ε k ) of a neutral molecule would represent KT-like ionization potential.One nonempirical DFT model for calculation of valence electron VIPs is unrestricted diffuse ionization (uDI) model, uDI(B88-P86)/cc-pVTZ, which is the third model and it is designated simply by uDI.The calculations were done with DeMon computational package.In uDI model, +0.5 charge is distributed only α-spin valence MOs.We used the exchange-correlation potential labeled as B88-P86, made from Becke's 1988 exchange functional 18 and Perdew's 1986 correlation functional 19 in the uDI calculations.We also calculated negative of k-th molecular orbital energy (-ε k ) of a neutral form of the nine small molecules by nonempirical DFT(B88-P86)/cc-pVTZ, which is our fourth model.It is designated as DFT(B88-P86)/cc-pVTZ(mKT). The other DFT model tested is -ε(SAOP)/TZP(mKT), which is the fifth model and it is designated simply by SAOP, using the ADF program with a triple-zeta polarized (TZP) basis set of Slater-type orbitals.In SAOP, electron density is calculated with the exchange-correlation potential V xc SAOP in SCF process, while energy is calculated with PW91x-PW91c functionals. 21In HF/cc-pVTZ (6 th model) and AM1 (7 th model) calculations with Gaussian 94 program, 38 KT was assumed.The use of the notation KT is restricted for HF and approximate HF methods (such as AM1 method) in this work, while mKT is reserved for the case of DFT, in which -ε of Kohn-Sham orbitals is used to approximate VIPs.The calculated results are compared to each other.Total energy (E) in DFT is exact in the sense that it includes "correlation energy", while total energy in HF theory does not contain the correlation energy.Because of this, ∆E method works well in DFT for calculation such as ionization energy of a system, in which ∆E is the difference of total energies between ionized and neutral molecules.
CEBEs of aniline, nitrobenzene, p-nitroaniline, and uracil were calculated by HAM/3 and DFT with the ∆E method.The ADF package was used for the DFT calculations.In all DFT calculations with ADF, we used the exchange-correlation potential labeled as PW86-PW91, made from Perdew-Wang 1986 exchange functional 20 and Perdew-Wang 1991 correlation functional. 21Basis set of triple zeta plus one polarization, TZP, was used.This model can be expressed as ∆E KS (PW86-W91)/TZP.Experimental CEBEs of uracil were observed in solid-state.CEBE (solid) and CEBE(gas) can be related by equation 1, CEBE(gas) = CEBE(solid) + WD (1)   where WD (W for work function and D for delta, other energies) is energy shift due to solid-state effects.We adopt molecular solids model.We estimate WD to be average absolute deviation between calculated CEBE (gas) and observed CEBE (solid).The results were compared with those observed and those calculated with nonempirical DFT and HAM/3.Experimental geometry was used for H 2 O and uracil.Molecular geometries were calculated for the remaining molecules with HF/6-31G*.

Results and Discussion
Table 1 summarizes VIPs of a total of 30 cases originating from the nine small molecules, calculated by the seven different models; (1) HAM/3(rDI), (2) HAM/3(-ε k ), ( 3) uDI(B88-P86)/cc-pVTZ, (4) DFT (B88-P86)/cc-pVTZ(-ε k ), (5) -ε(SAOP)/TZP, (6) HF/cc-pVTZ, and (7) AM1.In case of H 2 O, for example, all models, except the models (2) and ( 4), reproduce experimentally obtained vertical ionization energies fairly satisfactorily for the lowest three energies corresponding to 1b 1 , 3a 1 and 1b 2 .However, the error of the calculated VIP for the low-lying 2a 1 orbital by both HF and AM1 is more than 4 eV (see Figure 1).Errors of (1) HAM/3 (rDI) and the two DFT models, (3) uDI and (4) SAOP, for the same event are much smaller than 4 eV.In DFT, the error due to correlation effect is expected to be small.The fact that the large errors registered in VIP(2a 1 ) by HF and AM1 may be due to mainly electron correlation effect, since no correlation effect is taken into account in HF methods.In case of the two models, (2) HAM/3(-ε k ) and ( 4) DFT (B88-P86)/cc-pVTZ(-ε k ), the discrepancy between experimental ionization potentials and corresponding negative values of k-th molecular orbital energy (-ε k ) of the neutral form of H 2 O is very large.The -ε k values in the two models are uniformly smaller than the experimental ionization potentials by approximately 5 eV.The numerical -ε k values in the models (2) are fairly close to those in the model ( 4).Situations similar to those seen in H 2 O can be observed in Table 1.Ionization potentials (in eV) of nine small molecules calculated by seven different models; (1) HAM/3(rDI), (2) HAM/3(-ε k ), ( 3) uDI(B88-P86)/cc-pVTZ, ( 4) DFT (B88-P86)/cc-pVTZ(-ε k ) (5) -ε(SAOP)/TZP, (6) HF/cc-pVTZ, and ( 7 2) and (4), at the same time in the other pair of the two models, (1) and ( 3), HAM/3 can be considered as a semiempirical version of nonempirical DFT (B88-P86)/cc-pVTZ.On the other hand, AM1 semiempirical method reproduces VIPs of ab initio HF method well.This is because AM1 is a semiempirical version of HF method, and not that of DFT.The model (5) -ε(SAOP)/TZP resulted AAD of 0.65 eV, which is to be compared with AAD of 1.61 eV resulted from ( 6) ab initio HF/cc-pVTZ (KT).The model -ε(SAOP) with nonempirical DFT is superior to KT in ab initio HF.Since the model (2) HAM/3(-ε k ) and the model ( 4) DFT (B88-P86)/cc-pVTZ(-ε k ) resulted large systematic errors for calculating VIP's of the small molecules in Table 1, we shall not consider them any further and work only with the remaining five models for the rest of discussions.Table 2 compares 13 lowest VIPs of uracil calculated by the five different models.The experimental VIPs were observed in gas-phase. 39Only the five lowest ionization events (5a", 24a', 4a", 23a', and 3a") have been assigned.AAD, corresponding only to the five lowest ionization events, of HAM/3 (rDI) is ca.0.2 eV, while AADs of the three models, uDI, SAOP, and AM1, are in the neighborhood of 1 eV.AAD of HF is ca.1.5 eV.If the 13 ionization events are taken into account (based on the ordering of HAM/3 and SAOP), the AADs take the values of 0.19 eV (HAM/3), 0.80 eV (uDI), 0.94 eV (SAOP), 2.20 eV (HF) and 1.27 eV (AM1).AAD of HAM/3 (rDI) is the smallest; AADs of the  two DFT models are close to each other; and AAD of HF is more than twice to those of DFT.In other words, VIPs calculated with rDI in HAM/3 and the two nonempirical DFT models are closer to observed values than those calculated with KT in ab initio HF.The spectral patterns for the five lowest ionization events obtained by HAM/3 (rDI), uDI and SAOP reproduce that observed fairly well (see Figure 2).AM1 results parallel closely to those of HF, but neither reproduces the observed spectral pattern well.Shapley and Chong 40 calculated VIPs of C 2 F 4 with ∆E approach using seven different DFT exchange-correlation functionals.Most AADs were in the neighborhood of 1.5 eV.
The smallest AAD attained was 0.79 eV.We wanted to know if any one of our five methods can give better result than this.Table 3 lists VIPs of C 2 F 4 calculated with the five different methods.The value of AAD increases in the order HAM/3 (rDI) < SAOP < AM1< uDI < HF.The smallest AAD, 0.27 eV, is registered by the semiempirical HAM/3 (rDI) method.This may be because C 2 F 4 is one of the molecules that were used for parametrization of the method.Next, the SAOP model gave AAD of 0.36 eV, which is much smaller than 1.5 eV.In fact, Chong et al. 23 obtained AAD of 0.35 eV with SAOP model for calculated VIPs of ten perhalo molecules.Our result is just a confirmation of the success of the SAOP model.The semiempirical AM1 showed AAD of 0.77 eV.The AAD of nonempirical uDI and ab initio HF showed larger AAD values, more than 1.5 eV's.B88-P86 functional used in the uDI calculation was not appropriate for the type of molecule.Sharpley and Chong 40 obtained similar results as our uDI result using various other types of functionals.The -ε(SAOP)/TZP model proved to be a good choice for this type of molecules, as found earlier. 23rror due to correlation effect, together with the approximate nature of KT, is large in HF method for C 2 F 4 .
Table 4 lists CEBEs of substituted benzenes (X-C 6 H 4 -Y) such as benzene, aniline, nitrobenzene and p-nitroaniline calculated with HAM/3 and with ∆E KS (PW86-PW91)/TZP-ADF models.AAD obtained with HAM/3 is 3.25 eV, which is much larger than 0.14 eV, the AAD obtained with the nonempirical DFT calculations.The ∆E KS (PW86-PW91)/ TZP method works very well. 41AAD of 0.15 eV had also been obtained previously for exactly the same set of molecules using the uGTS model. 42The DeMon DFT program was employed in the uGTS calculations.We designate this as uGTS(DeMon) model.Both ∆E KS (PW86-PW91)/TZP(ADF) method and the uGTS(DeMon) model  resulted almost the same AAD.Although the quality of the two models is equivalent, one major difference is the presence/absence of a fortuitous partial cancellation between model error and functional error, making the ∆E KS method more reliable conceptually.On the other hand, the main difference between the two programs DeMon and ADF resides on efficiency.∆E KS (PW86-PW91)/TZP using ADF is almost five-to ten-fold more efficient than the uGTS(DeMon) model, or ∆E KS (PW86-PW91)/cc-pCVTZ(DeMon) method.It should be mentioned that ∆E KS (PW86-PW91)/cc-pCVTZ(DeMon) gives almost identical CEBEs as ∆E KS (PW86-PW91)/TZP(ADF), but much less efficiently.On the other hand, CEBEs calculated with HAM/3 are uniformly shifted toward smaller values with respect to the observed values by approximately 4 eV, in case of carbon atoms.
∆CEBE is the difference between CEBE of an atom in a substituted benzene and CEBE of the corresponding atom of unsubstituted benzene (the reference molecule).∆CEBE at C1 atom is positive in the two molecules with the values of ca.+0.9 eV in case of aniline and ca.+1.7 eV in nitrobenzene.These large positive shifts are due to the substituent, either -NH 2 or -NO 2 , at C1 position.The ∆CEBEs at C2, C3 and C4 atoms are negative in aniline because -NH 2 has strong electron donating nature.Electron density increases at C2, C3 and C4 in the phenyl ring that causes destabilization of the core electrons.Absolute value of shift increases in the order C3 < C2 < C4.This is due to the fact that the electron donating effect of -NH 2 in aniline increases in the order: C3 < C2 < C4.Nonempirical DFT results reproduce the observed results of ∆CEBEs and tendency very well (see Figure 3).HAM/3 resulted ∆CEBE for C1 position too large, by more than 1 eV, in comparison to the observed value.∆CEBEs of HAM/3 at C2, C3 and C4 atoms are approximately equal to those of observed values.However, the order of shift is C2 < C3 < C4, disagreeing with the observed order of C3 < C2 <C4 (Figure 3).HAM/3 did not reproduce correctly the order.In the case of nitrobenzene, the signs of ∆CEBEs at C2, C3 and C4 atoms are all positive.This is because -NO 2 group has a strong electron withdrawing character.Electron density in the phenyl ring decrease because of -NO 2. This stabilizes core electrons of the atoms in the ring causing the substantial increase of CEBEs.The magnitude of the shift As far as we know, HAM/3 is the only semiempirical SCF method that is capable of calculating molecular CEBEs by ∆E method.Consequently, calculations of ∆CEBEs of medium and large molecular systems can be done conveniently with HAM/3.Table 6 lists CEBEs of solid-state uracil calculated by equation 1 using CEBEs of gas-phase uracil calculated with ∆E KS (PW86-PW91)/TZP-ADF model and HAM/3.In principle, there is only one true value of WD for the system we study.However, we adopted two empirical values of WDs (WD1 and WD2) for the two different theoretical methods.This is a consequence of the method of approximation of WD that we adopted.Nonempirical DFT  The -ε(SAOP)/TZP model consistently gave results similar to (or better than) uDI(B88-P86)/cc-pVTZ model for calculation of VIPs.This indicates superiority of the exchange-correlation potential V xc = SAOP in comparison to B88-P86.The -ε(SAOP)/TZP model is like KT, but for Kohn-Sham orbitals, and so we call it meta-Koopmans' theorem, or mKT.In other words, if one has accurate functionals, then there is no need for uDI or uTS model.It is because B88-P86 functional and most other functionals are poor, especially in the large-r region, that we need to use uDI (or uTS) model to get better results.It appears that HAM/3 is approximating such an approximate procedure.Therefore, strictly speaking, HAM/3 is not an approximation for DFT, but rather an approximation for the restricted diffuse ionization (rDI) model, which is an approximation used in some DFT calculations of VIPs.We demonstrated numerically that values and trend of VIPs calculated by HAM/3 reproduces those obtained by uDI(B88-P86)/cc-pVTZ and -ε(SAOP)/TZP models.
In calculation of accurate CEBEs of the molecules, we confirmed high efficiency of the ∆E KS (PW86-PW91)/TZP method.HAM/3 also calculates CEBE by ∆E, with even greater efficiency, but with errors of more than 3 eV.Hence, the present study of CEBEs based on HAM/3 and accurate DFT procedure cannot answer the question in the title of this work.Fortunately, the results showed that ∆CEBEs (chemical shifts) obtained with HAM/3 are fairly close to the corresponding observed values.Hence, it may be useful for computational studies of relative chemical and/or biological activities of large systems.
) AM1.Average absolute deviations (AAD's) are listed in the last row of the table.rDI and uDI are the two transition state methods.-ε k signifies negative of orbital energy of a neutral molecule.KT signifies Koopmans's theorem the remaining eight small molecules in the Table1.The AADs obtained for the 30 cases in Table1of the seven different models,

Table 2 .
Ionization energies (in eV) of uracil calculated by HAM/3 (rDI), uDI, SAOP, HF, and AM1.Average absolute deviations (AAD's) are listed in the last row of the table a The experimental VIPs were observed in gas-phase.The values in parentheses are those read from photoelectron spectrum of uracil in the literature (Reference 39).

Table 3 .
Ionization potentials (in eV) of C 2 F 4 calculated by HAM/3 (rDI), uDI, SAOP, HF, and AM1.Average absolute deviations (AAD's) are listed in the last row of the table a Reference 49.

Table 5 .
Observed and calculated chemical shifts ∆CEBEs (in eV), of carbon atoms in benzene ring of substituted benzenes (C 6 H 5 -X).∆CEBEs in the table have been calculated as the difference between CEBE of carbon atom in the substituted benzene and CEBE of the carbon atom in benzene (equation 2) 43creases in the order C3 < C4 < C2 < C1 (observed).The nonempirical DFT calculations reproduce the observed results fairly well.The HAM/3 ∆CEBE at C1 atom are more than 2 eV in error.However, ∆CEBE values of C3, C4 and C2 are close to the observed ones.The order of shifts of C3 and C4 is inverted in HAM/3.Lindberg et al.43had shown that ∆CEBE correlate linearly to the Hammett sigma constants (σ) in substituted benzenes.Good agreements between Hammett σ's and ∆CEBEs at C3 (meta-) and C4 (para-) positions in C 6 H 5 -NO 2 can be seen in those observed, and calculated with DFT by ∆E KS (PW86-PW91)/TZP.Fairly good agreements are also obtained with HAM/3.This indicates a potential utility of HAM/3 in investigations of chemical reactivity and/or biological activity of series of molecules.When one works with a series of molecules studying relative chemical and/or biological activity, it is not necessary to have accurate CEBEs themselves, if reliable ∆CEBE values are available.