Heats of Formation of Alkanes Calculated by Density Functional

Heats of formation of organic molecules calculated by density functional theory: II. Alkanes. Jan Labanowski , Lawrence Schmitz , Kuo-Hsiang Chen , No...
1 downloads 0 Views 927KB Size
Download PDF
J. Phys. Chem. 1995,99, 9603-9610

9603

Heats of Formation of Alkanes Calculated by Density Functional Theory Norman L. Allinger* and Kazuhisa Sakakibara Computational Center for Molecular Structure and Design, Department of Chemistry, University of Georgia, Athens, Georgia 30602-2556

Jan Labanowski Ohio Supercomputer Center, 1224 Kinnear Road, Columbus, Ohio 43212-1163 Received: February 27, 1995@

The heats of formation for alkanes can be calculated using density functional theory with the accuracy comparable to that obtained by ab inifio Hartree-Fock and second-order many-body perturbation theory (MBPT[2] vel MP2) calculations. The density functional approach not only accounts for the correlation energy (though in approximate manner) but is much faster than comparable MBPT[2] calculations. The set of equivalents for C-H and C-C bonds and methyl, iso, and neo arrangements was derived from the experimental heats of formation for 12 molecules. These equivalents combined with the calculated total energy and simple statistical mechanical corrections yield values of heats of formation which are in excellent agreement with experimental data.

Introduction The heat of formation of a compound (Le., its enthalpy of formation from the elements in the standard state at 298.16 K and 1 atm) is an important quantity, and it is notoriously difficult to measure with high accuracy. They can be calculated by molecular mechanics if sufficiently good experimental data exist for the class of compounds of interest.'-4 Molecular mechanics calculations are ordinarily very good for this purpose but occasionally may yield a spurious result due to a glitch in the force field, resulting from our lack of understanding of some important point. Calculating heats of formation de novo is at the present time possible only for the smallest molecules with the most advanced ab initio methods, and even then, the results are not always satisfactory. In this case, the calculated quantity is ordinarily the atomization energy, i.e., the energy change upon dissociating the molecule from its equilibrium geometry into constituent atoms at 0 K. They can be converted (with some loss of accuracy) to heats of formation using thermochemical data on elements in their standard states (e.g., JANAF and NBS thermochemical t a b l e P ) and accounting for the vibrational motions (e.g., by computing them theoretically). Comparison of theoretically calculated atomization energies with experimental data is often used to assess the accuracy of the theoretical methods. Pople and co-workers introduced a database of 55 small molecules for which experimental atomization energies are known with errors of less than 1 kcal/mol as a benchmark for their Gaussian- 1 (Gl)738and Gaussian-2 (G2)9 procedures for molecular energies. While these two procedures, when used for de novo calculation of heats of atomization, provide an impressive agreement with experiment (average error of 1.6 and 1.2 kcaumol, for G1 and G2 procedures, respectively), they are complex and computer intensive and, therefore, can only be applied to very small molecules (the largest molecules in the series contain only three non-hydrogen atoms). They involve calculations at the correlated level (fourth-order perturbational theory, MP4; and the quadratic configuration interaction, QCI) and require basis sets of valence triple-l; quality with diffuse functions and extensive polarization. The same 55 molecules were also used to evaluate the potential of density functional @

Abstract published in Advance ACS Absrracfs, May 1, 1995.

0022-365419512099-9603$09.0010

theory (DFT) for de novo calculations of atomization energies.I0 The results for the basis-free DFT approach implemented in Becke's NUMOL program" were quite poor for the local spin density approximation (mean absolute error of 36.2 kcal/mol), but the gradient corrections dramatically improved the agreement with experiment (mean absolute error of 3.7 kcdmol). Similar accuracy (mean absolute error of 4.4 kcdmol) was obtained by the GTO-based DFT method with B-VWN functional.'* The DFT results can be further improved, however, not without cost. The hybrid method using exact Hartree-Fock (HF) exchange in conjunction with DFT functionals brings the mean absolute error down to 2.4 kcal/mol for the same series of 55 molecule^.'^-'^ According to Becke,I5 further improvement is unlikely unless exact-exchange information is considered. This, however, requires that the traditional, computer-intensive, HF calculation is performed. It is unlikely that the accurate de novo calculations of atomization energies will become feasible for much larger molecules anytime soon. The approximations used in routine calculations by traditional ab initio methods (basis set truncation and limited correlation treatment) and in the DFT approaches (limited grid density andor basis sets and approximate nature of energy functionals) lead to chemically unacceptable errors. However, these errors tend to be systematic for given chemical groups and bonds; Le., these errors can be removed by calibration over a small number of reliable experimental data. Extensive studies along these lines were reported previously'6-21 for small- and medium-size molecules (up to eight non-hydrogen atoms) using the ab initio Hartree-Fock method. The overall accuracy was better than 1 kcdmol. The Hartree-Fock calculations and molecular mechanics may use the same experimental database, but the Hartree-Fock calculations may suffer from different kinds of defects; for example, for a particular class of compounds, the importance of correlation energy may become apparent. There may also be special requirements for the basis sets. We found, not surprisingly, that when molecular mechanics and ab initio calculations give the same result, we can be rather confident that that result is in fact correct and an experimental value which differs significantly from that result is probably in error. 1995 American Chemical Society

9604 J. Phys. Chem., Vol. 99, No. 23, 1995

Allinger et al.

TABLE 1: Data for Alkanes energy, hartree compd

POP

TORS

T/R

methane ethane propane n-butane isobutane n-pentane isopentane neopentane 2,%-dimethylbutane 2,3-dimethylbutane cyclopentane cyclohexane

0.000 00 0.000 00

0.000 00 -0.000 67 0.OOO 00 0.OOO 67 0.000 00 0.001 34 0.000 67 0.000 00 0.OOO 67 0.000 67 ”0 00 0.000 00

0.003 82 0.003 82 0.003 82 0.003 82 0.003 82 0.003 82 0.003 82 0.003 82 0.003 82 0.003 82 0.003 82 0.003 82

0.000 00 0.000 48 0.000 00 0.001 04 0.000 10 0.000 00 0.000 00 0.000 11 0.000 00 0.000 00

AT,“ kcal mol-’ -17.89 -20.04 -25.02 -30.03 -32.07 -35.08 -36.85 -40.14 -44.48 -42.61 -18.44 -29.50

bondgroup count Me

is0

neo

c-c

C-H

0 2

0 0 0 0 1 0 1 0 0 2 0 0

0 0 0 0

0 1 2 3 3 4 4 4 5 5 5 6

4 6 8 10 10 12 12 12 14 14 10 12

2 2 3 2 3 4 4 4 0 0

0 0 0 1 1 0 0 0

Experimental heats of formation were taken from ref 48.

The primary limitation of the ab initio method is the long running times of the computations for larger molecules. We have previously used a 6-31G* basis set and optimized the geometry at that level. Molecules containing about eight firstrow atoms have running times on the order of many hours of Cray-YMP time and are the practical limit to molecular size at present. The hope is that DFT (which generally requires less computer time for small molecules and for which the formal increase in computation time with molecular size is less steep than for ab initio methods) will permit extension of the method to much larger molecules. Moreover, a major limitation of the HF method of unknown importance at present is that the results do not account for any correlation energy. The correlation corrections needed for the Hartree-Fock energies can in part be absorbed in the parametrization. Obviously, this is not a serious problem with alkanes, since we are able to fit experimental data quite well. And while it does not yet appear to be a serious problem with the various functional groups we have studied, it may become more serious in other kinds of molecules. DFT22-24accounts for the correlation energy explicitly (though in an approximate manner), so any correlation problems which arise in the Hartree-Fock calculations will hopefully be less serious in the density functional theory treatment. There are numerous reports (see, e.g., refs 11-15, 24-33, and 37 and references given there) comparing the accuracy of traditional ab initio methods and various DFT approaches for calculating ground-state molecular properties and reaction energies. Since there are still many areas unexplored by the DFT methodology and since new DFT approaches and functionals are actively being developed and refined, it is premature to make generalizations. However, it is widely believed that the state-of-the-art DFT approaches outperform in most cases the traditional Hartree-Fock calculations, as far as agreement with experiment is concerned. The DFT calculations are also less computer intensive. Depending on the problem at hand, the accuracy of the DFT calculations is comparable to MP2 or higher order correlated ab initio methods. The DFT methodology is still in the active developmental stage, and new ideas and approaches are proposed frequently. Moreover, the available codes are usually less well optimized and “tuned” compared to the more mature ab initio programs. However, as the current project shows, DFT has a potential to compete favorably with traditional ab initio methods in terms of computational speed, accuracy, and reliability of the results.

Methods We have examined a representative set of relatively small saturated hydrocarbon molecules, which are listed in Table 1. In the previous Hartree-Fock calculations, as in molecular

mechanics, we have found that six adjustable parameters are needed to fit the available experimental data on alkanes. These are the C-C and C-H bond energies and the structural features which take into account the fact that a carbon atom may be substituted by any number (from zero to four) of other carbon atoms. Thus, there are five different substitution pattems to consider. These require only three independent variables for their description, however. Thus, two of these pattems may be assigned relative energies of zero. We have chosen methane and the CH2 group as our standards, so parameters are assigned to methyl, iso, and neo carbon types; however, there are other equivalent ways in which this assignment can be made. The equation relating heat of formation, AHf, calculated energy, Etot, and bondlgroup parameters, a,, for the ith molecule is AH$] = Et,#]

+ POP[i] + TORS[i] + T/R[i] + 5

j= 1

POP[i] is the Boltzmann correction needed to account for the case when more than one low-energy conformation is significantly populated at room temperature. The TORS[i] term represents the component added to the enthalpy coming from low-lying (Le., torsional) vibrational levels beyond what is calculating using the harmonic approximation. The T/R[i] term represents the contribution from the three rotations and three translations of the (nonlinear) molecule as a whole, and an additional RT is needed to convert energy to enthalpy. The parameters aj are multiplied by the number of bonddgroups in the given molecule-n[Q. For further details on these statistical mechanical correlations, see refs 1 and 2. The bondlgroup parameters can be roughly estimated by trivial calculations. For example, the methane data can provide an estimate of the UC-H parameter as

as in this case, the terms TORS and POP are zero. Using UC-H so derived, the ac-c parameter may be estimated from the data for cyclohexane. It would yield 1

UC-C = ; ( w P [ C 6 H , 2 ]

- Etot[C6H12]- T/R - 12aC-H)

(3) since here TORS and POP are also negligible. The ac-c parameter from cyclohexane may be used to estimate the Me group parameter based on data for n-alkanes, etc. This simple

Heats of Formation of Alkanes

J. Phys. Chem., Vol. 99, No. 23, 1995 9605

TABLE 2: Atomic Orbital Basis Sets for the deMon Calculations hydrogen (41/1*) DZVP

carbon (621/41/1*) DZVP

shell

exponent

coeff

shell

exponent

s

50.999 180 0000 7.483 218 0000 1.777 468 0000 0.519 329 5000

0.009 660 4760 0.073 728 8600 0.295 858 1OOO 0.715 905 3000

s

s

0.154 1100000

1.0000000000

2808.064 OOO oo00 421.138 300 OOOO 95.586 620 OOOO 26.739 000 0000 8.432 827 0000 2.760 582 0000

0.002 017 8300 0.015 433 2000 0.075 581 5500 0.247 828 2000 0.479 372 5000 0.333 834 4000

p

0.7500000000

1.0000000000

s

P

5.447 004 0000 0.479 242 2000 0.146 156 5000 18.130 850 0000 4.099 883 0000 1.185 837 0000 0.368 597 4000

-0.077 840 7700 0.568 956 0000 1.000 000 0000 0.015 854 7300 0.095 682 7700 0.304 911 9000 0.493 501 7000

P d

0.109 720 0000 0.600 OOO 0000

1.000 OOO 0000 1.000 000 0000

S

procedure, however, tends to propagate errors in experimental heats of formation of individual molecules. By using multiple regression, one avoids the accumulation of errors. It is convenient to rearrange eq 1 for the regression analysis as

where n i l = AH&] - E,,,[i] - POP[i] - TORS[i] - T/R.The problem is suitable for the no-intercept multiple regression analysis, frequently called the regression through origin (see, e.g., the monograph of Smillie34). Parameters uj are calculated from the formula

where a is the column vector of coefficients, n represents the rectangular matrix whose ijth element represents the number of occurrences of the jth bondgroup in the ith molecule, Y is a vector of values Y, for each molecule as defined by the previous equation, and symbol T denotes matrix transposition. Table 1 lists experimental values of the heats of formation and other parameters used in these calculations. Calculation of the Total Energy by Quantum Methods. The total energies needed for the derivation of the heats of formation were calculated by several different quantum methods. The ab initio calculations were performed with the ACES I1 program.35 At the Hartree-Fock level, two sets of results were obtained: for a 6-31G* basis set with no polarization functions on hydrogen atoms, and for a 6-31G** basis set with polarization functions on all atoms. Full geometry optimization with analytical derivatives was performed in all cases. Full core MBPT[2] (frequently referred to as MP2) calculations were also performed with full geometry optimization for these molecules with the 6-31G** basis set. Geometry optimization was terminated when the largest component of the gradient was smaller than 0.0001 hartree/bob. The starting geometries for the IW6-3 lG* calculations were taken from molecular mechanics calculations with the MM3 program36 and, for some molecules, from our previous calculations with the same method.16 For the HF/6-31G** calculations, we used the final geometries from the 6-31G* calculations as input. The MBPT[2]/6-31G** input was prepared with the final coordinates from the HF/6-31G** results.

coeff

carbon (7111/411/1*) TZVP shell

exponent

coeff

s

5784.157 000 0000 869.303 500 0000 198.511 6000000 56.429 900 0000 18.285 460 OOOO 6.448 715 0000 2.341 860 0000

0.OOO 818 9950 0.006 293 4710 0.031 781 1800 0.1 17 273 4000 0.303 476 3000 0.453 521 4000 0.243 059 1000

s

5.459 533 0000 0.478 196 8000 0.145 730 1000 34.258 560 OOOO 7.863 895 OOOO 2.344 519 0000 0.796 171 5000 0.272 680 4000 0.089 260 4700 0.600 000 0000

1.000 000 0000 1.OOO 000 0000 1.000 000 0000 0.005 804 2900 0.040 640 3400 0.155 021 9000 0.353 144 4000 1.000 000 0000 1.OOO 000 000 1.OOO 000 0000

s

s p

p p d

For DFT calculations, d e M ~ nwas ~ ~used. The program deMon uses contracted Gaussians as atomic basis functions. However, the coefficients and exponents of contracted Gaussian functions optimized for ab initio calculations are not necessarily optimal for the DFT calculations. Better results are obtained if the basis sets are optimized by atomic DFT calculations; sets used in deMon were derived in this way.38 Two different atomic basis sets, of double-5 (DZVP) and triple-5 (TZVP) quality, were used in current calculations for carbon atoms. For hydrogen, only the DZVP basis set was used. These sets are listed in Table 2. For the DZVP carbon basis set, we also tried two selections controlling the grid density for integration/fitting of the exchange correlation contribution to energy and forces, namely, FINE (Fg) and EXTRA FINE (Xg). The deMon uses the Becke molecular grid^^^,^ approach. The FINE grid uses 832 grid points (32 radial shells of 26 angular points each) per atom while, the EXTRA FINE grid uses 1600 points (32 radial shells of 50 angular points) at the SCF level. The density of points is increased for the gradient calculations. The number of points per atom is 2968 for the FINE grid and 6214 for the EXTRA-FINE grid. However, in this case, the 32 radial shells differ in angular point density for the FINE grid with more grid points at the median shells. Additionally, deMon requires auxiliary basis functions for fitting the charge density, @(r), and exchange-correlation potential, v?(r). The expansion functions for these two quantities are structurally similar and consist of uncontracted Gaussian primitives. The first group is a set of s-type Gaussians (Le., Gaussians with angular portion corresponding to s-type orbitals-lone s-type Gaussians). The second group consists of one s-type Gaussian, three p-type Gaussians, and six d-type Cartesian Gaussians sharing the same value of the exponent (called the spd constrained set). The exponents of these primitive Gaussians are collected in Table 3. The calculations were performed for the VWM4' local spin density functional augmented with nonlocal gradient corrections, namely, Perdew correlation corrections4*and Becke exchange c0rrections.4~ The corrections were also incorporated during the calculations of the gradients required for geometry optimization. Starting coordinates were the same as for the HF calculations. The geometry optimization was carried out by the conjugate gradient method initially and then switched to the BFGS method close to the minimum energy geometry. The geometry optimization was terminated when the largest component of the gradient was smaller than 0.001 hartreehohr. It

Allinger et al.

9606 J. Phys. Chem., Vol. 99, No. 23, 1995 TABLE 3: Auxiliary Sets of Fitting Functions for the deMon Calculations charge density set shell

exuonent

exchange correlation set shell

exuonent ~

S S S

SPd

Hydrogen (3,1;3,1) A2 45.00 S 7.50 S 0.30 S 1.50 SPd Carbon (4,4;4,4) A2 1500.00 s 330.00 S 94.30 S 27.00 s 9.92 SPd 2.20 SPd 0.63 SPd 0.18 SPd

15.0 2.5 0.1 0.5 500.00 110.00 31.40 9.00 3.30 0.73 0.21 0.06

was, in our case, the limit of practical accuracy of the numerically calculated gradients. With larger molecules, we had problems with convergence close to the geometry of energy minimum. The expression for the gradient involves the expansion of exchange correlation potential in auxiliary functions; Le., it is no better than the quality of the auxiliary sets. Errors in gradients are especially awkward in the vicinity of minimum where gradients are small. These small errors manifest themselves as a lack of translational and rotational invariance. In such cases, we resorted to the VERSLUIS gradient optionu which performs more accurate, though computationally more expensive, gradient evaluations. The VERSLUIS approach tends to annihilate spurious one-center contributions by calculating gradients numerically. In cases where we still observed some oscillations close to the energy minimum, we simply accepted the lowest total energy result. For all those cases, the gradients were well below the 0.001 hartreehohr limit. The deMon program reports two different values of total energy. One, rather approximate, called “analytical” is calculated during cycles of the self-consistent field. For reasons of computational expediency, the corresponding exchange correlation potential is represented analytically by fitting on a relatively coarse grid. The other, which is calculated during geometry optimization, is called “numerical”. It is calculated on a denser grid (the extended grid is needed for the energy gradients), and here, the actual numerical values of exchange correlation densities are used, and the exchange correlation integrals are calculated numerically. While the analytical energy seems to be quite satisfactory in converging the density/spin orbitals during SCF iterations, its value was too noisy for the reliable fit of heats of formations. For example, for the DZVP basis set and the FINE grid option, the standard deviation between calculated and experimental heats of formation was 0.45 and 8.03 kcal/mol for numerical and analytical energies, respectively. For this reason, we report here only the more accurate numerical energies. Results and Discussion Values of bond and group equivalents derived from regression are collected in Table 4. Looking at this table, we see the same general patterns for the ab initio (HF and MBFT[2]) and DFT calculations. There are two large terms for the C-C and C-H bond energies. These contain not only the bonding energy but also the energy of putting together the atoms themselves from nuclei and electrons and are consequently very large numbers. The chain-branching parameters, a M e , aiso,and aneo, are much smaller and appear in the expected order. That is, branching

tends to stabilize the molecule, so aneois the most negative parameter, aiso is also negative but smaller, C I C H ~is zero by definition, and aMe is positive. Wiberg and Martin have shown that these chain-branching stabilizations are partly or completely due to electron correlati0n.4~ Therefore, we expect that the chain-branching parameters here would be much smaller in absolute value from DFT calculations than they were from HF calculations. The DFT branching parameters are indeed smaller than those from HF/ 6-31G* and HF/6-31G** calculations. The MBPT[2]/6-31G** calculations are still better in this respect, as they yield still smaller branching parameters. It is also desirable that these branching parameters be as small as possible for practical reasons. These terms have different values depending on the nature of the functions at the branch in functionalizedmolecules. The fewer such terms, or the smaller their values, the smaller the errors introduced by estimating them when experimental values are unavailable. The MBPT[2] results also yield the best root-mean-square (RMS) value when compared with experiment. Apart from the fact that experimental errors in heats of formation can approach 1 kcdmol, the excellent agreement for the MPBT[2] method may be partly fortuitous, since it is ~ e l l - k n o w n ~ ~ that this method (due to a fortunate cancellation of errors) sometimes gives better results with smaller basis sets (like 6-31G**) than with the more extensive and correlationconsistent basis sets. Heats of formation calculated by different methods together with the root-mean-square deviations between calculated and experimental values are listed in Table 5. We performed three sets of DFT calculations, different in the size of the basis set or the fineness of the fittingjintegration grid. Results from all of them are better than those from HF calculations, with the EXTRA FINE grid results being the best. The fact that in this case the quality of the grid weights more in the DFT calculations than the quality of the basis set seems to fall in place with our observation that we did not have problems with convergence during geometry optimization for the EXTRA FINE grid. However, in our case, for practical purposes, these three sets of calculations yield results of similar quality, and it seems that for production runs the slight improvement in accuracy cannot justify the additional expense of using the TZVP basis set or the EXTRA FINE grid. To compare the computational expense associated with each method, we recorded the average wall clock time for each geometry optimization iteration on the CRAY-YMP (all jobs were run in a single-processor mode). The time needed for the complete geometry optimization includes not only the solution of HF/MPBT[2] equations but also calculations of integral derivatives and construction of the Hessian matrix. On the other hand, the Hessian in the deMon program is not calculated analytically. These data are provided here for illustration only; the wall clock time was very dependent on the current load of the machine. Moreover, we did not attempt to vectorize the source code of deMon to run efficiently on the CRAY. We compiled the program with the standard cft77 compiler options and did not try to find the best compilation regime. Also, since we did not use a DIRECT mode in any of the methods, the wall clock times are biased with the disk access latencies which, again, are very load susceptible. Formally, the computation time in DFT, as implemented in deMon, increases approximately with N3 with the number of basis functions (assuming that the size of auxiliary functions follows the basis set size) as opposed to with N4 for the HF calculations and roughly fl for the MBPT[2] calculations. Of course, the exponents are sometimes made smaller by prescreening the density matrix elements and

Heats of Formation of Alkanes

J. Phys. Chem., Vol. 99, No. 23, 1995 9607

TABLE 4: Bond and Group Equivalents (in hartrees) Derived b-v Different Methods HF

DFT

group or bond

6-31G*

6-31G**

MBPT[2] 6-31G**

“Me “is0 aneo “C-c “C-H

0.002 3531 -0.005 6600 -0.01 1 9708 18.944 6560 10.040 71 10

0.002 291 1 -0.005 5697 -0.01 1 1877 18.944 4360 10.042 3440

0.001 5241 -0.001 0587 -0.001 6379 19.009 4670 10.084 3820

DZVP

+ Fg

0.001 7480 -0.004 0607 -0.009 8434 19.065 6830 10.123 3740

DZVP

+ Xg

0.001 7103 -0.004 0893 -0.009 6180 19.065 7010 10.123 3810

TZVP

+ Fg

0.001 7612 -0.004 0468 -0.009 6017 19.068 0920 10.124 8030

TABLE 5: Values of Heats of Formation [kcaYmol] Calculated from Equivalents Derived by Different Methods HF

DFT

molecule

6-31G*

6-31G**

MBPT[2] 6-31G**

methane ethane propane n-butane isobutane n-pentane isopentane neopentane 2,2-dimethylbutane 2,3-dimethylbutane cyclopentane cyclohexane RNS

-17.889 -20.066 -25.178 -29.904 -33.057 -34.535 -36.814 -41.039 -43.574 -42.128 - 18.430 -29.944 0.53035

-17.890 -20.082 -25.184 -29.900 -33.060 -34.519 -36.814 -41.036 -43.579 -42.129 -18.464 -29.925 0.52995

-17.889 -19.927 -25.115 -30.078 -32.204 -34.993 -36.819 -40.293 -44.3 24 -42.556 -17.562 -30.271 0.34962

TABLE 6: Approximate Regression Equations Relating the Wall Clock Time (t, s) Required per One Iteration of Geometry Optimization with the Number of Electrons (N) in the Molecule method

regression ea

HF/6-3 1G* HF/6-3 1G** MBPT[2]/6-31G** DFTDZVP Fg DFTDZVP Xg DFTRZVP Fg

t = 0.0004328N4 1 t = 0.0015748N4 88 t = 0.0000647fl+ 186 t = 0.007212fl+ 27 t = 0.02474fl 78 t = 0.0089377fl 48

+ +

~

+ + +

+

+

corr coeff 0.9997 0.9977 0.9985 0.9903 0.99 16 0.9930

integrals to be ~alculated.4~In fact, the cost of integral calculations in the HF method scales as N4 for small molecules only, while with integral prescreening, its asymptotic behavior for very large molecules is N‘,as stated in ref 12. Similarly, for very large molecules, the numerical integration in DFT would scale linearly with N for efficient implementations.I2 Moreover, different portions of the program may have different size dependences. However, for the sake of illustration, we fitted the iteration times to the simple equation t = anm b, where a and b are variable parameters, n is the number of electrons in the molecule (we assumed that n is roughly proportional to the number of basis functions for the given basis set), and m is the formal exponent of time/size dependence as described above. The equations arrived at for the different methods are summarized in Table 6. The results provide a ball-park estimate of the CRAY-YMP time required for one iteration of geometry optimization. The number of geometry iterations varied from molecule to molecule, and our experience suggests that the number of iterations required for reaching energy minimum is roughly proportional to the number of atoms in the molecule (provided that the starting geometries are of similar quality). The averaged values of the C-C and C-H bond lengths, the C-C-C, C-C-H, and H-C-H bond angles, and the standard deviation of dihedral angles relative to those calculated by MBPT[2]/6-31G** were picked for comparison of geometries and are listed in Table 7. As the standard values of the C-C-C, C-C-H, and H-C-H angles change due to the substitution pattems at the central carbon atom, these bonds angles were classified as one of three types: -CR2-, -CHR-,

+

DZVP

+ Fg

-17.891 -20.253 -25.381 -29.797 -32.637 -34.638 -36.509 -40.743 -43.883 -42.503 - 19.299 -29.187 0.44992

DZVP

+ Xg

-17.890 -20.266 -25.332 -29.867 -32.630 -34.545 -36.602 -40.793 -43.825 -42.452 - 19.135 -29.338 0.43082

TZVP

+ Fg

-17.890 -20.232 -25.383 -29.793 -32.652 -34.643 -36.512 -40.750 -43.872 -42.490 -19.199 -29.265 0.4335 1

and -CH2-. The geometries resulting from different methods are very similar for practical purposes. The C-C and C-H bonds lengths calculated by the DFT methods are about 0.01 8, longer than then those from HF or MBPT[2] calculations. It is thought that in general, MBPT[2] substantially improves geometries compared with the HF results. The improvement here is negligible, however, because the molecules we studied were alkanes. It should be noted that the DFT method in our case can give accurate geometries comparable to 6-31G** ab initio calculations, and with far smaller computational times. The geometries calculated herein are always some approximation to the equilibrium geometries. These differ from experimental (vibrationally averaged) geometries, but normally, one set may be approximately corrected to the other set by adding simple bond length corrections. The corrections to be added are about 0.01 8, smaller for the DFT (re)geometries than for the ab initio geometries. Thus, predicted experimental geometries from DFT and ab initio become essentially indistinguishable. The DFT calculations are faster by about a factor of 10 for alkane molecules having six carbon atoms (2,2-dimethylbutane, 2,3-dimethylbutane, and cyclohexane) over the corresponding HF ab initio calculations. This advantage in shorter computing times will surely be more important when we carry out calculations on larger molecules. For example, for naphthalene, a 10-carbon molecule with 68 electrons, we estimated (by using the equations in Table 6) that the wall clock time for a singleprocessor CRAY-YMP, would be 9.4 and 26.1 h for HF/631G** and MBPT[2]/6-31G** but only about 40 min per geometry optimization iteration for DFTDZPfFg. Of course, the HF and MBPT[2] calculations for larger molecules would be expected to take less time than projected by our simplistic equations, since with the larger molecular size, the values of many of the integrals involving basis functions located at distant atoms are negligible small. Modem programs are sophisticated enough to take advantage of this fact by prescreening integrals based on the overlap matrix. However, even with the 50% reduction in computational time due to integral prescreening, routine ab initio computations seem currently prohibitively costly for molecules of this size but feasible with the DFT method. Moreover, in popular implementations of ab initio

9608 J. Phys. Chem., Vol. 99, No. 23, 1995

Allinger et al.

TABLE 7: Comparison of Results Obtained with Different ab Initio and DFT Calculations" HF parameter

DFT

+ Fg

+ Xg

+ Fg

6-31G*

6-31G**

MBPT[2] 6-31G**

-40.195 172 5 1.0837 109.471

-40.201 705 27 1.0835 109.471

Methane -40.369 856 31 1.0843 109.471

-40.525 827 19 1.0892 109.471

-40.525 853 68 1.0982 109.471

-40.531 542 20 1.0964 109.471

-79.228 755 46 1.5273 1.0856 111.209 107.679 0.106

-79.238 235 217 1.527 1 1.0857 111.210 107.678 0.106

Ethane -79.553 713 249 1.5219 1.0878 111.191 107.698 0.000

-79.844 848 66 1.5378 1.1029 111.441 107.432 0.127

-79.844 854 238 1.5356 1.1030 111.491 107.378 0.127

-79.855 824 76 1.5339 1.1015 111.517 107.350 0.127

-118.263 651 195 1.5281 1.0865 112.799 109.387 111.179 106.285 107.711 0.186

-118.276 159 846 1.5281 1.0866 112.866 109.367 111.170 106.295 107.721 0.178

Propane -118.740 881 1033 1.5226 1.0893 112.421 109.478 111.052 106.322 107.844 0.000

-119.166 121 169 1.5353 1.1034 112.923 109.430 111.273 105.960 107.610 0.260

- 119.166 060 570 1.5353 1.1038 112.908 109.426 111.288 105.996 107.594 0.201

-194.182 400 194 1.5337 1.1026 112.988 109.435 111.313 105.865 107.567 0.275

-157.298 409 576 1.5286 1.0870 113.088 109.323 111.193 106.232 107.696 0.224

-157.313 949 2255 1.5283 1.0873 113.133 109.311 111.193 106.229 107.696 0.203

n-Butane -157.928 171 3046 1.5229 1.0902 112.925 109.352 111.051 106.289 107.845 0.000

-158.486 739 338 1.5370 1.1042 113.140 109.416 111.297 105.769 107.586 0.28 1

-158.486 899 1144 1.5361 1.1046 113.196 109.406 111.321 105.750 107.559 0.28 1

- 158.508 277

-157.298 978 580 1.5309 1.0867 111.024 107.869 111.153 107.738 0.303

-157.314 556 2384 1.5308 1.0869 111.066 107.825 111.146 107.746 0.339

Isobutane -157.930 875 3186 1.5241 1.0900 110.869 108.035 110.947 107.955 0.000

-158.487 802 327 1.5387 1.1040 110.897 108.005 111.158 107.732 0.141

-158.487 774 1073 1.5386 1.1040 110.918 107.982 111.158 107.732 0.144

-158.509 397 378 1.5376 1.1034 110.915 107.986 111.214 107.672 0.141

-196.333 097 1307 1.5290 1.0873 113.164 109.305 111.196 106.222 107.693 0.138

-196.351 664 4769 1.5286 1.0876 113.211 109.297 111.195 106.203 107.693 0.141

-197.808 115 569 1.5362 1.1044 113.296 109.318 111.320 106.011 107.560 0.160

-197.808 048 1926 1.5360 1.1049 113.326 109.319 111.328 105.975 107.553 0.188

-197.834 933 689 1.5344 1.1028 113.422 109.312 111.385 105.896 107.491 0.246

-196.331 812 1341 1.5325 1.0867 111.lo7 115.002 107.779 108.842 111.243 106.092 107.641 0.271

-196.350 432 5177 1.5324 1.0869 111.138 115.047 107.746 108.824 111.235 106.123 107.650 0.269

-197.807 175 575 1.5408 1.1045 111.075 114.796 107.815 108.951 111.374 105.856 107.503 0.397

-197.807 337 1933 1.5408 1.1044 111.117 114.834 107.770 108.941 111.338 105.856 107.541 0.412

-197.834 017 666 1.5372 1.1030 111.163 114.989 107.720 108.903 111.326 105.843 107.553 0.459

n-Pentane -197.115 465 1.5231 1.0908 113.093 109.312 111.050 106.272 107.847 0.000 Isopentane -197.117 230 1.5254 1.0904 111.060 114.505 107.830 108.931 111.032 106.288 107.865 0.000

DZVP

DZVP

TZVP

379 1.5340 1.1026 113.333 109.382 111.356 105.693 107.522 0.339

Heats of Formation of Alkanes

J. Phys. Chem., Vol. 99, No. 23, 1995 9609

TABLE 7 (Continued) HF parameter Got

t

dc-c

RMS,

Eta t_ dc-c

RMS,

DFT

6-3 lG*

6-31G**

-196.333 818 1358 1.5351 1.0865 109.471 111.146 107.746 0.000

-196.352 464 5236 1.5351 1.0865 109.471 111.146 107.746 0.000

-235.364 605 2698 1.5372 1.0862 109.462 117.111 108.301 111.261 105.986 107.620 0.208

-235.386 310 10203 1.5372 1.0862 109.462 117.111 108.300 111.261 105.986 107.620 0.207

-235.363 062 2677 1.5373 1.0862 111.257 107.617 111.315 107.563 0.283

-235.384 757 10040 1.5371 1.0864 1 11.284 107.588 111.304 107.574 0.278

-195.163 580 1140 1.5401 1.0857 104.620 11 1.299 107.032 1.676

-195.178 865 3867 1.5401 1.0857 104.620 1 11.299 107.032 1.676

-234.208 007 2350 1.5320 1.0880 111.460 109.670 106.569 1.410

-234.226 253 793 1 1.5320 1.0880 111.460 109.670 106.569 1.410

'

MBFT[2] 6-31G** Neopentane -197.122 942 1S273 1.0900 109.471 110.845 108.064 0.000 2,2-Dimethylbutane -236.308 266 1.5290 1.0900 109.465 116.352 108.452 111.005 106.219 107.892 O.OO0 2,3-Dimethylbutane -236.305 078 1.5292 1.0901 111.179 107.705 111.067 107.826 0.000 Cyclopentane -195.922 962 1.5358 1.0899 104.298 111.311 107.274 0.000 Cyclohexane -235.121 446 1.5265 1.0925 111.141 109.703 106.771 0.000

DZVP

+ Fg

DZVP

+ Xg

TZVP

+ Fg

-197.809 117 583 1.5446 1.1033 109.471 111.070 107.827 0.053

-197.809 430 1932 1.5445 1.1033 109.471 111.071 107.825 0.034

-197.836 206 614 1.5423 1.1017 109.471 111.120 107.773 0.042

-237.127 221 985 1.5463 1.1035 109.466 116.592 108.438 111.191 105.996 107.696 0.129

-237.127 391 3278 1.5460 1.1035 109.466 116.615 108.436 111.214 105.973 107.672 0.139

-237.159 550 1192 1.5450 1.1022 109.466 116.645 108.451 111.264 105.870 107.619 0.164

-237.126 854 957 1S442 1.1040 111.318 107.555 111.324 107.554 0.290

-237.126 753 3278 1.5443 1.1040 111.322 107.552 111.318 107.560 0.281

-237.158 965 1193 1.5431 1.1029 111.350 107.521 111.371 107.504 0.313

-196.596 729 430 1.5502 1.1034 104.936 111.240 106.965 3.243

- 196.596 629 1476 1S503 1.1031 104.938 111.239 106.964 3.248

- 196.622 905 570 1.5482 1.1021 104.943 111.272 106.833 3.248

-235.924 919 718 1.5406 1.1059 111.578 109.638 106.567 1.346

-235.925 351 248 1 1.5404 1.1060 111.573 109.650 106.521 1.353

-235.956 644 900 1.5402 1.1057 111.582 109.645 106.534 1.352

The following entries are provided for each molecule: total energy (Etot)in hartrees; averageowall clock time of single_processor of CRAYYMP per iteration of geometry optimization with (t) in seconds; average bond lengths (dc-x) in A; average bond angles (OA-C-B) with a central carbon atom belonging to group [-CXY-1, and root-mean-square (RMS,) between torsion angles calculated by the given method and MBFT[2]/ 6-3 1G**. For larger molecules, timditeration ( t ) for MBR[2]/6-3 1G** was not recorded.

correlated methods, the benefit of integral prescreening is essentially lost. For simple correlated methods, integral transformation from an atomic to a molecular orbital basis is a computational bottleneck, while for higher order correlated methods, computation and transformation of integrals is no longer a rate-limiting step. The cost advantage of the DFT method should continue to increase with further increasing the molecular size. Conclusions The DFT approach using moderate basis sets and grid density has been shown to give results of comparable quality to those from the MBPT[2]/6-31G** method in the case of estimating heats of formation but at the fraction of the computational

expense. At this point, we consider the DFT approach with double-5 basis sets augmented with polarization functions on all atoms and the FINE integratiodfitting grid as the best compromise for obtaining accurate heats of formation with the shortest computation times. Acknowledgment is made to the donors of the Petroleum Research Fund, administered by the American Chemical Society, for support of this research. This project was also supported by a grant of Cray-YMP computer time from Ohio Supercomputer Center. Supplementary Material Available: Atom numbering and tables containing final geometries of the 12 alkane molecules

9610 J. Phys. Chem., Vol. 99, No. 23, 1995 determined by DFT and ab initio methods used in this work (11 pages). Ordering information is given on any current masthead page. References and Notes (1) See, for example: Burkert, U.; Allinger, N. L. Molecular Mechanics; American Chemical Society: Washington, DC, 1982. (2) Allinger, N. L.;Schmitz, L. R.; Motoc, I.; Bender, C.; Labanowski, J. K. J . Am. Chem. SOC. 1992, 114, 2880. (3) Wiberg, K. B. J . Comput. Chem. 1984, 5, 197. (4) Ibrahim, M. R.; Schleyer, P. v. R. J . Comput. Chem. 1985, 6 , 157. (5) Chase, M. W., Jr.; Davies, C. A.; Downey, J. R., Jr.; Frurip, D. J.; McDonald, D. A,; Syverud, A. N. J. Phys. Chem. Re& Data 1985,14 (Suppl. 1).

(6) Wagman, D. D.; Evans, W. H.; Parker, V. B.; Schumm, R. H.; Halow. I.: Bailev. S. M.: Churnev. ' K. L.: Nuttall. R. L. J . Phvs. Chem. Ref. Dura 1982,jl (Suppl. 2). (7) Pede, J. A.: Head-Gordon, M.; Fox, D. J.; Raghavachari, K.; Curtiss, L.'A. J. Chem. Phys. 1989, 90, 5622. (8) Curtiss, L. A.; Jones, C.; Trucks, G. W.; Raghavachari, K.; Pople, J. A. J . Chem. Phys. 1990, 93, 2537. (9) Curtiss, L. A,; Raghavachari, K.; Trucks, G. W.; Pople, J. A. J . Chem. Phys. 1991, 94, 7221. (IO) Becke, A. D. J . Chem. Phys. 1992, 96, 2155. (1 1) Dickson, R. M.; Becke, A. D. J . Chem. Phys. 1993, 99, 3898. (12) Johnson, B. G.; Gill, P. M. W.; Pople, J. A. J . Chem. Phys. 1993, 98, 5612. (13) Becke, A. D. J . Chem. Phys. 1993, 98, 1372. (14) Becke, A. D. Int. J . Quantum Chem., Quantum Chem. Symp. 1994, 28, 625. (15) Becke, A. D. J . Chem. Phys. 1993, 98, 5648. (16) Allinger, N. L.; Schmitz, L. R.; Motoc, I.; Bender, C.; Labanowski, J. K. J . Phys. Org. Chem. 1990, 3, 732. (17) Schmitz, L. R.; Motoc, I.; Bender, C.; Labanowski, J. K.; Allinger, N. L. J . Phys. Org. Chem. 1992, 5 , 225. (18) Allinger, N. L.; Yang, L.; Motoc, I.; Bender, C.; Labanowski, J. K. Heteroatom Chem. 1992, 3, 69. (19) Schmitz, L. R.; Chen,'Y. R. J . Comput. Chem. 1994, 15, 1437. (20) Allinger, N. L.; Schmitz, L. R.; Motoc, I.; Bender, C.; Labanowski, J. K. J . Comput. Chem. 1992, 13, 838. (21) Liu, R.; Allinger, N. L. J . Phys. Org. Chem. 1993, 6 , 551. (22) Parr,R. G.; Yang, W. Density Functional Theory of Atoms and Molecules; Oxford University: New York, 1989. (23) Kryachko, E. S.; Ludena, E. V. Density Functional Theory of Many Electron Systems; Kluver: Dordrecht, 1990. (24) Ziegler, T. Chem. Rev. 1991, 91, 651. (25) Andzelm, J.; Wimmer, E. J . Chem. Phys. 1992, 96, 1280. (26) Dixon, D. A,; Andzelm, J.; Fitzgerald, G.; Wimmer, E. J . Phys. Chem. 1991, 95, 9197.

Allinger et al. (27) Fitzgerald, G.; Andzelm, J. J . Phys. Chem. 1991, 95, 10531. (28) Laming, G . J.; Handy, N. C.; Amos, R. D. Mol. Phys. 1993, 80, 1121. (29) Laming, G. J.; Termath, V.; Handy, N. C. J . Chem. Phys. 1993, 99, 8765. (30) Murray, C. W.; Laming, G. J.; Handy, N. C.; Amos, R. D. Chem. Phys. Lett. 1992, 199, 551. (31) Salahub, D. R.; Castro, M.; Proynov, E. I.; NATO ASI Ser., Ser. B 1994, 318, 411. (32) Sim, F.; St-Amant, A,; Papai, I.; Salahuz, D. R. J . Am. Chem. Soc. 1992, 114, 4391. (33) Ziegler, T. NATO AS1 Ser., Ser. C 1992, 367, 357. (34) Smillie, K. W. Introduction to Regression and Correlation; Academic: London, 1966. (35) Stanton, J. F.; Gauss, J.; Watts, J. D.; Lauderdale, W. J.; Bartlett, R. J. ACES2 quantum package; Quantum Theory Project of the University of Elorida at Gainesville. (36) The MM3 program is described by: Allinger, N. L.; Lii, J.-H.; Yuh, Y. H. J . Am. Chem. SOC.1989,111,8551,8566,8576. The program is available to all users from the Technical Utilization Corp., Inc., 235 Glen Village Ct, Powell, OH 43065; to commercial users only from Tripos Assoc., 1699 South Hanley Rd, St. Louis, MO 63144; and to academic users only from the Quantum Chemistry Program Exchange, University of Indiana, Bloomington, IN 47405. The current version is available to run on most types of computers, and interested parties should contact one of the distributors directly. (37) Salahub, D. R.; Fournier, R.; Mlynarski, P.; Papai, I.; St-Amant, A.; Ushio, J. In Density Functional Methods in Chemistry; Labanowski, J., Andzelm, J., Eds.; Springer-Verlag: New York, 1991; p 77. (38) Godbout, N.; Salahub, D. R.; Andzelm, J.; Wimmer, E. Can. J. Chem. 1992, 70, 560. (39) Becke, A. D. J. Chem. Phys. 1988, 88, 2547. (40) St-Amant, A. Ph.D. Thesis, DBpartementde chimie, Universitt de Montrtal, 1992. (41) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J . Phys. 1980, 58, 1200. (42) Perdew, J. P. Phys. Rev. B 1986, 33, 8822; (erratum) Ibid. 1986, 38, 7406. (43) Becke, A. D. Phys. Rev. A 1988, 38, 3098. (44) Versluis, L.; Ziegler, T. J . Chem. Phys. 1988, 88, 322. Versluis, L. Ph.D. Thesis, University of Calgary, Alberta, Canada, 1989. (45) Wiberg, K.; Martin, E. J. Am. Chem. SOC.1985, 107, 5035. (46) Del Bene, J. E.; Aue, D. H.; Shavitt, I. J . Am. Chem. SOC.1992, 114, 1631. (47) Gill, P. M. W.; Johnson, B. G.; Pople, J. A. Chem. Phys. Lett. 1994, 217, 65. (48) (a) Pedley, J. B.; Naylor, R. D.; Kirby, S. P. Thermochemical Data of Organic Compounds, 2nd ed.; Chapman and Hall: London, 1986. (b) Cox, J. D.; Pilcher, G. Thermochemistry of Organic and Organometallic Compounds; Academic: London, New York, 1970.

JF9506806