Comparison of ReaxFF, DFTB, and DFT for Phenolic Pyrolysis. 2

Oct 4, 2013 - Reaction paths for the loss of CO, H2, and H2O from atomistic models of phenolic resin are determined using the hybrid B3LYP approach...
0 downloads 0 Views 3MB Size
Download PDF
Article pubs.acs.org/JPCA

Comparison of ReaxFF, DFTB, and DFT for Phenolic Pyrolysis. 2. Elementary Reaction Paths Charles W. Bauschlicher, Jr.* Mail Stop 230-3, Entry Systems and Technology Division, NASA Ames Research Center, Moffett Field, California 94035, United States

Tingting Qi and Evan J. Reed Department of Materials Science and Engineering, Stanford University, Stanford, California 94305, United States

Antonin Lenfant Department of Materials Science and Engineering, Stanford University, Stanford, California 94305, United States and Department of Computer Science, Institut Supérieur d’Electronique de Paris, 28 Rue Notre-Dame des Champs, 75006 Paris, France

John W. Lawson Mail Stop 234, Thermal Protection Materials Branch, NASA Ames Research Center, Moffett Field, California 94035, United States

Tapan G. Desai Advanced Cooling Technologies, Inc., Lancaster, Pennsylvania 17601, United States ABSTRACT: Reaction paths for the loss of CO, H2, and H2O from atomistic models of phenolic resin are determined using the hybrid B3LYP approach. B3LYP energetics are confirmed using CCSD(T). The energetics along the B3LYP paths are also evaluated using the PW91 generalized gradient approximation (GGA), the more approximate self-consistent charge density functional tight binding (SCC-DFTB), and the reactive force field (ReaxFF). Compared with the CCSD(T)/cc-pVTZ level for bond and reaction energies and barrier heights, the B3LYP, PW91, DFTB(mio), DFTB(pbc), and ReaxFF have average absolute errors of 3.8, 5.1, 17.4, 13.2, and 19.6 kcal/mol, respectively. The PW91 is only slightly less accurate than the B3LYP approach, while the more approximate approaches yield somewhat larger errors. The SCC-DFTB paths are in better agreement with B3LYP than are those obtained with ReaxFF.

I. INTRODUCTION Phenolic resin has a long history of use in thermal protection systems (TPS). It was an important component of the heat shield material used for Apollo to protect the craft from the multithousand degree temperatures generated during highspeed entry of the Earth’s atmosphere and is used today as part of Phenolic Impregnated Carbon Ablator (PICA), the TPS used for the Mars Science Laboratory. Despite the many successes of the current TPS materials, there are some planetary entries that require materials that are beyond current capabilities. Improving the properties of phenolic resin offers one possible path for improved TPS materials. While much is known about the pyrolysis of phenolic resin, many critical factors are still unknown, which handicaps efforts to improve TPS materials. Molecular dynamics (MD) simulations are one © 2013 American Chemical Society

method that allows an atomic level understanding of the pyrolysis of phenolic resin. Two previous MD studies1,2 were performed using the reactive force field3 (ReaxFF) for the potential function. These studies were computationally intensive and therefore restricted to only the very early stages of pyrolysis and uncharacteristically high temperatures to accelerate the dynamics. The reliability of these MD simulations is limited to the accuracy of the underlying ReaxFF potential. While ReaxFF has been used successfully in many applications, it was not developed for phenolic pyrolysis, and therefore the accuracy of the MD simulations on phenolic Received: August 13, 2013 Revised: October 3, 2013 Published: October 4, 2013 11126

dx.doi.org/10.1021/jp408113w | J. Phys. Chem. A 2013, 117, 11126−11135

The Journal of Physical Chemistry A

Article

Figure 1. Three figures from the MD simulation showing the formation of water by the transfer of an H atom from one OH group to another.

harmonic frequencies were computed at the B3LYP level using analytic second derivatives. The vibrational frequencies were used to characterize the stationary points as minima or transition states. In general, B3LYP has proven to be quite accurate for bond energies and reaction barriers, but there are some systems where the B3LYP approach does not give highly accurate results. Therefore the accuracy of B3LYP results was tested by evaluating the path energetics using the coupled cluster singles and doubles approach,9 including the effect of connected triples determined using perturbation theory,10,11 CCSD(T) at the B3LYP geometries. For the open shell systems the spin unrestricted approach was used. In the CCSD(T) calculations, the carbon and oxygen 1s electrons were not correlated. In addition to the 6-31G** basis sets that are used in the B3LYP calculations, the correlation-consistent polarized valence tripleζ (cc-pVTZ) and quadruple-ζ (cc-pVQZ) sets of Dunning12 were used in the CCSD(T) calculations. The DFT and CCSD(T) calculations were performed using Gaussian0913 or earlier versions. As we show below, the CCSD(T) results confirm the use of B3LYP. We therefore evaluate the reaction path energetics using B3LYP geometries at several additional levels of theory. Functionals using the generalized gradient approximation (GGA) can have larger errors than hybrid methods like B3LYP, especially for reaction barriers. However GGAs are much faster than hybrid methods and become the method of choice for ab initio MD simulations, and this was the approach used in our simulations of phenolic pyrolysis. Therefore we also included the PW91 exchange-correlation functional14 (using the same 6-31G** basis set as the B3LYP) in our calibration study. Note this functional is called PW91PW91 in the Gaussian program. While the DFT based MD is expected to be reliable for the small systems studied in the companion manuscript, the more approximate DFTB approach allows longer times and larger systems to be used in the phenolic modeling and could be used to test ReaxFF in cases where the DFT becomes prohibitively

pyrolysis is not known. Because of the importance of phenolic resin in the TPS used by NASA, we have undertaken calibration studies of ReaxFF for the phenolic resin pyrolysis process.4 In order to make the problem tractable for higher levels of theory, a smaller model system was used in the calibrations studies than was used in the earlier MD studies. In addition to calibration of ReaxFF, it is worthwhile testing the self-consistent-charge DFT-tight binding5 (SCC-DFTB) approach. While generally more expensive and less parallel than ReaxFF, DFTB is much less expensive than DFT. If the DFTB is sufficiently accurate, it could be used to perform simulations to further test the ReaxFF approach for larger systems or for longer simulation time and then possible using DFT. Our calibration is composed of two parts: the first consists of MD studies at high temperatures similar to those reported previously,1,2 and the second considers the reaction paths for the formation of CO, H2, and H2O. These two studies are related, but they focus on different temperature regimes. However, we can use the high temperature MD results to guide our choice of paths for study. For example, an MD path where the H from one OH group is transferred to the OH group on the neighboring ring is shown in Figure 1. The C6 ring with the first OH group is already broken, and cross-linking is visible in this figure. However, this suggests a possible low temperature path where the H from one OH group is transferred to OH on the adjacent ring before any other reactions have occurred. This manuscript reports on the reaction paths, while the MD studies are described in the companion manuscript.

II. METHODS The paths studied were derived from a consideration of chemical intuition and paths observed in the MD simulations. While a significant effort was made to consider all paths, we cannot exclude the possibility that other paths exist. The reaction paths were fully characterized using DFT in conjunction with the 6-31G** basis set of Pople and coworkers.6 The hybrid7 B3LYP8 functional was used in all of the DFT calculations. The geometries were fully optimized, and the 11127

dx.doi.org/10.1021/jp408113w | J. Phys. Chem. A 2013, 117, 11126−11135

The Journal of Physical Chemistry A

Article

Table 1. C6H5OH Bond Energies (kcal/mol) at the DFT, CCSD(T), DFTB, and ReaxFF Levels of Theorya DFT

a

CCSD(T)

DFTB

bond

B3LYP +zpe

B3LYP

PW91

6-31G**

cc-pVTZ

cc-pVQZ

mio

pbc

ReaxFF

O−H C−OH C−H(1) C−H(2) C−H(3) C−H(4) C−H(5)

81.3 107.9 110.8 110.1 111.4 110.1 112.5

89.6 113.4 118.9 118.3 119.6 118.4 120.6

89.1 121.1 116.9 116.4 117.8 116.5 118.7

90.1 112.1 118.7 118.0 119.2 118.0 120.4

95.3 115.9 121.0 120.3 121.6 120.3 122.5

96.4 117.1 121.7 121.1 122.3 121.1 123.2

125.4 148.4 150.8 150.7 151.9 151.1 151.7

124.9 154.0 145.3 145.3 146.5 145.6 146.2

110.3 115.4 156.5 142.7 145.8 143.5 158.2

The C−H bonds are defined in Figure 2. The inclusion of zero-point energies is denoted with “+zpe”.

the tables, do not include zero-point energy. However, the vibrational frequencies are determined at the B3LYP level and confirm that all minima have no imaginary frequencies and the transition states have only one imaginary frequency. Those steps that have been fully characterized computationally are shown with a solid line. A dotted line shows a possible connection between two structures that has not been established by the determination of a transition state. The energetics at all levels of theory are summarized in Table 2 and shown graphically in Figure 6. The species in the figures are labeled as n.m, where n is the path and m is the step along the path. Path 1 is shown at the bottom of Figure 3. The first step is the transfer of the hydrogen from the oxygen to the neighboring carbon atom. This possibility was motivated by the weak O−H bond and by the fact that the simulations that lead to the formation of CO commonly showed no hydrogen on the oxygen before the loss of CO occurred. The product (1.2) is 17.6 kcal/mol above the starting system. The barrier to form (1.2) is sizable as the initial O−H bond is mostly broken before the C−H bond begins to form. However, hydrogen is a light atom, and tunneling is expected to reduce the barrier. The next step is the opening of the ring (1.4); this process is endothermic by 30 kcal/mol with only a small barrier in excess of the endothermicity. The final step is the loss of the CO and the closing of the ring to form a five-membered ring (1.6). This step has a high (73 kcal/mol) barrier because the CCO double bond is almost broken before the C−C bond begins to form. For path 1, the overall agreement between B3LYP and CCSD(T) is very good (see Table 2 and Figure 6). The PW91 results are also in good agreement with CCSD(T) but not quite as good as the B3LYP results. From the variation of the C6H5OH bond energies, the uncertainty in the CCSD(T) values is probably a few kcal/mol. Therefore it is difficult to conclude whether the B3LYP or PW91 results are more reliable, but it is clear that both B3LYP and PW91 supply an even handed description of this path, with the errors in the transition states being similar to those for the equilibrium structures. The DFTB results for both the mio and pbc parameter sets are in reasonable agreement with the DFT and CCSD(T) approaches, supplying a good overall description. ReaxFF gives a good overall reaction energy, but the barrier heights for (1.1) and (1.3) are too large. Path 2 is shown in Figure 4. In this path, the first step is the loss of the hydrogen from the oxygen. As with Path 1, the OH bond energy and MD simulations suggest this as a possible first step. Our MD simulations show that the loss of the H can be enabled by reactions with H or OH free radicals that have been

expensive. Thus we also calibrate the DFTB approach in this study. The DFTB calculations are performed utilizing the selfconsistent charge density functional tight binding method (SCC-DFTB)5 with the DFTB+ package,15 employing both “mio-1-1” and “pbc-0-3” parameter sets. The ReaxFF potential used the Chenoweth parameters3 including the C2 and CO correction terms. LAMMPS16 with the c implementation of ReaxFF was used evaluate the ReaxFF energies. The C2 value used is 70.

III. RESULTS AND DISCUSSION The energies required to break the C−OH, C−H, and O−H bonds in C6H5OH are given in Table 1. They do not include any zero-point energies (zpe), unless specifically noted in the table. We study the bonds in phenol because it is the fundamental unit of repetition in phenolic polymers. The C−H bonds are defined in Figure 2. Our most accurate value should

Figure 2. The structure of C6H5OH. The C−H bonds are numbered to aid in the discussion of bond energies.

be the CCSD(T) in the cc-pVQZ basis set. The small difference between the cc-pVTZ and cc-pVQZ shows that the large basis set calculations are close to the basis set limit. An inspection of the table shows that the B3LYP and CCSD(T) values are in good agreement. A significant result of these calculations is that the O−H bond is significantly weaker than the others, which all have relatively similar bond energies. The DFTB and ReaxFF values are larger than the CCSD(T) values but tend to follow the CCSD(T) and B3LYP trends. One exception is the ReaxFF value for the C−OH bond, which is more similar to the O−H bond than the other methods. Four paths for the formation of CO are shown in Figures 3−5. These paths are based on the MD studies and chemical intuition. The energies plotted in the figures, like those given in 11128

dx.doi.org/10.1021/jp408113w | J. Phys. Chem. A 2013, 117, 11126−11135

The Journal of Physical Chemistry A

Article

Figure 3. Paths 1 and 3 for CO formation at the B3LYP level. Solid lines indicate paths where the transition states have been computed, while dashed lines indicate paths where the transition states have not been computed. “T” signifies a triplet state.

Table 2. Energetics (kcal/mol) for CO Forming Reaction at Several Levels of Theorya DFT

a

CCSD(T)

reaction path step label

B3LYP +zpe

B3LYP

PW91

6-31G**

1.1 1.2 1.3 1.4 1.5 1.6 2.1 2.2 2.3 2.4 2.5 2.6 2.7 3.1 3.2 3.3 3.4 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10

65.1 16.8 51.6 44.7 114.2 32.2 81.3 132.8 129.7 137.7 120.1 124.4 110.7 143.0 111.6 116.8 101.4 −25.9 26.2 7.4 18.4 −5.7 26.6 −15.3 15.6 10.8 32.2

69.3 17.6 54.6 47.6 120.5 36.6 89.6 143.5 139.7 148.2 129.8 135.5 124.4 148.7 116.8 123.7 109.5 −31.4 23.8 1.9 14.9 −9.4 24.4 −19.7 14.2 10.2 36.6

61.0 17.4 50.3 47.9 124.1 43.0 89.1 135.7 133.1 139.5 127.8 133.7 130.3 149.3 119.1 125.8 119.3 −30.9 20.3 2.6 13.5 −6.3 24.1 −18.1 13.7 15.0 43.0

72.7 18.4 58.0 48.3 118.7 28.1 95.3 145.8 140.4 150.3 126.2 132.9 118.1 152.1 119.7 128.0 108.1 −26.7 30.2 3.9 20.4 −2.5 27.3 −20.4 16.7 8.7 24.3

DFTB

cc-pVTZ

mio

pbc

ReaxFF

−28.0 29.6 7.0 21.4 −0.8 28.9 −18.0 18.0 10.1 28.1

72.0 21.6 50.3 50.1 105.3 25.9 125.4 183.9 177.2 190.4 158.8 153.1 144.2 160.5 138.0 136.4 124.3 −49.3 12.9 −5.9 2.9 −13.6 14.3 −30.8 −5.1 0.2 26.0

80.6 26.6 58.4 58.3 115.9 34.2 124.9 176.7 174.7 184.7 158.7 149.4 147.1 162.4 145.7 146.7 134.9 −44.1 26.3 4.3 9.6 −0.1 27.6 −19.9 8.4 13.8 34.2

151.8 55.0 102.0 98.7 131.6 45.6 110.3 166.9 175.8 171.9 151.8 145.0 99.0 159.6 143.3 150.6 109.3 −45.9 17.9 14.3 19.4 9.6 41.5 1.0 26.8 23.5 45.5

The inclusion of zero-point energies is denoted with “+zpe”.

lost from the same molecule or another molecule. We believe that this step should not have any barrier; however, we find a small (about 4 kcal/mol) barrier as we stretch the O−H bond from equilibrium to infinite separation. We attribute this small

barrier to limitations of the DFT method and instead of including the barrier we connect the first two structures with a dotted line. This step is more endothermic than the first step in Path 1 (the H transfer), but, as we show below, this is smaller 11129

dx.doi.org/10.1021/jp408113w | J. Phys. Chem. A 2013, 117, 11126−11135

The Journal of Physical Chemistry A

Article

Figure 4. Path 2 for CO formation at the B3LYP level. Solid lines indicate species that have been computed, while dashed lines are suggested paths.

Figure 5. Path 4 for CO formation at the B3LYP level. Solid lines indicate species that have been computed, while dashed lines are suggested paths.

than the energy to open the ring as the first step. The next series of steps shifts a C−C bond, changing the initial sixmembered ring into a five-membered ring. The CO loss has only a very small barrier leading to final products. The highest barrier for this path is only slightly higher than that of Path 1, and actually lower if the initial abstraction occurs with a free radical that was created earlier in the process. As with Path 1, the B3LYP and CCSD(T) energetics are in good agreement. The agreement between PW91 and B3LYP is reasonable but not as good as found for path 1. For this path, ReaxFF and DFTB results are clustered together at higher energies, while the DFT and CCSD(T) results are clustered together at lower energies. For this path we also investigated an open structure, where the six-membered ring opens instead of converting to a five-

membered ring. With one fewer C−C bond, the open product (2.8) is less stable than the ring C5H5CO (2.5) product. In addition to the higher energy product, the barrier for this ringopening is rather high as π bonding is lost as the transition state is nonplanar. The open product can lose a CO and form either the same five member ring product (2.7) or form an open structure (2.9). Since the former product would involve a rotation about a C−C bond, we expect a high barrier as found when the ring opened. The linear product is quite high in energy. The open product path appears to be less favorable than the path that retains a ring structure throughout the process. Path 3 is the upper path in Figure 3. In this path, the first step is the opening of the ring. We were unable to find a singlet structure that did not form a ring, either returning to the 11130

dx.doi.org/10.1021/jp408113w | J. Phys. Chem. A 2013, 117, 11126−11135

The Journal of Physical Chemistry A

Article

Figure 6. Summary of the paths as a function of level of theory. The B3LYP geometries have been used for all levels of theory.

starting phenol or to some other ring, such as the one shown in Figure 3. We switched to the triplet surface to force an open structure, and as shown in Figure 3, this species is about 140 kcal/mol above the starting material and is higher than the highest point on path 1. We had a great deal of difficulty

determining subsequent reactions for this species. After breaking the O−H bond, we consistently found cyclic C6H5O (2.1), which is much easier to form by simply breaking the O− H bond than a process passing through a ring-opening step. We speculate that a H transfer reaction (3.1 to 3.2) would be a 11131

dx.doi.org/10.1021/jp408113w | J. Phys. Chem. A 2013, 117, 11126−11135

The Journal of Physical Chemistry A

Article

possible step, but we were unable to find the transition state that connects 3.1 to 3.2. If this step occurs, it is easy to lose the CO. While we did not find some of the transition states for the triplet surface, it seems safe to assume that the barriers will be even higher 3.1. Therefore, our overall impression is that path 3 is less favorable than either 1 or 2. The DFTB and ReaxFF results are somewhat larger than the B3LYP, but the agreement is better than that found for paths 1 and 2. Path 4 is shown in Figure 5 and appears to be the most favorable of the four CO forming reactions considered. This reaction was the most commonly observed in our MD studies. It starts by adding an H atom to C6H5OH, which is an exothermic process (4.1). This is followed by a hydrogen shift (4.3) and a ring-opening (4.5). The open structure closes to form a new five-membered ring (4.7). This is followed by the loss of HCO (4.9), which subsequently dissociates into H+CO (4.10). Since this reaction appears to be the most favorable, CCSD(T) calculations are performed using both the 6-31G** and cc-pVTZ basis sets. The CCSD(T) results differ by only a few kcal/mol, suggesting that these value are reliable. B3LYP and PW91 are in good agreement with the CCSD(T) results. The differences between the CCSD(T)/cc-pVTZ and the B3LYP results are smaller than 10 kcal/mol for all species. Excluding the final product where the difference is about 15 kcal/mol, the PW91 results are also within 10 kcal/mol of the CCSD(T)/cc-pVTZ values. As with the other reactions, the DFTB and ReaxFF are in reasonable agreement with the CCSD(T)/cc-pVTZ and B3LYP results. The differences exceed 20 kcal/mol for some species, while for others, the difference with CCSD(T) is only a few kcal/mol. While path 4 is the most favorable for the formation CO, we should note that (4.5) can undergo the loss of a hydrogen atom to form 1,3,5-hexatrien-1-one or 2,4,5-hexatrienal. These reactions are endothermic by 46.4 and 57.1 kcal/mol, respectively, at the B3LYP level. While these reaction energies are larger than the barrier for the conversion of (4.5) to the five-membered ring (4.7) at the same level of theory (33.8 kcal/ mol), it is expected that some loss of H will occur making path 4 somewhat less efficient than it appears. The next reaction we consider is the formation of H2. The path for this reaction is much simpler than those found for CO production and the reactants, transition state, and product are shown in Figure 7. We consider two models for this process: a single ring C6H4OHCH3 and a (C6H4OH)2CH2 model with two rings and a bridging CH2 group. We denote these as the CH3 and nonplanar models, respectively. The smaller model allows the use of a larger basis set in the CCSD(T) calculations. The energetics for this reaction are given in Table 3 and shown graphically in Figure 6. We first note that the two molecules yield similar results; that is replacing the second C6H5OH with a hydrogen atom does not have a significant effect. The agreement between CCSD(T) in the 6-31G** and cc-pVTZ basis sets is good as is the agreement between the B3LYP and CCSD(T) results. The PW91 results are in good agreement with CCSD(T) and B3LYP results. The DFTB transition states and products are too high in energy compared with the higher levels of theory. The ReaxFF transition state energy is even higher than that found using DFTB, while the ReaxFF treatment of the products is consistent with DFTB. The last product paths that we consider lead to the formation of H2O. The reactants, transition state, and products for an interchain reaction are shown in Figure 8, while the intrachain reactions are shown in Figures 8 and 9. The energetics of all

Figure 7. The reactants, transition state, and products for the formation of H2 using both the CH3 (top) and nonplanar (bottom) models. The geometries were obtained at the B3LYP level of theory.

H2O formation reactions are given in Table 4. For the intrachain reaction, we consider the same two models as used for the H2 formation reactions. For the nonplanar case, there are two paths, one involving a transfer of a hydrogen atom from the CH2 bridging group to the OH. This is the analogous hydrogen in the interchain and CH3 models. The second path is different, involving the transfer of the H from one OH group, not the CH2 bridging group, to the OH group on the neighboring ring. As with the H2 formation, the smaller system allows for the use of larger basis sets at the CCSD(T) level. Since H2O formation is believed to be one of the most important reactions in the pyrolysis process, we perform CCSD(T) calculations using the cc-pVQZ basis for the small model and the cc-pVTZ basis set for the larger model. The interchain reaction at the top of Figure 8 is probably our most qualitative path as the weak interaction between the two molecules results in many local minima, and we found several reactions, where the H from the OH on one molecule was transferred to a ring carbon on the other molecule. This is the only reaction we found where the OH from one molecule reacts with a hydrogen atom on the CH2 bridge of the other molecule to form H2O. As shown in Table 4 the reaction is nearly thermoneutral, being slightly endothermic at the DFT levels and slightly exothermic at the DFTB and ReaxFF levels. All methods have a sizable barrier. Overall, the methods are in reasonable agreement for this reaction. The intrachain reaction using the C6H4OHCH3 (denoted CH3) model is shown in the middle of Figure 8. It has a barrier that is a bit higher than the interchain model because the OH cannot approach the H as closely due to geometry constraints of being on the same molecule. In this case, the CCSD(T) results in the largest basis set support the DFT results. The reaction is very endothermic because no C−C bond forms as in the interchain case. Note that the product is above the transition state because we considered the product with the water completely removed instead of H2O complexed with the C6H4CH2 radical. Overall, the methods are in reasonable agreement for both the reaction energy and barrier. The CCSD(T) results suggest that both DFTB and ReaxFF underestimate the barrier. The intrachain reaction using the larger nonplanar model is show in Figure 9. This starts with the reaction of the OH and 11132

dx.doi.org/10.1021/jp408113w | J. Phys. Chem. A 2013, 117, 11126−11135

The Journal of Physical Chemistry A

Article

Table 3. Comparison of the H2 Reaction Path Energetics (kcal/mol) Computed Using DFT, CCSD(T), DFTB, and ReaxFFa DFT reaction path CH3 Model Barrier(TS) Products Nonplanar Model Barrier(TS) Products a

CCSD(T)

DFTB

B3LYP +zpe

B3LYP

PW91

6-31G**

cc-pVTZ

mio

pbc

ReaxFF

63.7 38.7

68.9 47.1

62.4 49.2

75.8 45.9

74.8 49.5

93.5 87.1

89.0 81.4

108.2 83.7

56.0 32.4

61.4 40.6

54.1 41.5

70.1 42.9

83.4 79.5

79.1 73.8

103.2 76.7

The B3LYP structures are shown in Figure 7. The inclusion of zero-point energies is denoted with “+zpe”.

H, as in the CH3 model. This leads to a water complex that is below the transition state. The loss of the water yields a product that is actually above the transition state. A rotation about a C− C bridging bond allows the radical site where the OH was lost to react with the other ring to form a five-membered ring. The strain associated with the five-membered ring results in a reaction which is more endothermic than the interchain case. As with the other H2O formation reactions, there is reasonable agreement between the methods; however, ReaxFF seems to underestimate the first barrier compared to the other methods. The formation of the third ring for the intrachain H2O formation reaction and the formation of the five-membered ring for the CO formation by path 2 are very interesting, because the MD studies show that during the pyrolysis of phenolic additional rings, both five- and six-membered form. The intrachain H2O formation reaction suggests that some ring formation could be occurring very early in the pyrolysis process, where as interchain reactions are expected to create bridges between the rings but not fused rings early in the pyrolysis and only have ring formation later in the pyrolysis process. The water formation involving two nearby OH groups is shown in Figure 8. This reaction leads to a six-membered ring containing an oxygen. The energetics are summarized in Table 4. This path is nearly thermoneutral like the interchain path. It has the lowest barrier of the H2O formation routes considered.

Figure 8. The reactants, transition state, and products for the formation of H2O using the interchain (top), CH3 (middle), and the nonplanar OH to OH transfer (bottom) models. The B3LYP path information given.

Figure 9. The path for the formation of H2O using the nonplanar model. The geometries were obtained at the B3LYP level of theory. 11133

dx.doi.org/10.1021/jp408113w | J. Phys. Chem. A 2013, 117, 11126−11135

The Journal of Physical Chemistry A

Article

Table 4. H2O Reaction Path Energetics (kcal/mol) at the DFT, CCSD(T), DFTB, and ReaxFF Levels of Theorya DFT B3LYP +zpe

reaction path step label

CCSD(T)

B3LYP

PW91

Interchain Model Barrier(TS) 85.4 87.1 Product 5.6 7.4 CH3 Model Barrier(TS) 95.4 97.8 Product 99.8 105.4 Nonplanar Model-Bridge H to OH Transfer 5.1 86.7 88.5 5.2 76.1 78.8 5.3 89.5 94.8 5.4 92.3 98.1 5.5 28.5 32.1 Nonplanar Model−OH to OH Transfer Barrier(TS) 64.6 67.9 Product 2.8 6.1 a

6-31G**

cc-pVTZ

DFTB cc-pVQZ

75.2 7.8 93.4 106.7

104.0 101.8

95.1 99.1

83.3 77.9 94.1 96.7 31.9

96.5 84.5 94.2 96.1 27.6

89.0 82.6 92.4 94.7 26.0

58.7 8.7

70.9 5.8

94.6 98.4

mio

pbc

ReaxFF

77.3 −3.1

75.4 −5.6

60.1 −13.2

82.2 98.4

78.3 98.3

53.8 112.4

75.0 80.2 85.2 94.8 25.5

71.3 80.4 85.1 94.7 22.9

35.6 106.1 100.6 103.6 43.1

64.9 1.4

65.3 1.3

85.8 44.7

The inclusion of zero-point energies is denoted with “+zpe”.

The B3LYP, PW91, and DFTB are in good agreement with each other and with the CCSD(T). The ReaxFF barrier is somewhat larger than the other methods and the reaction energy is much too endothermic. It appears the ReaxFF does not describe the C−O−C bridge as well as the other methods. It is possible that adding additional molecules to the ReaxFF training set could improve its description of this reaction. An overall assessment of the different methods is given in Table 5 and Figure 6. The highest level of theory is CCSD(T) using the cc-pVQZ basis set. Unfortunately, we have only nine values at this level, and none of them are for transition states.

The B3LYP and PW91 results are in good agreement with CCSD(T), while the more approximate methods show larger differences. There are more CCSD(T) results in the cc-pVTZ set including seven transition states. The agreement between the DFT, DFTB, and ReaxFF approaches and CCSD(T) in the cc-pVTZ basis set are similar to those obtained in the larger basis set. Interestingly, the differences with CCSD(T) for transition states versus bond and reaction energies do not show a clear pattern. Some of the errors are actually smaller for the transition state calculations. For the 6-31G** basis set calculations, the differences increase compared with the ccpVTZ and cc-pVQZ basis sets. Part of this probably arises from the fact that CCSD(T) is less reliable using the 6-31G** basis set and part arises from the different paths included in the 631G** study. Overall, B3LYP and PW91 show similar differences with CCSD(T). This supports the use of a nonhybrid functional, which is computationally less demanding. The two DFTB approaches show similar accuracy, while the errors in ReaxFF are slightly larger. It is interesting to note that the ReaxFF has a larger error for a few steps in the reactions studied. It is possible that the addition of some extra cases to the ReaxFF training set could give the ReaxFF a more even handed description of pyrolysis of phenolic resin.

Table 5. Average Absolute and Maximum Errors Compare with the CCSD(T) Results As a Function of Basis Set Used in the CCSD(T) Calculation basis set 6-31G** method no. of cases B3LYP PW91 DFTB(mio) DFTB(pbc) ReaxFF no. of cases B3LYP PW91 DFTB(mio) DFTB(pbc) ReaxFF no. of cases B3LYP PW91 DFTB(mio) DFTB(pbc) ReaxFF

avg

max

cc-pVTZ avg

max

All Values 47 26 3.3 12.3 3.8 8.6 5.8 18.7 5.1 14.9 18.2 41.2 17.4 37.6 16.3 41.9 13.2 38.1 23.9 79.1 19.6 53.4 Transition States 16 7 4.1 8.7 4.3 7.9 16.0 6.6 12.4 16.1 40.1 15.1 23.1 11.8 34.4 8.3 17.7 25.0 79.1 18.7 53.4 Bond and Reaction Energies 31 19 2.9 12.3 3.6 8.6 4.7 18.7 4.5 14.9 19.2 41.2 18.2 37.6 18.6 41.9 15.0 38.1 23.3 50.4 20.0 41.3

cc-pVQZ avg

max

9 3.8 4.9 24.4 22.4 23.1

7.0 8.3 31.3 36.9 40.8

IV. CONCLUSIONS We have considered the bond energies of C6H5OH and paths for the loss of CO. The loss of H2 and loss of H2O are studied using several models, namely C6H4OHCH3, which has a single ring, (C6H4OH)2CH2 with two rings and a CH2 bridge and a dimer of (C6H4OH)2CH2 molecules. CCSD(T) calculations in the 6-31G** and the cc-pVTZ basis sets agree reasonably well, except for the transition state for the formation of H2O, where there is a difference of about 10 kcal/mol. We find the average absolute differences for the B3LYP, PW91, DFTB(mio), DFTB(pbc), and ReaxFF with the CCSD(T)/6-31G** approach to be 3.3, 5.8, 18.2, 16.3, and 23.9 kcal/mol, respectively. The analogous values for the CCSD(T)/ccpVTZ approach are 3.8, 5.1, 17.4, 13.2, and 19.6 kcal/mol and for the CCSD(T)/cc-pVQZ approach 3.8, 4.9, 24.4, 22.4, and 23.1 kcal/mol. The B3LYP results are in the best agreement with the CCSD(T) calculations, which supports

0

9 4.9 24.4 22.4 23.1

31.3 36.9 40.8 11134

dx.doi.org/10.1021/jp408113w | J. Phys. Chem. A 2013, 117, 11126−11135

The Journal of Physical Chemistry A

Article

the use of the B3LYP geometries for the other levels of theory. The PW91 results differ only slightly more from the CCSD(T) than do the B3LYP results; this supports the use of the PW91 functional in ab initio MD simulations. The results obtained using DFTB are quite reasonable, although its differences with the CCSD(T) are somewhat larger than those obtained from the B3LYP and PW91 methods. From the difference with the CCSD(T), the pbc parametrization appears to be slightly better than the mio, but the differences are small. ReaxFF differs with the CCSD(T) by even more than the DFTB but incurs only a fraction of the computational cost of the DFTB. High temperature simulations using the DFTB and ReaxFF approaches are expected to give insight into phenolic pyrolysis. However, due to the exponential term in the Arrhenius equation, the errors in the energetics will reduce the reliability of the DFTB and ReaxFF rate constants considerably for lower temperatures. It is possible that reparameterization of ReaxFF using data generated in this work could improve the accuracy of ReaxFF for the study of phenolic pyrolysis.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS This work was partly supported by the Ames Center Innovation Fund. REFERENCES

(1) Jiang, D.; van Duin, A. C. T.; Goddard, W. A.; Dai, S. J. Phys. Chem. A 2009, 113, 6891. (2) Desai, T. G.; Lawson, J. W.; Keblinski, P. Polymer 2011, 52, 577. (3) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A., III J. Phys. Chem. A 2008, 112, 1040. (4) Qi, T.; Bauschlicher, C. W.; Lawson, J. W.; Desai, T. G.; Reed, E. J. J. Phys. Chem. A 2013, DOI: 10.1021/jp4081096. (5) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Phys. Rev. B 1998, 58, 7260. (6) Frisch, M. J.; Pople, J. A.; Binkley, J. S. J. Chem. Phys. 1984, 80, 3265 and references therein. (7) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (8) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623. (9) Bartlett, R. J. Annu. Rev. Phys. Chem. 1981, 32, 359. (10) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. Chem. Phys. Lett. 1989, 157, 479. (11) Watts, J. D.; Gauss, J.; Bartlett, R. J. J. Chem. Phys. 1993, 98, 8718. (12) Dunning, T. H. J. Chem. Phys. 1989, 90, 1007. (13) Frisch, M. J. et al. Gaussian 03, Revision D.02; Gaussian, Inc.: Pittsburgh, PA, 2003. (14) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. Rev. B 1992, 46, 6671− 87. (15) Aradi, B.; Hourahine, B.; Frauenheim, Th. J. Phys. Chem. 2007, 111, 5678. (16) Plimpton, S. J. Comput. Phys. 1995, 117, 1. Also see: http:// lammps.sandia.gov (accessed Oct 8, 2013).

11135

dx.doi.org/10.1021/jp408113w | J. Phys. Chem. A 2013, 117, 11126−11135