Subscriber access provided by CMU Libraries - http://library.cmich.edu
Article
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 AlvarezIdaboy, Junwei (Lucas) Bao, and Donald G. Truhlar J. Phys. Chem. A, Just Accepted Manuscript • Publication Date (Web): 17 Sep 2015 Downloaded from http://pubs.acs.org on September 22, 2015
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 25
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
Sep. 13, 2015. Prepared for J. Phys. Chem. A
Hydrogen Abstraction Reactions from Phenolic Compounds by Peroxyl Radicals: Multireference Character and Density Functional Theory Rate Constants Annia Galano,1 Leonardo Muñoz-Rugeles,2 Juan Raul Alvarez-Idaboy,2,∗ Junwei Lucas Bao,3 and Donald G. Truhlar3,∗ 1
Departamento de Química. Universidad Autónoma Metropolitana-Iztapalapa. San Rafael Atlixco 186, Col. Vicentina. Iztapalapa. C. P. 09340. México D. F. México. 2
Departamento de Física y Química Teórica, Facultad de Química, Universidad Nacional Autónoma de México, México DF 04510, México. 3 Department of Chemistry, Chemical Theory Center, and Supercomputing Institute, University of Minnesota, Minneapolis, Minnesota 55455-0431, USA.
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 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, by the size of the phenolic compound, and by 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. Keywords: rate constants, static correlation, free radical scavengers, transition states, activation energies
* To whom correspondence should be addressed. E-mail:
[email protected], Phone: +52 55 5622 3776;
[email protected], Phone: +1 612 624 7555. *
1
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 2 of 25
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 exchangecorrelation 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 is usually lower for systems with inherently multiconfigurational wave functions than 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 Since an antioxidant is defined as “any 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 to be 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 2
ACS Paragon Plus Environment
Page 3 of 25
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
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 specific 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 pointed out 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 on 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 coworkers 13 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 ]. 3
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 4 of 25
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 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), and so we will consider the B1 diagnostic. The B1 diagnostic is defined as follows:14 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 exchange-correlation 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, and 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, and so it is of great interest to study it. 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 4
ACS Paragon Plus Environment
Page 5 of 25
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
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, since 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 only carried out 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 kcal/mol. 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 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)//M052X/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-COP+π space includes 5
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 6 of 25
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-hydroxylbenzaldehyde, 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 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 since 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 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 multi–structural variational transition state theory30 (MS–VTST) calculations for the reaction of 1,4-benzenediol with HOO• in aqueous solution. The small-curvature tunneling approximation 31 (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 multi–structural rovibrational partition functions with torsional anharmonicity32 were carried out at the same M05/6– 311+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 6
ACS Paragon Plus Environment
Page 7 of 25
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
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. 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 lastnamed is the peroxyl radical generated from azobisisobutyronitrile. 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 QM-ORSA protocol25 to quantify antioxidant activities.
7
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 8 of 25
Scheme 1. Set of phenolic compounds studied in this work.
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 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 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 exchangecorrelation functional increases the calculated barrier height, and 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 8
ACS Paragon Plus Environment
Page 9 of 25
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
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 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), and thus the magnitude–averaged |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 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, and thus this reaction is classified as single–reference. In this case, B1 is consistent with T1 since |B1| < 10 kcal/mol. Table 1. Diagnostic values obtained with different basis sets for the phenol + CH3OO• reaction in the gas phase. T1
6-311+G(d,p) 6-311++G(2df,2p)
0.0445 0.0437
M 0.105 0.102
B1
GB1
-1.24 -1.23
-11.33 -11.36
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 4 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, provided that the currently accepted thresholds are used. 9
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 10 of 25
This is only partly surprising since the unreliability of T1 diagnostics have 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 is higher in the gas phase than 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.
10
ACS Paragon Plus Environment
Page 11 of 25
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
Table 2. Multireference diagnostics. Phenol
ROO•
Medium
T1
M
B1
GB1
Ph1a Ph1a Ph1a Ph1a Ph1a
R= H R= H R= CH3 R= CH3 R= t-butyl
gas phase water gas phase water gas
0.041 0.038 0.045 0.037 0.044
0.105 0.104 0.110 0.105 0.113
-1.00 -1.47 -1.24 -1.75 -1.49
-12.11 -9.56 -11.33 -8.89 -10.83
Ph1b Ph1b Ph1b Ph1b
R= H R= H R= CH3 R= CH3
gas phase water gas phase water
0.038 0.038 0.043 0.038
0.105 0.103 0.107 0.104
-1.07 -1.61 -1.31 -1.89
-12.12 -9.45 -11.32 -8.47
Ph1c
R= AIBN
gas
0.032
0.111
-2.08
-9.75
Ph1d Ph1d
R= H R= H
gas phase water
0.041 0.038
0.113 0.111
-1.02 -1.34
-12.34 -10.89
Ph2a Ph2a Ph2a Ph2a
R= H R= H R= CH3 R= CH3
gas phase water gas phase water
0.039 0.032 0.036 0.032
0.102 0.101 0.102 0.101
-1.78 -2.05 -2.03 -2.33
-14.21 -11.99 -13.27 -12.06
Ph2a(-H)
R= H
water
0.025
0.064
-1.69
-8.30
Ph2b
R= H
water
0.033
0.102
-1.74
-8.18
Ph2c
R= H
water
0.031
0.100
-1.90
-8.05
Ph2d Ph2d Ph2d Ph2d
R= H R= H R= CH3 R= CH3
gas phase water gas phase water
0.038 0.030 0.035 0.029
0.101 0.067 0.101 0.067
-1.53 -1.90 -1.77 -2.18
-14.16 -8.52 -13.02 -8.34
Ph2e Ph2e
R= H R= H
gas phase water
0.038 0.029
0.102 0.069
-1.32 -1.87
-11.28 -8.51
Ph3a Ph3a
R= H R= H
gas phase water
0.037 0.036
0.098 0.074
-2.04 -2.18
-8.31 -8.02
Ph3b
R= H
water
0.032
0.102
-1.96
-12.04
Ph3c
R=H
water
0.027
0.097
-2.27
-9.10
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 11
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 12 of 25
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 is higher for the reactions involving CH3OO 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 vs. 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 solvent. Therefore the latter is
12
ACS Paragon Plus Environment
Page 13 of 25
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
proposed as responsible for reducing the multireference character of the studied transition states. Table 3. T1 values obtained with various geometries (SMD vs. gas phase). Ph1a
Ph1a
T1
geometry
HOO
gas gas water water
gas water water gas
0.041 0.0385 0.038 0.035
CH3OO 0.0445 0.037 0.037 0.033
Ph1b
Ph2a
HOO
0.038 0.038 0.0375 0.032
0.0385 0.039 0.032 0.033
HOO
General conclusions on diagnostics. Analyzing the gathered data 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 gas phase are all bigger than in 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 single-determinant-referenced methods, the HOO• radical might be preferred since it leads to transition states with lower multireference character than larger radicals such as the CH3OO•. B1 diagnostics for all the cases studied here are far smaller than 10 kcal/mol in magnitude, and therefore B1 diagnostics lead to the conclusion that all of these hydrogen 13
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 14 of 25
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: (I)
Ph1a + t-butyl- OO•, in isopentane solution
(II)
Ph1c + AIBN-OO•, in chlorobenzene solution.
(III)
Ph2b + HOO•, aqueous solution (pH = 0.5-1.5)
(IV)
Ph2c + HOO•, aqueous solution (pH = 0.4-3.5)
(V)
Ph3b + HOO•, aqueous solution (pH = 0.5-1.5)
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 M05-2X.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 exchange-correlation functionals. First we consider the calculations by conventional transition state theory with Eckartapproximation 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 to V are reported in Table 4. The corresponding reaction barriers in terms of conventional TST 14
ACS Paragon Plus Environment
Page 15 of 25
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
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 4. Experimental and conventional TST rate constants (M-1 s-1) in solution at 298.15 K.
Reaction I II III IV V
Conventional TST M05-2X M05 2.9E+00 7.4E+01 1.1E+02 1.3E+03 4.4E+01 4.3E+03 2.4E+03 2.8E+05 6.8E+00 1.2E+04
Exp.
Ref.
3.2E+01 9.4E+03 4.1E+03 8.5E+03 2.3E+03
53 54 55 56 55
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). Reaction I II III IV V
M05-2X
M05
ΔG‡,°
κ ΔG‡,° κ ω‡ ω‡ 19.2 2765i 58 16.0 1455i 6 18.1 2837i 322 15.0 1897i 19 19.9 2858i 1493 15.9 2412i 146 16.0 2317i 106 12.3 2506i 25 21.6 2884i 2383 16.9 3395i 1477 , • ‡ ΔG ° = G°(conventional transition state) – G°(phenol) – G°(ROO )
In all cases except reaction II, the two theoretical values in Table 4 bracket the experimental ones. Reaction I has the transition state with the highest multireference 15
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 16 of 25
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. Since 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 exchange-correlation functionals. Next, we investigate the effects including of variational effects, anharmonicity, and multidimensional tunneling. Since most of the rate constants computed by M05 in Table 4 agrees 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 multi–structural variational transition state theory (MS–VTST) with the small–curvature tunneling approximation (SCT) based on a minimum energy path (MEP) computed with M05/6–311+G(2df,2p). The rate constant is computed as MS-T k CVT/SCT k MS-CVT/SCT = Fact
(1)
where k CVT/SCT is the canonical variational transition state theory (CVT) rate constant (with small curvature tunneling approximation) based on the lowest–energy reaction path; MS–T the multi–structural and torsional anharmonicity are included in the Fact factor30,38 by
the MS–T method; and the vibrational 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 coefficients κ SCT , 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 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 16
ACS Paragon Plus Environment
Page 17 of 25
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
value 8.5×103 M-1 s-1; the computed and experimental phenomenological Gibbs energies of ° activation ΔGact (see below for definition) are 10.9 and 12.1 kcal/mol at 298.15 K, which
are also in good agreement. In order to investigate whether the differences in geometries of transition states optimized with M05 and M05-2X functionals influence the multireference character of these species, the corresponding T1 and M diagnostics were compared for reaction IV. Their values were found to be 0.031 (0.025) and 0.100 (0.100) for the M05-2X (M05) TSs. Once again we see that the behavior depends on the chicer of diagnostic. In particular, the geometry differences between the TSs optimized with these two functionals lead to only minor differences in the M diagnostic, but the T1 value is significantly lower for the M05 TS. We also carried out some calcuations with the M08-HX exchange-correlation functional. In the gas phase, the M08–HX functional usually shows better performance for reaction barrier heights than M05,45 but for the present reaction it seriously underestimated the rate constant. It is not clear whether this is due to errors in the implicit solvation model we used in the calculation or due to M08-HX, which – like M05-2X and M06-2X – has a high fraction of Hartree-Fock exchange, giving an inaccurate barrier height due to multireference character in this reaction. The classical barrier height computed by M08– HX/6–311+G(2df,2p) is 10.9 kcal/mol, which may be compared to 5.9 kcal/mol by M05/6– 311+G(2df,2p). MS–T Table 6. Computed Fact factors, tunneling transmission 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. T/K
MS–T Fact
273.15
1.49
298.15
Γ CVT
k MS-CVT/SCT
27.0
0.98
4.3E+04
1.49
16.6
0.98
5.9E+04
300.00
1.49
16.0
0.98
6.0E+04
350.00
1.49
8.0
0.97
1.1E+05
κ SCT
17
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 18 of 25
373.15
1.48
6.3
0.97
1.4E+05
Since kinetics are often discussed in terms of the free energy of activation, we have also analyzed the rate constants in those terms for reactions I to V. In particular, the ° phenomenological Gibbs energies of activation ( ΔGact ) were derived from the rate
constants (calculated and experimental) are58
"k T % ° ΔGact = RT ln $$ B '' # hkC 0 & where k is the calculated or experimental rate constant, R is the gas constant, T is the temperature (in K), kB and h are the Boltzmann and Planck constants, and for bimolecular reactions C 0 is the concentration in the standard state, taken here as 1 molar in solution. The values obtained this way are in Table 7, and they make the trends easier to see and statistically quantify than the trends in the rate constants. For the tested reactions the MUE values were found to be 1.0 and 2.2 kcal/mol for the M05 and M05-2X calculations, respectively, which is within the uncertainties expected45 for barrier height calculations with these approximate exchange-correlation functionals as well as being a within a margin of 2 kcal/mol (twice the widely used value of “chemical accuracy”) 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 (
