A Theoretical Approach for Accurate Predictions of the Enthalpies of

Nov 3, 2010 - Donald W. Rogers , Andreas A. Zavitsas , Nikita Matsunaga. Wiley Interdisciplinary Reviews: Computational Molecular Science 2013 3 (1), ...
0 downloads 0 Views 3MB Size
Download PDF
12334

J. Phys. Chem. A 2010, 114, 12334–12344

A Theoretical Approach for Accurate Predictions of the Enthalpies of Formation of Carotenes Jenn-Huei Lii, Fu-Xing Liao, Hong-Yi Hsieh, and Ching-Han Hu* Department of Chemistry, National Changhua UniVersity of Education, Changhua 50058, Taiwan ReceiVed: August 23, 2010; ReVised Manuscript ReceiVed: October 19, 2010

A computational approach has been designed for accurately determining enthalpies of formation (∆Hf) for the carotene species. This approach, named correlation corrected atomization (CCAZ), is based on the concept of bond and group additivity, and is applied along with density functional theory (DFT). Corrections to the deficiencies in DFT were divided into 1,2-, 1,3-, and 1,4- atomic interactions, which were determined by comparisons with the G3 data of the training set. When comparing predictions from CCAZ combined with two different DFT methods (B3LYP and MPWB1K), fairly accurate prediction is expected. In contrast, DFT using the atomization and isodesmic schemes resulted in poor predictions of ∆Hf. The equivalent methods, atomic equivalent (AEQ) and group equivalent (GEQ) provide improved predictions; however, the accuracies are lower than that of CCAZ. 1. Introduction Carotenes play important roles in biological systems. They act as antenna in the light-harvesting units in chlorophyll and bacteriochlorophyll, and they also protect biological tissues from the deleterious effects of singlet oxygen and triplet bacteriochlorophyll.1,2 Dietary carotenes have been shown to offer antioxidation effects.3,4 For example, Di Mascio et al. demonstrated lycopene to be most effective in quenching singlet oxygen.5 Some carotenes such as β-carotene (Scheme 1) are converted to retinol (vitamin A). Thus, carotenes are vital for vision, skin and membrane protection, and anticancer activity.6 Natural carotenes consist of the carbon skeletal framework and one or two of the seven end groups shown in Scheme 1. For example, β-carotene consists of two β end groups, while lycopene consists of two ψ end groups (Scheme 1). With few exceptions, most carotenes consist of a C40 skeleton.7 Thermodynamic data, particularly the enthalpy of formation (∆Hf) of molecules, is a fundamental property in chemistry. Using this, the enthalpies of a reaction can be obtained before performing an experiment. However, no such data exists for carotenes. In fact, the ∆Hf for hydrocarbons of more than three conjugated double-bonds has not been reported. The extensively delocalized polyene system of carotenes discourages the determination of their ∆Hf by experimental or computational means. In the past decade, accurate predictions of the enthalpies of formation for gas-phase molecules have been achieved via highlevel quantum chemistry methods. Several theoretical approaches have been developed for the purpose of accurate predictions for the ∆Hf of medium-sized molecules, including radicals and ionic species.8-16 Among them, the GX theories and their modified versions have received much attention among chemists.10-14 For neutral, closed-shell molecules consisting of first and second period atoms, chemical accuracy ((1 kcal/mol) can be expected via the G3 theory. G3 theory with atomization scheme was developed via best fit with experimental data and has been demonstrated as being superior to G2 theory.11 G3 is also more computationally efficient than G2, as it reduces the size of the basis set at the QCISD(T) level of theory.11 With * To whom correspondence should be addressed.

SCHEME 1: Structures of β-Carotene and Lycopene, and the Seven End Groups Found in Natural Carotenes

the G2/97 test set (148 molecules), the average absolute deviation (AAD) of G3 predicted ∆Hf at 298 K is only 0.94 kcal/mol, a significant improvement compared with the G2 result of 1.56 kcal/mol.11 For hydrocarbons, the AAD of G3 theory is only 0.68 kcal/mol, compared with the AAD of 1.29 kcal/mol for G2 theory.11 Several variants of the G3 model have been proposed to reduce the computational costs. The G3MP2 model significantly reduced computational cost and disk storage by avoiding MP4 single-point energies.14 The G3B3 model replaced the MP2 geometries and scaled Hartree-Fock frequencies of G3 with those obtained using density functional theory (DFT), namely B3LYP.15 Similarly, the G3MP2B3 model replaced the MP2 geometries and scaled Hartree-Fock frequencies of G3MP2 with those obtained using B3LYP.15 Another line of endeavor for accurate energy computations has been the complete basis set (CBS) approach of Peterson et

10.1021/jp107973s  2010 American Chemical Society Published on Web 11/03/2010

Accurate ∆Hf Predictions for Carotenes al.8 In this approach, extrapolation to the CBS limit is achieved using asymptotically converging pair natural orbital expansions. The CBS model is comparable to the G2 model.8 The accuracy of CBS was further improved using B3LYP optimized geometries and vibrational frequencies. The variant, CBS-QB3, is more accurate than CBS.17 The size of carotenes forbids applications of ab initio methods. DFT is the natural alternative; however, previous DFT studies pose warnings for their applications in the computations of the ∆Hf for carotenes. For n-alkanes, Curtiss et al. demonstrated that B3LYP involves an accumulation of errors as the molecular size increases.18 Systematic error was shown by Schleyer et al. to be the inadequacy in describing long-range nonbonded interactions, such as the protobranching stabilization effects of DFT.19 For n-alkanes, DFT underestimates the stabilization energy of extended alkane chains. This underestimation is even more so for branched alkanes. Grimme demonstrated that DFT presents serious difficulties in describing the stereoelectronic alkylation effects.20 Despite the problems with DFT, it remains a valuable alternative to ab initio methods for its computational effectiveness. Thus, in these studies, isodesmic or homodesmic reactions were applied, and improved results were obtained for extended alkanes.18 In a study for the predictions of the ∆Hf of C40H56 carotenes, we have demonstrated that predictions using the isodesmic and homodesmic reactions are more accurate than those obtained using atomization reactions.21 However, noticeable accumulations of error were observed. Bond and group additivity approaches based on an empirical experimental data set have been widely used in estimating enthalpies of formation.22-28 The concepts of atomic equivalents (AEQs) and bond and group equivalents (BGEs) are among the many attempts to predict the thermodynamics properties of molecules. Both originated from a similar idea, and both have been explored by several research groups.29-32 For example, Jorgensen et al. used a large set of training molecules to make accurate predictions with semiempirical methods that included AM1, PM3, and MNDO.30 Liu et al. used the AEQ scheme and converted DFT energies to enthalpies of formation to achieve an AAD for hydrocarbons of 0.81 kcal/mol.33 These approaches utilize a set of molecules with known enthalpies of formation. The AEQs of elements and BGEs of functional groups were obtained through least-squares fit between experimental enthalpies and computed energies. With a few fitted parameters, the enthalpies of formation could be obtained using the less accurately computed energies of molecules. The size of a carotene forbids high-level quantum chemistry computations such as G3. The application of AEQ and BGE, along with DFT or semiempirical methods may offer a remedy to this situation. In addition to the AEQ approach, a group equivalent (GEQ) method was proposed in the study of C40H56 carotene systems.21 The GEQ approach resembles AEQ, in which the atomization reaction is replaced by an isodesmic reaction, while the AEQs are replaced by GEQs (Vide infra).21 AEQ and GEQ predictions are more accurate than those obtained from isodesmic or homodesmic reactions; however, the accuracy is still less than satisfactory.21 The drawback of the AEQ and GEQ approaches is that these methods translate the computed enthalpies into predicted enthalpies of formation via the optimized equivalents. Differences in computed enthalpies of the isomers are retained. For this reason, the notable deficiencies of DFT in describing

J. Phys. Chem. A, Vol. 114, No. 46, 2010 12335 the protobranching effects and the ring to open chain differences are uncorrectable. Despite the documented artifacts of DFT, it is still a reasonable choice for the predictions of the ∆Hf for carotenes. One important functional feature of carotene is its long conjugated chain, which is not adequately described using the semiempirical methods. For this reason, the correlation corrected atomization (CCAZ) method was developed. The CCAZ approach is based on a concept similar to BGE, in which possible errors in a computational theory (DFT in this case) are corrected by adding 1,2-, 1,3-, and 1,4- functional interaction terms. 2. Computational Approach The ∆Hf of the seven end groups were calculated using atomization and isodesmic34 (bond separation) reactions. G2,10 G3,11 G3MP2B3,15 and CBS-QB38 theories were applied. G3MP2B3 is a variation of G3MP2,14 in which B3LYP/631G(d) geometries and zero-point energies (ZPEs) were used to replace MP2/6-31G(d) geometries and HF/6-31G(d) ZPE. The computations were performed using the Gaussian 03 series of programs.35 In addition to the ab initio approaches, several DFT functionals were applied and compared with those obtained using ab initio theories. The tested functionals include pure, generalized gradient approximation (GGA) methods, such as BLYP,36,37 BPW91,36,38 and PBEPBE.39 The accuracy of hybrid GGAs (B3LYP37,40 and MPW1PW91),38,41 and the kinetic energy density-dependent meta-GGA functionals (BB95,42 MPW1K,41,43 MPW1B95,41,42 and MPWB1K)44 were also ascertained. The DFT theories were applied with the 6-31G(d) basis set to obtain optimized geometries and harmonic vibration frequencies. Single-point calculations were made using the 6-311+G(3df,2p) basis set. The vibrational frequencies were scaled by 0.96, and were then used to obtain ZPE and thermal corrections to 298 K. A scale factor of 0.96 was chosen, as has been used in the G3MP2B3 approach.15 The DFT-predicted ∆Hf values were compared with those obtained by high-level theories. The molecular formulas of the end groups are C9H16 (ψ, β, ε, and γ), C9H18 (κ), and C9H12 (φ and χ). When using the isodesmic reactions, equations that retain the same number of single and double C-C bonds were used. For example, two equations were used in the computations for the C9H16 isomers consisting of two CdC bonds (ψ end group, eq 1) and one CdC bond (β, ε, and γ end groups, eq 2).

C9H16 + 7 CH4 f 2 C2H4 + 6 C2H6

(1)

C9H16 + 9 CH4 f 1 C2H4 + 8 C2H6

(2)

AEQ, GEQ, and CCAZ schemes were applied to a larger training set of molecules. The optimized parameters obtained from these schemes were later used to compute the enthalpies of formation for the carotene species. The enlarged set of molecules (Scheme 2) includes all-trans polyenes (Pn), substituted alkenes, and ring compounds. The AEQ approach bears a similarity to atomization, in which the predicted ∆Hf of CnHm is obtained from the computed enthalpies [E(CnHm), E(C), and E(H)] and the experimental enthalpies of formation for C and H atoms [∆Hf (C) and ∆Hf (H)] in eq 3.

12336

J. Phys. Chem. A, Vol. 114, No. 46, 2010

Lii et al.

∆Hf(CnHm) ) E(CnHm) - n[E(C) - ∆Hf(C)] - m[E(H) - ∆Hf(H)] ≈ E(CnHm) - nεC - mεH

∆Hf from the reference is corrected by structural terms, which account for the energy difference.

(3) The AEQ approach replaces the terms in the brackets with equivalents εC and εH, which are least-squares fitted to reproduce ∆Hf (CnHm) from the computed enthalpies E(CnHm) for a test set of molecules. Due to the lack of experimental data, the ∆Hf from the atomization scheme of G3 was used as a reference. Using the isodesmic reaction of the ψ end group as an example, the equation is ∆Hf(C9H16) ) E(C9H16) + 7[E(CH4) - ∆Hf(CH4)] 2[E(C2H4) - ∆Hf(C2H4)] - 6[E(C2H6) - ∆Hf(C2H6)] ≈ E(C9H16) + 7εCH4 - 2εC2H4 - 6εC2H6

(4) Similar to the AEQ approach, in the GEQ approach the equivalents (εCH4, εC2H4, and εC2H6) are least-squares fitted to reproduce the G3 enthalpies of formation from E(C9H16). As in the isodesmic scheme, the equations are expressed to retain the same number of C-C and CdC bonds. CCAZ originates from a different perspective from the AEQ and GEQ methods. For a DFT functional, the deviation of the SCHEME 2: Compounds Included in the CCAZ Parameterization

∆Hf(CnHm) ≈ E(CnHm) - n[E(C) - ∆Hf(C)] m[E(H) - ∆Hf(H)] +

∑ aiHi1,2 + ∑ bjHj1,3 + ∑ ckH1,4 k i

j

k

(5) In eq 5, ai, bj, and ck are the frequencies of 1,2, 1,3-, and 1,41,3 1,4 interactions H1,2 i , Hj , and Hk , respectively. The equation resembles the atomization reaction, while ∆Hf predicted by DFT is corrected by the inclusion of 1,2-, 1,3-, and 1,4-interactions. The CCAZ method resembles the BGE methods and is very similar to Allinger’s approach,45-47 which is commonly used in MM2/MM3/MM4 enthalpy of formation calculations. Both the CCAZ and Allinger’s methods adapt Benson’s “bond energy” scheme.48 The CCAZ method includes atomization energies in the calculations, while Allinger’s method implicitly spreads atomization energies over all parameters, especially bond enthalpy terms. Furthermore, MM4 uses conformational energies instead of total energies in the calculations. Therefore, MM4 parameters, especially bond enthalpy terms, contain energies to compensate for the difference between the total energy and the conformational energy of the system. For comparison, the MM4 molecular mechanics approach was applied.47 In this study, 14 CCAZ interactions (Scheme 3) were included. Four 1,2-interactions (C2dC2, C2-C3, C3-C3, and C-H) were introduced to describe chemical bonds in the molecules. The notions 2 and 3 represent sp2 and sp3 carbons, respectively. Nine 1,3-interactions are considered, six of which (C3A-C3F) consist of sp3 carbon centers and three consist of sp2 carbon centers (C2A-C2C). The 1,4-interaction (H14) term is introduced to account for the CdC-tert-butyl functional group. A computer program (HGC) was written to count the number of various CCAZ interactions. 3. Results and Discussion Relative Enthalpies of Formation of the Seven End Groups. Table 1 summarizes the ∆Hf of end groups predicted using G2, G3, G3MP2B3, and CBS-QB3 theories with the atomization and isodesmic schemes. Of the seven end group models, only two experimental values are available, and these values are rather small in magnitude. For the purpose of comparisons among theoretical methods, predicted ∆Hf of C2H4 (ethylene), C4H6 (trans-butadiene) and C6H8 (trans-hexatriene) are included in the table. For these five molecules with experimental ∆Hf, G3 predictions using atomization and isodesmic schemes are very close, and are in good agreements with the experimental values. G2 theory within the atomization scheme provides higher (more positive) values than G3. In contrast, G2 along with the isodesmic schemes tend to provide more negative values than G3. G3MP2B3 values are more negative with the atomization scheme, and are more positive with the isodesmic scheme compared with those of G3. CBS-QB3 predictions using atomization tend to overestimate. With the isodesmic scheme however, CBS-QB3 predictions are in good agreement with the experiment and with the G3 results. We chose ∆Hf obtained from the G3 atomization scheme as a reference in this study. Table 2 summarizes the deviations in the ∆Hf predictions of the end groups obtained with various DFT functionals. The table

Accurate ∆Hf Predictions for Carotenes

J. Phys. Chem. A, Vol. 114, No. 46, 2010 12337

SCHEME 3: Correlation Correction Interactions Included in This Study

presents data as deviations from G3 atomization results. With the atomization scheme, inferior results were obtained using the chosen functionals. The most serious overestimations were observed in BLYP, while the most serious underestimations were observed in PBEPBE. Among the four C9H16 end groups, the errors of the ring-compounds: β, ε, and γ are similar in magnitude; however, these errors are more positive than that of their open-chain isomer, the ψ end group. Nevertheless, none of the tested functionals encourages the application of this approach in the prediction of the ∆Hf of the end groups.

With the isodesmic scheme, the deviations were significantly reduced; however, the errors are still large. Among the predictions, those of the aryl (φ and χ) end groups are the most accurate. Comparatively, the predictions of the C9H16 end groups (ψ, β, ε, and γ) were noticeably overestimated. All tested functionals offer more positive ∆Hf values for the C9H16 ring species (β, ε, and γ) than for their open-chain isomer ψ. Thus, the disadvantage of applying DFT in the computations of the ∆Hf is partly attributed to its deficiency in describing the relative energies of the ring structures with respect to their open-chain isomers. Since the equivalent schemes are based on the concept of translating total computed enthalpies into the ∆Hf, the embedded error of DFT in describing the relative enthalpies of isomers is retained in the AEQ and GEQ schemes. Among the tested functionals, the errors of the meta-GGA functionals MPW1K, MPW1B95, and MPWB1K are the most consistent. For the AEQ and GEQ approaches, these functionals are more correctable than others are. (A systematic comparison of relative enthalpies of end groups with DFT is summarized in a table in the Supporting Information). For this reason, we used B3LYP and MPWB1K in the AEQ, GEQ, and CCAZ schemes. MM4, AEQ, GEQ, and CCAZ Approaches. In addition to the end groups, we added 30 molecules into our training set. These include (Scheme 2) all-trans polyenes up to P9, five- and six-membered rings, and relevant functional groups that are involved in the structure of carotenes. Within this training set, parameters involved in AEQ, GEQ, and CCAZ schemes are least-squares fitted to the ∆Hf from G3 atomization, for which B3LYP and MPWB1K functionals were applied. Table 3 summarizes the predicted ∆Hf. The table also includes the ∆Hf computed using the molecular mechanics approach MM4.47 Table 4 summarizes the values of the equivalents in the AEQ and GEQ schemes, while Table 5 summarizes the values of CCAZ interactions. As seen in the AADs in Table 3, AEQ of B3LYP exhibits the poorest results, with an AAD of 4.80 kcal/mol. The aforementioned overestimations for the ∆Hf of the ring species (β, ε, and γ) with respect to the open-chain C9H16 end group (ψ) were still observed. Two noteworthy overestimations are those for PI and PJ, which are 11.27 and 11.79 kcal/mol higher than the reference. The ∆Hf valies of polyenes are ∼4 kcal/ mol underestimated with AEQ-B3LYP. By contrast, with MPWB1K, the AAD of the AEQ scheme was reduced to 1.25 kcal/mol. The improvement is attributed to the fact that errors from MPWB1K predictions are more consistent (see Table 2) and, thus, are relatively more correctable via parametrization. The GEQ scheme significantly improves the accuracy of

TABLE 1: ∆Hf Values of End Groups (in kcal/mol) Predicted Using Atomization and Isodesmic Schemes atomization ψ β ε γ κ φ χ C 2 H4 C 4 H6 C6H8 a

isodesmic

exp.a

G3

G3MP2B3

G2

CBS-QB3

G3

G3MP2B3

G2

CBS-QB3

-3.33 ( 0.27 -2.29 ( 0.30 12.5 ( 0.1 26.00 ( 0.19 40.1 ( 0.6

-2.50 -24.52 -25.22 -23.11 -46.32 -2.92 -1.77 12.50 26.98 40.54

-3.34 -24.86 -25.56 -23.43 -45.93 -4.18 -3.18 11.96 25.63 38.46

-0.86 -22.82 -23.49 -21.63 -45.46 0.69 1.85 12.91 28.30 42.80

0.60 -21.54 -22.17 -19.90 -43.36 -1.00 -0.01 13.35 28.58 42.90

-2.83 -24.90 -25.60 -23.50 -46.56 -3.59 -2.43

-1.94 -24.02 -24.72 -22.59 -46.02 -1.47 -0.47

-3.77 -25.84 -26.50 -24.64 -47.72 -3.83 -2.66

-2.79 -24.69 -25.33 -23.06 -46.15 -4.89 -3.90

26.79 40.16

27.08 40.83

26.65 39.90

26.74 40.07

Experimental ∆Hf are obtained from the NIST (National Institute of Standards and Technology) database.52

12338

J. Phys. Chem. A, Vol. 114, No. 46, 2010

Lii et al.

TABLE 2: ∆Hf Values (in kcal/mol) of End Groups Predicted Using Atomization and Isodesmic Schemes at Various Levels of DFTsa B3LYP

BLYP

BPW91

MPWPW91

ψ β ε γ κ φ χ

15.35 23.54 23.99 24.93 30.09 13.98 14.78

31.23 42.86 43.67 44.48 54.30 24.66 26.55

3.58 9.77 10.31 11.51 21.02 -9.51 -8.82

-22.05 -16.53 -16.15 -15.00 -6.87 -33.24 -32.71

ψ β ε γ κ φ χ

7.45 12.98 13.42 14.37 19.08 4.31 5.12

8.42 13.90 14.71 15.52 20.61 5.16 6.05

6.53 11.30 11.86 13.05 18.05 1.06 1.74

4.91 10.34 10.72 11.87 16.41 0.78 1.31

a

PBEPBE

BB95

MPW1K

MPW1B95

MPWB1K

atomization -47.31 -43.31 -42.95 -41.94 -34.56 -58.82 -58.43

-14.16 -7.22 -7.11 -5.93 1.17 -25.96 -26.34

-7.45 -8.27 -8.15 -7.18 -5.89 -10.72 -10.30

-24.38 -22.55 -22.90 -21.80 -21.22 -28.89 -29.61

-22.33 -22.53 -22.90 -21.94 -22.83 -24.32 -24.29

isodesmic 3.84 9.27 9.63 10.65 14.85 0.11 0.50

5.39 8.91 9.03 10.21 12.62 -0.45 -0.83

5.13 9.45 9.57 10.53 14.77 1.12 1.54

3.66 7.16 6.81 7.91 9.36 -0.92 -1.64

3.75 6.91 6.55 7.51 9.49 -0.61 -χ0.59

Values are in deviations from G3 with the atomization scheme.

TABLE 3: Enthalpies of Formation (in kcal/mol) Predicted Using Various Methodsa B3LYP exp.b AAD signed average ψ β ε γ κ φ χ H P1 P2 P3 P4 P5 P6 P7 P8 P9 PA PB PC PD PE PG PH PI PJ PM PS PT PV PW PX PY PZ FA FB FG a

G3

MM4

∆

1.46 -1.37

-3.33 ( 0.27 -2.29 ( 0.30 12.5 ( 0.1 26.00 ( 0.19 40.1 ( 0.6

18.11 ( 0.16 -9.92 ( 0.21 -15.96 20.4 ( 0.4 18.09 ( 0.24 19.77 ( 0.22 -43.26 -40.14 ( 0.15 -32.07 ( 0.15 -36.99 ( 0.25 19.80 ( 0.20 36.00 ( 2.00 35.11 ( 0.24

-2.50 -24.52 -25.22 -23.11 -46.32 -2.92 -1.77 17.62 12.51 26.98 40.54 53.88 67.06 80.24 93.51 106.77 119.90 19.27 -10.58 11.83 -16.10 10.02 38.34 20.57 -5.26 -6.40 18.67 20.63 -43.72 -40.24 -31.97 3.32 19.43 -36.46 20.86 35.46 36.54

-4.36 -24.76 -25.67 -23.41 -45.46 -2.82 -2.13 16.28 12.52 26.20 38.60 50.95 63.27 75.59 87.92 100.24 112.56 18.57 -10.38 10.97 -16.02 9.21 36.42 19.63 -6.24 -8.09 18.40 19.09 -43.42 -40.67 -32.36 2.63 19.46 -37.09 19.88 35.40 35.23

AEQ

∆

GEQ

4.80 -0.64 -1.86 -0.24 -0.45 -0.30 0.86 0.10 -0.36 -1.34 0.01 -0.78 -1.94 -2.93 -3.79 -4.65 -5.59 -6.53 -7.34 -0.70 0.20 -0.86 0.08 -0.81 -1.92 -0.94 -1.12 -1.69 -0.27 -1.54 0.30 -0.43 -0.39 -0.69 0.03 -0.63 -0.98 -0.06 -1.31

-5.59 -19.45 -19.71 -16.66 -38.65 0.48 1.92 12.27 7.30 22.32 36.29 49.88 63.25 76.57 89.79 103.04 116.21 13.66 -15.94 6.83 -20.66 5.52 37.01 16.58 5.60 5.39 14.83 15.28 -40.60 -45.80 -39.59 5.39 17.12 -36.76 22.13 47.25 40.38

∆

1.88 0.01 -3.09 5.07 5.51 6.45 7.67 3.40 3.69 -5.35 -5.21 -4.66 -4.25 -4.00 -3.81 -3.67 -3.72 -3.73 -3.69 -5.61 -5.36 -5.00 -4.56 -4.50 -1.33 -3.99 11.27 11.79 -3.84 -5.35 3.12 -5.56 -7.62 2.07 -2.31 -0.30 1.27 11.79 3.84

-3.63 -25.62 -25.88 -14.69 -45.72 -3.89 -2.44 15.14 13.87 28.21 41.52 54.44 67.14 79.79 92.34 104.93 117.42 18.77 -11.74 11.16 -17.24 7.61 39.99 20.90 -1.24 -1.45 19.94 20.39 -46.89 -42.49 -35.50 8.14 21.45 -42.27 20.13 36.67 37.71

MPWB1K CCAZ

∆

0.36 0.01 -1.13 -1.10 -0.66 8.42 0.60 -0.97 -0.67 -2.48 1.36 1.23 0.98 0.56 0.08 -0.45 -1.17 -1.84 -2.48 -0.50 -1.16 -0.67 -1.14 -2.41 1.65 0.33 4.43 4.95 1.27 -0.24 -3.17 -2.25 -3.53 4.82 2.02 -5.81 -0.73 1.21 1.17

-3.01 -24.52 -25.79 -23.25 -46.15 -3.36 -1.92 17.36 12.35 27.12 40.88 54.24 67.36 80.48 93.50 106.51 119.44 19.17 -10.42 11.43 -14.47 9.65 39.26 20.31 -5.24 -5.59 18.59 20.79 -43.65 -40.77 -31.59 2.80 18.62 -36.84 21.45 35.36 36.74

AEQ

∆

1.25 -0.07 -0.51 0.00 -0.57 -0.14 0.17 -0.44 -0.15 -0.26 -0.16 0.14 0.34 0.36 0.30 0.24 -0.01 -0.26 -0.46 -0.10 0.16 -0.40 1.63 -0.37 0.92 -0.26 0.43 0.81 -0.08 0.16 0.07 -0.53 0.38 -0.52 -0.81 -0.38 0.59 -0.10 0.20

-2.79 -25.01 -26.07 -23.00 -46.82 -5.75 -4.57 16.45 13.64 28.15 41.76 55.06 68.23 81.36 94.45 107.53 120.59 19.31 -10.78 11.58 -15.69 8.73 39.35 21.39 -3.48 -3.90 19.32 20.70 -45.44 -40.36 -33.42 6.28 21.09 -39.71 18.05 33.03 34.60

GEQ

∆

1.05 -0.16 -0.29 -0.49 -0.85 0.11 -0.50 -2.83 -2.80 -1.17 1.13 1.17 1.22 1.18 1.17 1.12 0.94 0.76 0.69 0.04 -0.20 -0.25 0.41 -1.29 1.01 0.82 2.28 2.50 0.65 0.07 -1.72 -0.12 -1.45 2.96 1.66 -3.25 -2.81 -2.43 -1.94

-3.07 -24.13 -25.19 -23.28 -45.81 -5.12 -3.94 16.04 12.70 27.30 41.01 54.41 67.68 80.90 94.08 107.26 120.42 18.58 -11.38 10.96 -16.17 8.43 38.92 20.77 -2.50 -2.92 18.59 19.97 -44.54 -40.84 -34.01 5.88 20.47 -38.92 18.33 34.54 34.98

CCAZ

∆

0.41 -0.07 -0.57 0.39 0.03 -0.17 0.51 -2.20 -2.17 -1.58 0.19 0.32 0.47 0.53 0.62 0.66 0.57 0.49 0.52 -0.69 -0.80 -0.87 -0.07 -1.59 0.58 0.20 3.26 3.48 -0.08 -0.66 -0.82 -0.60 -2.04 2.56 1.04 -2.46 -2.53 -0.92 -1.56

-2.86 -24.52 -25.63 -23.31 -46.34 -3.33 -2.16 17.45 12.11 26.75 40.53 53.96 67.30 80.57 93.79 107.04 120.25 19.10 -10.07 11.25 -14.94 9.63 39.58 20.36 -5.34 -5.62 18.60 20.49 -43.78 -40.32 -31.88 2.91 18.65 -36.55 18.49 36.00 35.46

-0.36 0.00 -0.41 -0.20 -0.02 -0.41 -0.39 -0.17 -0.40 -0.23 -0.01 0.08 0.24 0.33 0.28 0.27 0.35 -0.17 0.51 -0.58 1.16 -0.39 1.24 -0.21 0.42 0.78 -0.07 -0.14 -0.06 -0.08 0.09 -0.41 -0.78 -0.09 -2.37 0.54 -1.08

∆ values are deviations from ∆Hf obtained from G3 atomization. b Experimental ∆Hf values are obtained from the NIST database.52

predictions: the AADs for GEQ with B3LYP and MPWB1K are 1.88 and 1.05 kcal/mol, respectively.

The MM4 predictions for the training set (AAD 1.46 kcal/ mol) are not as good as those predicted for conjugated alkenes

Accurate ∆Hf Predictions for Carotenes

J. Phys. Chem. A, Vol. 114, No. 46, 2010 12339

TABLE 4: The Equivalents (in hartrees) of the AEQ and GEQ Schemes AEQ εC εH GEQ εCH4 εC2H4 εC2H6

B3LYP

MPWB1K

-38.130404 -0.579678

-38.108412 -0.577484

-40.459394 -78.589984 -79.747902

-40.416577 -78.524962 -79.679985

TABLE 5: Optimized Parameters of the Correlation Correction Terms of the CCAZ Scheme B3LYP

C2dC2 C2-C3 C3-C3 C-H C3A C3B C3C C3D C3E C3F C2A C2B C2C H14

MPWB1K

in hartrees

in kcal/mol

in hartrees

in kcal/mol

-0.001317 -0.004013 -0.004359 0.000419 0.000470 -0.004759 -0.008145 0.001190 -0.009768 -0.003537 -0.002614 -0.002805 -0.004373 -0.006029

-0.826 -2.518 -2.735 0.263 0.295 -2.987 -5.111 0.747 -6.129 -2.220 -1.640 -1.760 -2.744 -3.783

0.001784 0.001913 0.002766 -0.001545 0.002070 -0.003773 -0.005497 0.001514 -0.005566 -0.001011 -0.002040 -0.000820 0.000427 -0.003123

1.119 1.200 1.735 -0.969 1.299 -2.368 -3.450 0.950 -3.493 -0.634 -1.280 -0.515 0.268 -1.960

(AAD 0.6 kcal/mol).49 The signed average of -1.37 kcal/mol reveals that MM4 has a tendency to underestimate the ∆Hf. The most serious deviations were seen in the polyene system, in which the magnitude of underestimation is proportional to the polyene size. This feature of MM4 would affect its applications for carotenes (Vide infra). However, the agreement among G3, MM4, and the 19 experimentally reported values within the training set is good (see Figure 1). The AADs for G3 and MM4 against experimental values are 0.50 and 0.39 kcal/mol, respectively. The CCAZ scheme provides the best estimations with either B3LYP or MPWB1K functionals. The AADs are 0.36 and 0.41 kcal/mol with B3LYP and MPWB1K, respectively. The good accuracy of CCAZ within the training set is attributed to its

Figure 1. ∆Hf predicted using G3 and MM4 compared with available experimental data (diagonal line) in the training set. The AADs for G3 and MM4 against experimental values are 0.50 and 0.39 kcal/mol, respectively.

much larger number of parameters (14), compared with those in AEQ (2 parameters) and GEQ (3 parameters). With CCAZB3LYP, only one data item (PD) exceeds the reference by more than 1 kcal/mol, while for CCAZ-MPWB1K there are four (PD, PG, FA, and FG). For the polyenes and end groups, the CCAZ predictions are in good agreement with the reference data and the experiments. On the other hand, the CCAZ scheme resolves the difficulty in describing n- and branched alkanes using DFT.20 With CCAZB3LYP, the ∆Hf of neopentane and n-pentane are -40.77 and -34.03 kcal/mol, respectively; with CCAZ-MPWB1K, the predictions for them are -40.32 and -35.75 kcal/mol, respectively. These predictions are in good agreement with the experimental results (-40.14 and -35.08 kcal/mol for neopentane and n-pentane), both in absolute and relative values. Another example demonstrates the much improved accuracy of CCAZ. The compounds PI and PJ are extended from β and ε end groups to possess one more CdC bond. The ∆Hf of PI and PJ predicted using AEQ and GEQ schemes, either with B3LYP or MPWB1K, were overestimated by several kilocalories per mole. MM4 predictions are more accurate in these two cases. However, CCAZ accurately predicts the ∆Hf of PI and PJ. The G3 prediction of a ∆Hf of PJ that is 1.14 kcal/mol lower than that of PI was surprising since PI consists of a delocalized π-system. Two features within these species might explain this. First, the ∆Hf of end group ε is 0.70 kcal/mol lower in energy than that of β. Second, the two CdC bonds of PI adopt a torsional angle (54.4° with B3LYP) caused by steric repulsion. None of the AEQ or GEQ schemes provides a satisfactory prediction for the absolute and relative enthalpies for these species. In contrast, the CCAZ scheme combined with either B3LYP or MPWB1K provides satisfactory predictions. Among the 14 correlation correction terms for B3LYP (Table 5), 11 of them are stabilizing (negative value). The particularly stabilizing interactions include the C-C-H interaction C3B, C-C-C interactions C3C and C3E, and the 1,4-interaction H14. In contrast, the correlation corrections for MPWB1K are relatively less stabilizing and are smaller in magnitude compared with those for B3LYP. Similar to B3LYP, in MPWB1K the C3C and C3E terms are the most stabilizing. These terms correspond to the notorious deficiency of DFT in describing the stereoelectronic or protobranching effects.19,20 Both the H-C-H interactions C3A and C3B are destabilizing for B3LYP and MPWB1K; however, the magnitudes are smaller than those of the C-C-C interactions. In general, the correlation correction is smaller with MPWB1K. ∆Hf of Carotenes. Table 6 summarizes the ∆Hf of 22 carotenes predicted using MM4, AEQ, GEQ, and CCAZ. The smallest investigated carotene is septapreno-β-carotene (C35H50), and the largest is dodecapreno-β-carotene (C60H80). Seven C40H56 isomers, including lycopene and β-carotene, were investigated. When performing geometry optimizations, MM4 was applied to generate the lowest energy isomers. DFT optimizations were performed on these candidate species, and only the lowestenergy isomer was considered in this study. For example, the lowest-energy isomers of β-carotene and lycopene adopt a ∼Ci symmetry, and are marginally lower in energy than their ∼C2 counterparts. For β-CAR, our predicted bond angles and bond distances agree well with those obtained from theory and experiment.50,51 The torsional angle between the conjugated chain and the β-ionone ring is in good agreement with the theoretical prediction of 47.2°50 and the experimental value of 42.5°.51 Figure 2 shows the structures of four of the carotenes.

12340

J. Phys. Chem. A, Vol. 114, No. 46, 2010

TABLE 6: Enthalpies of Formation (in kcal/mol) of Carotenes Predicted Using Various Methods

Lii et al.

Accurate ∆Hf Predictions for Carotenes

J. Phys. Chem. A, Vol. 114, No. 46, 2010 12341

TABLE 6: Continued

The other carotenes are illustrated in the Supporting Information. Overall, the bond distances and bond angles predicted using B3LYP and MPWB1K (in parentheses) are similar, although B3LYP tends to predict longer bond distances than MPWB1K. Table 6 reveals that the ∆Hf values predicted using AEQB3LYP are significantly larger than those obtained using AEQMPWB1K. In marked contrast, ∆Hf values predicted using GEQ and CCAZ schemes are more consistent within the two DFT functionals. The AAD between GEQ-B3LYP and GEQMPWB1K is 2.07 kcal/mol, and the AAD between CCAZB3LYP and CCAZ-MPWB1K is 1.92 kcal/mol. The comparison of CCAZ-B3LYP data against CCAZ-MPWB1K is shown in Figure 3. The predictions of AEQ and GEQ are plotted against CCAZ results in Figure 4. As demonstrated in Figure 4, GEQ results are much closer to the CCAZ results than AEQ. The results from the AEQ scheme for B3LYP are more scattered and overestimated than those obtained from the GEQ scheme. The AEQ results for MPWB1K are slightly underestimated, while GEQ results for MPWB1K are in close agreement with CCAZ. In a previous study, ∆Hf values of the four C40H56 isomers;lycopene, prolycopene, R-carotene, and β-carotene;were predicted using

B3LYP with the GEQ scheme in a previous study.21 However, the predictions for the ∆Hf in this previous study now seem significantly overestimated when compared with the CE and CCAZ values obtained in this study. As shown in the figure, MM4 predictions are noticeably smaller than those from CCAZB3LYP and CCAZ-MPWB1K. This trend has been observed in the training set, that MM4 tends to underestimate ∆Hf of the more extended polyene systems. Overall, the CCAZ predictions for the ∆Hf of 22 carotenes from B3LYP and MPWB1K functionals are close. The largest deviation occurrs in dodecapreno-β-carotene (C60H80), for which B3LYP (135.63 kcal/mol) is 4.36 kcal/mol smaller than that predicted using MPWB1K (139.99 kcal/mol). The rather small deviation in the ∆Hf of carotenes resulting from different exchange-correlation functionals speaks to the reliability of the CCAZ approach. Extension of the CCAZ method to predict ∆Hf of oxygen-containing carotenoids (xanthophylls) is in progress. 4. Conclusion For the predictions of the ∆Hf of carotenes, we proposed a method, CCAZ, which corrects for the various bonding interac-

12342

J. Phys. Chem. A, Vol. 114, No. 46, 2010

Lii et al.

Figure 2. Optimized structure and selected geometrical parameters (bond lengths in Å and bond angles in degrees) of (a) β-carotene (C40H56), (b) lycopene (C40H56), (c) renieratene (C40H48), and (d) chlorobactene (C40H52). The B3LYP data are shown, with MPWB1K data in parentheses. Structures of the other carotenes are illustrated in the Supporting Information.

tions of DFT computations. Within the training set of 37 molecules, CCAZ with B3LYP has an AAD of 0.36 kcal/mol from the G3 data. In contrast, CCAZ with MPWB1K has an AAD of 0.41 kcal/mol. The CCAZ approach is recommended

when computing the enthalpies of formation for carotene species. As expected, DFT, using the atomization and isodesmic schemes, achieves poor results when computing the ∆Hf. The equivalent methods, AEQ and GEQ, provide improved predic-

Accurate ∆Hf Predictions for Carotenes

J. Phys. Chem. A, Vol. 114, No. 46, 2010 12343 scheme is suitable not only for carotenes, but also for hydrocarbons including acyclic and cyclic alkanes, alkenes, polyenes, and species consisting of aryl groups. Acknowledgment. The authors acknowledge that the National Science Council of Taiwan, Republic of China, supported this work. We also thank the National Center for HighPerformance Computing for computer time and facilities. Supporting Information Available: Optimized structures of the end groups and carotene species. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes

Figure 3. ∆Hf of the carotenes predicted using CCAZ-B3LYP compared with those obtained CCAZ-MPWB1K (diagonal line). The AAD is 1.92 kcal/mol.

Figure 4. ∆Hf of the carotenes predicted using MM4, AEQ, and GEQ compared with those obtained using CCAZ (diagonal line) along with (a) B3LYP and (b) MPWB1K.

tions; however, the accuracies are lower than those found for the CCAZ approach. Our DFT-based bond and group correction

(1) Koyama, Y.; Fujii, R. The Photochemistry of Carotenoids; Kluwer Academic: Dordrecht, 1999. (2) Telfer, A. Philos. Trans. R. Soc. London, B 2002, 357, 1431–1440. (3) Edge, R.; McGarvey, D. J.; Truscott, T. G. J. Photochem. Photobiol. B: Biol. 1997, 41, 189–200. (4) Sies, H.; Stahl, W. Am. J. Clin. Nutr. 1995, 62, 1315S–1321S. (5) Di Mascio, P.; Murphy, M. E.; Sies, H. Arch. Biochem. Biophys. 1989, 274, 532–538. (6) Rock, C. L. Carotenoids in Health and Disease; Mercer Dekker: New York, 2004. (7) Britton, G. FASEB J. 1995, 9, 1551–1558. (8) Ochterski, J. W.; Petersson, G. A.; Montgomery, J. A., Jr. J. Chem. Phys. 1996, 104, 2598–2619. (9) Mayer, P. M.; Parkinson, C. J.; Smith, D. M.; Radom, L. J. Chem. Phys. 1998, 108, 604–615. (10) Curtiss, L. A.; Raghavachari, K.; Trucks, G. W.; Pople, J. A. J. Chem. Phys. 1991, 94, 7221–7230. (11) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Rassolov, V.; Pople, J. A. J. Chem. Phys. 1998, 109, 7764–7776. (12) Pople, J. A.; Head-Gordon, M.; Fox, D. J.; Curtiss, L. A. J. Chem. Phys. 1989, 90, 5622–5629. (13) Curtiss, L. A.; Raghavachari, K.; Pople, J. A. J. Chem. Phys. 1993, 98, 1293–1298. (14) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K.; Rassolov, V.; Pople, J. A. J. Chem. Phys. 1999, 110, 4703–4709. (15) Baboul, A. G.; Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. J. Chem. Phys. 1999, 110, 7650–7657. (16) Henry, D. J.; Sullivan, M. B.; Radom, L. J. Chem. Phys. 2003, 118, 4849–4860. (17) Montgomery, J. A., Jr.; Frisch, M. J.; Ochterski, J. W.; Petersson, G. A. J. Chem. Phys. 1999, 110, 2822–2827. (18) Redfern, P. C.; Zapol, P.; Curtiss, L. A.; Raghavachari, K. J. Phys. Chem. A 2000, 104, 5850–5854. (19) Wodrich, M. D.; Corminboeuf, C.; Schleyer, P. v. R. Org. Lett. 2006, 8, 3631–3634. (20) Grimme, S. Ang. Chem. Int. Ed. 2006, 45, 4460–4464. (21) Tu, C.-Y.; Guo, W.-H.; Hu, C.-H. J. Phys. Chem. A 2008, 112, 117–124. (22) Wiberg, K. B. J. Comput. Chem. 1984, 5, 197–199. (23) Cohen, N.; Benson, S. W. Chem. ReV. 1993, 93, 2419–2438. (24) Ibrahim, M. R.; Schleyer, P. v. R. J. Chem. Phys. 1985, 6, 157– 167. (25) Dewar, M. J. S.; Storch, D. M. J. Am. Chem. Soc. 1985, 107, 3898– 3902. (26) Leal, J. P. J. Phys. Chem. Ref. Data 2006, 35, 55–76. (27) Gronert, S. J. Org. Chem. 2006, 71, 1209–1219. (28) Wodrich, M. D.; Schleyer, P. v. R. Org. Lett. 2006, 8, 2135–2138. (29) Allinger, N. L.; Schmitz, L. R.; Motoc, I.; Bender, C. J. Am. Chem. Soc. 1992, 114, 2880–2883. (30) Repasky, M. P.; Chandrasekhar, J.; Jorgensen, W. L. J. Comput. Chem. 2002, 23, 498–510. (31) Winget, P.; Clark, T. J. Comput. Chem. 2004, 25, 725–733. (32) Cioslowski, J.; Liu, G.; Piskorz, P. J. Phys. Chem. A 1998, 102, 9890–9900. (33) Mole, S. J.; Zhou, X.; Liu, R. J. Phys. Chem. 1996, 100, 14665– 14671. (34) Hehre, W. J.; Ditchfield, R.; Radom, L.; Pople, J. A. J. Am. Chem. Soc. 1970, 92, 4796–4801. (35) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery Jr., J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Lyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Adamo, C.; Jaramillo, J.;

12344

J. Phys. Chem. A, Vol. 114, No. 46, 2010

Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A.; Gaussian 03, revision B.05; Gaussian, Inc.: Pittsburgh, PA, 2003. (36) Becke, A. D. Phys. ReV. A 1988, 38, 3098–3100. (37) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785–789. (38) Perdew, J. P.; Burke, K.; Wang, Y. Phys. ReV. B 1996, 54, 16533– 16539. (39) Perdew, J. P.; K. Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865–3868. (40) Becke, A. D. J. Chem. Phys. 1993, 98, 5648–5652. (41) Adamo, C.; Barone, V. J. Chem. Phys. 1998, 108, 664–675. (42) Becke, A. D. J. Chem. Phys. 1996, 104, 1040–1046. (43) Lynch, B. J.; Fast, P. L.; Harris, M.; Truhlar, D. G. J. Phys. Chem. A 2000, 104, 4811–4815.

Lii et al. (44) Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2004, 108, 6908– 6918. (45) Labanowski, J.; Schmitz, L.; Chen, K.-H.; Allinger, N. L. J. Comput. Chem. 1998, 19, 1421–1430. (46) Lii, J.-H.; Allinger, N. L. J. Mex. Chem. Soc. 2009, 53, 95–106. (47) Allinger, N. L.; Chen, K.; Lii, J.-H. J. Comput. Chem. 1996, 17, 642–668. (48) Bensen, S. W. Thermochemical Kinetics: Methods for the Estimation of Thermochemical Data and Rate Parameters; John Wiley & Sons: New York, 1976. (49) Nevins, N.; Lii, J.-H.; Allinger, N. L. J. Comput. Chem. 1996, 17, 695–729. (50) Requena, A.; Cero´n-Carrasco, J. P.; Bastida, A.; Zu´n˜iga, J.; Miguel, B. J. Phys. Chem. A 2008, 112, 4815–4825. (51) Senge, M. O.; Hope, H.; Smith, K. M. Z. Naturforsch., C: J. Biosci. 1992, 47, 474–476. (52) Afeefy, H. Y.; Liebman, J. F.; Stein, S. E. Neutral Thermochemical Data. In NIST Chemistry WebBook, NIST Standard Reference Database Number 69; Linstrom, P. J., Mallard, W. G., Eds.; National Institute of Standards and Technology: Gaithersburg, MD, 2005 (http://webbook.nist.gov).

JP107973S