Volume 208, number $6
CHEMICAL PHYSICS LETTERS
18 June 1993
Density functional theory applied to proton-transfer systems. A numerical test Claude Mijoule Luboratoirede Dynamiquedes InteractionsMokculaires, Universitt! Pierreet Marie Curie, T22, 75252 ParisCedex, France
Zdzislaw Latajka ’ and Daniel Borgis Laboratoirede PhysiqueThkoriquedes LiquiZes,Universitk Pierreet Marie Curie, T16. 752St ParisCedex, France Received 18 February 1993; in fial form 7 April 1993
The structure, interaction energy and proton-transfer features of the (Hs02) + complex are determined by using a nonlocal density functional theory method developed recently by Salahub and co-workers, The results are compared to those of a highorder Moller-PIesset electron-correlation calculation. The DFT calculation predicts correctly a symmetrical configuration of the complex at the equilibrium O-O distance, and yields a correct value of the H-bonding interaction energy. For a larger 0-O separation, the electronic correlation effects on the double-well proton potential curve are overestimated in the DFT approach, whereas they are well reproduced by the second-order Mpller-Plesset expansion.
ments of the DFT in a chemical context, see also ref.
1. Intmduction
[IllThere is growing interest in the quantum-chemical community for computational methods based on density functional theory (DFT) mainly because of their extensive use and success in solid-state electronic calculations [ 11, their potential computational advantages and their possible extension to the simulated annealing or the ab initio molecular dynamics simulation of rather complex systems [ 2,3]. DFI methods have been applied to atomic solids, liquids and clusters #*, to metal surfaces [ 51 and to isolated atoms and molecules [ 1,6,7], In the latter case, the agreement with experimental data and/or extensive CI calculations was found to be satisfactory when nonlocal correlation and exchange effects are incorporated in the expression of the energy functional [S-lo]. For some other recent develop-
t Permanent address: Institute of Chemistry, University of Wroclaw, Poland. st For recent applications to condensed atomic phases, see e.g. ref [4].
364
H-bonded complexes are important typical systems in which intermolecular interactions play a dominant role [ 12 1. From a theoretical point of view, long experience has shown that high-level electroncorrelation methods, associated with basis sets of sufficient quality, are required to describe the structure and interaction energy of the H-bond properly [ 13 1. Among H-bonded systems, [ H20H-OH*] + is a simple prototype for the study of proton-transfer reactions. It is also an example of a strong H-bond, an intermediate step between strong intramolecular bonds and the generally weak intermolecular interactions. Many computational insights and reliable results at the electron-correlation level are available for this complex [ 14,15 1. The purpose of this Letter is to test the up-to-date DFT methods for the (H502)+ complex. The corresponding results will be compared to those of a high-level electron-correlation calculation, in order to see to what extent the density functional approach is able to cope with a strong intermolecular interaction. A related study has been presented recently
0009-2614/93/$06.00
@ 1993 Elsevier Science Publishers B.V. All rights reserved.
Volume 208, number 54
CHEMICAL PHYSICS LETTERS
in the case of weak inter- and intra-molecular hydrogen-bonding interactions, including the weak Hbond in the water dimer [ 161.
2. Computational details For the DFT calculations, we have used the LCGTO-LSD (linear combination of Gaussian-type orbitals-local spin density) program deMON developed by Salahub and co-workers [ 8 1. This program goes beyond the local density approximation by incorporating nonlocal, density-gradient terms in the expression of the exchange and correlation energy functional, according to the parameterizations proposed by Perdew for exchange/correlation terms [8,17,18] and Becke for the exchange term only [191The electron-correlation calculations were performed by the Mailer-Plesset expansion method up to the fourth order [ 201. We used the GAUSSIAN 88 package [ 2 1] _ We chose the standard 6-3 1G( d, p) basis set, which was shown to be flexible enough to describe correctly the proton-transfer features in strongly H-bonded complexes [ 14,15 1. The Dm calculations were performed with the same 6-3 1G (d, p) basis set for comparison. We have also used an equivalent basis set specially optimized for DFT by Salahub and co-workers [ 221. It consists of a double zeta basis set with polarization (DZVP) for the oxygens, and a slightly augmented DZVPZ basis set, including four s-type primitive Cartesian Gaussians and one p-type polarization function, for the hydrogens. Previous studies of proton-transfer systems using conventional correlation methods have shown that the refinement of the H representation is important for the proper description of the transferring hydrogen, Hhb in fig. l_ For each calculation, we have tried different possible parameterizations of the
Fig. 1. Schematic picture of the complex and definition of the structural parameters.
18 June 1993
nonlocal exchange and correlation energy functional, using either the Perdew-Perdew (PP) [ 17,18 ] or Becke-Perdew (BP) [ 191 expressions, respectively. The various levels of DFT calculation are noted below according to the basis set and nonlocal energy functional employed, i.e. for example 6-3 1G (d, p ) / PP or DZVP/BP (note that 6.31G is in essence a DZVP basis set too. The difference in notation is only introduced for a convenient discrimination between basis sets),
3. Results 3.1. Comparison energies
of optimized structures and
Structures were determined in both methods by full geometry optimization. The results for the isolated subunits (H,O and H30+ ) are presented in table 1. In the case of water, for which the experimental structure is known accurately, the best agreement regarding the O-H bond length is obtained with the Msller-Plesset approach; the results are identical at the MP2 and MP4 level. The HOH angle is slightly too small. Nonlocal density functional theory gives a slightly too long OH bond length but an essentially correct HOH bond angle. The results are improved by choosing the DFT-optimized basis set. For protonated water, no experimental data are available, but the discrepancy between nonlocal DFT and MP2/MP4 is of the same order as for the pure water molecule. The protonation energy is found to depend significantly on the DFT method used, although it is rather independent of the order of the MP expansion in the electron-correlation approach. The results go from 180.7 kcal/mol for MP4 to 171.4 kcal/mol for DZVP/PP (see table 1). A schematic picture of the complex is given in fig. 1. Various results concerning bond lengths, bond angles, and interaction energies for the optimized geometry are collected in table 2. Except for the Hartree-Fock approximation, all the calculations, and particularly all the variants of the DFI’, predict a symmetrical configuration in which the proton is located at the midpoint of the O-O bond, this is the mark of a single-minimum potential in the proton 365
Volume 208, number $6
CHEMICAL PHYSICS LETTERS
18 June 1993
Table 1 Calculated geometries for isolated subunits (O-H in A and HOH in deg) at the different Meller-Plesset levels and for different DFT selections. A& is the protonation energy, in kcal/mol HF
MP2
MP3
MP4
6-3lG(d,p)/PP
DZVP/PP
DZVP/BP
Exp.
H20
O-H HOH
0.943 105.97
0.96 1 103.84
0.959 104.27
0.962 103.84
0.979 102.91
0.975 105.0
0.974 104.4
H,O+ O-H HOH A-%
0.961 114.14 - 179.9
0.979 112.59 - 179.8
0.976 113.06 - 180.7
0.978 112.13 - 180.7
0.989 115.3 - 178.8
0.992 112.9 - 171.4
0.990 113.9 -174.0
0.957 104.43
-
Table 2 Calculated geometries and interaction energies of the ( Hs02) + system. See fig. 1 for the definition of the different parameters. Bond lengths are in A, angles in deg, and energies in kcal/mol. In the case of the DF;T, no symmetry requirement was imposed a priori and a slight symmetry breakdown is observed in the results, due to a rather loose convergence criterion for the gradient ( 10m4au) HF O--O O-H, o-wb o-Hi,b H.OH. H,OHu 6, e, AE
2.409 0.953 0.949 1.078 111.7 109.2 144.7 160.8 - 33.06
MP2 2.396 0.969 0.969 1.198 108.4 108.4 141.0 141.0 - 38.89
MP3 2.389 0.966 0.966 1.189 108.6 108.6 141.2 141.2 -37.14
coordinate. A double-well topology is expected when the O-H-O complex is stretched beyond its equilibrium distance, or at equilibrium distance in the case of systems with a weaker H-bond, like [NH-N] + complexes [ 231. The DFT calculations and the MP2/MP4 ones lead to O-O equilibrium distances which are in essential agreement, around 2.40 A, although, again, DFT leads to slightly longer distances. On the other hand, we find that the interaction energy, defined with respect to the isolated subunits Hz0 and HsO+, depends rather strongly on the level of DFT calculation used. The best agreement with the MP4 value AE= - 37.9 kcal/mol is reached with theDZVP/PP selection (AE= -35.6 kcal/mol). The DZVP/BP selection yields a value of the interaction energy which is = 7 kcal/mol smaller than the MP4 result (A,!?= -31.6 kcal/mol), and is even x 1.5 kcal/mol smaller than the HF result. This appears as 366
6-3lG(d,p)/PP
MP4 2.395 0.969 0.969 1.198 108.5 108.5 141.0 141.0 -37.87
2.449 0.985 0.985 1.224 107.8 107.8 140.2 140.2 -44.69
DZVP/PP
DZVP/BP
2.445 0.983 0.983 1.223
2.427 0.980 0.980 1.213 108.1 108.1 137.9 137.9
108.0 108.0 147.0 147.0 - 35.6
-31.6
a deficiency of the BP parameterization for H-bonded systems, since the same trend was observed for the water dimer (see ref. [ 161). In this case, the PP interaction energy, A,??=-5.99 kcal/mol, is in satisfactory agreement with the experimental value, AE= -5.4 kcal/mol, and with the MP4 result of Frisch et al., AE= - 5.34 kcal/mol 1241. Besides, the BP parameterization with the same basis set leads to a value of -4.5 1 kcal/mol, which is * 20% smaller than the MP4/experimental one. This discrepancy is consistent with that observed here for (H50,)+. It should be noted also that the 6-31G(d, p)/PP selection, in which the basis set is optimized for standard correlation methods rather than for DFI’, gives a binding energy which is ~9 kcal/mol larger than that obtained with a reoptimized basis set, i.e., with the DZVP/PP selection, This is another manifestation of the extreme sensitivity of the system that
CHEMICAL PHYSICS LETTERS
Volume 208, number 5,6
we study on the level of approximation used and, within that, on the quality of the basis set. 3.2. Correlated proton-potential curve It is known that the proton-potential curve at a given O-O distance, including position of the minima and height of the barrier when the latter exists, is also sensitive to the basis set and level of approximation employed. It is useful therefore to test the DFI method for this particular problem. In order to explore a situation in which a substantial proton barrier is expected, we have chosen to fix the O-O distance at 2.74 A [ 151. The proton curves corresponding to the various methods are plotted in fig. 2. As already discussed in previous work [ 15 1, the barrier height decreases when going from the HF to the MP2 approximation. The higher-order MP4 result is close to the MP2 one. It was shown also that more advanced electron-correlation methods like coupled cluster with single and triple excitation corrections lead to potential curves in close agreement with the MP2/MP4 ones [ 15 J. The nonlocal DFT calculation gives a much lower barrier height: 1.7 kcal/mol for DZVP/PP instead of 5.1 kcal/mol for MP4. The latter number agrees well also with the 4.9 kcal/mol
lSJune1993
obtained by using the much larger TZZP basis set [ 151. It appears therefore that the DFT method overestimates the correlation effects in the proton potential.
4. Conclusion The H-bonded complex studied in this work is a good model system for testing ab initio methods. In particular, the proton curve appears to result from the compensation of fine effects and it is sensitive to the level of approximation used in the calculation. Past experience with electron-correlation methods, which is confirmed in this work, has shown that the Msller-Plesset expansion is well converged at the second order for this type of intermolecular complex. Accurate structural data, interaction energies, and proton potentials can be obtained in this way; they are in close accordance with the results of more elaborate methods. Referring to those results, we have found that up-to-date density functional methods, incorporating nonlocal correlation and exchange effects, may lead to accurate values of the H-bonding interaction energy. Nevertheless, the agreement with the MP2 / MP4 results is dependent on which parameterization of the correlation/exchange functional is used. The double-well shape of the proton potential, which is a characteristic feature of the system at large O-O separations, is respected but with too small a barrier height and too long an OH bond. The corrections brought by the electronic correlation are found to be too strong. More work on related systems will be needed to understand the ability of the new DFT methods to describe properly the intermolecular interactions and, among them, the H-bonding. This is the price to pay before being able to use them for the dynamical simulation of chemical events in H-bonded media, a potential future application which is in many minds.
Acknowledgement Fig. 2. Proton-potential curve for for various levels of approximation sion: (- - - ) SCF, (-- -) MP2; (-) represents the DFT results with the
an O-O distance of 2.74 A of the Meller-Plesset expanMP4. The upper solid curve DZVP/PP selection.
ZL acknowledges a visitor grant from the French Ministry of Research and Technology. This work was also partially supported (ZL) by the State Committee for scientific Research of Poland (207449101). 361
Volume 208, number 5,6
CHEMICAL PHYSICS LETTERS
References [ 1] R.O. Jones andO. Gunnarsson, Rev. Mod. Phys. 6 1 ( 1989) 689. [2]R. Car and M. Parrinello, Phys. Rev. Letters 55 (1985) 2471. [3] D. Remler and P. Madden, Mol. Phys. 70 ( 1990) 921. [ 41 R. Car and M. Paninello, Phys. Rev. Letters 60 ( 1988) 204; P. Ballone, W. Andreoni, R. Car and M. Parrinello, Phys. Rev. Letters60 (1988) 271; G. Galli, R.M. Martin, R. Car and M. Parrinello, Phys. Rev. Letters 62 (1989) 555; 63 (1989) 988. [5] V. Russier, D.R. Salahub and C. Mijoule, Phys. Rev. B 42 ( 1990) 5046; C. Mijoule and V. Russier, Surface Sci. 254 ( 1991) 329. [6] J.K. Labanowksi and J.A. Andzelm, Density functional methods in chemistry (Springer, Berlin, 1991) . [7]T. Ziegler, Chem. Rev. 91 (1991) 651. [8]A. St-Amant and D.R. Salahub, Chem. Phys. Letters 169 ( 1990) 387, and references therein; P. Mlynarski and D.R. Salahub, Phys. Rev. B 43 (1991) 1399. [9] J.A. Pople, P.M.W. Gill and B.G. Johnson, Chem. Phys. Letters 197 (1992) 499; 199 (1992) 157; B.G. Johnson, P.M.W. Gill and J.A. Pople, J. Chem. Phys. 97 (1992) 7846. [lo] C.W. Murray, G.J. Laming, N.C. Handy and R.D. Amos, Chem. Phys. Letters 199 (1992) 551. [ 111 P.M. Boerrigter, G. TeVelde and E.J. Baerends, Intern. I. Quantum Chem. 33 ( 1988) 87; B.J. Delley, J. Chem. Phys. 92 ( 1990) 508;
368
18Junel993
AD. Becke and R.M. Dickson, I. Chem. Phys. 92 (1990) 3610. [ 121 H. Ratajcrak and W.J. Orville-Thomas, eds., Molecular interactions (Wiley, New York, 1980). [ 13] S. Scheiner, in: Theoretical models of chemical bonding, ed. Z.B. Maksic (Springer, Berlin, 1991). [ 141 S. Scheiner, Accounts Chem. Res. 18 (1985) 174; K. Luth and S. Schemer, J. Chem. Phys. 97 (1992) 7507, 7519. [ 151 Z. Latajka and S. Schemer, J. Mol. Stmct. THEOCHEM 234 (1991) 373. [ 161 F. Sim, A. St-Amant, I. Papai and D.R. Salahub, J. Am. Chem. Sot. 114 (1992) 4391. [ 171J.P. Perdew, Phys. Rev. B 33 (1986) 8822. [ 181 J.P. Perdew, Phys. Rev. Letters 55 (1985) 1665; J.P. Perdew and Y. Wang, Phys. Rev. B 33 ( 1986) 8800. [ 191A.D. Becke, Phys. Rev. A 38 (1988) 3098. [ 201 R. Krishnan and J.A. Pople, Intern. J. Quantum Chem. 14 (1978) 91; R. Krishnan, M.J. Frisch and J.A. Pople, J. Chem. Phys. 72 ( 1980) 4244. [ 2 1 ] M.J. Friseh, J.S. Binkley, H.B. Schlegel, K. Raghavachari, C.F. Melius, J.J.P. Stewart, F.W. Bobrowicz, C.M. Rohlfing, L.R. Kahn, D.J. DeFrees, R. Seeger, R.A. Whiteside, D.J. Fox, E.M. Fleuder and J.A. Pople, GAUSSIAN 88 (Gaussian Inc., Pittsburgh, PA, 1988). [22] F. Sim, D.R. Salahub, S. Chin and M. Dupuis, J. Chem. Phys. 95 (1990) 4317. [23] L. Jaroszewski, B. Lesyng, 5.3. Tanner and J.A. McCammon, Chem. Phys. Letters 175 (1990) 282. [24] M.J. Frisch, J.E. de1 Bene, J.S. Binkley and H.F. Schaefer III, J. Chem. Phys. 84 ( 1986) 2279.