Density functional study of a highly correlated molecule, oxygen

J . Phys. Chem. 1991, 95,9197-9202

9197

Density Functional Study of a Highly Correlated Molecule, FOOF David A. Dixon,* Central Research & Development,t Du Pont, P.O. Box 80328, Experimental Station, Wilmington, Delaware 19880-0328

Jan Andzelm, George Fitzgerald, and Erich Wimmer

Downloaded by KTH ROYAL INST OF TECHNOLOGY on August 24, 2015 | http://pubs.acs.org Publication Date: November 1, 1991 | doi: 10.1021/j100176a031

Industry, Science and Technology Department, Cray Research, Inc., Eagan, Minnesota 551 21 (Received: October 9, 1990)

The molecular and electronic structure of FOOF has been calculated in the density functional framework with both numerical and Gaussian basis sets. The geometry is well predicted by the local density functional method with the 0-0bond length calculated to better than 0.01 A and the 0-F bond length to better than 0.015 A. The predicted vibrational frequencies are in reasonable agreement with experiment with the frequencies all greater than the experimental values. The largest error in the calculated frequencies is found for the symmetric 0-0-F bend, which is predicted to be about 20% too large. An analysis of the force constants in terms of simple internal coordinates reveals large couplings between the stretches. The rotational barriers are calculated and shown to be larger than the dissociation energy leading to loss of F. The computed dissociation energies for loss of one and two F atoms are found to be larger than the experimental values. Whereas geometries and frequencies are well described within the local density approximation, inclusion of nonlocal corrections is required for the dissociation energies to have even qualitative agreement with experiment.

Introduction

One of the most difficult systems to describe in terms of molecular orbital methods is the molecule F0OF.l This molecule has a very short 0-0 bond like that in 02,and it is 0.26 A shorter than the "comparable" 0-0 bond in HOOH.2 FOOF also has very long 0 - F bonds which are 0.17 A longer than the bonds in The molecule has a nonplanar structure with an FOOF dihedral angle near 90°. FOOF is a very powerful fluorinating agent which can convert compounds containing actinides such as U to volatile fluorides at room temperature. Such fluorides can only be generated at this low temperature by reaction with FOOF or KrF2, and thus FOOF has potential use in the processing of nuclear materiakg A simple model for the electronic structure of FOOF in generalized valence bond terms4 can be given as follows considering only the valence p orbitals. The orbitals on 0 are like the orbitals in 02,

U

U

0

0

U 0

y\.Jq& .. F

0 2

FOOF

The molecule has a half-occupied p orbital on each 0 atom, and these orbitals are orthogonal to each other. If the 0-0 axis is the z axis, then there is a single electron in a px orbital on 0 atom 1 and a single electron in a pu orbital on 0 atom 2. The remaining p electrons are distributed as the 0-0 u bond and as a doubly occupied p,, on 0 atom 1 (denoted by a single bracket) and a px pair on 0 atom 2. The F atoms then form weak bonds (denoted by a dashed bracket) with the available single electrons but do not perturb the 0-0 bonding. They, of course, do change the overall molecular spin of the molecule as O2is a triplet and FOOF 'Contribution no. 5639.

0022-3654/9 1 /2095-9197$02.50/0

is a singlet. In the model of L i p s c ~ m b , lthe - ~ result of these interactions is the formation of two three-center, two-electron bonds which gives a significantly stronger 0-0 bond than in HOOH and correspondingly weaker 0-F bonds. There have been a large number of theoretical studies612 aimed at elucidating the bonding in this unique molecule. Although the Hartree-Fock (HF) description usually gives excellent geometries,13 it completely fails to give the proper description of the bonding in this molecule. This was a surprise as one of the early successes of H F calculations with polarized double-{ basis sets was the correct prediction of the electronic structure of HOOH.I4 As shown in Table I, the 0-0bond is almost 0.1 A too long and the 0-F bond is 0.2 8, too short at the S C F level with adequate basis sets. Smaller basis sets give even worse results. This is in contrast to the usual result of a bond distance that is 0.01-0.02 8, too short when DZ+P basis sets are used. Traditional methods for treating the correlation problem do little to resolve the differences between the calculations and experiment. Configuration interaction (CI) methods8 give similar results to the HartreeFock methods. Coupled cluster methods8 lead to a significant increase in the 0-F bond but still predict a value almost 0.1 8, too short as compared to experiment. The MP-2 methodI2 apparently gives a qualitatively correct structure with an 0-F bond longer than experiment and an 0-0bond shorter than experiment. However, (1) Jackson, R. H. J. Chem. SOC.1962,4585; unpublished work by W. N. Lipscomb quoted therein. (2) Harmony, M. D.; Laurie, V . W.;Kuczkowski, R. L.; Schwendeman, R. H.; Ramsay, D. A.; Lovas, F. J.; Lafferty, W. J.; Maki, A. G.J . Phys. Chem. Ref.Data 1979, 8, 619; see p 640. (3) (a) Streng, A. Chem. Reu. 1%3,63,607. (b) Asprey, L. B.; Kinkead, S. A.; Eller, P. G. Nucl. Technol. 1986, 73, 69. (c) Asprey, L. B.; Eller, P. G.;Kinkead, S. A. Inorg. Chem. 1986, 25, 670. (4) (a) Goddard, W. A., 111; Dunning, T. H., Jr.; Hunt, W. J.; Hay, P. J. Acc. Chem. Res. 1973,6, 368. (b) Moss, B. J.; Bobrowicz, F. W.; Goddard, W. A., Ill. J. Chem. Phys. 1975, 63, 4632. ( 5 ) Cotton, F. A.; Wilkinson, G. Aduanced Inorganic Chemistry, 4th ed.; Wiley-Interscience: New York, 1980. (6) Lucchese, R. R.; Schaefer, H. F., Ill; Rodwell, W. R.; Radom, L.J. Chem. Phys. 1978, 68, 2507. (7) Clabo, D. A., Jr.; Schaefer, H. F., 111. Int. J . Quantum Chem. 1987, 31, 429. (8) Lee, T. J.; Rice, J. E.; Scuseria, G. E.; Schaefer, H. F., 111. Theor. Chim. Acta 1989, 75, 81. (9) Ahlrichs, R.; Taylor, P. R. Chem. Phys. 1982, 72, 287. (IO) Rohlfing, C. M.; Hay, P. J. J. Chem. Phys. 1987,86, 4519. ( I I ) Raghavachari, K.; Trucks, G. W. Chem. Phys. Lett. 1989,162, 511. (12) Mack, H.-G.; Oberhammer, H. Chem. Phys. Lett. 1988, 145, 121. (13) Hehre, W. J.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio Molecular Orbital Theory; John Wiley & Sons: New York, 1986. (14) Dunning, T. H.; Winter, N. W. Chem. Phys. Lett. 1971, 1 1 , 194.

0 199 1 American Chemical Society


9198

The Journal of Physical Chemistry, Vol. 9.5, No. 23, 1991

careful studies8 have shown that the MP-2 results are basis set dependent, and an improvement in the basis set leads to worse agreement until a very large TZ+ZP basis set augmented with f functions is reached. A CEPA9 treatment with a polarized basis set gives reasonable agreement with experiment with an 0-F bond distance that is too long by 0.06 A. The best calculated structures are those derived from a multiconfiguration S C F plus CI calculationlo and a calculation including higher excitations, QCISD(T)." However, quite small basis sets were used for these latter two calculations. In only two cases has the force field been calculated for FOOF.7," These are at the SCF level with the 6-31G*," DZ+P,' and DZ+P + diffuse7 basis sets and at the QCISD(T) level with the 6-31G* basis set." Since the H F calculations cannot predict the geometry correctly, they do an even worse job on predicting the vibrational frequencies. The 0-0 stretch is predicted to be slightly lower than experiment15 and of comparable magnitude to the 0-F stretches, whereas experimentally the 0-F stretches are about half the frequency of the 0-0 stretch. The predicted frequencies are improved at the QCISD(T) level, but the 0-0 stretch is now 200 cm-I below the experimental value. Clearly, this is a difficult structure to predict. Because computational resources are not growing with the same scaling factor as in traditional molecular orbital calculations, especially for methods that include correlation in the wave function ( Pwith m > 4,where N is the number of basis functions), we are interested in the behavior of new computational methods that do not scale with such a large exponent but are based on rigorous theory with no empirical parameters.I6 One promising method is based on density functional theory (DFT) which has been used with great success in solid-state p h y s i ~ s , ' especially ~-~~ within the local density approximation. It has the advantage of scaling as N 3 and, as in the case of the Hartree-Fock method, may have a lower exponent depending on implementation. Because of its form," it does include correlation effects in its calculation of the energy. We are interested in determining whether this method can be used successfully in molecular calculations and are testing its ability to predict molecular properties for a wide range of systems.2',22 Because FOOF is such a difficult molecule to treat with the usual theoretical methods, we have calculated its geometry and force field using the local density functional (LDF) method.

Methods The calculations described below were done with the program systems D M o I ~and ~ DGau~s.*~ DFI' theory is based on a theorem of Hohenberg and Kohn2s which states that the total energy E , ( 1 5 ) Kim, K. C.; Campbell, G. M. J . Mol. Struct. 1985, 129, 263. (16) (a) Dixon, D. A. In Science and Engineering on Cray Supercompu-

ters; Proceedings of the Third International Symposium, Cray Research: Minneapolis, MN, 1987; p 169. (b) Dixon, D. A.; Capobianco, P. J.; Mertz, J . E.; Wimmer, E. Science and Engineering on Cray Supercomputers; Proceedings of the Fourth International Symposium, Cray Research: Minneapolis, MN, 1988; p 189. (c) Dixon, D. A,; Van-Catledge, F. A. f n t . J . Supercomputer Appl. 1988, 2 (2), 62. (17) Parr, R. G.; Yang, W. Density Functional Theory of Atoms and Molecules; Oxford University Press: New York, 1989. (18) Salahub, D. R. I n A b Initio Methods in Quantum Chemistry-If; Lawley, K. P., Ed.; Wiley & Sons: New York, 1987; p 447. (19) Wimmer, E.; Freeman, A. J.; Fu, C.-L.; Cao, P.-L.; Chou, S.-H.; Delley, B. In Supercomputer Research in Chemistry and Chemical Engineering Jensen, K . F., Truhlar, D. G., U s . ; ACS Symposium Series No. 353; American Chemical Society: Washington, DC, 1987; p 49. (20) Jones, R. 0.; Gunnarsson, 0. Reu. Mod. Phys. 1989, 61, 689. (21) Dixon, D. A.; Andzelm, J.; Fitzgerald, G.; Wimmer, E.; Delley, B. In Science and Engineering on Cray Supercomputers; Pitcher, E. J., Ed.; Computational Mechanics Publications: Southampton, England, 1990: p 285. (22) Dixon, D. A.; Andzelm, J.; Fitzgerald. G.; Wimmer, E.; Jasien. P. I n Density Functional Methods in Chemistry; Labanowski, J.. Andzelm, J., Eds.; Springer-Verlag: New York, 1991; Chapter 3, p 33. (23) Delley, B. J . Chem. Phys. 1990, 92, 508. DMol is available commercially from BIOSYM Technologies, San Diego, CA. (24) (a) Andzelm, J.; Wimmer, E.; Salahub, D. R. In The Challenge oj d and f Electrons: Theory and Computation; Salahub. D. R.. Zerner. M. C Eds.; ACS Symposium Series 394; American Chemical Society: Washington, DC, 1989; p 228. (b) Andzelm, J. W. I n Density Functional Methods in Chemistry; Labanowski, J., Andzelm, J . W., Eds.; Springer-Verlag: New York, 1991; Chapter I I , p 155. (25) Hohenberg, P.; Kohn, W. Phys. Rec. E 1964. 136, 184.

Dixon et al. is a function of the charge density, p , and can be written as where T is the kinetic energy of the noninteracting electrons of density p , U is the classical Coulomb electrostatic energy, and E,, includes all of the many-body contributions to the energy. For a molecular system, the most important terms in E,, are the exchange energy and the correlation energy. The density is given as a sum over the squares of the occupied molecular orbitals as in standard molecular orbital treatments. The first term in eq 1 is the kinetic energy of the electrons. The second term can be written as (3) where the first term in eq 3 corresponds to the nuclear attraction energy, the second term to the electron-electron repulsion, and the third term to the nuclear-nuclear repulsion. These quantities and the kinetic energy term can all be evaluated by use of straightforward techniques. Up to this point, the solution of eq 1 is exact, and it is in the last term of eq 1 that the local density approximation is introduced. It has been found that a good approximation for the final term can be taken from the exchange-correlation energy of the uniform electron gas.26s27The assumption is that the charge density varies slowly on the scale of exchange and correlation effects, and one can obtain E,, by integrating the uniform electron gas result: (4)

with exc being the exchange-correlation energy per particle in the uniform electron gas. In the calculations described below, the form derived by von Barth and Hedin2* is used in the DMol program and that of Vosko, Wilk, and N ~ s a i in r ~the ~ DGauss program. The value of E, is determined by optimizing the variations in E, with respect to variations on p subject to the usual orthonormality constraints for molecular orbitals. This leads to a set of coupled differential equations which must be solved iteratively until the density and the energy are converged just as in standard molecular orbital treatments. Equation 4 is the first term of the more general formula valid in the case of slowly varying density.25 It contains additional gradient density terms as shown:

substantial literature exists on nonlocal (gradient) corrections for both the exchange and correlation energy DMol employs numerical functions for the atomic basis sets. The atomic basis functions are given numerically on an atomcentered, spherical-polar mesh. The radial functions from the numerical solution of the atomic LDF equations are stored as sets of cubic spline coefficients so that the radial functions are piecewise analytic, a necessity for the evaluation of gradients. The use of exact spherical atom results offers some advantages. The molecule will dissociate exactly to its atoms within the LDF framework, although this does not guarantee correct dissociation energies. Because of the quality of the atomic basis sets, basis set superposition effects should be minimized and correct behavior at the nucleus is obtained. Since the basis sets are numerical, the various integrals arising from the expression for the energy equation need to be evaluated Hedin, L.; Lundquist, B. I . J . Phys. C 1971, 4, 2064. Ceperley, D. M.; Alder, B. J . Phys. Reu. Lett. 1980, 45, 566. van Barth, U.; Hedin, L. J . Phys. C 1972, 5 , 1629. Vosko, S. H.; Wilk, L.; Nusair, M. Can. J . Phys. 1980, 58, 1200. (30) (a) Becke, A. D. In The Challenge of d and f Electrons: Theory and Computation; Salahub, D. R., Zerner, M. C., Eds.; ACS Symposium Series 394; American Chemical Society: Washington, DC, 1989; p 166. (b) Becke, A. D. Int. J Quantum Chem. Symp. 1989, 23, 599. (31) Perdew, J . P.Phys. Reu. B 1986, 33, 8822. (32) Becke. A. D. J . Chem. Phys. 1986,84,4524. Stoll, H.; Pavlidou, C . M E : Preuss, H Theor. Chim. Acta 1978, 49, 143. (26) (27) (28) (29)

The Journal of Physical Chemistry, Vol. 95, No. 23, 1991 9199

Density Functional Study FOOF TABLE I: Molecular Geometry for FOOP method R(F-0) R(O-O) 0 1.217 109.5 expt 1.575


LDF LDF LDF LDF LDF LDF LDF LDF LDF HF CISD CCSD MP2 CEPA CAS+CI+D QCISD(T)

1.563 1.563 1.561 1.561 1.560 1.561 1.560 1.563 1.564 1.363 1.407 1.482 1.604 1.63 1.600

1.610

*Bond distances in in the optimization.

1.220 1.218 1.212 1.211 1.218 1.218 1.218 1.210 1.206 1.310 1.301 1.278 1.178 1.22 1.233 1.234

ref

7

87.5 110.3 88.7 110.4 88.6 110.5 88.3 110.5 88.3 110.5 88.2 110.4 88.2 110.5 88.2 110.5 88.3 110.6 88.3 106.2 85.1 106.2 85.9 109.5’ 87.5’ 110.0 88.2 109.5’ 87.5’ 109.5’ 87.5’ 107.8’ 86.8’

1 DMol DMol, DG-I DG-2 DG-3 DG-4 DG-5 DG-6 DG-7 7 8 8 8

L,,, = L

9

IO 11

A. Bond angles in deg. ’Constrained parameter

over a grid. The integration points are generated in terms of angular functions and spherical harmonics. The number of radial points NR is given as where Z is the atomic number. The maximum distance for any function is 12 au. The angular integration points No are generated at the N R radial points to form shells around each nucleus. The value of No ranges from 14 to 302 depending on the behavior of the density.” The Coulomb potential corresponding to the electron repulsion term could be solved by evaluation of integrals. However, since the method is based on the density, it was found to be more appropriate to solve Poisson’s equation for the charge density directly: -VZVe(r) = 4i7e2p(r)

(7) All of the DMol calculations were done with a double numerical (DN) basis set augmented by polarization functions. This can be thought of in terms of size as a polarized double-( basis set. The quality of the DN numerical basis sets is such that addition of one set of polarization functions is usually all that is required to give good predictions of energies and geometries for most molecules.34 However, because we are using exact numerical solutions for the atom, this basis set is of significantly higher quality than a normal molecular orbital polarized double-( basis set. The fitting functions contain an angular momentum number 1 greater than that of the polarization function; thus a value of L = 3 was used for the fitting function as d polarization functions were included. An efficient analytical approach to variationally fit the density has been developed p r e v i o ~ s l y . ~The ~ - ~analytical ~ implementation of the LDF method requires that Gaussian basis sets be used. Additionally, this allows one to build on the wealth of experience gained from traditional molecular orbital calculations. In the DGauss program, Cartesian Gaussians are used as primitives and are contracted in the same way as in Hartree-Fock methods.36 The Gaussian basis sets employed in the DGauss calculations are summarized in Table I1 together with the grid information. Although the experience about basis sets gained from HartreeFock molecular calculations is invaluable in terms of defining the basis set size and the need for additional polarization and diffuse functions, Hartree-Fock basis sets used in an LDF calculation (33) This grid can be obtained by using the FINE parameter in DMol. (34) Delley. B.; Freeman, A. J.; Weinert, M.; Wimmer, E. Phys. Reu. B 1983, 27, 6509. (35) Dunlap, B.;Connolly, J.; Sabin, J. J. Chem. Phys. 1979, 71, 3396. (36) (a) Huzinaga, S.; Andzelm, J.; Klobukowski, M.; Radzio, E.; Sakai,

Y.; Tatasaki, H. Gaussian Basis Sets f o r Molecular Calculations; Elsevier: Amsterdam, 1984. (b) Dunning, T. H.,Jr.; Hay, P. J. In Methods of Elecfronic Sfrucfure Theory; Schaefer, H.F., 111, Ed.; Plenum Press: New York, 1977; Chapter 1.

TABLE II: Basis Sets for DGauss Calculations basis size contraction fitting set DG-I DG-2 DG-3 DG-4 DG-5 DG-6 DG-7

9/5/1

91511 10/6/1

10/6/1 10/6/1 10/6/2 10/6/2

621/41/1 621/41/1 721/51/1 721/51/1 721/51/1 721/51/2 721/51/11

grida

TZ TZ TZ

M F

TZ

F F

QZ QZ QZ

M F F

“ M = medium, F = fine.

may exhibit large basis set superposition errors (BSSE). In order to minimize BSSE, LDF optimized basis sets were developed,” and these are the ones presented below. The basis sets are all contracted from a larger primitive set to a [3,2,1] set following the ideas of H ~ z i n a g a . ~ ~ ~ The approach used in DGauss is an analytical implementation of the LDF method as opposed to the purely numerical approach taken in DMol. In the analytical approach, the density can be variationally fit leading to exact Coulomb The remaining exchange-correlation energy term is a smooth function of the density and can be accurately fit on an adaptive set of grid points. It is this numerical evaluation of the exchange-correlation energy term that leads to the differences in the geometries obtained from the minimum in the energy and the one with zero gradient as discussed below. An additional auxiliary basis set is required to fit the electron density as obtained from the molecular orbitals and the exchange-correlation potential and energy. This additional basis set is constructed from a set of atom-centered Gaussian-type basis functions with s, p, and d character. It is found that even-tempered expansions of triple-{ (TZ) and quadruple-( (QZ) quality are adequate to reproduce the electron density. The coefficients in this density expansion are determined from a variational condition that requires that the residual (second-order) Coulomb energy term arising from the difference between the exact and the fitted density be minimized. In practice, this leads to analytic expressions for the fitting coefficients involving Coulomb-type three-index, two-electron integrals. The exchange-correlation potential is also expanded in Gaussian-type basis functions. However, the determination of the expansion coefficients is carried out numerically by using a grid similar to the one described above for the DMol program. The grid selection in DGauss is derived from an adaptive procedure based on the approximation of the exchange-correlation energy in the angular shell around each The radial distribution of points is accomplished by using the method of B e ~ k e .Once ~~ the Gaussian expansions of the exchange-correlation potential are determined, the matrix elements involving the three-index, oneelectron integrals are calculated. To this end, the DGauss approach leads to an method, which in practice can be made close to i V by using sparse, direct matrix algorithms for the integral calculations, numerical integration, and density synthesis. This can be compared to the four-index, two-electron integrals required for the Hartree-Fock method which formally scale as N4. All three-index, two-electron integrals are calculated by using an algorithm based on the work of Obara and Saika.40 Geometries were determined by optimization using analytic gradient First derivatives in the LDF framework ~~

(37) 520. (38) 6371. (39) (40) (41)

~

Andzelm, J.; Radzio, E.; Salahub, D. R. J . Comput. Chem. 1985.6, Fournier, R.; Andzelm, J.; Salahub, D. R. J. Chem. Phys. 1989,90,

Becke, A. J. Chem. Phys. 1988.88, 2547. Obara, S.; Saika, A. J. Chem. Phys. 1986,84, 3963. (a) Komornicki, A.; Ishida, K.; Morokuma, K.;Ditchfield, R.; Conrad, M. Chem. Phys. Lett. 1977.45, 595. (b) Pulay, P. In Applications of Electronic Structure Theory; Schaefer, H. F., 111, Ed.; Plenum Press: New York, 1977; p 153. (42) Delley, B. In Density Functional Methods in Chemistry; Labanowski, J. K., Andzelm, J. W., Eds.; Springer-Verlag: New York, 1991; Chapter 7, p 101.

9200 The Journal of Physical Chemistry, Vol. 95, No. 23, 1991 TABLE 111: Vibrational Frequencies (cm-I) for FOOF assign 0-0 s 0-F


HF/6-31GZ” QC ISD( T) I I LDF DMol L D F DG-I L D F DG-2 L D F DG-3 L D F DG-4 L D F DG-5 L D F DG-6 L D F DG-7 exptIs

I161 I010 1334 1384 1389 1327 1374 1445 1419 1464 1210

s

Dixon et al.

F-0-0 s

1144 756 749 738 734 733 734 747 735 725 630

can be calculated efficiently and only take on the order of 3-4 S C F iterations or 10-30% of the calculation of the energy. This is in contrast to the evaluation of the derivatives in traditional molecular orbital methods which usually takes at least a factor of 1-2 times the evaluation of the energy. There are two problems with evaluating gradients in the LDF framework which are due to the numerical methods that are used. The first is that the energy minimum does not necessarily correspond exactly to the point with a zero derivative. The second is that the sum of the gradients may not always be zero as required for translational invariance. These tend to introduce errors on the order of 0.001 8, in the calculation of the coordinates if both a reasonable grid and basis set are used. This gives bond lengths and angles with reasonable error limits. The frequencies were determined by using numerical differentiation of the gradients. A two-point difference formula was used, and a displacement of 0.01 au was used for the DMol calculations and 0.02 au for the DGauss calculations. It is important to note that analytic second derivatives are available for H a r t r e e - F ~ c kmethods ~~ and some correlated methods whereas analytic second derivatives have not yet been implemented for density functional methods.45 This means that, at present, second derivatives in the traditional molecular orbital methods may be done more quickly than in the LDF methods, and this discrepancy increases with molecular size.

Results Geometry. The geometry of FOOF predicted by the LDF method is given in Table I. The numerical LDF geometry is in excellent agreement with the experimental structure. The calculated F-0 bond length is 0.014 8, too short, and the 0-0bond length is 0.003 8, too long. The result of the calculated values being the same or less than the experimental values is in contrast to our previous restricted set of results for LDF where the bond distances tend to be too long. This will affect the frequencies as described below. The LDF method also predicts both the bond and torsion angles correctly. Variation of the quality of the mesh to a smaller number of mesh points did not affect the geometry. The use of a value of L = 2 for the fitting functions instead of L = 3 gave essentially the same geometry as with L = 3 for both the standard and fine mesh sets. The only real variation is in the 0-0 distance which becomes 0.002 8, shorter. The results from the Gaussian LDF calculations (see Table 11 for a description of the basis sets) also show in essence the same level of agreement with experiment although a sensitivity to the quality of the basis set is found. The 0 - F bond distances are slightly short just as found by the numerical LDF. The 0-0 bond distance is also too short in contrast to the numerical LDF which produces a slightly long 0-0 bond. The angles at the Gaussian LDF level show the same good agreement with experiment as found with the numerical LDF. As found previously, the geometry is not particularly sensitive to the grid or the fitting functions at the Gaussian LDF (43) Versluis, 1.; Ziegler, T.J . Chem. Phys. 1988, 88, 3322. (44) King, H . F.; Komornicki, A. In Geomerrical Derioatioes of Energy Surfaces and Molecular Properties; Jargensen, P., Simons, J., Eds.; NATO AS1 Series C, Vol. 166; D. Reidel: Dordrecht, 1986; p 207. Pople, J. A,; Krishnan, R.: Schlegel. H. B.; Binkley, J. S. Inr. J . Quantum Chem. Symp. 1979, 13. 225. (45) Expressions for second derivatives for Gaussian-based LDF have recently been presented. Fournier, R. J . Chem. Phys. 1990, 90, 5422.

556 307 459 457 459 454 456 448 460 413 355

F-0-0-F 210 164 20 1 202 208 200 213 197 189 224 202

s

F-0 a 1135 688 715 668 719 686 692 678 692 652 614

F4-0 a 709 510 529 528 528 520 524 522 533 535 470

level. At the level of a double-[valence basis set augmented by polarization functions (DZVP), the results are converged to 0.008 8, in the bond distances and 0.3’ in the angles. Within these limits, the results are not sensitive to the choice of grid, fitting set size, or number of contracted coefficients in the core (if six or more functions are used). The results found for DG-7 which has two uncontracted d functions indicate that care may need to be taken in the design of more uncontracted sets and that better optimization of the d exponents may be needed. Thus the LDF method provides a good description of the molecular geometry of FOOF. Vibrational Spectrum. The calculated vibrational spectrum provides an even more stringent test of the method, and the results are summarized in Table 111. The theoretical vibrational transitions should be greater than the experimental valuesI5 as the calculated spectrum is harmonic whereas the experimental values include anharmonicity. This is exactly what is found by both the numerical and Gaussian LDF calculations with the exception of the torsion which should have a small anharmonicity correction. The 0-0 stretch is about 7% too large with the numerical LDF and is 13-15% too large with the Gaussian LDF. The 0 - F stretches are about 15-1 8% too large at the LDF level. These differences are consistent with the geometry predictions which show that LDF predicts the 0-0 bond to be almost the experimental value, whereas the 0-F bonds are predicted to be too short. The biggest discrepancy between theory and experiment is for the symmetric 0-0-F bend which is predicted to be 22% too high by the LDF method except for one basis set which halves the error. This result is probably consistent with large coupling of the 0-0 and symmetric 0-F stretches through the bend as will be discussed below. The QCISD(T) results” show an 0-0 stretch that is too small by 17% and an 0 - F symmetric stretch that is similar to the LDF values. The QCISD(T) asymmetric stretch shows better agreement with experiment than do the LDF values, but the splitting of the symmetric and asymmetric 0 - F stretches is too large. The LDF method predicts a much smaller splitting which is similar to the experimental value. As might be expected, the QCISD(T) value for the symmetric 0-0-F bend is too low and the value for the asymmetric 0-0-F bend is comparable to the LDF results. The torsion frequency at the QCISD(T) level is surprisingly too low by about 40 cm-I. The H F results7,” as noted above cannot account for any of the vibrational transitions except for the torsion. The comparison of the results with different basis sets, grids, and fitting functions for the Gaussian LDF calculations provides further information about the sensitivity of the potential energy surface to small changes in the electronic character of the density. The 0-0 stretch is sensitive to the fit as a better representation of the core is employed. It shows an increase with both improved fitting function and more accurate grid as shown by comparing the results of DG-3, DG-4, and DG-5. Improving the form of the d polarization function leads to a decrease in the 0-0 stretch (compare DG-6 with DG-5) and also a decrease in the predicted torsion frequency. Splitting the d polarization function (DG-7) leads to an increase in the 0-0stretch and in the torsion but leads to a decrease in the 0 - F asymmetric stretch. In the DGauss calculations of the vibrational frequencies, the rotations and vibrations have not been transformed out, which could slightly influence the lower frequencies. The DGauss frequencies display a slight sensitivity to the convergence criteria for both the SCF

The Journal of Physical Chemistry, Vol. 95, No. 23, 1991 9201

Density Functional Study FOOF TABLE IV: Force Constants (mdyn/A) for FOOF in Terms of Internal Coordinates

internals 1 2


3

FI 1 FI2 F14 F16 F22 F23 F24 F25 F26 F44 F45 F46 F66

internals

0-0 stretch

4 5 6

0-F stretch 0 - F stretch LDF, DMol 8.99 1.18 0.28 -0.01 2.50 0.09 -0.12 0.09 0.03 1.27 -0.07 0.02 0.33

TABLE V: Geometrical Parameters for Planar FOOF and Rotation Barriers

QCISD(T)I1 5.87 1.07 0.46 -0.50 1.98 -0.12 0.07 -0.25 -0.02 1.31 -0.08

OQFbend 04-Fbend F-0-O-F torsion expt" expt" expt4' 10.25 1.27 1.71

0-0-F

7.00

GeOmetriesO trans

DMol 1.198 1.736 114.9

DG-5 1.194 1.737 114.9

1.70

1.85

0.14

0.122

0.18

LDF DMol LDF DG-5 LDF DG-5 a T

1.19

1.11

1.05

0.47

0.014 0.27

0.11 0.24

cis

cis

DMol 1.202 1.736 114.9

DG-5 1.200 1.724 1 13.9

Energies (kcal/mol) method

1.50

-0.1 1

0.27

7.14

parameter 0-0 0-F

trans

+ B-P

trans 40.2 40.3 36.2

Bond distances in A. Bond angles in deg.

= Oo for cis.

cis 41.5 41.8 39.1 T

= 1 80° for trans and

is 1.1 kcal/m01.4~ It is of interest to examine the rotation barrier if FOOF especially since the experimental estimate of the barrier is quite high, 30 k~al/mol.l,~'This estimate is based on fitting the torsion potential to a cosine expansion. The high barrier is anomalous since the dissociation energy for loss of F (FOOF FOO' + F')is only 20 k c a l / m ~ l . ~ ~

-

and the gradient. The DGauss gradients are good to no better than lO-5 au, and this may slightly affect the vibration frequencies. The largest error is expected for the lower frequencies and should not be greater than 3%. In order to gain more insight into the predicted vibrational spectra and the unique electronic structure of FOOF, we determined the force constants for FOOF at the numerical LDF level in terms of internal coordinates. These force constants are given in Table IV. The diagonal force constants for the LDF and QCISD(T)" calculations are similar in form except that the LDF stretching force constants are larger. What is interesting to note are the magnitudes of the off-diagonal elements which lead to coupling of the internal modes. As noted previously, there are significant couplings of the 0-F stretches to the 0-0 stretch. The off-diagonal couplings of the 0-F stretches to the 0-0 stretch are about 50% of the 0-F diagonal force constant, showing the very strong coupling between the stretches. The coupling of the bends to the 0-0stretch are also not inconsequential and are about 22% of the diagonal bending force constant. The QCISD(T)" calculation predicts that this off-diagonal element is about 35% of the diagonal bending force constant. The numerical LDF predicts very little coupling of the 0-0stretch to the torsion in contrast to the QCISD(T) results which predict a large coupling. The remaining off-diagonal couplings are small with the largest being the coupling of the 0-F stretch with the bend; this value is about 50% of the off-diagonal element coupling the bend to the 0-0 stretch. Clearly, the force constants show how closely coupled the various internal coordinates are and how a small change in one geometric parameter will affect the other ones. The experimental assignments for FOOF were difficult to make as there were no gas-phase measurements until recently, and early measurements were made in a low-temperature matrix.Thus there was considerable ambiguity as to the value of the 0-0 stretch. This was resolved in the gas-phase experiment^.'^ This has also made a force constant analysis difficult. As shown in Table IV, our force constants are in reasonable agreement with the values given by experiment except for the two calculated stretching force constants which we know are based on the frequencies that are too large. It is worth noting that the a b initio force constants are a complete harmonic set whereas the experimental values are only a partial set. This accounts for the difference in signs for F24. Rotation Barrier. One of the early successes of Hartree-Fock theory was the prediction of the rotation barrier in HOOH.I4 The cis barrier is 7.6 kcal/mol from experiment, and the trans barrier

Thus it is unlikely that FOOF will rotate without dissociating. A simple orbital model shows what occurs. As the structure is rotated, the O2orbitals will go into one component of the IA state. Thus the rotation barrier should be at least the order of the 32 to * A transition in O2 which is 0.98 eV (or 23 kcal/mol).s' However, the possibility of the 02 orbitals giving a 1~ coupling may outweigh the 0-F bond interactions, and the moleeule can dissociate with loss of a fluorine atom. We calculated the rotation barrier to both the cis and trans structures. The structures and barriers are summarized in Table V. In order to obtain convergence at the two planar structures, it was necessary to allow the HOMO and LUMO to mix. This yielded a wave function in which about half an electron is transferred to the LUMO from the HOMO. For 7 = 0", the HOMO is at 8.7 (8.6) eV with 1.41 (1.38) e and the "LUMO" is at 8.1 (8.0) eV with 0.59 (0.62) e whereas a t T = 180°, the HOMO is a t 9.1 (9.0) eV with 1.48 (1.46) e and the "LUMO" is at 8.4 (8.3) eV with 0.52 (0.54) e with the DGauss values in parentheses. This result is not surprising as the HOMO of FOOF is a t 9.5 eV and the LUMO is at 6.0 eV. Thus destabilization

(46) Loos. K. R.; Goetschel, C. T.; Campanile, V. A. J . Chem. Phys. 1970. 52. 4418. (47) Burdett, J. K.; Gardiner, D. J.; Turner, J. J.; Spratley, R. D.; Tchir, P. J . Chem. SOC., Dalron Trans. 1973, 1928. (48) Jacox, M . E. J. Mol.Spectrosc. 1980, 84,74.

(49) Hunt, R. H.; Leacock, R. A.; Petewrs, C. W.; Hecht, K. T. J . Chem. Phys. 1965, 42, 1931. (50) Lyman, J. L. J . Phys. Chem. Re/. Data 1989, 18,799. (51) Huber, K. P.; Herzberg, G. Consranis ofDiaromic Molecules;Van Nostrand Reinhold: New York. 1979.

F

0

0 ;

trans-FOOF

9202 The Journal of Physical Chemistry, Vol. 95, No. 23, 1991 TABLE VI: Bond Dissociation Energies (kcal/mol) method FOOF ---L FOO + F FOO


L D F DG-5 LDF DG-5 expt

+ B-P

52 30 19.9

4

F + 0,

49 28 13.5

due to rotation would lead to the HOMO and LUMO approaching each other, and mixing of the orbitals can occur. This result is consistent with the fact that at the planar structures the molecule is probably trying to dissociate to lose an F. The calculated energetics are consistent with this model with a barrier of about 40 kcal/mol at the LDF leveL5* The Gaussian LDF results were corrected for nonlocal effects by using the Becke correction30with an additional term due to Perdew3' (B-P correction). These nonlocal corrections were applied after a self-consistent energy was obtained at the optimum LDF geomtries. The Becke correction is to the exchange energy, and the Perdew correction is to the correlation energy. The energy lowering due to the correction term is small, lowering the rotation barrier by 3-4 kcal/mol. The important result is that the calculated barriers are very high, consistent with the experimental estimates, and are greater than the experimental 0-F dissociation energy. The geometries for the planar structures are also consistent with the energetics and the model for dissociation. The 0 - F bond length at the numerical LDF level has lengthened by 0.17 A and the 0-0 bond distance has shortened to less than the value in O2(1.208 A).51 Thus the molecule is clearly destabilized and partially dissociated. The Gaussian LDF results differ somewhat as the 0-0distance is predicted to lengthen slightly and the 0-F bonds are not quite as long. Dissociation Energy. The dissociation energy of FOOF is an interesting theoretical target because of the weak 0 - F bond.50 The geometry of FOO was optimized in the ZA" state with DG-5 and gave the following parameters: r(O-0) = 1.212 A, r(O-F) = I .594A, and (0-0-F) = 11 1 . 7 O . These values can be compared to the experimental values of 1.200 A, 1.649A, and 1 1 I .2O? The agreement between theory and experiment for the 0-0bond and the 0-0-F angle is good, but the 0 - F bond is calculated to be 0.055 too long. The dissociation energy is calculated to be 52 kcal/mol, which is too high as expected from previous LDF Inclusion of the nonlocal Becke-Perdew correction gives 30 kcal/mol. This is qualitatively in agreement with the experimental value of 20 kcal/mol. Inclusion of a different nonlocal correction term based on a less rigorous theoretical gives 21 kcal/mol for the dissociation energy FOOF FOO' F'. This result suggests that better models need to be developed for these nonlocal terms. More importantly, the calculations predict a bond dissociation energy less than the torsion barrier when nonlocal effects are included. Thus the conclusions reached above that the molecule will probably dissociate on rotation are consistent with this dissociation energy calculation. One can also calculate the dissociation energy of FOO if the energy of O2is known. The optimized bond distance in O2with DG-5 is 1.228 A. The bond dissociation energies in Table VI give results that are similar to those found for dissociating FOOF with the dissociation of FOO requiring slightly less energy as is found

-

+

Conclusions

The present results demonstrate that local density functional theory provides a good description of the molecular geometry and ( 5 2 ) If one electron is placed in the HOMO and one in the LUMO, the barrier via the trans is 48 kcal/mol and via the cis is 49 kcal/mol. (53) Yamada, C.; Hirota, E. J . Chem. Phys. 1984, 80, 4694. (54) A nonspherical (broken symmetry) solution for the F atoms is employed.

Dixon et al. vibrations in the highly correlated FOOF molecule. Bond distances are predicted to within 0.015 A of experiment, bond angles within about l o , and vibrational frequencies within about 20%. Such deviations in predicted geometries are typical for LDF theory and are also found in many other molecular and solid-state systems. This consistent behavior of LDF theory is in contrast to Hartree-Fock (HF) theory which fails for FOOF whereas typically H F predicts geometries of noncorrelated molecules within a few hundredths of an angstrom and within 1-2' for angles. In the present investigation, the LDF equations have been solved by two different computational techniques, one using a numerical basis set with numerical integrations for the matrix elements and the other employing analytic Gaussian-type orbitals with analytic integration methods. Both approaches lead to results of similar agreement with experiment. The bond lengths calculated with the Gaussian basis sets are converged to within 0.01 A at the DZVP level, but there is still a need to do a systematic study on the role of more uncontracted basis sets with additional polarization functions. The most sensitive bond, as might be expected, is the 0-0 bond. An important question to at least address qualitatively is why do the LDF calculations yield such a good description of the potential energy surface near the minimum for FOOF. LDF theory can be expected to work well when (i) there is a fairly high electron density and (ii) this electron density does not show steep gradients. Furthermore, as a third criterion, LDF theory underestimates the ionic character in a bond and should be expected to work well when electronegativity differences between bonded atoms are small. If these criteria are met, LDF theory correctly accounts for all many-body effects due to the electron-electron interactions; Le., the theory correctly describes correlation effects. The FOOF molecule should meet these three criteria well as the following: (i) both 0 and F have a high number of valence electrons; (ii) the electron density should be fairly evenly distributed over the entire molecule; and (iii) the electronegativity difference between 0 and F is not large, are valid. Inclusion of nonlocal corrections to the exchange-correlation energies appears to be required in order to obtain a reasonable description of the dissociation energies. Those corrections are necessary, since the local density approximation overestimates binding energies, as is well-known from a number of cases. In summary, the present LDF results for the geometry and vibrational frequencies of FOOF are substantially better than those obtained from Hartree-Fock theory on the single-determinant level and compare in quality with those achieved with a traditional correlated ab initio method such as the QCISD(T) method. The LDF results on FOOF are consistent with those on other small molecules where a correlated treatment is required.20,24bs30-55 Because of their scaling in the computational efforts, those traditional correlated methods are limited to small molecules. In contrast, LDF calculations with the same quality in basis sets and numerical grids as used here scale by less than the third power in the number of basis functions. Thus, the LDF approach is applicable to molecular systems containing 100 atoms and more using present-day large-scale computing technology. This suggests that the LDF method will be very useful in predicting the electronic structure of large molecular systems which must be treated at a correlated level at the same accuracy as achieved for the four-atom system in this study.

Acknowledgment. We thank Nathalie Godbout for help in developing the basis sets. Registry No. FOOF,7783-44-0. ( 5 5 ) Salahub, D. R.; Fournier, R.; Mxynarski, P.; Papai, 1.; St-Amant, A.; Ushio, J. In Density Functiond Methods in Chemistry; Labanowski, J. K., Andzelm, J. W., Eds.; Springer-Verlag: New York, 1991; Chapter 6, p 77.

Density functional study of a highly correlated molecule, oxygen

Recommend Documents