ARTICLE pubs.acs.org/JPCA
Theoretical Investigation of the Isomerization of trans-HCOH to H2CO: An Example of a Water-Catalyzed Reaction Phillip S. Peters,†,‡ Denis Duflot,† Alexandre Faure,‡ Claudine Kahane,‡ Cecilia Ceccarelli,‡ Laurent Wiesenfeld,*,‡ and Celine Toubin*,† †
Laboratoire de Physique des Lasers, Atomes et Molecules (PhLAM), UMR CNRS 8523, Universite Lille 1, Villeneuve d’Ascq Cedex 59655, France ‡ Universite Joseph Fourier-Grenoble 1/CNRS-INSU, Institut de Planetologie et d’Astrophysique de Grenoble (IPAG) UMR 5274, Grenoble F-38041, France
bS Supporting Information ABSTRACT:
A concerted hydrogen atom transfer mechanism has been elucidated for the isomerization of trans-HCOH to H2CO using a variety of ab initio and density functional theory methods. This work places specific emphasis on the role water molecules can play as a catalyst for this reaction and the mechanism by which this is achieved. This is of particular importance in the context of molecular ices in the interstellar medium because the presence of water in this reaction reduces the activation energy by at least 80%, which is accompanied by a significant enhancement of the reaction rate, at e300 K.
’ INTRODUCTION Formaldehyde is one of the most abundant and ubiquitous molecules observed in the interstellar medium (ISM).1,2 In addition, formaldehyde is believed to be the first step of a process leading to the more complex, organic molecules observed in the ISM and, in particular, in star-forming regions.3,4 Although formaldehyde was detected many years ago and despite its important role in interstellar chemistry, there is still considerable debate about how formaldehyde is formed in the ISM. At present, the formation process that is favored is the hydrogenation of CO on the surfaces of interstellar grains.5 7 The idea is that CO freezes out onto the interstellar dust grains in the dense and cold (e30 K) regions of the ISM and that H atoms landing on the grain surfaces react with CO to form HCO, H2CO, and up to CH3OH in successive H addition steps. Although laboratory experiments8 10 support the validity of such a scheme, there is not a universally accepted theory explaining how this may happen. Having such a theory, on the other hand, is extremely important as it may improve our understanding of the chemistry on the interstellar grains. It is currently believed that the synthesis reactions are catalytically enhanced by many orders of magnitude when they occur on the surface of icy grain mantles. A previous theoretical study of the reaction of H with CO has been conducted in the presence of water by Woon,11 where the water being added in a physisorption-type manner had no effect on the activation r 2011 American Chemical Society
energy. By contrast, recent studies on the HNC HCN isomerization12 and aminoacetonitrile synthesis13,14 provided evidence of the catalytic role that water can play through a concerted exchange mechanism of hydrogen atoms or protons. The aim of our study is to see if the low-temperature ice mantle can significantly alter the energetics and rate of formation of formaldehyde under interstellar conditions. If this should be the case, a detailed analysis of the mechanisms involved will thus be conducted. Two recent theoretical studies have fully characterized the potential energy surfaces (PESs) for the photodissociation of H2CO15 and for the reaction of carbon atoms with water, leading to H2CO formation.16 One mechanism proposed by Schreiner16 begins with the complexation of C atoms to the lone pair of the water oxygen followed by rearrangement to HCOH. From here, two pathways are possible, either dissociation to H2 and CO or internal rearrangement to H2CO, though it should be noted that for this reaction only the singlet pathway, the reaction of singlet carbon atoms with water, may be possible under interstellar conditions because the reaction of triplet carbon atoms with water possesses a barrier that would require a large flux of UV photons or high temperatures to be surmountable. Received: March 3, 2011 Revised: July 8, 2011 Published: July 12, 2011 8983
dx.doi.org/10.1021/jp202052h | J. Phys. Chem. A 2011, 115, 8983–8989
The Journal of Physical Chemistry A
ARTICLE
Scheme 1. Sequential Hydrogenation of CO to H2CO and the Possible Isomers That May Also Be Formeda
a Activation energies for CO + H and the isomerization COH to HCO taken from ref 17. Relative Energy, Re, of trans-HCOH relative to H2CO taken from ref 18. At present, it is not known if the addition of H to HCO to form HCOH possesses an activation barrier for either the singlet or triplet case.
During the hydrogenation of CO (see Scheme 1), addition of hydrogen to both the carbon and oxygen atoms is possible. Normally, only the direct addition to carbon is considered within this sequence. However, during the second hydrogenation, the addition of hydrogen to the oxygen atom of the formyl radical, HCO, would lead to the formation of HCOH. In such a case, there would be competition between trans-HCOH and H2CO formation. Should trans-HCOH be formed, an internal rearrangement would then have to occur to yield H2CO. Depending on the rate of reaction for this internal rearrangement, one may reasonably see one of the following scenarios. First, the time frame for rearrangement is short, so trans-HCOH is readily converted to H2CO. Alternatively, the time frame for the rearrangement is long, and trans-HCOH could then be considered to compete with H2CO in the subsequent hydrogenation reactions. In the matrix isolation experiments by Schreiner et al.,18 HCOH is observed to isomerize in a matter of hours to H2CO. Because the activation energy of this process in the gas phase is clearly too high to be surmounted under the conditions of the experiment, it has been postulated that this process occurs via hydrogen tunneling. In this study, tunneling is considered using the Eckart19 and Wentzel Kramers Brillouin (WKB)20,21 methods. Our goal is to elucidate an alternative route from trans-HCOH to H2CO when water molecules are allowed to actively participate in the isomerization.
’ COMPUTATIONAL METHODS Molecular Structure. The geometries of the reactants, transition states, and products were computed using the aug-cc-pVTZ basis set of Dunning.22,23 Each structure was fully optimized and verified to be a stationary point or a transition state by a normal mode analysis. Intrinsic reaction coordinate (IRC)24 calculations were used to confirm the transition state (TS) geometries. Harmonic zero-point energies (ZPEs) were computed from the relevant frequencies.
As in previous computational studies, we have used the Becke three-parameter hybrid functional and the exchange functional of Lee, Yang, and Parr (B3LYP).25 To validate the accuracy of this functional, Mø ller Plesset perturbation theory at second order (MP2) and coupled cluster [CCSD and CCSD(T)] calculations were also performed. For comparison of the DFT method with ab initio methods, we also used the M06 family of hybrid functionals parametrized by Truhlar.26 This was done because some of the variations of M06 are parametrized to better describe hydrogen bonding interactions compared to B3LYP.26 Note that for the complexes that include water molecules, all geometries were first optimized using the B3LYP/6-31+G** basis set before a second optimization with the aug-cc-pVTZ basis set and the various investigated methods. The impact of basis set superposition error (BSSE) is expected to be small by cancellation because it affects all species equally. It is thus neglected in this study. The MP2 and B3LYP electronic structure calculations were performed with the Gaussian suite of programs,27,28 while for the coupled cluster calculations, the CFOUR29 program package was also used. For comparison, the gas phase transHCOH to H2CO isomerization reaction was also studied using multiconfigurational methods, namely, CASSCF/MRCI (with Davidson’s correction)30 33 and CASSCF/RS2C34,35 levels. For these calculations, the active space consisted of 12 electrons and the 10 valence molecular orbitals (MO’s). However, for HCOH, the final CASSCF resulted in permutation of one valence MO and the oxygen 1s MO, which became active, giving incorrect MRCI and RS2C energies. The phenomenon occurred with the aug-cc-pVTZ basis set and not with the 6-31+G** basis set and was reproduced by both Molpro36 and Gamess.37,38 The only way to solve this problem was to include the 1s orbitals within the active space, which increased the computational cost. The properties of the electronic density were studied in the framework of the Atoms In Molecules (AIM) theory derived by Bader.39 Topological analysis was performed with AIM200040,41 on the CCSD/aug-cc-pVTZ wave functions for each of the stationary points in the gas phase and for complexes including one water molecule. 8984
dx.doi.org/10.1021/jp202052h |J. Phys. Chem. A 2011, 115, 8983–8989
The Journal of Physical Chemistry A
ARTICLE
Table 1. Activation Energies (Ea, in kilocalories per mole) with and without the ZPE Correction for the Cases in Which N = 0 and 1 N=0 method B3LYP M06-HF M06 M06-L M06-2X MP2 CCSD CCSD(T) CASSCF/RS2C
basis 6-31+G**
Ea
Ea + ZPE
Ea
Ea + ZPE 4.94
34.21
37.06
8.14
aug-cc-pVTZ 34.27
30.54
8.93
5.69
6-31+G**
35.33
6.69
3.76
39.24
aug-cc-pVTZ 38.88
34.79
7.83
5.06
6-31+G** 33.15 aug-cc-pVTZ 32.84
29.47 29.11
10.15 10.45
7.05 7.34
6-31+G**
31.87
28.08
9.89
6.58
aug-cc-pVTZ 31.29
27.57
10.01
6.46
6-31+G**
35.64
31.84
9.02
5.90
aug-cc-pVTZ 35.27
31.45
9.84
6.66
6-31+G**
35.49
9.48
6.27
39.34
aug-cc-pVTZ 30.46
26.61
6.84
4.47
6-31+G** 38.51 aug-cc-pVTZ 37.39
34.41 33.37
15.30 13.74
12.16 10.62
6-31+G**
35.80
31.82
13.17
9.87
aug-cc-pVTZ 34.31
30.29
11.10
7.93
6-31+G**
34.37
30.54
aug-cc-pVTZ 30.59
26.31
CASSCF/MRCI(+Q) 6-31+G**
a
N=1
34.89
30.81
aug-cc-pVTZ 33.48
29.42
CCSD(T)a CCSD(T)b
aug-cc-pVTZ 34.75 cc-VTZ
30.5
CCSDT(Q)c
CBS limit
29.7
CCSD(T)//QCISDd
aug-cc-pVTZ
30.3
From ref 15. b From ref 16. c From ref 18. d From ref 42.
Tunneling. Following the calculation of the activation energies, the tunneling probability and rate of reaction were computed using the Eckart19 and Wentzel Kramers Brillouin (WKB)20,21 methods. For details of these models, see the Supporting Information. In both cases, the allowed energies were taken to be those of the harmonic oscillator and the temperature dependence of the rate constant was considered by averaging over these energy levels using a Boltzman distribution. Such a model offers some consideration of the difference between the activation energy of the process and the enthalpy of activation above 0 K. We also note that the Eckart19 potential fits only the three stationary points of the IRC curve, the reactant, TS, and product, and is as such not very accurate, particularly in the reactant region of the potential. However, it has been shown that the Eckart model produces results above 250 K very similar to those of the WKB method, which is based upon the full IRC curve.42 By contrast in the low-temperature regime, the WKB rate is typically 1 order of magnitude higher than the Eckart rate. It should also be noted that formerly the Eckart equations form an exact analytical solution to the Schr€odinger equation for an unsymmetrical potential, whereas the WKB method is a semiclassical approximation using the IRC potential.
’ RESULTS AND DISCUSSION Quantum Chemistry Results. In the study by Zhang et al.,15
the isomerization reaction of hydroxymethylene to formaldehyde
has been considered as part of the photodissociation of formaldehyde producing H2 and CO. The activation energy computed at the CCSD(T)/aug-cc-pVTZ level is reported to be equal to 34.75 kcal mol 1. Another computational study16 considers this isomerization in a series of reactions starting from ground state carbon atoms associated with water. The zero-point corrected activation energy calculated at the CCSD(T)/cc-pVTZ level is calculated to be 30.5 kcal mol 1 in their study. In addition to these studies, the isomerization was also considered at the CCSDT(Q)/CBS + ZPE level by Schreiner et al.,18 where the barrier is reported to be 29.7 kcal mol 1. In our study, we find that the activation energies differ according to the methods employed within the range of 30.46 38.88 kcal mol 1 and between 26.31 and 34.79 kcal mol 1 with the ZPE correction. The B3LYP energies agree very well with the values previously computed at the CCSD(T)/cc-pVTZ level by Schreiner et al.,16 whereas the noncorrected one is in closer agreement with that of Zhang et al.15 With respect to the CCSDT(Q) results,18 the CCSD energies calculated with the augmented basis set overpredict the barrier height by almost 4 kcal mol 1 (Table 1), whereas the MP2 energies underpredict the barrier by almost 3 kcal mol 1. We see that the M06, M06-L, and CASSCF/MRCI(+Q) results are in good agreement with the CCSDT(Q) results,18 suggesting that the M06 functional may perform reasonably well, and given that the close agreement between the CCSDT(Q) results18 and the explicit multireference method suggests that the wave function is described well by a single determinant. We also note that within the M06 family of functionals we see a degree a variation, with the M06-HF functional performing the worst in comparison with the CCSDT(Q) results.18 This is significant because in the calculation of reaction rates, a difference of 1 kcal mol 1 in the activation energy can lead to a difference of up to 1 order of magnitude in the rate. In contrast with the energetics, the geometries do not vary greatly with the method used. While the variation in the angles is much greater than that in the bond lengths, the actual arrangement of the atoms does not change significantly. Additionally in comparing the MP2/aug-cc-pVTZ geometries with the CCSD(T)/ cc-pCVQZ geometries of Schreiner et al.,18 we see only minor variations. For illustrative purposes, the gas phase transition state (TS) is shown in Figure 1a; for comparison, the bond lengths computed by Schreiner et al.18 are shown in parentheses within the figure. As shown for the HNC HCN isomerization13 and for some other reactions that have unfavorably high activation energies in the gas phase, water can play the role of a catalyst. We decided to add one explicit water molecule that can actively participate in the reaction. In the TS structure, shown in Figure 1b, the dissociating hydrogen of HCOH is transferred to the water molecule and one of the hydrogen atoms of the water molecule is subsequently transferred to the carbon atom. Analysis of the normal modes shows that within this TS structure the H atom transfer corresponds to only one imaginary mode, making it consistent with a concerted mechanism. With respect to the different methods, the absolute geometrical parameters show some variations, more especially for the lengths of the two hydrogen bonds present in the TS, whereas the angles and lengths of the covalent bonds have similar values. Consequently, the five-membered ring geometry is conserved regardless of the method chosen. By contrast with the results for the hydrogenation reaction,11 the addition of one water molecule leads to a significant 8985
dx.doi.org/10.1021/jp202052h |J. Phys. Chem. A 2011, 115, 8983–8989
The Journal of Physical Chemistry A
ARTICLE
Figure 1. Geometries of the transition states computed at the MP2/aug-cc-pVTZ level with N water molecules (N = 0, 1, 2, or 3). Angles are given in degrees and lengths in angstroms. Values in parentheses were computed at the CCSDT(Q)/CBS + ZPE level by Schreiner et al.18
reduction in the activation energy, Ea, to values ranging between 6.84 and 13.74 kcal mol 1 (Table 1), by comparison with the gas phase values in the range of 30.46 38.88 kcal mol 1. As already observed for the gas phase, the MP2 noncorrected barrier is calculated to be the lowest, but after the zero-point energy correction, the MP2 energy differs by only 0.48 kcal mol 1 from the B3LYP value. However, the activation energies predicted by both the MP2 and B3LYP methods are approximately half of the CCSD barrier. This suggests that these methods overestimate the stabilization effect caused by the interaction with the water molecule in the transition state. Additionally, we find that the M06 functionals give activation energies that are intermediate with respect to those of MP2 and CCSD(T), with the M06 functional reproducing the CCSD(T) energy to within a few tenths of a kilocalorie per mole. For the gas phase reaction, we note the similarity in values (Table 1) between the MRCI(+Q) results and the results of the single-reference methods. Because these values are so close, one may reasonably conclude that none of the species involved with the reaction possess any significant multireference character; this supports the previous study by Matus et al.,43 who found the singlet triplet gap to be 30 kcal mol 1. We also notice the somewhat eratic comparison of the energies computed with the various methods and the two basis sets; this is typical for these cases because of the nonvariational nature of MP2, DFT, and CCSD(T). While these results clearly do not converge to yield a quantitative value for the activation energy in the presence of one water molecule, they suggest, at least in a qualitative manner, that there is a significant reduction in the activation energy with the participation of one water molecule. From this point, two questions
arise. (i) How does the water molecule lead to such a stabilization? (ii) Does the addition of more water molecules further affect the activation energy? We attempt to give an explanation of the first point using the AIM analysis of the electronic density, and we then discuss the second point. Electronic Density Analysis. Using AIM analysis,40,41 the molecular wave function is analyzed to identify critical points within the electron density. Such analysis provides information about the bonding properties of the complex. The analysis is performed using AIM2000. During the analysis, the CCSD/ aug-cc-pVTZ wave functions for the transition states of the gas phase reaction and the reaction in the presence of a water molecule are employed. Within the Atoms in Molecules (AIM) theory,40 an interaction between two atoms is characterized by a bond critical point along the bond path connecting the nuclear critical points. Additionally, a molecule or molecular complex can be considered to be a ring system when a ring critical point is identified within the density. Bonding interactions may be further classified by type on the basis of the value of the electronic density, F, and the value of the Laplacian L = 32F of the density at the critical point. Typically, a large value of F and a negative value of L are indicative of a covalent or shared electron interaction, whereas a low value of F and a positive value of L correspond to a closed shell interaction such as a hydrogen or ionic bond. In the gas phase transition state without any water molecules present, three bond critical points are evidenced (Figure 2a). For each of them, the values of F and L are typical of covalent bonds. The C O bond is not altered much in the transition state by 8986
dx.doi.org/10.1021/jp202052h |J. Phys. Chem. A 2011, 115, 8983–8989
The Journal of Physical Chemistry A
ARTICLE
Figure 2. Molecular graphs of the transition states when N = 0 (a) and N = 1 (b). The molecular graphs were obtained from the CCSD/aug-cc-pVTZ wave function. Within the molecular graphs, nuclear critical points are colored black, bond critical points are colored red, and the ring critical point is colored blue.
Figure 3. Trend in activation energy, with the zero-point energy correction, as a function of the number, N, of water molecules. For all the methods of computation shown, the aug-cc-pVTZ basis set was used.
Table 2. Activation Energies (Ea, in kilocalories per mole) with and without the ZPE Correction for Cases in Which N = 2 and N = 3 N=2 method
Ea
N=3
Ea + ZPE
Ea
Ea + ZPE
B3LYP
5.08
0.93
5.28
0.84
MP2 CCSD
4.41 9.35
0.36 5.31
4.80 9.87
0.54 5.62
CCSD(T)
7.33
3.28
7.87
3.62
comparison with HCOH and keeps a strong covalent character. Moreover, no interaction is observed between the transferring hydrogen atom and the carbon. This suggests that the transition state is a HCOH derivative with a large degree of angle strain in ^ angle. the COH In comparison with the values of F and L for the TS with one water present, we note that the H C bond of the formyl backbone and the nonparticipating O H bond retain their covalent nature. More interestingly, there is a change in the character of the C O bond between the gas phase TS and the TS H2O complex. This is evidenced by a change in sign of L: in
the gas phase, L is negative, whereas in the TS H2O complex, it is positive. In spite of this change in L, we see that F remains large at the bond critical point. Such behavior suggests that the density is being accumulated outside of the bonding region. This can be rationalized by considering the interaction between the water hydrogen, Hw1, and C; the values of F and L at the critical point are consistent with density being accumulated between the nuclei. These data suggests that the carbon atom is transferring density away from the C O bond to form a covalent bond with the hydrogen atom being donated by the water molecule. We see that the other interactions within the ring show a degree of covalency with the interaction between Hf2 and the water oxygen, Ow, being more typical of a hydrogen bonding interaction.41 This analysis suggests that the stabilization of the transition state is caused by relief of angle strain on going from the gas phase to the one-water molecule complex. This can be rationalized by the change in geometry between the two TS structures, with the latter five-membered ring structure being much more sterically favored. Further stabilization is caused by a significant electron density rearrangement provided by the H2O molecule bringing the TS structure closer to the product complex. The addition of one and two more water molecules is now considered to quantify the effect of further solvation on the barrier. Effect of Further Solvation. When two water molecules are present, the corresponding transition state, where the water molecules may exchange hydrogen atoms, is a seven-membered ring structure, as shown in Figure 1c. The concerted transfer of the hydrogen atoms between the water molecules and transHCOH leads to the formation of formaldehyde complexed with two water molecules. Considering the geometry of the TS 2H2O complex (Figure 1c), all of the hydrogen atoms are shared between the molecules rather than belonging to distinct units, and as observed with a single water molecule present, all of the hydrogen atoms are transferred in a single step. This feature is also observed for the case with three water molecules (Figure 1d), but in this case, the atoms enclose a nine-membered ring system that shows greater flexibility. For the energetics, the overall trend in the presence of water molecules is reported in Figure 3, and values for the corresponding activation energies are listed in Table 2. The major reduction of the activation energy, by a factor of 3.5 7 depending on the method, occurs upon addition of the first molecule. The second 8987
dx.doi.org/10.1021/jp202052h |J. Phys. Chem. A 2011, 115, 8983–8989
The Journal of Physical Chemistry A
Figure 4. Rate of reaction for the isomerization of trans-HCOH to H2CO plotted as a function of temperature and number, N, of water molecules. Eckart rates computed using the optimized CCSD(T)/augcc-pVTZ geometries and WKB rates determined by the interpolation of the CCSD(T)/aug-cc-pVTZ//MP2/aug-cc-pVTZ IRC.
Figure 5. Comparison of the interpolated CCSD(T)/aug-cc-pVTZ// MP2/aug-cc-pVTZ IRC and Eckart potentials in the region of the transition state for the N = 1 case. E0, E4, and E8 are the corresponding energy levels of the harmonic oscillator [En = hν0(n + 1/2)], with n values of 0, 4, and 8, respectively.
water molecule tends to stabilize the TS by 3 5 kcal mol 1 depending on the method, giving a ZPE corrected activation energy lower than 1 kcal mol 1, for the MP2 and B3LYP results, by comparison with the 30 40 kcal mol 1 barrier calculated for the gas phase. This is likely due to a further decrease in the steric strain within the system. The addition of the third water molecule has no significant effect on the barrier, the value of the activation energy varying by only a few tenths of a kilocalorie per mole. This result suggests that the inclusion of two explicit H2O molecules provides optimal conditions for the concerted hydrogen exchange. Addition of extra water molecules can be then considered to mimic the solvation of the complexes in the ice mantle. From a quantum chemistry perspective, even if the methods agree fairly well on the global trend, the absolute values of the computed activation energies vary depending on the method used; in some cases, differences of up to 30% are observed. The B3LYP functional, despite its efficiency, is known to be unsatisfactory in predicting barrier heights and noncovalent interactions, such as hydrogen bonds, that are present in the investigated mechanism. This may in part explain the deviation between the CCSD(T) result and the B3LYP result, though it should be noted that the B3LYP results do not perform with any less
ARTICLE
efficiency than the MP2 results for which no such deficiency in describing hydrogen bonding interactions is known. For the cases in which N > 1, no results are reported using the M06 family of functionals, because when N = 0 and N = 1 we see that they provide a wide range of values for the activation energy, with only the M06 functional performing considerably better than the B3LYP functional in comparison with the CCSD(T) results. Furthermore, there were numerical difficulties in converging the product geometries with these functionals, which could not be resolved by modifying the convergence criteria or the integration grid. We now consider the role water plays in altering the rate of reaction for this system. Tunneling Rates. Using the Eckart and WKB methods,19 21 which have been discussed previously, we have computed the tunneling rate for this isomerization. Figure 4 shows how the rate varies as a function of temperature and with the number of water molecules present in the reaction complex. We see in this case with no water present that the time scale of H atom transfer is on the order of a few thousand hours in the low-temperature regime. In contrast, the WKB method gives a rate between 10 and 100 h. In comparison with the results of Kiselev et al.,42 our rates are 2 orders of magnitude lower. This corresponds to the fact that our activation energy at the CCSD(T)/aug-cc-pVTZ level is 2 kcal mol 1 higher than their CCSD(T)/aug-cc-pVTZ// QCISD/aug-cc-pVDZ barrier (Table 1). We note that the zero-temperature limit of the WKB rate computed by Kiselev et al. corresponds to a half-life of 2 h, which is in perfect agreement with the results of the matrix isolation experiments.18 In contrast, this time scale is significantly reduced to
