11596
J. Phys. Chem. 1996, 100, 11596-11599
Calculation of Proton Affinities Using Density Functional Procedures: A Critical Study A. K. Chandra and Annick Goursot* Ecole Nationale Superieure de Chimie de Montpellier, URA 418 CNRS, 8 Rue de l’Ecole Normale, 34053 Montpellier, Cedex-1, France ReceiVed: February 7, 1996; In Final Form: April 24, 1996X
Downloaded via UNIV OF CALGARY on June 20, 2018 at 06:25:43 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.
Density functional theory has been used with different combinations of exchange and correlation functionals to study the proton affinities of six organic molecules, namely H2CO, CH3CHO, CH3OH, C2H5OH, HCOOH, and CH3COOH. Complete geometry optimizations have been carried out for both the neutral and protonated species with all combinations of functionals. The proton affinity values are then compared with the corresponding experimental values. It has been observed that combination of Perdew’s and Becke’s exchanges with Proynov’s correlation functional is the most effective in reproducing the proton affinity.
1. Introduction Density functional theory (DFT) has achieved considerable degree of success in the computation of molecular properties.1 The Gaussian-2 (G2) method,2 introduced by Pople and coworkers, has been found2,3 to provide estimate of proton affinity (PA) within a target accuracy of 2 kcal/mol. Smith and Radom4 also showed that the more economical G2(MP2) method5 produces proton affinities within this target accuracy. It is wellknown that computational demand is always much less in DFT than in traditional ab initio molecular orbital (MO) calculations when density is used as the basic variable. DFT-based studies have also been devoted to the predictions of thermochemical properties including proton affinity.4,6-8 Very recently, Schmiedekamp et al.9 calculated the PA of some nitrogen- and oxygencontaining molecules with a density functional scheme, at different levels of sophistication. All these studies have shown that a satisfactory prediction of PA values necessitates the use of gradient-corrected exchange and correlation functionals, the local density approximations yielding underestimated values. There are a few DFT studies10-14 on hydrogen-bonded systems and on dispersion-dominated van der Waals interactions and the results are diverse, showing that improvement of exchange and correlation energy functional (henceforth referred to as XC functionals) is still needed in order to achieve a better description of weakly bound systems. Recently, Proynov et al.15 proposed a kinetic energy density-dependent correlation functional, developed in the context of the adiabatic connection procedure. Extensive tests have shown that this functional gives very reliable results on H-bonded and weakly interacting systems. Hence we felt it would be interesting to explore various combinations of exchange and correlation potential functionals with special attention on Proynov’s correlation functional and find out the most effective one for the study of proton affinity (PA). Six organic molecules have been chosen to this end. There are recent proton affinity calculations16,17 by HartreeFock (HF) based methods for some of the molecules chosen for our study. Five combinations of exchange and correlation potentials have been chosen for the present study. The first one corresponds to the local level of theory using the DiracSlater exchange term18 and the Vosko-Wilk-Nusair (VWN) parametrization19 for the correlation energy. In the second procedure we used Perdew’s exchange and correlation functionals20 (henceforth designated as PD86). Becke’s exchange21 and Perdew’s correlation functional have been taken in the third X
Abstract published in AdVance ACS Abstracts, June 1, 1996.
S0022-3654(96)00375-9 CCC: $12.00
procedure. The fourth set consists of Perdew’s exchange functional and a correlation functional, recently proposed by Proynov et al.15 (hereafter referred to as LAP). In the last procedure we used Becke’s exchange and Proynov’s correlation functionals. The five procedures mentioned above will be designated henceforth as VWN, PP, BP, PLAP, and BLAP, respectively. Although this study does not present exhaustive results for exchange and correlation functionals, we hope it will provide basic information about the calculation of PA of molecules with DFT and contribute to an on going extensive body of work on the efficiency of DFT calculations in reproducing experimental data.1,22 2. Method Standard DFT calculations have been performed with different exchange and correlational functionals using the deMon program.23-25 As already mentioned in the Introduction, five procedures, VWN, PP, BP, PLAP, and BLAP, have been used. In the last four procedures, the calculations correspond to the so-called nonlocal level of approximation where nonlocal XC functionals are used for energy and potential. All different procedures are placed in the same numerical environment (grid, basis set, etc.). Grids for numerical integration correspond to 64 radial points and a number of angular points varying from 50 to 194 with the distance from the nucleus. For C and O (6311/311/1) and H (41/1) orbital basis sets with corresponding (5,2:5,2) and (5,1:5,1) auxiliary bases26 were employed. The (6311/311/1) basis set is equivalent to the well-known 6-31+G* basis used in ab initio MO calculations.27 The geometries of both the neutral and protonated molecules were fully optimized at the nonlocal level. 3. Proton Affinity The proton affinity of a molecule, B, is defined as negative of the change in enthalpy for the reaction
B + H+ f BH+
(1)
For evaluating the absolute proton affinity, the following expression is used:16,28-30
PA(B) ) -∆H(eq 1) ) -∆Eo - ∆ZPE + ∆Ev(T) + C (2) The first term represents the difference in electronic energies between protonated and unprotonated species, and the second © 1996 American Chemical Society
Calculation of Proton Affinities
J. Phys. Chem., Vol. 100, No. 28, 1996 11597
TABLE 1: Total Energies (in au) of the Stationary Points for Neutral and Protonated Molecules (Bottom Line for Each Molecule) VWN
PP
BP
PLAP
BLAP
H2CO -113.599 40 -114.633 30 -114.518 79 -114.617 05 -114.502 23 -113.869 29 -114.907 92 -114.796 37 -114.898 25 -114.786 34 CH3CHO -152.564 62 -154.016 56 -153.853 90 -153.980 03 -153.816 90 -152.862 48 -154.317 34 -154.157 95 -154.285 61 -154.125 57 CH3OH -114.796 39 -115.853 13 -115.737 68 -115.825 41 -115.709 35 -115.081 45 -116.143 66 -116.030 57 -116.122 99 -116.009 34 C2H5OH -153.749 80 -155.225 56 -155.061 46 -155.177 70 -155.012 94 -154.043 98 -155.525 96 -155.364 67 -155.485 83 -155.323 07 HCOOH -188.354 12 -189.972 99 -189.793 10 -189.942 81 -189.762 58 -188.637 16 -190.259 47 -190.082 93 -190.234 79 -190.057 86 CH3COOH -227.314 32 -229.351 07 -229.123 37 -229.301 08 -229.072 49 -227.616 24 -229.655 94 -229.431 51 -229.610 51 -229.385 92
TABLE 2: Absolute Proton Affinity Values (PA298 in kcal/mol) Obtained from Our DFT Calculations Using Eq 2 VWN
PP
BP
PLAP
BLAP
expta
H2CO
162.42
165.39
167.25
169.52
171.34
HCOOH CH3OH CH3CHO
171.28 172.33 180.51
173.44 175.76 182.34
175.54 177.26 184.39
176.90 180.20 185.36
178.96 181.72 187.29
C2H5OH CH3COOH
178.24 183.26
182.15 185.11
183.91 187.16
187.00 187.97
188.25 190.50
171.7 170.4b 178.8 181.7 186.6 183.6b 188.3 190.2
system
a
Reference 32. b Reference 17.
term is the difference in zero-point energies (ZPE) of BH+ and B; the third term appears due to change in population of vibrational levels with the change in temperature, and the last term introduces the correction for translational and rotational energy changes assuming classical behavior and a ∆nRT term required to convert an energy to enthalpy assuming ideal gas behavior. The first term can easily be evaluated from the energy difference between the protonated and unprotonated species. Table 1 shows the total energies of the molecules and the corresponding protonated species at their respective optimized geometries for all the five procedures. The second and third terms are calculated from the frequencies of the normal modes of vibration by performing a vibrational analysis of both the neutral and protonated species. Let us mention that vibrational analysis was performed only for the PP procedure, the same correction being applied to the other procedures. We checked the validity of this approximation for H2CO, where the ZPE corrections obtained with PP and PLAP were 8.38 and 8.51 kcal/mol, respectively. We believe that such a small difference in zero-point energy corrections in no way upset our final conclusions. The fourth term is taken as 1.5 kcal/mol16 at 298 K for all the molecules included in the present study. In Table 2, we have tabulated the absolute proton affinity (PA298) values obtained through eq 2. Here the molecules are arranged in increasing order of their experimental PA values in order to observe, additionally, whether our DFT calculations also predict the same order of PA. From the table, it is obvious that the PA value of every molecule increases on going from VWN to BLAP through PP, BP, and PLAP. However, the results obtained with PLAP and BLAP, particularly BLAP, are in very good agreement with experiment. BSSE corrections
TABLE 3: Difference between Experimental and Calculated Proton Affinity Values (in kcal/mol), and the Average Absolute Error for Comparison system a
H2CO HCOOH CH3OH CH3CHOd C2H5OH CH3COOH error
VWN
PP
BP
PLAP
BLAP
7.98 7.52 9.37 6.09 10.06 6.94 7.99
5.01 5.36 5.94 4.26 6.15 5.09 5.30
3.15 3.26 4.44 2.21 4.39 3.04 3.41
0.88 1.90 1.50 1.24 1.30 2.23 1.51
-0.94 -0.16 -0.02 -0.69 0.05 -0.30 0.36
a Experimental values are taken as 170.4 and 186.6 kcal/mol for H2CO and CH3CHO, respectively.
have not been performed for all systems, because this correction, which is expected to be very small, has indeed been found negligible for H2CO (0.04 kcal in the case of PLAP). In order to illustrate the uncertainty in PA calculations the differences between experimental and calculated PA values are presented in Table 3. The PA values obtained at the local level are far off from the experimental values and the average error is around 8 kcal/mol. At the nonlocal level, the PP procedure is also not in close agreement with experiment (mean error is 5.3 kcal/ mol) and always underestimates the PA values. But dramatic improvement occurs when Proynov’s correlation functional is used instead of Perdew’s. The absolute PA values obtained with PLAP and BLAP are in very good agreement with experiment. The average absolute error goes down to 1.51 and 0.36 kcal/mol for PLAP and BLAP, respectively. The agreement is very good, particularly in view of the target accuracy in ab initio MO calculations. The PA values obtained through BP are also consistent with the average error is 3.41 kcal. From Table 3, it is also obvious that, generally, PLAP slightly under estimates the PA values whereas the reverse is true for BLAP. It is interesting to know that a computationally efficient DFT procedure can provide PA values comparable with those obtained through very elaborate ab initio MO calculations. From these results, it is seen that the correlation functional plays a very critical role in this evaluation. Indeed, the use of the LAP correlation functional instead of the Perdew functional increases systematically the PP and BP values by 2-5 kcal. However, the incidence of the exchange functionals (Perdew to Becke) is not negligible, since there is a systematic difference of 2 kcal/ mol between the results obtained with Perdew and Becke exchange functionals. Before leaving this section, it will be useful to note some points about the geometries of the neutral and protonated molecules. The geometries of the neutral molecules used in the present investigation are well established both experimentally and theoretically. Thus, instead of presenting the individual structure, we calculated the average absolute errors in the evaluation of bond lengths and bond angles. For bond length and bond angle (in bracket) the errors are 0.014 (1.1), 0.008 (1.2), and 0.007 (1.2) for PP, PLAP, and BLAP, respectively. This clearly indicates that the procedures used here are good in reproducing the equilibrium geometries of the neutral molecules. For the protonated species, we obtained31 the same type of structures as those obtained from Hartree-Fock SCF calculations with 6-31G** basis functions.16 Moreover, the geometrical parameters change very little (not more than 0.01 Å for bond length and 1° for bond angle) with the change in XC functionals. At the same time, substantial geometrical rearrangements take place after protonation. For example in PLAP, the CdO bond in H2CO increases from 1.214 to 1.257 Å upon protonation. Similar effects have also been observed for the other molecules.
11598 J. Phys. Chem., Vol. 100, No. 28, 1996
Chandra and Goursot
TABLE 4: Net Negative Charge on Oxygen Atom of Different Molecules at Which Proton Is Getting Attached (Obtained from Mulliken Population Analysis) molecule
PP
BP
PLAP
BLAP
H2CO HCOOH CH3OH CH3CHO C2H5OH CH3COOH
0.370 0.410 0.377 0.363 0.352 0.422
0.355 0.404 0.366 0.352 0.333 0.417
0.397 0.424 0.381 0.382 0.358 0.444
0.382 0.420 0.370 0.371 0.342 0.443
TABLE 5: Minimum Molecular Electrostatic Potential Values (kcal/mol) Near the Oyxgen Atom at Which Proton Attachment Takes Placea molecule
PP
BP
PLAP
BLAP
H2CO HCOOH CH3OH CH3CHO C2H5OH CH3COOH
37.4 36.8 46.6 39.0 46.3 40.4
37.8 37.8 47.1 40.6 47.3 41.9
39.0 38.3 49.3 43.0 48.4 42.7
39.4 39.3 49.7 44.4 49.1 42.8
a
Figure 1. Plot of differences in electron density (electrons/Å3) distribution in the molecular plane (XZ plane with CdO bond along the X-axis) between PLAP and PP for the H2CO molecule. The densities have been suitably scaled for easy visualization. The two hydrogen atoms connected to carbon atom of H2CO (not marked in the figure) are symmetrically placed about the O-C bond with C2V symmetry.
The potentials are always negative.
4. Role of Exchange and Correlation Functionals From Table 3, it is seen that going from PP to PLAP (the exchange potential remaining the same (PD86)), the PA values increase by 3-5 kcal, which means that the use of LAP, as a correlation functional, introduces a large improvement in the PA values. From PLAP to BLAP or from PP to BP, the exchange functional is changed keeping the same correlation functional (PD86 or LAP). In both cases, the Becke’s exchange functional leads to an increased binding energy for the additional proton by about 2 kcal and this increment is not very sensitive on the correlation functional. PA is the measure of the stabilization undergone by a molecule after proton attachment. When the proton binds to the host atom, the positions of the nuclei and the electron density rearrange. Relative values of proton affinity must then be considered from a geometric and electronic point of view. As mentioned above, the geometries of both neutral and protonated molecules obtained with the different exchange and correlation functionals are not significantly different. The increase in PA values going from VWN to BLAP for every molecule is thus related essentially to electronic redistribution. The success of PLAP (or BLAP) is then attributed to a more realistic description of the electron density distribution over the molecule. In order to compare the electron densities obtained with the different functionals, one can simply check the charge on the oxygen atom at which the proton is getting attached. Qualitatively, it could be expected, for a given molecule, that the greater the negative charge on the oxygen atom the larger will be the proton affinity. In Table 4, we report the net charges on oxygen atoms obtained from the Mulliken population analysis. One observes that, generally, the net negative charge on oxygen atom at which proton attachment takes place increases from PP to PLAP or from BP to BLAP as expected from the PA values. As the total number of electrons in the system is constant, it indicates that there must be depletion of charge density at some other center. Indeed, analysis of the Mulliken charges shows that going from PP to PLAP (or BP to BLAP), the net positive charge on the next carbon atom increases for all the molecules studied here. This clearly indicates a shift of more electron density from carbon to oxygen atoms when we use LAP as correlation functional. However, the net charges on oxygen atoms reported in Table 4 do not exhibit the same dependence on XC functionals as PA values for a given molecule.
Figure 2. Plot of differences in electron density (electrons/Å3) distribution on the molecular plane (XZ plane with CdO bond along the X-axis) between BP and PP for the H2CO molecule. The densities have been suitably scaled for easy visualization. The hydrogen atoms connected to carbon atom of H2CO (not marked in the figure) are symmetrically placed about the O-C bond with C2V symmetry.
Moreover, it is difficult to explain the effect of exchange functional from Mulliken charges alone. It is well-known that the molecular electrostatic potential (MEP) is a much more reliable property for the description of the electron distribution among the atoms of a molecule. It is thus expected that the minimum MEP values found near the oxygen atom where proton attaches itself will correlate with the changes in PA values for a given molecule. Table 5 shows the minimum MEP values obtained with the different functionals. It is gratifying to note that the MEP values increase in the same sequence as PA for a particular molecule. Since DFT deals with the electron density, maps of electron density differences also give clear indication of electronic redistributions. In order to illustrate these changes associated with different XC functionals, differences in total electron densities of the H2CO molecule are also reported. Figure 1 illustrates the difference, drawn in the molecular plane, between the total electron density obtained from a PLAP calculation and that obtained from a PP calculation. Figure 2 shows the same difference between BP and PP calculations. Changing the correlation functional (Figure 1) or the exchange functional (Figure 2) leads to a clear redistribution of electron density near the carbon and oxygen atoms. Very similar features are obtained in a plane perpendicular to the sigma plane and containing the CdO bond, although the differences are less pronounced. Both pictures exhibit a shift of electron density in the C region, which could
Calculation of Proton Affinities correspond to some transfer from the carbon atom to the C-O bond. The LAP correlation functional always brings more electron density near the oxygen atom. As mentioned above, the nonlocal calculations reported here have been performed using both nonlocal XC functionals for the energy and the potential. The various XC potential functionals are thus acting differently on the density, leading to different electron distributions. We would like to stress that, indeed, PA is a very sensitive property which allows to probe electron density redistribution. A careful analysis of the additive contributions to the total energy for the different functionals (kinetic, nuclearelectron and nuclear-nuclear interactions, Coulomb and XC terms) shows that the changes in XC energy, although not negligible (particularly the correlation energy which is always 2-4 kcal/mol larger with LAP), are not the most substantial ones. Both Coulomb and electron-nuclear energies are very sensitive to the XC functionals. This result indicates that the action of the XC potential functional is preponderant, for a realistic distribution of electron density in space. PA results then from a delicate balance of the energy terms. It is worth noting that, if the MEP values reflect, in some sense, the variation of electron density with the XC functionals for a given molecule, their evolution for the different molecules of Table 5 does not follow the evolution of PA values. For example, the MEP values are larger for CH3OH than for CH3COOH, whereas the PA values are in reverse order. This result can be explained by the consideration of geometric rearrangements. Indeed, the MEP values are evaluated with the geometry of the neutral molecule. When the proton binds to the oxygen atom, there are simultaneous geometric and electronic rearrangements, which are implicitly taken into account in the PA but not in the MEP values. This effect is qualitatively illustrated by the following remarks. If we remain in the same “family” of compounds, as H2CO and CH3CHO or HCOOH and CH3COOH, the geometric rearrangement upon protonation is comparable and both MEP and PA values increase. For CH3OH and C2H5OH the MEP and PA values do not follow the same sequence, which may be associated with different torsional rearrangements around the C-O-H+ unit due to the presence of an additional methyl group in C2H5OH. This leads to a different structure for the protonated C2H5OH molecule. 5. Conclusion Proton affinity can be calculated quantitatively from density functional procedures provided right combination of exchange and correlation functionals is used. In this context, we observe that the combination of Perdew’s or Becke’s exchange and Proynov’s correlation functionals is very fruitful in reproducing the proton affinity values within a target accuracy of 2 kcal/ mol. We also like to mention that we are now carrying out PA calculations of nitrogen-containing molecules with PLAP and the results are very encouraging (absolute errors are less than 1 kcal). Moreover, this study has helped to clarify some aspects of the dependence of calculated properties on XC functionals.
J. Phys. Chem., Vol. 100, No. 28, 1996 11599 The PA results from a very delicate balance of all energetic contributions, due to a combination of geometric and electronic rearrangements. It is thus a very sensitive probe for testing XC functionals for energy and potential. Acknowledgment. The authors thank Drs. Vladislav Vassiliev and G. Gil for their help during computation. Thanks are also due to the referees for their constructive suggestions. References and Notes (1) Parr, R. G.; Yang, W. Density Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1989. (2) Curtiss, L. A.; Raghavachari, K.; Trucks, G. W.; Pople, J. A. J. Chem. Phys. 1991, 94, 7221. (3) Smith, B. J.; Radom, L. J. Am. Chem. Soc. 1993, 115, 4885. (4) Smith, B. J.; Radom, L. Chem. Phys. Lett. 1994, 231, 345. (5) Curtiss, L. A.; Raghavachari, K.; Pople, J. A. J. Chem. Phys. 1993, 98, 1293. (6) Fitzgerald, G.; Andzelm, J. J. Phys. Chem. 1991, 95, 10531. (7) Becke, A. D. J. Chem. Phys. 1992, 97, 9173; 1993, 98, 1372, 5648. (8) Gill, B. G.; Johnson, P. M. W.; Pople, J. A.; Frisch, M. J. Chem. Phys. Lett. 1992, 194, 499. (9) Schmiedekamp, A. M.; Topol, I. A.; Michejda, C. J. Theor. Chim Acta 1995, 92, 83. (10) Sim, F.; St-Amant, A.; Papai, I.; Salahub, D. R. J. Am. Chem. Soc. 1992, 114, 4391. (11) Kim, K.; Jordan, K. D. J. Phys. Chem. 1994, 98, 10089. (12) Lacks, D. J.; Gordon, R. G. Phys. ReV. A 1993, 47, 4681. (13) Kristyan, S.; Pulay, P. Chem. Phys. Lett. 1994, 229, 175. (14) Perez-Jorda, J. M.; Becke, A. D. Chem. Phys. Lett. 1995, 233, 134. (15) Proynov, E. I.; Ruiz, E.; Vela, A.; Salahub, D. R. Int. J. Quantum Chem., Quantum Chem. Symp. 1995, 29, 61. (16) von Nagy-Felsobuki, E. I.; Kimura, K. J. Phys. Chem. 1990, 94, 8041. (17) Smith, B. J.; Radom, L. J. Phys. Chem. 1995, 99, 6468. (18) Slater, J. C. Quantum Theory of Molecules and Solids; McGrawHill: New York, 1994; Vol. 4. (19) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200. (20) Perdew, J. P. Phys. ReV. B 1986, 33, 8822; erratum in: Phys. ReV. B 1986, 34, 7406. (21) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (22) Salahub, D. R.; Fournier, R.; Mlynarski, P.; Papai, I.; St-Amant, A.; Ushio, J. In Theory and Applications of Density Functional Approaches in Chemistry; Labanowski, J., Andzelm, J., Eds.; Springer-Verlag: Berlin, 1991. (23) St-Amant, A.; Salahub, D. R. Chem. Phys. Lett. 1990, 169, 387. (24) Daul, C.; Goursot, A.; Salahub, D. R. In Numerical Grid Methods and Their Application to Schrodinger’s Equation; Cerjan, C., Ed.; NATO ASI Series Vol. 412; Kluwer Academic Press: Dordrecht, 1993; p 153. (25) St-Amant, A. Ph.D. Thesis, Universite de Montreal, 1991. (26) For a description about basis sets, see: Papai, I.; Goursot, A.; Fajula, F.; Weber, J. J. Phys. Chem. 1994, 98, 4654. (27) Frisch, M. J.; Pople, J. A.; Binkley, J. S. J. Chem. Phys. 1984, 80, 3265. (28) Curtiss, L. A.; Pople, J. A. J. Phys. Chem. 1988, 92, 894. (29) Dixon, D. A.; Gole, J. L.; Komornicki, A. J. Phys. Chem. 1988, 92, 2134. (30) Dixon, D. A.; Lias, S. G. In Molecular Structure and Energetics; Liebman, J. F., Greenberg, A., Eds.; VCH: New York, 1987; Vol. 2. (31) Optimized structures of both protonated and neutral molecules can be obtained on request from the authors as supplementary material. (32) Lias, S. G.; Liebman, J. F.; Levin, R. D. J. Phys. Chem. Ref. Data 1984, 13, 695. Lias, S. G.; Bartmess, J. E.; Liebman, J. F.; Holmes, J. L.; Levin, R. D.; Mallard, W. G. J. Phys. Chem. Ref. Data 1988, Suppl. 17, No. 1.
JP9603750
