10988
J. Phys. Chem. A 2010, 114, 10988–10996
A Quantum Mechanical Study of Methacrylate Free-Radical Polymerizations Matthew D. Miller† and Andrew J. Holder*,‡ AstraZeneca, Waltham, Massachusetts 02451, and UniVersity of Missouri-Kansas City, Kansas City, Missouri 64110 ReceiVed: May 8, 2010; ReVised Manuscript ReceiVed: August 30, 2010
Theoretical calculations at the semiempirical and ab initio levels of theory have been completed for a series of methacrylate compounds reported in the literature as measured by pulsed-laser initiated polymerization in conjunction with size-exclusion chromatography (PLP-SEC). Modeling includes calculation of the Gibb’s free energies (∆G‡) and activation energies (Ea). These results were then compared to experimental results. Semiempirical ∆G‡ using AM1-CI calculations successfully predicted relative activation energies (R2 ) 0.89). HF and DFT methods more accurately predicted absolute activation energies, but the relative values were less reliable. Accurate quantitative structure property relationship (QSAR) models for propagation rate coefficients, kp, were developed using AM1-CI and DFT. The semiempirical model included two charge descriptors, partial negatively charged surface area (PNSA), and minimum net atomic charge for oxygen (R2 ) 0.959). The DFT information, which included two quantum chemical descriptors (1-electron reactivity index for carbon and point charge component of the molecular dipole) calculated from the ground state structure, had improved statistics (R2 ) 0.979). A second DFT model is reported for 10 hydrocarbon methacrylate structures based on the 1-electron reactivity index for carbon (R2 ) 0.979). Theoretical results were also analyzed to provide an explanation for the unexpectedly large experimental kp values observed in the case of larger methacrylate monomers. 1. Introduction Radical reactions involving methacrylates are important in a variety of products including dental composites, and paints. In the form of hydrogels, these materials serve in biomedical and biotechnological applications such as drug delivery systems, contact lenses, and food technology.1 The reaction mechanism is understood, but the influence of substituents on the reactive group continues to defy a consistent explanation. This information is obviously critical for a systematic molecular design process. Therefore, a quantitative structure activity relationship (QSAR) model focusing on improving propagation rate constants would be valuable in screening new and untested compounds prior to laboratory synthesis. Such a model has not been developed in the past, likely due to the broad variability in the experimentally reported propagation rates (kp). This is apparent from kp values for the same compounds in such references as the Polymer Handbook,2 even under seemingly identical experimental conditions. Recently, a new procedure has been developed that enables reproducible rate constant measurements by combining pulsedlaser initiated polymerization and size exclusion chromatography (PLP-SEC).3 Briefly, photopolymerization is initiated by a laser pulse. The large number of free-radicals produced by each pulse may be studied under systematic variation of the pulse cycle time and the number of propagation steps. Further details of this method can be obtained elsewhere.4-7 Advantages of the method include reproducibility, reduced complications from side reactions,7,8 and a number of self-consistency checks ensuring accurate and precise measurements. The reliability of PLP-SEC * To whom correspondence should be addressed. Phone: 816-235-2293; Fax: 816-235-6543;
[email protected]. † AstraZeneca. ‡ University of Missouri-Kansas City.
for linear polymer systems is further increased by combining the results with FT-NIR.9 PLP-SEC is now the method of choice for kp measurement as designated by IUPAC.10,11 PLP-SEC data for a set of 13 methacrylates have been reported at 50 °C, including activation energies (Ea).5 Experimental data were obtained for bulk homopolymerization reactions performed at ambient pressure. Additional details may be obtained from the original reference. Theoretical values from semiempirical and ab initio quantum mechanical models are comparable with one another and with the PLP-SEC experimental data, providing an indication of the computational methods’ quantitative and qualitative reliability. In addition, potential relationships between ∆G‡, Ea, and kp have been investigated here. The reliability of PLP-SEC kp measurements also now provides data of sufficient quantity and quality for deriving QSAR models, elucidating the relationship between monomer structures and kp results. Although there are published QSAR models on free-radical polymerization chain-transfer constants12 and glass transition temperatures for acrylates and methacrylates,13 the present study is the first to predict kp for methacrylates. The models presented here are based on descriptors computed using information from semiempirical and ab initio quantum mechanical calculations. The quantum mechanical approach also made possible a comparison of monomer structures and kp values, aimed at determiners of reactivity, such as exothermicity, radical stability, and ionization energies/electron affinities. 2. Experimental Methods Semiempirical calculations were performed with the AM114 method using AMPAC 8.16 with Graphical User Interface.15 We modeled reactions for the dimerization of each compound producing an initial radical structure by neutral protonation of
10.1021/jp104198p 2010 American Chemical Society Published on Web 09/27/2010
Methacrylate Free-Radical Polymerizations the terminal alkene of the reactive group with doublet spin. Developing similar results for larger oligomers (trimers or tetramers) will probably not impact the calculated energies, so this was not attempted. Coote’s16 results for free-radical acryzhylonitrile and vinyl chloride polymerization support this procedure, since only modest effects on the activation energy (less than 0.5 kcal/mol) using B3-LYP were observed. The procedure to locate a transition state involved a reaction path with a bond length constrained to incremental steps between the product and reactants while all other bond lengths remained free to optimize. Upon location of a plausible transition state (TS) geometry, a gradient minimization was performed on the completely unconstrained system. The potential activated complex was characterized by carrying out a frequency calculation. Connection of the reactant with product via this TS was confirmed by the presence of a single negative eigenvalue along an eigenvector motion showing an appropriate trajectory. The proposed transition state was allowed to follow the negative eigenvector in both directions until minima were reached. The proposed product was then fully optimized with an energy minimization and frequency calculation to ensure a true local minimum with no negative eigenvalues.17 The same procedure was followed for the reactant complex for semiempirical calculations. This procedure was repeated multiple times for methyl methacrylate, with the lowest transition state energy assumed to be the most likely transition state structure. The syndiotactic mode of attack was followed throughout the entire series of methacrylates, as it is the likely mechanism.18 Conformational analysis was performed on each individual monomer and radical structure to identify the global minimum geometry by systematic dihedral rotation about rotatable bonds. For the series of analogous reactions, appropriate substituents replaced the methyl group on the methyl methacrylate transition state in a similar orientation to the global minimum of the individual structure. In cases where the transition state for the modeled reaction did not lead to a single negative eigenvalue and appropriate eigenvector upon gradient minimization, the procedure was repeated for the full methacrylate and radical monomer in a similar fashion to that used to locate the lowest energy transition state structure for methyl methacrylate. As above, calculations were performed on the TS for each compound, confirming correct products and reactants. Ab initio calculations were performed with Gaussian03.19 Beginning with AM1 structures for each reaction as starting geometries, TSs and products were optimized at the corresponding levels of theory. Frequency calculations were performed as before to ensure a single negative eigenvalue for the TS geometry and zero negative eigenvalues for the product. Default convergence criteria were maintained for these calculations. Thermal corrections and ZPEs were provided by frequency results for each method. These were scaled for the 6-31G(d) basis set20 using 0.913521 for UHF calculations and 0.980422 for UB3-LYP23,24 functionals. This DFT method with 6-31G(d) has been relatively successful for radical reactions25 and was expected to provide reasonable prediction of energies. The reaction rate is directly related to the free energy of reactions shown in eq 1, where kB is the Boltzmann’s constant, T is the relevant temperature, h is Planck’s constant, ∆G‡ is the free energy of activation (the difference between the TS and the sum of the individual reactants), and R is the universal gas constant. Application of Maxwell’s relationship (eq 2) separates the rate into entropic (∆S‡) and enthalpic (∆H‡) portions. This expression is commonly condensed to the Arrhenius
J. Phys. Chem. A, Vol. 114, No. 41, 2010 10989 equation, eq 3, where Ea is calculated from eq 4 and m is the molecularity of the reaction.26
)
(1)
∆G‡ ) ∆H‡ + T∆S‡
(2)
k)
(
kB T -∆G‡ exp h RT
( )
k ) A exp
-Ea RT
Ea ) ∆H‡ + mRT
(3) (4)
Experimental ∆G‡ values were calculated from eq 1 using reported experimental rate values at 323 K. Semiempirical and ab initio ∆G‡ were calculated from eq 2. For ab initio methods, theoretical Ea were calculated at T ) 323 K from eq 4 with m ) 1. Semiempirical Ea were calculated without the correction term mRT, as this is included in the parametrization. Enthalpies of reaction (∆Hrxn) were calculated from the difference in energies between the products and reactants. Reactant energies were calculated as the sum of the energies for each individual reactant. This was not possible for the AM1CI method, as the method requires a priori selection of the CI active orbitals. For consistency, these orbitals included only the highest doubly occupied molecular orbital (DOMO), the singly occupied molecular orbital (SOMO), and the lowest unoccupied molecular orbital (LUMO). Inclusion of additional virtual orbitals has a substantial effect on both energy and wave function, and as the ground state monomer has no SOMO orbitals, we defined the reactants as a molecular complex for the purpose of calculating activation barriers and reaction enthalpies rather than the sum of energies for the individual reactants as before. Descriptors for QSAR development were calculated using the program CODESSA27 on previously described structures. A heuristic algorithm28 was used to derive several correlations from these descriptors resulting in a final set of multilinear regression equations. For ab initio calculations, Gaussian03 required the keywords “density ) current” and “pop ) (regular, nboread)” with “bndidx” included in the NBO section to generate data for subsequent QSAR analysis. The “CODESSA” keyword was used in AMPAC for the same purpose. Overall correlations were evaluated using R2, adjusted R2adj, cross-validated R2CV (leave one out method), the F-test, and t-test. Variance inflation factors (VIF < 5) and significance (p < 0.01) were also calculated for each descriptor using SPSS 11.5 to ensure orthogonality and non-colinearity of the descriptors with one another. In our approach, correlations are developed using a forward stepwise method,29 beginning with a single descriptor and adding descriptors as long as criteria of chemical sensibility, inclusion of additional variance, and orthogonality of the descriptor are satisfied. 3. Results and Discussion 3.1. Activation Energies. Several comparisons of theoretical versus experimental activation energies for free-radical polymerization reactions have already been published.16,30-33 These show that a higher level of theory, such as G3(MP2)-RAD reduces error in the calculated activation energies when compared with lower levels of theory like B3-LYP.16,30,32 However, these computationally intensive studies were limited to smaller data sets. In this case relatively large molecules (such as dodecyl
10990
J. Phys. Chem. A, Vol. 114, No. 41, 2010
Miller and Holder
TABLE 1: Calculated Gibb’s Free Energies of Activation (∆G‡, kcal/mol) for Methacrylates AM1 CI
BOMA (isobornyl) BUMA (butyl) BZMA (benzyl) CHMA (cyclohexyl) DMA (dodecyl) EHMA (2-ethylhexyl) EMA (ethyl) GMA (glycidyl) HEMA (hydroxyethyl) HPMA (2-hydroxypropyl) IBUMA (isobutyl) IDMA (isodecyl) MMA (methyl) mean error
21.77 21.88 24.89 21.94 21.17 20.74 22.04 21.42 21.57 21.21 21.45 21.95 21.99 7.35
AM1 exp UHF UHFa UB3LYPa ∆G‡b
exp Eac
16.12 17.96 19.60 17.05 18.23 16.27 18.07 17.73 19.11 19.55 17.86 17.80 21.32 3.71
5.38 5.47 5.14 5.33 5.02 4.88 5.59 5.23 5.23 4.97 5.21 4.97 5.35
26.82 28.80 28.06 28.55 28.79 28.56 28.84 27.93 28.57 28.40 28.65 28.80 28.00 13.87
23.17 22.19 22.38 21.92 23.30 23.33 22.90 22.37 22.93 21.92 22.57 24.56 22.32 8.26
14.80 14.78 14.70 14.53 14.67 14.55 14.56 14.40 14.37 14.39 14.52 13.92 14.26
a Calculated with 6-31G(d) basis set. b Calculated from ref 5 based on eq 1. c Reported values based on ref 5. In the case of CHMA, GMA, and BZMA, multiple values are reported. Data listed came from the experimental lab that determined kp based on Arrhenius fit covering lower temperatures to reduce potential of depropagation. According to ref 5, depropagation should be negligible at T < 90 °C.
Figure 1. Relationship between the AM1 Gibb’s free energy of activation (∆G‡) and the experimental Ea.
methacrylate) are being considered, and it is impractical to locate a transition state for that system at the G3(MP2)-RAD level of theory. The Gibb’s energies of activation, ∆G‡, are reported in Table 1 for AM1-CI, AM1-UHF, and the ab initio methods. AM1UHF provides ∆G‡ barriers that are closer to the absolute experimental ∆G‡ than other methods, with a mean absolute error of 3.7 kcal/mol. The AM1-CI calculated ∆G‡ values are slightly better than the DFT results and are substantially improved versus UHF 6-31G(d). The relationship between experimental and calculated ∆G‡ values was generally poor in relative terms, although a reasonable trend exists between AM1CI ∆G‡ and experimental Ea. Following the removal of two outliers, a correlation was located where R2 ) 0.89 (Figure 1). As previously mentioned, only the highest DOMO, SOMO, and LUMO were included in the CI manifold. It is not surprising that BZMA was poorly predicted since the four highest DOMOs and four lowest unoccupied MOs are less than 0.1 eV apart. For comparison, nonaromatic compounds often have >0.6 eV between levels. The second outlier, IDMA, had a vibrational
TABLE 2: Calculated Activation Energies (Ea, kcal/mol) for Semiempirical and ab Initio Methodsa Ea BOMA BUMA BZMA CHMA DMA EHMA EMA GMA HEMA HPMA IBUMA IDMA MMA mean error
AM1 CI
AM1 UHF
14.57 14.94 16.24 15.80 14.77 14.74 14.63 15.38 14.88 16.90 16.17 14.94 14.19 10.03
11.57 11.69 12.85 11.69 12.35 11.84 11.67 11.66 13.30 14.83 11.66 11.70 12.99 7.08
c
UHF
11.59 12.97 12.14 12.87 12.96 13.06 13.03 12.53 12.95 12.69 13.02 12.96 12.49 7.50
c
UB3LYP 7.08 7.28 6.91 7.36 7.14 7.60 7.41 7.18 7.18 6.35 7.28 7.28 7.16 1.96
b
Exp.
5.38 5.47 5.14 5.33 5.02 4.88 5.59 5.23 5.23 4.97 5.21 4.97 5.35
a Ab initio methods were calculated based on thermally corrected structure energies. b Reference 5. c Calculated with the 6-31G(d) basis set.
frequency that is much lower than other reactant complexes (1.49 cm-1), including DMA. This is likely due to a flat potential energy surface near the TS. Although the absolute AM1-CI ∆G‡ values are uniformly poor, they appear to be useful for Ea trends. The relationship between the AM1-CI ∆G‡ and experimental Ea values is counterintuitive. According to eq 4, Ea should correlate better with ∆H‡ rather than ∆G‡. It is possible that the inclusion of entropy in the AM1 calculations corrects for isolated methacrylate monomers versus a bulk, condensed phase to some extent, although this likely is pure coincidence. However, the relationship between ∆G‡ and experiment is worth noting. This was the only reasonable trend found between any of the calculated ∆G‡ barriers and experimental ∆G‡ or Ea data. Calculated activation energies (Ea) are listed in Table 2 for the semiempirical and ab initio methods. It was previously pointed out that the size of linear substituents has little effect on Ea.34 The distance between the substituents and the reactive site reasonably explains this observation. Experimentally, Ea is relatively invariant with respect to solvent for alkyl methacrylates MMA and BOMA.35 Given the general lack of solvent influence on Ea for alkyl methacrylates, gas phase calculations for radical reactions at sufficient levels of theory will likely provide good estimates of the TS in the condensed phase.36 Semiempirical and UHF 6-31G(d) Ea overestimated the magnitude of Ea significantly with no discernible trend. DFT values are much closer to experiment, with an error of only 1.96 kcal/ mol. This is similar to other results in the literature demonstrating overestimation of Ea by UHF versus DFT for radical reactions.18 Although the mean error for all methods exceeds the entire range of the experimental data (0.71 kcal/mol), with the relative trends apparently masked by error, the DFT results provide acceptable Ea as compared to experiment. 3.2. QSAR. 3.2.1. Semiempirical QSAR Model. There is a normal distribution of kp data covering 648-1504 L/mol · s for 12 of the 13 data points addressed here (experimental values listed in Table 3). The thirteenth data point, HEMA, had a value of 2563 L/mol · s, and we did not use it in the training set to prevent artificial inflation of statistics. Descriptors were calculated from: (1) the ground state form of each compound, (2) the radical form (based on AM1-CI), (3) the transition state, and (4) reactant complex structures. The descriptors were then collected to create a pool of descriptors from all four forms. A
Methacrylate Free-Radical Polymerizations
J. Phys. Chem. A, Vol. 114, No. 41, 2010 10991
TABLE 3: Calculated Enthalpies of Reaction (∆Hrxn, kcal/mol) and Reaction Rates (L/mol · s) for Methacrylate Homopolymerizations BOMA BUMA BZMA CHMA DMA EHMA EMA GMA HEMA HPMA IBUMA IDMA MMA range
AM1 CI
UHFa
UB3LYPa
rateb
-12.98 -12.43 -16.44 -11.67 -12.59 -12.81 -12.45 -12.06 -12.48 -10.65 -11.34 -12.43 -12.84 c 2.33
-13.63 -11.59 -12.30 -11.53 -11.60 -11.64 -12.07 -12.07 -11.57 -11.71 -11.42 -11.60 -12.05 2.22
-12.63 -12.44 -13.12 -12.29 -12.49 -12.21 -12.60 -12.60 -12.48 -13.21 -12.40 -12.31 -12.57 1.00
1002 757 1224 1212 995 944 671 1273 2563 1504 794 959 648 4-fold
a Calculated with the 6-31G(d) basis set. b Reported values based on ref 5. c Range excludes BZMA.
successful model (R2 ) 0.959) was derived including a descriptor calculated from the reactants complexes (RCT) and transition state (TS) structures (eq 5). Both descriptors are based on electrostatics and have negative coefficients and descriptor values. The negative descriptor values and coefficients lead to a net increase in kp for an increase in the descriptor value. The first descriptor, RCT-PNSA2Zefirov37 or partial negatively charged surface area, is calculated from the product of the sum of all negative charges of atoms with the sum of the solvent accessible surface areas for these atoms (eq 6). Although PNSA has been found to correlate with electrophilicity,38 in this instance it more likely describes the size of the substituents. As more carbon atoms are introduced into each monomer, the total sum of the negative charges of the molecule increases slightly while the solvent-accessible surface area of the molecules increases substantially. The second descriptor is the minimum (most negative) electrostatic potential charge for an oxygen atom for the TS complex. This descriptor is a quantum chemical descriptor that is usually centered in the methacrylate portion of the molecule. The descriptor has lower values for CHMA and HPMA transition states and this change may explain the inability for the first descriptor to predict the kp if PNSA2Zefirov does describe the size of the molecule.
kp ) -4.9(RCT-PNSA2Zefirov) 4975.7(TS-ESP min net atomic charge for O)-1163.3 (5)
PNSA2Zefirov )
( ∑ surface areai)Qtotal
(6)
Although the ratio of the number of variables (descriptors) to observations (test compounds) is reasonable for QSAR development (>1:5),39 the large number of descriptors increases the risk of pure chance for what appears to be a statistically valid model.40,41 Unfortunately, the small number of compounds available to us did not allow the creation of separate training and test sets. To further test the validity of the model, the experimental data was scrambled (using a perl script) which required the new value to be different from the reported measurement. This process was repeated six
Figure 2. Semiempirical model for propagation constant kp including two descriptors. Statistics include R2 ) 0.959, R2cv ) 0.911, R2adj ) 0.950, F ) 105, t-test > 6.7, p-significance < 0.001, RMSE ) 59, VIF ) 1.0, N ) 12.
times, half as many as the total number of data points. No statistically valid model was found following this technique with the entire pool of descriptors (listed in Supporting Information). This result shows that the likelihood of chance correlations is low. Prediction of the property value for the small molecule HEMA (not included in the training set) using eq 5 resulted in a calculated kp value of 1260 L/mol · s, far from the reported value of 2563. Although this deviation is unacceptable, the nonmethacrylate portion of the compound contains a hydroxyl group in close proximity to the reactivity site. As this type of structure is represented by only one compound in the training set (HPMA), it is not surprising that HEMA would not be predicted well. Beuermann and Nelke suggested hydrogen bonding may be responsible for the enhanced rate of these methacrylates,42 and with only two alcohol compounds in the data set, such a descriptor was not included in the model and poor predictive ability is not surprising. 3.2.2. DFT QSAR Model. Similar to the procedure used in the development of a semiempirical model, descriptors were calculated for the ground state, radical, transition state, and reactant complex using the DFT level of theory. For 12 structures (without HEMA), a statistically valid 2-descriptor model (R2 ) 0.979) was developed from the ground state monomers (eq 7). The first descriptor (point charge component of the molecular dipole, µpoint-charge) is charge related and describes the polarizability of the compounds. The second descriptor, the average one-electron reactivity index for carbon (1ERIC) is based on the atomic orbital coefficients of the HOMO and LUMO of each carbon atom. The largest reactivity indices are located on the alkene carbon atoms of the methacrylate group, and other carbon atoms in the molecule are generally at least 2 orders of magnitude smaller. This descriptor is the average of all reactivity indices of carbon atoms and is directly related to the size of the molecule. As the descriptor value decreases with more carbon atoms, lower kp values are predicted due to the negative coefficient. The coefficient is large due to the magnitude of the descriptor values ( 6.2, p-significance < 0.001, RMSE ) 42, VIF ) 1.1, N ) 12.
Figure 4. DFT QSAR model for propagation constant kp including only hydrocarbon substituent methacrylates with one descriptor. Statistics include R2 ) 0.963, R2cv ) 0.934, R2adj ) 0.958, F ) 206, t-test > 14.3, p-significance < 0.001, RMSE ) 42, VIF ) 1.0, N ) 10.
descriptor was possible (R2 ) 0.979) when the training set included only hydrocarbon methacrylates (eq 8).
kp ) 399.6(µpoint-charge) 157730(avg one-electron rct index for C) + 306.9 (7) kp ) -168028(avg one-electron rct index for C) + 1116.6 (8) As in the semiempirical correlation, randomization of the property value was performed six times to test the validity of each model. No statistically valid model was found for either set of compounds with the scrambled propagation rate values. External testing of HEMA resulted in a predicted rate constant of 1482 based on eq 7, the poor prediction of which has been discussed for the semiempirical model. Both computational methods predict propagation rates with reasonable accuracy. As in the experimental data, there is little variation for cyclic ester groups, with a broader distribution for linear alkyl groups.6 The advantage of rapidly calculating
Miller and Holder geometries for AM1 is substantial versus DFT, although the results provided by the DFT model are closer to experiment. 3.3. Calculation Results to Explain Propagation Rates. Although the QSAR models provide predictions, the model variables describe the structures in the model rather than explaining the increase in propagation rate found with larger monomers. A variety of causes influence either the reactivity and/or activation energy of a particular compound in a homologous series of structures such as the methacrylates in this study. Unfortunately, the DFT calculated and experimental activation energies for the 13 methacrylates in this study varied by only 1.25 and 0.71 kcal/mol, respectively. We had hoped that a comparison between the polymerization rate constants and calculated activation energies would provide insight into the much larger (4 fold) range of kp rate constants for these reactions. Reaction rates can be affected by steric effects, charge transfer, resonance stabilization of the radical reactants or products, and enthalpic effects. The importance of each varies depending on the particular reaction species.18,25,26,43-45 Of these effects, reaction enthalpy by nature is actually a function of the other three and may therefore be a combination of their collective influence.43 Steric hindrance of the reactive site has obvious implications, with bulkier substituents affecting the pre-exponential factor (A) in eq 3 and reducing the reaction rate, all other factors being equal. For methacrylate polymerization reactions, a consensus has not been reached as to what portion of the reaction rate can be attributed to either the activation energy or to the preexponential factor.7,46,47 More noteworthy, the size of the substituent on methacrylate (as well as acrylate) compounds corresponds to an increase in the experimentally determined reaction rate.48 This is not only counterintuitive, but is opposite of the trend observed for other radical polymerization reactions, such as itaconate compounds.49 A proposed explanation is that smaller methacrylates have increased polar interactions near the charge centers of growing polymer chains and prevent them from reacting by effectively lowering the bulk monomer concentration. In contrast, larger monomers such as dodecyl methacrylate with longer alkyl chains do not show this dilution, effectively shielding the apparent reduction in local monomer concentration.34 Thus, steric hindrance may exercise an indirect correlation to the rate and effectively increase the local monomer concentration. These were minimized in this study by constraining the transition states through syndiotactic attack, in which substituents were on opposite sides of the active site, well separated from intramolecular interactions. Huang et. al18 showed that syndiotactic attack matches better with experimental results via this mechanism. Many studies have found a relationship between activation energies and reaction enthalpies for radical reactions.25,26,50 This is based on the early work for proton transfer reactions by Bell51 and on the interrelationships between reaction rates, exothermicities, and activation energies of Evans and Polanyi.52 The effect of reaction enthalpy on activation energy should be more prominent when the TS shows a strong resemblance to the product;18 that is, the reaction is significantly exothermic and the distance between the atoms of the bond formed resembles the distance of the same atoms at the TS. However, the DFTpredicted bond length formed in the dimer product was 1.57 Å, whereas the distance between the same two atoms at the TS was 2.26-2.27 Å. The lack of similarity between TS and products would indicate little resemblance between the two and any relationship between reaction enthalpies and reaction rate seems unlikely. The DFT calculated exothermicity of MMA
Methacrylate Free-Radical Polymerizations (-12.6 kcal/mol) is reasonably close to experiment (-13.0 kcal/ mol)43 and similar experimental exothermicities for the other methacrylates seems reasonable for such a homologous series of reactions. The calculated ∆Hrxn for all methods have only a small range of values (Table 3), further supporting a minor enthalpic effect on the reaction rate for these systems. This is different from results reported by Wong et al. for methyl radical attack on alkenes where significant correlation between exothermicity and reactivity were found.50,53,54 In these experiments, strong electronic effects were found near the reactive sites and would be expected to have significant affect on ∆Hrxn. In the current study, lack of strong electronic effects may be inferred. Radical stability also influences the reaction rate as well as activation energy of each homopolymerization reaction. Such an approach was taken by Tedder to understand the relationship between regioselectivity and resonance stabilization of the radicals.55 One method to describe radical stability is the Wiberg bond order between the carbonyl and carbon radical. A greater bond order indicates enhanced stability of the radical due to resonance considerations, resulting in decreased monomer reactivity. For the entire set of methacrylates, no apparent relationship between reaction rates and bond orders of the carbon radical - carbonyl bonds was found for either the isolated radicals or TS structures. Likewise, changes in the bond order from isolated radical to TS showed no apparent trend. A third measure of radical stability may be inferred from the calculated bond length or the radical-carbonyl, which similarly showed no relationship. Another consideration for reactivity of radical reactions is charge transfer, CT. The so-called curve-crossing model56,57 has been used to explain reactivity by a number of researchers for radical reactions.25,50,53,54,58-60 Similar to the work of Wong et al., D+ and D- represent the respective ionization energies and electron affinities of the radical species, whereas A+ and Arepresent the corresponding ionization energies and electron affinities of the monomers. Combinations are calculated from the difference in energies (D+ minus A- or A+ minus D-). As the radical and monomers approach each other to form a dimer, the CT configurations mix into the transition state electronic configuration and overall wave function. Both cases of extreme CT are shown in Figure 5. The relative importance of each configuration may be estimated by comparison of the calculated ionization energies and electron affinities at infinite separation. The lower energy form between D+A- and D-A+ represents the dominant charge transfer species qualitatively. For example, if the D+A- species is substantially lower in energy (more stable), then the radical imparts nucleophilic character at the transition state, whereas in the reverse case the radical exhibits electrophilic character. If a minor difference between the two forms was present, neither would dominate59 and the radical may impart both nucleophilic and electrophilic character.61 The calculated energy for the methacrylates in this study is: Rad- > Rad0 > Rad+ > Mon0 > Mon- > Mon+ where “Rad” is the radical and “Mon” is the monomer. The ionization energy is the energy of the cation less the energy of the ground state of the appropriate species (D+ and A+). Because the electron affinity is the energy of the ground state less that of the anion, only A- is negative based on the relative stability (Table 4). The calculated charge transfer energies are listed in Table 4, with data sorted by the difference between the D+A- and D-A+ values. The expected outcome from this treatment is a net lowering of the activation energy for each reaction would be observed if either form was significantly lower in energy.25,45,50,53,54,58 Thus, the small range of the calculated activation energies for the reac-
J. Phys. Chem. A, Vol. 114, No. 41, 2010 10993
Figure 5. Charge transfer (CT) configurations of the D+A- and D-A+ species of radical (D) attack on a monomer (A).
TABLE 4: Adiabatic Electron Affinities and Ionization Energies of Radical (D) and Monomer (A) Methacrylates in eV GMA CHMA EHMA MMA IBUMA EMA BUMA HEMA HPMA IDMA DMA BOMA BZMA
D+
A-
D-
A+
7.469 7.213 7.418 7.610 7.455 7.505 7.455 7.562 7.270 7.429 7.429 7.282 7.357
-0.526 -0.568 -0.635 -0.763 -0.688 -0.753 -0.739 -0.677 -0.286 -0.729 -0.732 -0.544 -0.507
0.909 0.865 0.785 0.647 0.749 0.699 0.729 0.938 1.139 0.743 0.740 0.896 0.923
9.117 8.817 8.943 9.074 8.912 8.963 8.911 9.135 8.524 8.712 8.630 8.263 8.328
D+A- D-A+ difference 7.995 7.780 8.053 8.373 8.143 8.257 8.194 8.239 7.556 8.159 8.161 7.826 7.864
8.208 7.952 8.158 8.427 8.163 8.264 8.182 8.197 7.385 7.968 7.890 7.367 7.404
-0.214 -0.172 -0.104 -0.054 -0.020 -0.007 0.011 0.042 0.171 0.190 0.271 0.459 0.460
a
CT
0.0030 0.0039 0.0070 0.0044 0.0078 0.0070 0.0066 0.0106 0.0071 0.0069 0.0067 0.0037 0.0055
a NBO charge transfer from the monomer to the radical calculated from the sum of atomic charges of the monomer at the TS.
tions in this study indicate that CT should play only a minor role in the propagation rates. The largest difference between the two forms is 0.46 eV, which is much smaller than the generally accepted value of 1 eV to be interpreted as significant CT for an effect on the activation energy.54,59,60 As a result, this study is consistent with previous interpretation attributing only a minor CT effect on the activation barrier for this series of methacrylates. Studies to identify CT also compare net charge transfer from the radical to the monomer (sum of the atomic charges of the monomer) in the TS structure to support the presence or absence of CT. A similar comparison has been made here using charges from the natural bond orbital (NBO) population analysis method.62 Unlike previous studies, the TS net charge transfer in the current study always flows away from the monomer to the radical (CT is positive), rather than a shift from negative to positive as might be expected from the shift found with the differences between D+A- and D-A+. Further analysis of the NBO atomic charges of each corresponding methacrylate atom shows only minor differences regardless of the methacrylate (
