Hydrogen Abstraction Reactions from Phenolic ... - ACS Publications

0 downloads 0 Views 770KB Size Report
Sep 17, 2015 - Peroxyl Radicals: Multireference Character and Density Functional. Theory Rate ... substance that when present at low concentrations compared to that of an ..... is only partly surprising because the unreliability of T1 diagnostics has .... predict rate constants of the reactions between phenolic compounds ...
Article pubs.acs.org/JPCA

Hydrogen Abstraction Reactions from Phenolic Compounds by Peroxyl Radicals: Multireference Character and Density Functional Theory Rate Constants Annia Galano,† Leonardo Muñoz-Rugeles,‡ Juan Raul Alvarez-Idaboy,*,‡ Junwei Lucas Bao,§ and Donald G. Truhlar*,§ †

Departamento de Química, Universidad Autónoma Metropolitana-Iztapalapa, San Rafael Atlixco 186, Col. Vicentina. Iztapalapa, 09340 México D.F., México ‡ Departamento de Física y Química Teórica, Facultad de Química, Universidad Nacional Autónoma de México, México D.F. 04510, México § Department of Chemistry, Chemical Theory Center, and Supercomputing Institute, University of Minnesota, Minneapolis, Minnesota 55455-0431, United States S Supporting Information *

ABSTRACT: An assessment of multireference character in transition states is considered to be an important component in establishing the expected reliability of various electronic structure methods. In the present work, the multireference characters of the transition states and the forming and breaking of bonds for a large set of hydrogen abstraction reactions from phenolic compounds by peroxyl radicals have been analyzed using the T1, M, B1, and GB1 diagnostics. The extent of multireference character depends on the system and on the conditions under which the reaction takes place, and some systematic trends are observed. In particular, the multireference character is found to be reduced by solvation, the size of the phenolic compound, and deprotonation in aqueous solution. However, the deviations of calculated rate constants from experimental ones are not correlated with the extent of multireference character. The performance of single-determinant density functional theory was investigated for the kinetics of these reactions by comparing calculated rate constants to experimental data; the results from these analyses showed that the M05 functional performs well for the task at hand.



INTRODUCTION Phenolic compounds have become a subject of great interest in recent decades. One of the reasons is their strong antioxidant activity, which has been mainly related to two different reaction mechanisms: electron transfer and hydrogen transfer (HT).1 While the first mechanism becomes particularly important when the phenolic OH group is deprotonated, HT can be important regardless of the protonation states of the reacting species. Numerous theoretical investigations have been devoted to such systems, which, due to the relatively large size of naturally occurring and synthetic phenolic compounds, are usually performed using Kohn−Sham density functional theory (KS-DFT) with approximate exchange-correlation functionals. The reliability of KS-DFT depends strongly on the exchangecorrelation functional that is used, but all KS-DFT calculations are based on representing the density with a single Slater determinant, and just as for wave function-based methods with single-determinantal reference functions, its reliability with currently available functionals for systems with inherently multiconfigurational wave functions is usually lower than that © 2015 American Chemical Society

for systems whose wave function is represented qualitatively correctly by a single Slater determinant. We will follow the usual convention of labeling inherently multiconfigurational systems as multireference systems (they are also sometimes called “strongly correlated”). Peroxyl radicals (ROO•, where R denotes an organic fragment) are among the biologically relevant radicals that can be effectively scavenged by phenolic compounds to retard oxidative stress.2 In fact, it has been proposed that trapping peroxyl radicals is the main antioxidant function of phenolic compounds. 3,4 This is because the half-lives of ROO • compounds are not too short for the efficient interception by phenolic compounds.5 Because an antioxidant is defined as “any Special Issue: Piergiorgio Casavecchia and Antonio Lagana Festschrift Received: August 8, 2015 Revised: September 16, 2015 Published: September 17, 2015 4634

DOI: 10.1021/acs.jpca.5b07662 J. Phys. Chem. A 2016, 120, 4634−4642

Article

The Journal of Physical Chemistry A

diagnostic is that A25% is based on the total atomization energy (TAE), and B1 is based on a bond energy (BE). For the purposes of considering HT reactions, the multireference character of the O−H bonds seems more relevant than a diagnostic that also considers the multireference character of the whole molecule (including that due to spectator bonds or portions of the molecule that are not strongly coupled to the reaction of interest); therefore, we will consider the B1 diagnostic. The B1 diagnostic is defined as follows:14

substance that when present at low concentrations compared to that of an oxidizable substrate would significantly delay or prevent oxidation of that substrate”,6,7 it should react faster than the species being protected. Therefore, rate constants (k) for reactions of phenolic compounds with peroxyl radicals are important in assessing antioxidant protection. It has recently been demonstrated that the transition state involved in the HT reaction from the OH group in phenol, when reacting with the CH3OO• radical, has significant multireference character.8 If this were also the case for other reactions between phenolic compounds and peroxyl radicals (ROO•), it would possibly cast doubt on the reliability of some methods for calculating reaction barriers and consequently on some estimations of the rate constants. Consequently, to assess the reliability of computed rate constants, it becomes relevant to further investigate the possible multireference character of transition states involved in HT reactions between phenolic compounds and peroxyl radicals. In particular, although some exchange-correlation functionals are known for their good general performance for calculating rate constants, it is crucial to establish their reliability when used for the specif ic case of phenolic reactions with peroxyl radicals. As discussed above, this reliability may be influenced by the multireference character of the transition states. Unfortunately, judging multireference character is not a settled issue. There are several diagnostics available for that purpose.9 Probably the most frequently used is the T1 diagnostic, which is based on analyzing the Euclidean norm of the t1 amplitudes optimized in coupled-cluster calculations with single and double excitations (CCSD). This diagnostic is expected to allow qualitative assessments of the relevance of static correlation. The larger the T1 value, the less reliable the results obtained from CC, or other single-reference, calculations. It has been proposed that T1 values ≥0.02 and 0.045 indicate a significant multireference character for closed-shell and open-shell systems, respectively.10,11 However, it has been reported that these threshold values are not rigorously defined and can be system- and method-dependent.9 An alternative diagnostic has been recently proposed, which consists of quantitatively evaluating the multireference character in terms of the occupation numbers of a completeactive-space self-consistent-field (CASSCF) wave function.12 It is known as the M diagnostic, and has been used for measuring the multireference character in open-shell reagents and transition states. The proposed threshold values for the M diagnostic are 0.05 and 0.10, defining three ranges (M < 0.05, 0.05 ≤ M ≤ 0.10, and M > 0.10), which correspond to small, moderate, and large multireference characters, respectively.8 In a recent study, Martin and co-workers13 employed some very expensive wave function diagnostics of multireference character and found that the best of these is the contribution, which is usually unaffordable to calculate, of connected quadruple and connected pentuple excitations to the atomization energy; they call this %TAE[T4 + T5 ]. Then they considered affordable diagnostics and reported that “of the wave function-based diagnostics, the only one that has a reasonable correlation with %TAE[T4 + T5 ]” is the M diagnostic. They also recommended a very simple diagnostic that they call A25%, which is very similar to the earlier B114 diagnostic. Aside from normalization and the choice of generalized gradient approximation for the exchange-correlation density functional (to which the A25% diagnostic is not sensitive), the difference of the A25% diagnostic from the B1

B1 = EBE(BLYP) − EBE(B1LYP//BLYP)

where EBE denotes equilibrium bond energy (i.e., De, which is a difference of potential energies, rather than D0), and where the first bond energy is calculated with the BLYP15 exchangecorrelation functional and the second bond energy is calculated with the B1LYP16 functional at the geometry optimized by BLYP. The B1LYP functional differs from the BLYP functional simply in that 25% of the density-based exchange functional is replaced by Hartree−Fock exchange. A bond with a B1 diagnostic greater in magnitude than 10 kcal/mol is considered to be a multireference bond. One must, however, issue a cautionary note about the B1 diagnostic (and this would apply to the A25% diagnostic as well). In particular, the B1 diagnostic is more indirect than the T1 and M diagnostics, which are based on coefficients of configuration state functions and hence directly address the question of multireference character. The B1 diagnostic is based on an observed correlation; thus, it may sometimes be large for reasons other than multireference character. However, the B1 diagnostic has the very significant advantage of being much less expensive than other diagnostics; therefore, it is of great interest to study this diagnostic. It is also interesting to test a generalized B1 diagnostic. In this diagnostic, we calculate the classical barrier height (potential energy at the lowest-energy saddle point minus the potential energy at the equilibrium structure of reactants) by both BLYP and B1LYP//BLYP. The difference will be called the GB1 diagnostic. If the GB1 diagnostic exceeds 10 kcal/mol in magnitude, the reaction will be considered to be a multireference one. In the present work we have used the T1, M, B1, and GB1 diagnostics to investigate the amount of multireference character in a series of reactions between phenolic compounds and peroxyl radicals. In addition, because the ultimate goal is to obtain reliable rate constants, we also present direct comparisons of the calculated rate constants with the available experimental data. Multireference diagnostics were calculated both in the gas phase and in aqueous solution, but the rate constant calculations were carried out only in aqueous and organic solvents because no experimental rate constants are available in the gas phase.



COMPUTATIONAL DETAILS In applying the T1 and M diagnostics we calculate the diagnostic for the transition state of the reaction. In calculating the B1 diagnostic, we calculate the average of the B1 diagnostics for the O−H bond that is broken and for the O−H bond that is formed (that is, the B1 diagnostic values reported in this work are averages of the B1 of the phenolic O−H bond and the ROO−H bond). The GB1 diagnostic is explained above. T1 and M are unitless; B1 and GB1 are given in units of kilocalories per mole. 4635

DOI: 10.1021/acs.jpca.5b07662 J. Phys. Chem. A 2016, 120, 4634−4642

Article

The Journal of Physical Chemistry A Geometries. Full geometry optimizations, frequency calculations, and the B1 and T1 diagnostic calculations were carried out with Gaussian 09 package of programs,17 while the CASSCF calculations needed for the M diagnostics were performed with the ORCA program.18 The M05-2X functional and the 6-311+G(d,p) basis set were used for geometry optimizations and frequency calculations for the T1 and M diagnostics; the BLYP functional and 6-311+G(d,p) basis set were used for geometries for the B1 and GB1 diagnostics. Local minima and transition states were identified by the number of imaginary vibrational frequencies (0 or 1, respectively). Diagnostics. The T1 diagnostic was calculated at the CCSD/6-311+G(d,p)//M05-2X/6-311+G(d,p) level of theory. For the M diagnostic, the CASSCF wave functions were constructed using the correlated participating orbitals scheme, in particular the active space denoted as nom-CPO+π.12,19 For the reactions studied here, the nom-CPO+π space includes the σXH, σ*XH, and SOMOY orbitals and the π orbitals of the aromatic system. The level of theory for these calculations was CASSCF(9,9)/6-311+G(d,p), except for 4-hydroxyl-benzaldehyde, in which case a larger active space, namely (11,11), was used because of its extended conjugated system. B1 and GB1 diagnostics were performed with the 6-311+G(d,p) basis with ultrafine integration grids and a tight optimization criterion. Kinetics. Rate constants (k) in liquid solution were calculated by two methods. The first method is conventional transition-state theory (TST)20−22 with harmonic vibrational frequencies and Eckart tunneling.23,24 The Eckart barrier was fit to reproduce the gas-phase zero-point-inclusive energies of reactants, saddle point, and product. The computational details for the conventional TST calculations with Eckart tunneling are in line with the protocol of the quantum-mechanics-based test for overall free radical scavenging activity (QM-ORSA);25 this has been validated by comparison with experimental results, and its uncertainties were found to be no larger than those arising from experiments.25 These rate constant calculations were carried out with both the M05 and the M05-2X exchangecorrelation functionals because the latter is overall more accurate for main-group chemistry, but it was previously reported that its accuracy is lower than that of the former for reactions involving transition states with significant multireference character.8 Nevertheless, the TST+Eckart method is not the most quantitative method available for calculating rate constants because it neglects variational effects (i.e., the effect of placing the transition-state dividing surface at the variational transition state rather than at the saddle point),26 anharmonicity,27 the reaction-specific shape of the effective potential and reduced mass for tunneling, and corner-cutting tunneling.28,29 The second method was chosen to improve on these aspects. In particular, we employed multistructural variational transition-state theory30 (MS-VTST) calculations for the reaction of 1,4-benzenediol with HOO• in aqueous solution. The smallcurvature tunneling approximation31 (SCT) was employed for computing the multidimensional tunneling transmission coefficients. The minimum-energy path was computed with M05/6-311+G(2df,2p). Conformational searches and calculations of multistructural rovibrational partition functions with torsional anharmonicity32 were carried out at the same M05/6311+G(2df,2p) level of theory with the MSTor program;33 canonical variational transition-state theory calculations were performed with the Gaussrate34 and Polyrate35 programs. All of these calculations were performed with a density functional

integration grid of 99 radial shells and 974 angular points per shell. Solvation. Solvation was treated by continuum solvent models. All kinetics calculations included solvent effects by the solvation model based on density (SMD).36 The multireference diagnostics for the solution phase were based on the SMD model for the T1, B1, and GB1 diagnostics carried out with Gaussian 09 or the conductor-like screening model (COSMO)37 for the M diagnostics carried out with ORCA. All solvation free energies are standard-state values with a standard state of 1 M.



RESULTS AND DISCUSSION Diagnostics in the Gas Phase. The set of phenolic compounds investigated in this work are shown in Scheme 1. Scheme 1. Set of Phenolic Compounds Studied in This Work

This set comprises both monophenols and polyphenols, with OH groups at various positions. In addition, other substituents with electron donor or electron acceptor character are included, as well as one anionic species. Such structural variation is designed to include several of the structural features that can be present in phenolic compounds with antioxidant activity. In addition, it may help identifying if there is any relationship between the structure of the phenolic compound and the multireference character. Also, with that purpose in mind, four different peroxyl radicals were included in this investigation. They are HOO•, CH3OO•, t-butyl-OO•, and AIBN-OO•, where the last-named is the peroxyl radical generated from azobis(isobutyronitrile). CH3OO• was chosen because it is the one involved in a previous case study with multireference character,8 and the others because of the experimental kinetic data that is available for comparisons. In addition, HOO• is the free radical most frequently used in conjunction with the QMORSA protocol25 to quantify antioxidant activities. The T1 and M values used to assess the multireference character were obtained using the 6-311+G(d,p) basis set. Using a larger basis set would not be feasible at a reasonable computational cost for some of the largest systems investigated in this work. However, the effect of increasing the basis set was tested to ascertain the magnitude of the effect. The phenol (Ph1a) + CH3OO• reaction was chosen for this test, and the 4636

DOI: 10.1021/acs.jpca.5b07662 J. Phys. Chem. A 2016, 120, 4634−4642

Article

The Journal of Physical Chemistry A

Recall that the B1 diagnostic values reported in this work are averages of the B1 of phenolic O−H bond and ROO−H bond. When we look at the individual bonds, we find that the B1 values for all of these bonds are negative; therefore, the magnitudes of the averages are the same as the average of the magnitudes. However, in the case of the H2 + HO2 reaction, B1(H−H) = 0.51 kcal/mol while B1(OO−H) = −0.08 kcal/mol (computed with the 6-311+G(d,p) basis); thus, the magnitudeaveraged |B1| is 0.30 kcal/mol while the direct averaged B1 is 0.22 kcal/mol. Both of these B1 diagnostic values for H2 + HO2 reaction are positive, and this is because De of H−H bond predicted by BLYP (109.2 kcal/mol) is higher than that of B1LYP//BLYP (108.7 kcal/mol). The previously reported39 T1 diagnostic value for the transition-state structure of the H2 + HO2 reaction is 0.024, which is smaller than the 0.045 threshold for open-shell molecules; thus, this reaction is classified as single-reference. In this case, B1 is consistent with T1 because |B1| < 10 kcal/mol. The calculated diagnostics of the transition states of the investigated reactions are reported in Table 2. The first thing that stands out is that for most of the studied cases the results from the T1 and M diagnostics are in contradiction. While all the T1 values are below the 0.045 threshold for open-shell systems, the M values are larger than 0.10 for all but four cases. These four cases are the transition states of Ph2d with HOO• and CH3OO• in aqueous solution and those of Ph3a with HOO• in both gas phase and water. This means that assessing the extent of multireference character would depend on the diagnostic chosen for most of the investigated systems,

results are shown in Table 1, which shows that increasing the basis set from 6-311+G(d,p) to 6-311++G(2df,2p) has only a Table 1. Diagnostic Values Obtained with Different Basis Sets for the Phenol + CH3OO• Reaction in the Gas Phase 6-311+G(d,p) 6-311++G(2df,2p)

T1

M

B1

GB1

0.0445 0.0437

0.105 0.102

−1.24 −1.23

−11.33 −11.36

minor effect on the T1, M, B1, and GB1 values. Moreover, the values obtained from these diagnostics slightly decrease with the larger basis set. Therefore, the values reported in this work are reliable and do not appear to systematically underestimate the T1 and M values. As is usually the case, including more Hartree−Fock exchange in the exchange-correlation functional increases the calculated barrier height; thus, the hybrid B1LYP functional gives higher barrier heights than the local functional BLYP, and this results in negative GB1 diagnostic values for all cases in this article. The B1 diagnostics are not so easily rationalized. For hydrogen abstraction reactions by HO2 radical, the bond dissociation energy, De, correlates very well with the barrier heights;38 and De predicted by B1LYP tends to be higher than the one predicted by BLYP. Consequently we find that all B1 diagnostic values are also negative. This contrasts with previous work on B1 diagnostic values for transition-metal compounds, which were found to be positive.14 Table 2. Multireference Diagnostics phenol Ph1a Ph1a Ph1a Ph1a Ph1a Ph1b Ph1b Ph1b Ph1b Ph1c Ph1d Ph1d Ph2a Ph2a Ph2a Ph2a Ph2a(-H) Ph2b Ph2c Ph2d Ph2d Ph2d Ph2d Ph2e Ph2e Ph3a Ph3a Ph3b Ph3c

R R R R R R R R R R R R R R R R R R R R R R R R R R R R R

ROO•

medium

T1

M

B1

GB1

= = = = = = = = = = = = = = = = = = = = = = = = = = = = =

gas phase water gas phase water gas gas phase water gas phase water gas gas phase water gas phase water gas phase water water water water gas phase water gas phase water gas phase water gas phase water water water

0.041 0.038 0.045 0.037 0.044 0.038 0.038 0.043 0.038 0.032 0.041 0.038 0.039 0.032 0.036 0.032 0.025 0.033 0.031 0.038 0.030 0.035 0.029 0.038 0.029 0.037 0.036 0.032 0.027

0.105 0.104 0.110 0.105 0.113 0.105 0.103 0.107 0.104 0.111 0.113 0.111 0.102 0.101 0.102 0.101 0.064 0.102 0.100 0.101 0.067 0.101 0.067 0.102 0.069 0.098 0.074 0.102 0.097

−1.00 −1.47 −1.24 −1.75 −1.49 −1.07 −1.61 −1.31 −1.89 −2.08 −1.02 −1.34 −1.78 −2.05 −2.03 −2.33 −1.69 −1.74 −1.90 −1.53 −1.90 −1.77 −2.18 −1.32 −1.87 −2.04 −2.18 −1.96 −2.27

−12.11 −9.56 −11.33 −8.89 −10.83 −12.12 −9.45 −11.32 −8.47 −9.75 −12.34 −10.89 −14.21 −11.99 −13.27 −12.06 −8.30 −8.18 −8.05 −14.16 −8.52 −13.02 −8.34 −11.28 −8.51 −8.31 −8.02 −12.04 −9.10

H H CH3 CH3 t-butyl H H CH3 CH3 AIBN H H H H CH3 CH3 H H H H H CH3 CH3 H H H H H H

4637

DOI: 10.1021/acs.jpca.5b07662 J. Phys. Chem. A 2016, 120, 4634−4642

Article

The Journal of Physical Chemistry A

Table 3. T1 Values Obtained with Various Geometries (SMD vs Gas Phase)

provided that the currently accepted thresholds are used. This is only partly surprising because the unreliability of T1 diagnostics has been noticed by many workers in the past. Despite these differences, both diagnostics agree on the direction of the solvent effects, i.e., including continuum solvation effects in the calculations systematically decreases the values of both T1 and M diagnostics. In particular, we find that the multireference character in the gas phase is higher than that in solution for every combination of a phenolic molecule with an ROO• radical. The largest reduction in multireference character was found for methylated dihydroxybenzenes (Ph2d and Ph2e), regardless of the investigated isomer. In general, solvation of polar solutes is dominated by electrostatic polarization, so this is an electrostatic polarization effect. The average decrease due to solvation is 14% for the T1 diagnostic and 15% for the M diagnostic. This is relevant for the studies on the antioxidant activity of phenolic compounds because the reactions of interest take place in aqueous solution. Solvation is not the only aspect of a trend in the diagnostics that is the same for the T1 and the M diagnostics. The regiochemistry of substituent effects is another such aspect. For example, the more substituted phenols tend to have lower multireference character. This is in line with previous findings for a phenolic compound, canolol, with a larger number of substituent groups than those presented here.40 In that case it was found that the T1 value, corresponding to the transition state involved in the H abstraction from the phenolic moiety by HOO•, is lower than any of those reported in Table 2 (0.023). This is consistent with the present observation that the multireference character decreases as the number of substituent groups in the phenolic ring increases. In addition, for the less substituted phenols the multireference character of the transition state for the reactions involving CH3OO• is higher than for those involving HOO•. For the more substituted phenols, an opposite trend is found for T1, but the difference is small. All these trends support that the phenol (Ph1a) + CH3OO• reaction is among those within the studied set with higher multireference character, in line with previous findings.8 The other transition states with relatively high T1 and M values are those corresponding to Ph1b + CH3OO• and Ph1d + HOO•. The latter case is interesting when compared with the reaction between Ph1b and the same radical. Both compounds are monophenols; the difference is the nature of the substituent at the para position, electron donor in the case of Ph1b and electron acceptor in Ph1d. This indicates that the multireference character increases with the electron-withdrawing ability of the substituent, although the effect is only moderate. Another interesting comparison is that between the acid and base forms of catechol (Ph2a). In this case the differences in multireference character are rather large, with deprotonation significantly reducing this behavior. Diagnostics in Aqueous Solution. In an attempt to identify if the decrease in the multireference character when the calculations are performed using continuum solvent models is caused by the geometrical differences arising from optimization in vacuum versus optimization in aqueous solution, additional T1 calculations were performed for selected systems. They consist of using gas-phase geometries for SMD calculations, and vice versa (using SMD geometries for gas-phase calculations). The results of these tests are in Table 3, and they suggest that the changes in geometry are of minor importance, compared to the role of electrostatic polarization interactions with the

Ph1a

Ph1a

Ph1b

Ph2a

T1

geometry

HOO•

CH3OO•

HOO•

HOO•

gas gas water water

gas water water gas

0.041 0.0385 0.038 0.035

0.0445 0.037 0.037 0.033

0.038 0.038 0.0375 0.032

0.0385 0.039 0.032 0.033

solvent. Therefore, the latter is proposed as being responsible for reducing the multireference character of the studied transition states. General Conclusions on Diagnostics. When the gathered data is analyzed altogether, it becomes evident that assessing the multireference character for reactions between phenolic compounds and peroxyl radicals as a whole set, based on established diagnostics, is a rather difficult task. The amount of multireference character depends on the particular system under study and on the conditions, i.e., solvent and pH, under which the reaction takes place. Nevertheless, it is possible to make some generalizations. The most straightforward generalizations are that the multireference character is reduced by the presence of the solvent, by increasing the size of the phenolic compound, and by deprotonation in aqueous solution. For example, for the same system, M values in the gas phase are all greater than in the water phase, and the magnitudes of GB1 values have the same trend. The phenolic compounds with higher biological importance as antioxidants are substituted phenols, and the relevant conditions for such activity involve the presence of solvents. Therefore, the above-mentioned trends suggest that the multireference character of the transition states involved in antioxidants HT reactions with peroxyl radicals is not expected to be particularly high. In assessing the peroxyl radical scavenging activity of phenolic compounds using singledeterminant-referenced methods, the HOO• radical might be preferred because it leads to transition states with multireference character lower than that of larger radicals such as the CH3OO•. B1 diagnostics for all the cases studied here are far smaller than 10 kcal/mol in magnitude; therefore, B1 diagnostics lead to the conclusion that all of these hydrogen abstraction reactions are single-reference systems. When we compare the GB1 and M diagnostics, if we do not consider borderline cases (i.e., those where the M diagnostic value is close to 0.10 or the GB1 diagnostic value is close to 10 kcal/mol), the GB1 and M diagnostics are consistent, and they lead to the same conclusions. Kinetics. We have seen above that using the threshold values currently accepted for classifying the multireference character of transition states may lead to some contradictory assessments. Therefore, we address the original question behind this research more directly, that is, when is it safe to use KS-DFT with given exchange-correlation functionals to predict rate constants of the reactions between phenolic compounds and peroxyl radicals? We will compare calculated rate constants to experimental ones to answer this question. Unfortunately, experimental kinetic data on reactions of phenolic compounds with ROO• are rather scarce. For the studied compounds we could not find any gas-phase reaction rates, but we found data for five reactions in solution: 4638

DOI: 10.1021/acs.jpca.5b07662 J. Phys. Chem. A 2016, 120, 4634−4642

Article

The Journal of Physical Chemistry A Ph1a + t ‐butyl‐OO• , in isopentane solution

In all cases except reaction II, the two theoretical values in Table 4 bracket the experimental values. Reaction I has the transition state with the highest multireference character, according to both T1 and M diagnostics, and this is the only case where both estimates agree with experiment within a factor of 11; furthermore, in this case M05-2X leads to a rate constant that is closer to the experimental value than M05. Because the deviations from experiment are not largest for the cases where the multireference character is highest, we tentatively conclude that the multireference character in the transition states is not the dominant source of error in the rate constants. The conclusion is tentative because of the neglect in Table 4 of variational effects, anharmonicity, and multidimensional tunneling, as well as the approximate nature of the two exchangecorrelation functionals. Next, we investigate the effects of including variational effects, anharmonicity, and multidimensional tunneling. Because most of the rate constants computed by M05 in Table 4 agree very well with experiments (within a factor of 10) except for reaction IV (which is 33 times larger than the experimental rate constant), we selected reaction IV for these higher-level kinetics calculations. We used multistructural variational transition-state theory with the small-curvature tunneling approximation based on a minimum energy path computed with M05/6-311+G(2df,2p). The rate constant is computed as

(I)

Ph1c + AIBN‐OO• , in chlorobenzene solution

(II)

Ph2b + HOO• , aqueous solution (pH 0.5−1.5)

(III)

Ph2c + HOO• , aqueous solution (pH 0.4−3.5)

(IV)

Ph3b + HOO• , aqueous solution (pH 0.5−1.5)

(V)

It has been proposed that for systems with significant multireference character, the M06 and M05 exchange-correlational functionals perform better than their partners with higher proportions of exact exchange, that is, than M06-2X and M052X.8,41,42 On the other hand, some benchmark studies testing the performance of different functionals for kinetics have demonstrated that for chemical reactions in general the trend is the opposite, i.e., the 2X variants have a significantly better performance for kinetics,41−52 provided that the studied systems do not involve transition metals. Accordingly, the kinetics data in this work has been obtained from calculations performed with both the M05 and the M05-2X exchangecorrelation functionals. First we consider the calculations by conventional transitionstate theory with Eckart-approximation tunneling. The rate constants (k) calculated in this work using M05 and M05-2X reaction barriers, as well as those experimentally obtained, for reactions I−V are reported in Table 4. The corresponding

MS‐T CVT/SCT k MS‐CVT/SCT = Fact k

(1)

CVT/SCT

where k is the canonical variational transition-state theory (CVT) rate constant (with the small-curvature tunneling approximation) based on the lowest-energy reaction path; the multistructural and torsional anharmonicity are included in the FMS‑T factor30,38 by the MS-T method; and the vibrational act anharmonicity is included by scaling all the normal-mode frequencies by a factor of 0.977.57 These calculations were performed in aqueous solution by using the SMD solvation model. MS‑T The computed Fact factors, tunneling transmission SCT coefficients κ , variational transmission coefficients ΓCVT, and MS-VTST/SCT rate constants at various temperatures for reaction IV are shown in Table 6. The variational effect is not

Table 4. Experimental and Conventional TST Rate Constants (M−1 s−1) in Solution at 298.15 K conventional TST reaction

M05-2X

M05

I II III IV V

× × × × ×

× × × × ×

2.9 1.1 4.4 2.4 6.8

0

10 102 101 103 100

7.4 1.3 4.3 2.8 1.2

exptl 1

10 103 103 105 104

3.2 9.4 4.1 8.5 2.3

× × × × ×

ref 1

10 103 103 103 103

53 54 55 56 55

reaction barriers in terms of conventional TST standard-state free energies of activation (ΔG‡), the imaginary frequencies of the transition state (ω‡), and the tunneling transmission coefficients (κ) are reported in Table 5. For all the studied cases except for reaction II the rate constants calculated using M05 are overestimated.

Table 6. Computed FMS‑T Factors, Tunneling Transmission act Coefficients κSCT, Variational Transmission Coefficients ΓCVT, and MS-VTST/SCT Rate Constants (in M−1 s−1) at Various Temperatures for Reaction IV (Ph2c + HOO• in Aqueous Solution) with the M05/6-311+G(2df,2p) Exchange-Correlation Functional

Table 5. Details of the Calculations in Table 4: Conventional TST Standard-State Quasiclassical Free Energy of Activation (ΔG‡,° in kcal/mol, at 298.15 K), Imaginary Frequency (ω‡ in cm−1), and Eckart-Approximation Tunneling Transmission Coefficient (Unitless κ at 298.15 K)a M05-2X ‡,



M05 ‡,

reaction

ΔG °

ω

κ

ΔG °

ω‡

κ

I II III IV V

19.2 18.1 19.9 16.0 21.6

2765i 2837i 2858i 2317i 2884i

58 322 1493 106 2383

16.0 15.0 15.9 12.3 16.9

1455i 1897i 2412i 2506i 3395i

6 19 146 25 1477

T (K)

FMS−T act

κSCT

ΓCVT

kCVT/SCT

273.15 298.15 300.00 350.00 373.15

1.49 1.49 1.49 1.49 1.48

27.0 16.6 16.0 8.0 6.3

0.98 0.98 0.98 0.97 0.97

4.3 5.9 6.0 1.1 1.4

× × × × ×

104 104 104 105 105

important (less than 3%) for this system; and the SCT tunneling transmission coefficient is close to the one computed based on Eckart tunneling (17 vs 25). The computed rate constant at 298.15 K is 5.87 × 104 M−1 s−1, which agrees within a factor of 6.9 with the experimental value 8.5 × 103 M−1 s−1; the computed and experimental phenomenological Gibbs energies of activation ΔGact ° (see below for definition) are

ΔG‡,° = G°(conventional transition state) − G°(phenol) − G°(ROO•)

a

4639

DOI: 10.1021/acs.jpca.5b07662 J. Phys. Chem. A 2016, 120, 4634−4642

Article

The Journal of Physical Chemistry A

that has previously been proposed to judge the accuracy of phenomenological Gibbs energies of activation for reactions in solution.51 The results support the better performance of the M05 functional for the studied reactions between phenolic compounds and peroxyl radicals. The MUE corresponding to the M05 calculations is much lower than the 2 kcal/mol criterion, while that of the M05-2X data is close to it. For the individual reactions, all the M05 deviations with respect to the experimental value are less than or equal to 2.1 kcal/mol, with the highest discrepancy corresponding to reaction IV (2.1 kcal/ mol) and the lowest to reaction II (