Quantum Chemical Study of the Dissociative Adsorption of OH and

Dec 14, 2006 - Horizontal and vertical bulk effects have been considered by variation of the model system size and number of graphene layers. Potentia...
11 downloads 13 Views 1MB Size
Download PDF
J. Phys. Chem. C 2007, 111, 1355-1365

1355

Quantum Chemical Study of the Dissociative Adsorption of OH and H2O on Pristine and Defective Graphite (0001) Surfaces: Reaction Mechanisms and Kinetics S. C. Xu, S. Irle,* D. G. Musaev,* and M. C. Lin* Cherry L. Emerson Center for Scientific Computation and Department of Chemistry, Emory UniVersity, Atlanta, Georgia 30322 ReceiVed: September 19, 2006; In Final Form: October 20, 2006

Quantum chemical studies of the dissociative adsorption mechanisms and kinetics of OHx (x ) 1, 2) on graphite (0001) surface models were carried out with emphasis on the influence of surface defects. Finitesize model systems for Stone-Wales type (SW), mono- (1V), and divacancy (2V) defects are studied in addition to the defect-free (0D) surface models. Horizontal and vertical bulk effects have been considered by variation of the model system size and number of graphene layers. Potential energy surfaces for dissocative adsorption of one and two OH radicals on pristine graphite and for H2O adsorption on 1V, 2V, and SW defects are presented at the dispersion-augmented density functional tight binding (DFTB-D), integrated ONIOM(B3LYP/6-31+G(d):DFTB-D), and density functional theory B3LYP/6-31+G(d) levels of theory. OH is found to oxidatively erode the graphite surface by producing gaseous CO and hydrogenated vacancy defects (L1H1V), while water is found to react only with monovacancy defects. The rate constants for the dissociative adsorption have been predicted by RRKM theory. Quantum chemical molecular dynamical (QM/ MD) simulations at high temperatures of adsorption products of OH and H2O on graphite surface models have been carried out using the DFTB-D method.

1. Introduction Graphite is an important surface lining material for systems operating under high temperature and high pressure (high-T/ high-P) and has become popular as surface material for rocket nozzles1 as well as plasma facing components.2 Despite this popularity, very little is known about the high-T/high-P processes causing graphite erosion due to reactions with oxidizing agents from fuel combustion, most importantly H2O and CO2.1 Only very recently, numerical simulations have given insight into the dynamics of reactive turbulent combustions3 and have stressed the importance of OxHy (x ) 1, 2; y ) 0-2) species. While some ideas are currently being developed regarding the graphite erosion process in nuclear fusion reactions,4 still nearly nothing is known about the oxidative processes related to fuel combustion, causing rocket nozzle material to fail after prolonged exposure, and how such erosion could be mitigated by improving the surface material. It is obvious that with most modern fuels, OHx, COx, and NOx are important exhaust species that may induce the graphite oxidation. However, if and how graphite defect formation caused by the interaction of these species with graphite lining material plays a role at the high-T/high-P conditions is still unclear. We have recently developed a methodology that allows us to investigate a priori high-T/high-P dissociative adsorption processes on graphite surfaces.5 This methodology is based on a description of the surface graphite (0001) layer using mono-, bi-, and trilayers of dicircumcoronene C96H24 as a model of graphite, employing density functional tight binding6,7 augmented with London dispersion (DFTB-D)8,9 quantum chemical molecular dynamics (QM/MD) of canonical ensembles, ONIOM(DFT:DFTB-D) calculations for the characterization * Corresponding authors. E-Mail: (S.I.) [email protected]; (D.G.M.) [email protected]; (M.C.L.) [email protected].

of intermediate surface species and irreversible adsorption products, and reaction rate constant predictions based on RRKM.5 Applying this methodology, we predicted high-T/ high-P dissociative adsorption processes and their reaction kinetics for H2O,5 COx,10,11 and NOx11 (x ) 1, 2) species. In the investigation of the H2O-graphite system, we found that irreversible oxidation of a pristine graphite (0001) surface by the H2O f H + OH reaction is virtually impossible due to associated high O-H bond activation barriers and very small reverse barriers, which stems from rapid self-healing of an intact graphene surface, leading to immediate ejection of any chemisorbed adsorbate.5 In the investigation of dissociative COx and NOx adsorption on pristine graphite (0001),11 we found CO210 as well as the radical species NO and NO2 causing irreversible oxide defects on the graphite surface11 with CO leaving as a thermodynamically highly stable, volatile species. We also noticed the incorporation of N and subsequent C/N replacement that occurs in the case of NO.11 In this work, we extend our efforts on the studies of the dissociative addition mechanisms and kinetics on graphite (0001) to include a Stone-Wales type 5/7-defect (SW), mono- (1V) and divacancy (2V) graphite defects in addition to the pristine (defect-free, 0D) graphite surface. While the Stone-Wales defect has been known for a long time12 and its energetics are studied in great detail,13-15 vacancy defect formations occur frequently in HR-TEM observations of graphite and carbon nanotubes as a consequence of the high-energy electron beam,16,17 and their structures, energetics,13,18-27 and surface diffusion dynamics have been studied both experimentally28 and theoretically29 in detail before using quantum chemical potentials. HR-TEM experiments have provided direct evidence for the existence of 5/9 and 5/8/5 1V and 2V defects and their migrations,16 and even larger defects such as tri- and quadrovacancies have been recently studied theoretically.27 The 1V

10.1021/jp066142i CCC: $37.00 © 2007 American Chemical Society Published on Web 12/14/2006

1356 J. Phys. Chem. C, Vol. 111, No. 3, 2007 contains a divalent, very reactive carbon reaction center, and we have included a hydrogenated radical defect, 1H1V, with a single H-atom attached to the carbene center, saturating the σ-bond network but creating an unpaired π-electron in our simulations. The same system had been investigated by Lu et al. before using a tight-binding approach in the case of 1Vdefective nanotubes.22 We assume that these vacancy defects are likely to be introduced in graphite lining material under the turbulent conditions of propellant combustion, and we study their reactions with OHx adsorbates with x ) 1, 2. Until now, only one experimental study for the reaction of OH radicals with graphite at 298 K has been reported30 in which the authors have shown that the OH radicals (a) react rapidly with graphite at 298 K and produce approximately equal amounts of CO and CO2, and (b) are much more reactive with graphite than free oxygen atoms under comparable conditions. Previously, several theoretical studies for the interactions of O31-34 and H34-38 atoms on pristine graphite surface have been reported. For instance, Jelea et al.34 used PBE DFT and DFTbased QM/MD to study the addition energetics and surface diffusion of atomic H and O on a one-layer model of graphite. In these studies, they have shown that O atoms can bind to the graphite surface to form an epoxide-like structure, and that hydrogen addition to the adsorbed oxygen centers occurs readily, resulting in the formation of surface-attached hydroxyl (OH) radicals that are bound by only 12 kcal/mol to the otherwise undisturbed graphite surface.34 Further H addition to these surface OH groups is reportedly barrierless and results in loss of water, while the graphite sheet reassumes a honeycomb lattice structure.34 Earlier, Incze et al. had investigated the interaction of atomic O with graphite surfaces, in particular, on graphene edges.33 A work similar in scope was carried out by Sorescu et al.;32 in all cases, the epoxide structure resulting from a [2 + 1] cycloaddition appeared to be the most prominent and stable on-surface oxide-defect species. To the best of our knowledge, no attempt has been made to investigate high-T/high-P dissociative adsorption reactions of OHx species on pristine 0D as well as SW, 1V, 1H1V, and 2V vacancy-defective graphite surface models, and we will present the results of our investigation below in Sections 3-5 with separate sections for potential energy surface (PES) explorative studies, kinetic reaction rate constant calculations, and QM/MD simulations of the reaction products, respectively. 2. Computational Methodology The exploration of the OHx (x ) 1, 2) reactions with pristine and defective graphite PESs was performed using the hybrid density functional B3LYP/6-31+G(d), the integrated ONIOM(B3LYP/6-31+G(d):DFTB-D),39-41 and self-consistent charge density functional tight binding6,7,42 plus London dispersion8 (DFTB-D) levels of theory on finite-size molecular systems, which are shown in Figure 1 and will be described in detail in Section 2.2. 2.1. Methods. As in our previous investigations, we made use of the ONIOM approach to include the horizontal bulk effect (mainly π-conjugation) of graphite, while the monolevel, costeffective DFTB-D method was employed in the description of a bilayer graphene system to qualitatively include the graphite vertical bulk effect also. Whenever radical species were encountered, a spin-unrestricted DFT calculation was employed, while DFTB calculations were left spin-unpolarized because no exchange interaction term exists in the selected DFTB formalism. This means that the highest occupied molecular orbital (HOMO) is treated as a singly occupied molecular orbital

Xu et al. (SOMO) in DFTB. Depending on the model system, we have optimized equilibrium geometries and transition states (TSs) using either monolevel B3LYP/6-31+G(d), ONIOM(B3LYP/ 6-31+G(d):DFTB-D), and monolevel DFTB-D approaches and employed the latter approach to reoptimize structures and calculate numerically vibrational frequencies and zero-point energy (ZPE) corrections. The quantum chemical methods (in particular, the combination of B3LYP/6-31+G(d) and DFTB-D levels of theory) and most of the model systems employed (presented below) have been successfully tested by us in previous extensive benchmarks on the graphite-water reaction system in comparisons with higher level DFT and MP2 calculations.5 A comparison of B3LYP/6-31+G(d) versus B3LYP/6-311+G(3df,2p) in this work showed that the systematic error introduced by the use of the smaller basis set is on the order of -0.02 to -0.05 Å regarding intermolecular hydrogen bonds in water clusters and far less for bond distances involving water and graphite or intramolecular bonds, thus the use of the smaller basis set in the present studies is appropriate. We used the GAUSSIAN 03 revision C.143 implementation of ONIOM in all calculations, using standard convergence criteria and the “external” keyword to perform DFTB-D energy and gradient calculations with a standalone code. QM/MD simulations at the DFTB-D level were carried out for the dissociation products of OHx on the surface of a graphite monolayer model, using a Verlet algorithm with a 0.12 fs time step. This time interval was checked to conserve total energy to within 3 kcal/ mol accuracy during microcanonical dynamics with temperatureadjusted initial velocities for 5000 K. A target temperature of 5000 K was maintained by using the scaling of velocities approach with a 20% overall scaling probability. Although at very high temperatures electronic excitations should be abundant, internal conversion and thereby quenching of electronic excitations due to vibrational excitations is very efficient. We therefore consider it justified to consider only the electronic ground state in our QM/MD simulations. The reaction rate constants for the adsorption and dissociation reactions have been determined using the ChemRate program,44 based on ONIOM energetics and DFTB-D vibrational frequencies. 2.2. Description of Model Systems and Benchmarks for Defect Energies. Because different defects have different spatial extensions, we have selected different finite-size model systems for their description. We are investigating the reactions of OHx (x ) 1, 2) species with molecular models for 0D, SW, 1V, 1H1V, and 2V defective graphite (0001) surfaces. We are classifying the sizes of different graphene models roughly based on a coronene C24H12 core: structures related to C24H12 with less than a full hexagon ring surrounding them are denoted by the prefix S (small), structures related to circumcoronene C54H1845 with less than a full hexagon ring surrounding them are denoted by the prefix M (medium), and structures immediately related to dicircumcoronene C96H2445 are denoted by the prefix L (large). In Figure 1, low level ONIOM carbon atoms (treated at the DFTB-D level) are indicated by unfilled circles, and high-level ONIOM carbon atoms (treated at the B3LYP/ 6-31+G(d) level) are indicated by gray-filled circles. In monolevel calculations, the respective circle notations are used to indicate which level of theory was used. Black atoms mark the atoms corresponding to the defect structure itself and are always part of the ONIOM high-level if ONIOM was used in the calculation. The pristine dicircumcoronene L0D C96H24 (large 0D graphite model) structure possesses D6h symmetry and was only treated at the ONIOM(DFT:DFTB-D) level of theory. The 1V and

Quantum Chemical Studies on Graphite (0001) Surfaces

J. Phys. Chem. C, Vol. 111, No. 3, 2007 1357

Figure 1. Model structures for pristine and defective graphite (0001) surfaces. Low-level ONIOM carbon atoms (DFTB-D) are indicated by unfilled circles, and high-level ONIOM carbon atoms (B3LYP/6-31+G(d)) are indicated by gray-filled circles. In monolevel calculations, respective colors are used to the level of theory. Black atoms mark the atoms corresponding to the defect structure itself and are always part of the ONIOM high-level if applicable. Shown from left to right and top to bottom are: dicircumcoronene L0D C96H24, L1V C95H24, L1H1V HC95H24, M2V C64H20, S2V C30H14, SSW C32H14, and M2V-M0D C130H40.

1H1V defects also were studied using the large models, namely L1V C95H24 and L1H1V HC95H24, respectively. In these cases, the D6h symmetry is obviously broken, and in fact the defect is positioned off-center with respect to the formal D6h boundary of the graphene molecule. Figure 1 shows their respective positions in the ONIOM high-level treatment in which the divalent carbon is located either inside (L1H1V) or at the border (L1V) of the ONIOM model system. Finding that the energetics of the 1H1V defects are mainly dictated by the defects themselves and less influenced by the horizontal bulk influence,

we decided to investigate M2V C64H20 5/8/5 defect using both DFTB-D and ONIOM(B3LYP/6-31+G(d):DFTB-D). We have investigated smaller models for S0D C32H14, S1V C31H14, and S2V C30H14 using both DFTB-D and B3LYP/6-31+G(d). As for the Stone-Wales SW defect, we have only investigated the smallest member of the series, namely SSW C32H14, and treated it at the single-level B3LYP/6-31+G(d) for reasons which we will explain later in the text. Finally, we have estimated the effect of vertical graphite bulk in the case of the divacancy defect using an A-B stacked bilayer arrangement

1358 J. Phys. Chem. C, Vol. 111, No. 3, 2007

Xu et al.

TABLE 1: Defect Creation Energies ∆Ed for Singlet State Monovacancy S1V and L1V, and Singlet State Divacancy S2V and L2V Defects by the Removal of a Single and Two Individual Carbon Atoms in their 3P Electronic Ground Statea model

method

C atom

0D

1V

2V

L-Model

DFTB-D ONIOM DFTB-D B3LYP

0.00 0.00 0.00 0.00

0.00 0.00 0.00 0.00

381.7 345.9 341.9 307.7

568.5 486.0 533.4 452.1

S-Model

a The defect creation energy is defined as ∆Ed ) E(nV) + nE(C) E(0D) and is given in kcal/mol at both DFTB-D and ONIOM(B3LYP/ 6-31+G(d):DFTB-D) levels of theory for larger graphite L-models and DFTB-D and monolevel B3LYP/6-31+G(d) for smaller graphite S-models.

of M2V C64H20 and an M0D C66H20 species, yielding the M2V-M0D structure C130H40. Because of the high computational costs associated with this structure, we have only performed DFTB-D calculations on this bilayer species. To validate the use of the DFTB-D and ONIOM approaches, we have computed the energies for defect creation of L1V and L2V structures with respect to the removal of a single and two isolated carbon atoms in the 3P electronic ground state. It was found previously that 1V defects requiring approximately 170 kcal/mol formation energy at the GGA DFT level of theory for C removal20 are associated with a 3-fold degenerate Jahn-Teller relaxation from a D3h symmetric structure to pentagon/nonogon (5/9 ring combination) fused rings, containing one divalent carbon atom in the outer nonogon (9-membered ring, 9MR) boundary opposite to the 5/9 fusing bond,19 that their electronic ground state corresponds to a singlet (a triplet being higher in energy by about 12 kcal/mol),19 and that surfacediffusion of 1V defects costs about 23 kcal/mol per C atom hopping. The 5/8/5 2V defects do not possess a divalent center and require a total formation energy of less than that of two monovacancies by about 180 kcal/mol.29 This means it would be more favorable to eliminate two carbon atoms directly and create a 5/8/5 divacancy defect rather than two separate 5/9 monovacancy defects, but the probability for this to happen in experiment is rather low. Rather, monovacancies are expected to coalesce to divacancy defects in an exothermic process.29 Table 1 lists our results for defect creation energies ∆Ed ) E(nV) + nE(C atom) - E(0D), which use the energy of a carbon atom rather than the energy of a carbon atom in bulk graphite, therefore our numbers are shifted by the atomization energy of graphite (about 170 kcal/mol) when comparing to solidstate values. It can be seen that DFTB-D, ONIOM(B3LYP/631+G(d):DFTB-D), and monolevel B3LYP/6-31+G(d) approaches agree very well with the results of El-Barbary et al.19 for the creation of 1V and 2V defects when the graphite atomization energy is added, and that the energy lowering of the 2V defect with respect to the creation of two isolated monovacancies by about 206 kcal/mol at the ONIOM level for the large model system is comparable to the values reported elsewhere. We note that larger model systems increase the defect creation energy due to the greater loss of π-conjugation. Therefore, we are confident that our treatment of the graphite defect structures with the methods employed in this work is well justified. 3. Results and Discussions For reasons of simplicity, in the text we will discuss geometrical parameters and energetic data at the highest available levels of theory only. All computed structural parameters and energetics are given in the Supporting Information.

3.1. Dissociative Adsorption of OH Radical Attacks on Pristine Graphite. Figure 2 displays the PES and most important structural information for the reaction of the first hydroxyl radical with the L0D dicircumcoronene system. Reaction energies are defined here as ∆E ) Etot(compound) [Etot(L0D) + Etot(OH)]. The first step of the radical attack is naturally an edge-on coordination to form L0D-OH. This structure is characterized by the relatively long C-O bond with 1.46 Å (the B3LYP/6-31G(d) bond length for the C-O bond in methanol is 1.41 Å) and two significantly elongated C-C bonds of about 1.51 Å connecting the attacked carbon atom with the remainder of the honeycomb lattice. The associated ∆E of -4.9 kcal/mol is rather small at the ONIOM level of theory. In comparison, Jelea et al. had predicted a reaction energy of -8 kcal/mol for the OH radical on a pristine graphite surface,34 but they have used a nonhybrid PBE density functional prone to overbinding. The exothermicity for the formation of the initial L0D-OH complex is comparable to that of the NO2 reaction with L0D in which we predicted the corresponding ∆E value to be around -8 kcal/mol.11 Thus, both OH and NO2 can coordinate to a pristine graphite surface at low temperatures; however, in agreement with the results of Jelea et al.,34 the OH group will readily detach upon increasing the temperature. All further rearrangements of L0D-OH start with the formation of H-surface bond, which is required to overcome a 57.7 kcal/mol barrier at the transition state TS1. The alternative pathway, direct H dissociation from the coordinated OH fragment of L0D-OH, would require an even greater energy of about 85 kcal/mol, which is the bond dissociation energy of an O-H bond attached to an aromatic system.46 Therefore, below we will focus on a hypothetical dissociative adsorption pathway only, although it is clear that OH would immediately detach with high probability at high temperatures. As seen in Figure 2, the calculated O-H bond distance in the four-membered ring transition state TS1 is 1.591 Å, which is a significantly longer than 0.98 Å of the free OH radical. The resulting intermediate complex P1 lies 45.0 kcal/mol higher than the prereaction complex L0D-OH and possesses a low reverse-barrier of only 12.7 kcal/mol. In complex P1, the fourmembered ring structure of TS1 has changed into an epoxide structure and an adjacent hydrogenated center, separated by a 1.529 Å C-C single bond. The O-H distance is 2.362 Å in P1. From intermediate P1, the reaction may proceed via two different pathways: isomerization and vacancy defect formation. For the isomerization pathway, we investigated the feasibility of further O-H dissociation. Such a hypothetical process proceeds via several intermediate and transition state structures and leads to a product P3 in which H and O atoms are located in the 1,4-position of the affected hexagon with the epoxide group extending to the C-C bond of an adjacent hexagon. The entire process, P1 f P2 f P3 is found to be nearly thermoneutral with only 2.3 kcal/mol endothermicity and proceeds with a 37.1 kcal/mol rate-determining barrier at transition state TS2. Therefore, the relative position of the dissociated O and H defects on pristine graphite does not matter at all; the recombination is a fairly local event, and no long-distance effects of such adatom defects are observed. As shown in Figure 2, the initial barrier for vacancy defect formation calculated from the intermediate P1 is 23.6 kcal/mol at TS4, which is 13.5 kcal/mol lower than the barrier for isomerization. Analysis of the structure of TS4 shows that this barrier corresponds to C-C bond cleavage: the C-H bond of P1 migrates on the graphite surface and breaks the C-C bond to which O and H are attached. The product of this process is

Quantum Chemical Studies on Graphite (0001) Surfaces

J. Phys. Chem. C, Vol. 111, No. 3, 2007 1359

Figure 2. PES and main structural features of the dissociative adsorption reaction of OH on the defect-free graphite L0D C96H24 model. Relative energies are given at the ONIOM(B3LYP/6-31+G(d):DFTB-D) + ZPE(DFTB-D) level.

the DP1 intermediate with a very long C-C distance of 2.441 Å. The P1 f DP1 transformation is found to be exothermic by 5.8 kcal/mol. Somewhat surprisingly, this C-C cleaved structure is the second lowest energy structure on the potential energy surface of the reaction of OH with L0D and the lowest energetic structure of all investigated reaction pathways originating from the L0D-OH. In the next step, the O adatom can interact with the C center of a neighboring hexagon and thereby lead to the formation of the DP2 ketyl intermediate, which lies 67.0 kcal/ mol higher in energy than DP1 with ∆E(DP2) ) 100.3 kcal/ mol. The DP1fDP2 transformation occurs with a 67.9 kcal/ mol barrier at TS5. The calculated ∆E value of TS5 is 102.2 kcal/ mol; thus a reverse barrier for DP1 r DP2 is only 1.9 kcal/ mol. CO dissociation from the DP2 is predicted to occur with a modest 24.0 kcal/mol barrier at the transition state TS6 and is exothermic by 6.6 kcal/mol with respect to DP2; the parting CO fragment would thereby lead irreversibly to the formation of an L1H1V. Thus, the reaction of OH with pristine graphite can create an 1H1V defect, in particular under high-T/high-P conditions in which the chances for overcoming the high barriers are high. Previously, we reported the irreversible defect creation in pristine graphite by the reaction of NO in which CO left as a gaseous fragment and C/N replacement occurred on the surface.11 For comparison, the NO-induced defect formation occurs with a barrier that is 13 kcal/mol higher in ∆E value than the highest barrier for the OH-induced 1H1V defect creation at TS6. As can be seen in Figure 3, the dissociative addition of a second OH radical to the L1H1V hydrogenated monovacancy product leads to a second CO gaseous fragment and an irreversible L1H2V structure. This process is initiated by 35.6 kcal/mol drop in energy leading to a monohydrogenated, single hydroxy containing monovacancy defect with all dangling bonds saturated in the structure L1H1V-OH. This intermediate also could be a product of the reaction L1V + H2O. The mechanism of the L1V + H2O reaction will be discussed in the next section.

Subsequent H migration from the hydroxy radical in L1H1V-OH to the only left dangling bond in this structure leads to a highly exothermic dihydrogenated, keto group containing monovacancy defect DG-P1 with ∆E of -53.0 kcal/ mol relative to L1V + H2O. The formation of this nice looking (in terms of the presence of a dominant local Lewis-structure corresponding to cyclopentadiene-1-one) species occurs with a barrier of 31.1 kcal/mol at DG-TS2 with a ∆E of 5.5 kcal/mol relative to L1H1V + OH. A full incorporation of the O atom into the carbon framework via formation of a six-membered C5O ring in DG-P2 is energetically less favorable by 28.7 kcal/ mol and occurs with a 43.8 kcal/mol barrier at the transition state DG-TS3. The calculated energy destabilization is due to the re-emergence of a divalent carbon atom in the carbon framework, which is obviously less desirable. However, subsequent CO elimination to yield L2H2V is an irreversible step leaving a structure in which two carbon atoms were oxidized to CO by two OH molecules. The calculated reaction barrier of 72.4 kcal/mol at the transition state DG-TS4 can easily be overcome at high-T/high-P conditions. The overall CO elimination reaction by the OH from L1H1V is calculated to be exothermic by 16.2 kcal/mol. The sum of observations described above indicates that a once-formed 1H1V defect is very prone to further oxidation by OH radicals. 3.2. Reactions of H2O with Vacancy-Defective Graphite Surfaces. As mentioned above, previously we have shown that the reaction of water molecules (or water clusters) with defectfree graphite cannot irreversibly create defects on a graphite surface.5 This finding is in agreement with the dynamics simulations of Jelea et al. who found that OH and H fragments on pristine graphite surfaces associatively detach without a barrier.34 However, the reaction of a water molecule with a defective graphite surface promises to be more attractive, as theoretical calculations have shown of reduced SW barriers in the presence of hydrogen atoms,47 as well as exothermic addition of H2O to SW defects in a (7,0) carbon nanotube.48 Below, we

1360 J. Phys. Chem. C, Vol. 111, No. 3, 2007

Xu et al.

Figure 3. PESs and main structural features of the isostoichiometric dissociative adsorption reactions of OH on the hydrogenated monovacancy 5/9 defective graphite model L1H1V HC95H24 and of H2O on the monovacancy 5/9 defective graphite model L1V. Relative energies are given at the ONIOM(B3LYP/6-31+G(d):DFTB-D) + ZPE(DFTB-D) level.

discuss our findings on the reaction of a water molecule with 1V, 2V, and SW defective graphite surface models. 3.2.1. Reaction of H2O with L1V. As shown in Figure 3, the first step of the reaction of H2O with a monovacancy 5/9 defective graphite L1V is the H-OH bond cleavage, which proceeds with 18.2 kcal/mol barrier at the transition state DG-TS1 and leads to an exothermic L1H1V-OH intermediate, the same species we already encountered in the attack of an OH radical to L1H1V. This process is an exothermic process by 23.4 kcal/mol. Subsequent reactions follow the pathways as outlined in Section 3.1. The energy of the entire dihydrogenated, divacancy defect formation is found to be slightly exothermic with a ∆E of -4.0 kcal/mol. The barrier for adsorptive dissociation of H2O on the L1V at DG-TS1 is 64.2 kcal/mol lower than that on the L0D studied in our previous work,5 clearly indicating that the defect creation by water on OHand H-introduced defects is feasible under high-T/high-P conditions. 3.2.2. Reaction of H2O with Mono- and Bilayer 2V Defects. Figure 4 summarizes the potential energy surface of the reaction of water with a 2V 5/8/5 defective graphite surface using monolayer models M2V C64H20 at the integrated ONIOM(B3LYP/6-31+G(d):DFTB-D) and pure DFTB-D levels of theory and S2V C30H14 at the B3LYP/6-31+G(d) level. To illustrate the influence of vertical bulk π-stacking effects, in this Figure we also show single-level DFTB-D energetics for a bilayer model M2V C64H20 π-stacked with an M0D C66H20 species, yielding the combined M2V-M0D structure of C130H40. As is apparent from a general comparison of the energetics obtained using different models and methods, the ONIOMcalculated barriers for the reaction of M2V with water are 513 kcal/mol higher than the B3LYP-calculated barriers for the reaction of S2V with water. Because in the ONIOM studies the high-level model of M2V is identical in structure and theoretical treatment with the S2V model, this finding indicates that the 5/8/5-ring divacancy defect may have a strong influence from the surrounding carbon atoms, stabilizing the intermediates

and transition states by extended π-conjugation. On the other hand, the DFTB-D-calculated barriers and intermediate structures for the reaction of the bilayer M2V-M0D are up to 14 kcal/mol higher than those for the monolayer M2V system. This clearly indicates the existence of appreciable destabilization of outer-layer defects due to the second layer of graphite. In summary, we find in agreement with chemical intuition that horizontal bulk effect stabilizes defects due to π-conjugation, while vertical bulk effect destabilizes defects due to π-stacking. Generally, the reactions involving of 1,4-adducts (2V-RC2) require slightly lower energy barriers and lead to slightly less endothermic intermediates 2V-P2 than those involving 1,2adducts (2V-RC1) in which water attacks the 5/9 fused interring bond. Although the calculated barriers are large (50-67 kcal/mol) and the corresponding endothermicity of the reactions is high (30-60 kcal/mol) for the energetically most favorable 1,4-adduct pathways, these values are 17-30 kcal/mol lower than those for the reaction of H2O on defect-free graphite.5 The recombination and detachment of H2O, however, should still be more competitive, as further reconstruction pathways leading to irreversible defect formation are too complicated. This is in noticeable disagreement with the reactions of water with 1V defects, which are expected to readily proceed to further defect formation. In our calculations, 2V defects seem to act more like pristine graphite because of the absence of divalent carbon with dangling bonds as in 1V defects, thereby eliminating an important point of oxidative attack. Although this finding is in itself interesting, for practical applications the density of 2V defects is presumably rather low. As already mentioned, under experimental conditions their concentration should be dominated by the surface diffusion rate of 1V defects, and the widely dispersed monovacancy defects themselves are probably attacked by water and other reactive species before they can recombine to divacancy defects. 3.2.3. Reaction of H2O with SW DefectiVe Graphite. Figure 5 shows the B3LYP/6-31+G(d)-calculated PES and important structural features of the reaction of H2O with SSW C32H14.

Quantum Chemical Studies on Graphite (0001) Surfaces

J. Phys. Chem. C, Vol. 111, No. 3, 2007 1361

Figure 4. PES and main structural features of the dissociative adsorption reaction of H2O on divacancy 5/8/5 defective graphite model systems: monolayer 2V defects using ONIOM(DFT:DFTB-D) for compound M2V C64H20 (plain numbers), monolevel B3LYP/6-31+G(d) for compound S2V C30H14 (numbers in parenthesis), and DFTB-D for compound M2V C64H20 (numbers in square brackets); bilayer 2V defect M2V-M0D C130H40 (numbers in curly brackets) in which the relative energies are corrected for zero-point vibrational energy contributions at the ZPE(DFTBD) level.

Figure 5. PES and main structural features of the dissociative adsorption reaction of H2O on the SSW. Relative energies are given at the ONIOM(B3LYP/6-31+G(d):DFTB-D) + ZPE(DFTB-D) level.

The first step of the reaction is a formation of the 1,2- and 1,4-adducts, SSW-RC1 and SSW-RC2, respectively. The processes (O-H bond activation and subsequent isomerization) starting from both initial complexes require very large overall barriers with 93.4 kcal/mol and with 98.2 kcal/mol from SSWRC1 and SSW-RC2 adducts, respectively. The overall reactions H2O + SSW f SSW-P3 and H2O + SSW f SSW-P4 are found to be endothermic with ∆E values of 73.3 kcal/mole and 67.9 kcal/mol, respectively. Although the calculated reaction

barriers for dissociative adsorption of H2O on SSW are 2-7 kcal/mol lower than those for H2O on defect-free graphite,5 they are still very large, and dissociative adsorption of H2O on SW defects in graphite is therefore unlikely to occur. For this reason, we did not study larger model systems, which would merely lead to further destabilization of adduct structures. This result stands in remarkable contrast to the B3LYP/STO-3G calculations of Ellison et al. who predicted a 10 kcal/mol exothermicity for the reaction of H2O with SW defects on a (7,0) zigzag single-

1362 J. Phys. Chem. C, Vol. 111, No. 3, 2007

Xu et al.

walled carbon nanotube.48 In comparing our B3LYP/6-31+G(d) graphene result and their BLYP/STO-3G nanotube result, one has to consider two factors: (a) small-diameter carbon nanotubes such as (7,0) are much more reactive than largediameter tubes or graphite,49 and (b) the small STO-3G basis set consistently overestimates the binding energies.50

The rate constants (ki) given above are defined by51

d[X]surf/dt ) ki (θ/As) [X]g

4. Reaction Rate Constants Predicted by RRKM Theory

which has the unit of a flux, molecule cm-2 s-1. In this equation, θ represents the fraction of available surface sites, As is the surface area, and [X]g is the gas-phase concentration of OH in molecules/cm3.

We also have predicted rate constants of the above-presented processes. The following nine elementary reactions are involved for the prediction of rate constants of these processes:

5. QM/MD Simulations of OH Dissociative Adsorption Products on Graphite

k2 ) 1.42 × 10-12 exp(-35 400/T)

The DFTB-D based QM/MD simulations reported here were carried out at the target temperature of 5000 K, similar in spirit to our previous simulations of the high-T dynamics of the reaction of water5 and COx/NOx.10,11 with pristine graphite. The simulations trace the way “back” to initial conditions if feasible, and if not, we confirm the formed defects were formed irreversibly. Simulation times did not exceed 40 fs because processes occur very rapidly at the target temperature. The QM/MD simulations for the L0D-OH, P1, P2, DP1, and DP2 products of the reaction of OH with a pristine graphite sheet L0D C96H24 are shown in Figure 6a. The undissociated OH group in L0D-OH rapidly detaches and leaves the surface rapidly as free the OH radical, leaving behind a wildly bending but otherwise unmodified graphene sheet with the exception of edge defects caused by the finite size of the graphene model. For the products P1, P2, DP1, and DP2 with O and H atoms, we find that the attached individual addend atoms can migrate on the surface but remain on the surface. It is interesting to note that we observed reattachment of a temporarily recombined OH radical in the dynamics of the structure P1 in which OH seems to leave the graphene sheet at 30 fs with a C-O bond length longer than 2.5 Å, but dissociative adsorption follows at 40 fs. Generally, epoxide structures are seen to be unstable at the target temperature of 5000 K, and the associated three-membered rings break to create a wildly dangling ketyl C-O fragment on the surface. The QM/MD simulations for the L1H1V-OH, DG-P1, and DG-P2 products of second OH attack on L0D are shown in Figure 6b. The O-H bond in L1H1V-OH cleaves, and the H atom recombines to H2 with the H located on the 5/9 hydrogenated 1V defect edge, leaving a purely oxidized surface. In the simulations from both the DG-P1 and DG-P2 we observe instantaneous CO-formation, while the irreversibly created 2H2V structures are not annealed to their equilibrium 5/8/5 geometries yet. As usual, we observe in the QM/MD simulations of graphite at 5000 K that within 40 fs some C-C bonds of the graphene edges also have been broken, causing either the formation of dangling polyacetylene chains or SW transformations.

k3 ) 4.43 × 10-10 exp(-68 400/T)

6. Conclusions

k1: OH + L0D f P1

(1)

k2: OH + L0D f DP1

(2)

k3: OH + L0D f CO + L1H1V

(3)

k4: H2O + L1V f CO + L2H2V

(4)

k5: OH + L1H1V f CO + L2H2V

(5)

k6: H2O + SSW f SSW-P1

(6)

k7: H2O + SSW f SSW-P2

(7)

k8: H2O + 2V f 2V-P1

(8)

k9: H2O + 2V f 2V-P2

(9)

in which reactions (1), (2), and (3) can be treated as monowell association reactions by transition states TS1, TS4, and TS6, respectively. However, reactions (4) and (5) can be treated as two-well association reactions: DG-TS1

DG-TS2

H2O + L1V 98 L1H1V - OH 98 DG-TS4

DG-P1 98 CO + L2H2V and Barrierless

OH + L1H1V 98 DG-TS2

DG-TS4

L1H1V - OH 98 DG-P1 98 CO + L2H2V respectively. The predicted rate constants for reactions 1-9 in the temperature range from 1000 to 5000 K are (in units of cm3/s):

k1 ) 1.46 × 10-12 exp(-30 000/T)

k4 ) 9.37 × 10-12 exp(-28 000/T) k5 ) 2.95 × 10-12 exp(-20 900/T) k6 ) 1.12 × 10-12 exp(-46 000/T) k7 ) 1.09 × 10-13 exp(-53 400/T) k8 ) 2.09 × 10-12 exp(-38 500/T) k9 ) 9.18 × 10-14 exp(-35 000/T)

From the above discussion, we can draw the following conclusions: 1. The OH radical, which is abundant in propellant exhaust at high-T/high-P, can react with L0D and irreversibly lead to the formation of L1H1V and gas CO molecules. The calculated overall barrier of the reaction L0D + OH f L1H1V + CO is ca. 125 kcal/mol, which is 13 kcal/mol lower than that reported for the reaction L0D + NO f LN1V + CO.11 2. The dissociative adsorption of OH on the defective graphite model L1H1V has a much lower barrier than on the L0D: The reaction L1H1V + OH f L2H2V + CO is calculated to be

Quantum Chemical Studies on Graphite (0001) Surfaces exothermic by 16.2 kcal/mol, and proceeds by only ca. 36 kcal/ mol barrier. This indicates that once a defect is created on a pristine graphite (0001) surface, oxidatiVe erosion proceeds in this location at a higher speed. 3. By comparing our previous results of H2O reactions with a pristine graphite (0001) surface with the those of H2O attack on L1V, L2V, and SW defects, we find that the dissociative adsorption of H2O is only feasible on a L1V, which proceeds by only 48 kcal/mol barrier and produces gas CO molecules and L2H2V defect. The 2V and SW defects possessing only sp2 carbon ring members are much less prone toward attack by closed-shell species such as water. 4. The increase of π-conjugation (horizontal bulk effect) stabilizes both intermediates and transition state structures and

J. Phys. Chem. C, Vol. 111, No. 3, 2007 1363 reduces barriers by 5-13 kcal/mol, while the inclusion of the second layer into calculations (π-stacking or vertical bulk effect) increases the barriers up to 14 kcal/mol. Hence, in both cases the energy changes are rather moderate given the combustion conditions, and we conclude that a relatively small model of graphite surface is sufficient to elucidate the major reaction pathways of dissociative adsorption of combustion products on pristine and defective surfaces in high-T/high-P environments. 5. The magnitude of k1 for the OH dissociative adsorption on graphite is considerably greater than those of analogous reactions involving H2O, COx, and NOx. The magnitude of k3 for the dissociation of CO following the OH addition process on graphite is considerably greater than that of analogous reaction involving NO. The magnitudes of k6-9 for the H2O

1364 J. Phys. Chem. C, Vol. 111, No. 3, 2007

Xu et al.

Figure 6. (a) The dynamics structures for the QM/MD simulations of L0D-OH, P1, P2, DP1, and DP2 on the surface of the pristine L0D C96H24 graphite model at 5000 K using the DFTB-D method. (b) The dynamics structures for the QM/MD simulations of L1H1V-OH, DG-P1, DG-P2 on the surface of the L1H1V HC95H24 defective graphite model at 5000 K for producing gas CO molecules using the DFTB-D method.

dissociative adsorption on the defective graphite models are considerably greater than those of H2O dissociative adsorption on the defect-free graphite. 6. The high-T QM/MD simulations on L0D-OH, P1, P2, DP1, and DP2 intermediate structures of the reaction L0D + OH have shown that nondissociated OH radicals immediately leave the surface at temperatures of 5000 K, while the O and H

atoms of the dissociated OH remain attached to the surface. In these cases, CO formation is frequently observed, leading to the formation of larger vacancy defects, in the both cases of single and double attack of OH to pristine surfaces. Similarly, H2O attack on monovacancy sites leads to CO formation and divacancy formation, hence oxidative erosion of already defective graphite (0001) surfaces.

Quantum Chemical Studies on Graphite (0001) Surfaces In summary, we have shown that OH radicals can easily create monovacancy defects, which in turn are very reactive toward further attack by OH, H2O, COx, NOx, and other combustion products. We find that divacancy and SW defects are less reactive, as they do not contain dangling bonds. The studies of PESs and the dynamics of larger defect creation by high concentrations of OH radicals are one of the interesting future directions and will be the subject of future investigations. Acknowledgment. We gratefully acknowledge financial support from the Office of Naval Research under a MURI grant (BAA#03-012 MURI #12). M.C.L. acknowledges support from the National Science Council of Taiwan for a Distinguished Visiting Professorship at the National Chiao Tung University in Hsinchu, Taiwan. S.C.X. thanks the Cherry L. Emerson Center of Emory University for a Visiting Fellowship, and the authors acknowledge the Emerson Center for the use of its resources. Supporting Information Available: Table S1 lists names and total energies using ONIOM, B3LYP, and DFTB-D methods and imaginary frequencies using the DFTB-D method for all structures of the reaction pathways, and Table S2 lists the corresponding Cartesian coordinates. Figures S1 to S4 display all relevant structural parameters of these structures. This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) Keswani, S. T.; Andiroglu, E.; Campbell, J. D.; Kuo, K. K. J. Spacecr. Rockets 1985, 22 (4), 396-397. (2) Federici, G.; Skinner, C. H.; Brooks, J. N.; Coad, J. P.; Grisolia, C.; Haasz, A. A.; Hassanein, A.; Philipps, V.; Pitcher, C. S.; Roth, J.; Wampler, W. R.; Whyte, D. G. Nucl. Fusion 2001, 41 (12), 1967-2137. (3) Hawkes, E. R.; Sankaran, R.; Sutherland, J. C.; Chen, J. H. J. Phys. Conf. Ser. 2006, 16, 65-79. (4) Ferro, Y.; Marinelli, F.; Allouche, A.; Brosset, C. J. Nucl. Mat. 2003, 321, 294-304. (5) Xu, S.; Irle, S.; Musaev, D. G.; Lin, M. C. J. Phys. Chem. A 2005, 109, 9563-9572. (6) Porezag, D.; Frauenheim, T.; Koehler, T.; Seifert, G.; Kaschner, R. Phys. ReV. B 1995, 51, 12947-57. (7) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Phys. ReV. B 1998, 58, 7260-7268. (8) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. J. Chem. Phys. 2001, 114 (12), 5149-5155. (9) Kumar, A.; Elstner, M.; Suhai, S. Int. J. Quantum Chem. 2003, 95, 44-59. (10) Xu, S.; Irle, S.; Musaev, D. G.; Lin, M. C. The JANNAF 2005, 2005-0244.AE, 1-10. (11) Xu, S. C.; Irle, S.; Musaev, D. G.; Lin, M. C. J. Phys. Chem. B 2006, 110, 21135-21144. (12) Stone, A. J.; Wales, D. J. Chem. Phys. Lett. 1986, 128 (5,6), 501503. (13) Kaxiras, E.; Pandcy, K. C. Phys. ReV. Lett. 1988, 61 (23), 26932696. (14) Eggen, B. R.; Heggie, M. I.; Jungnickel, G.; Latham, C.; Jones, R.; Briddon, P. R., Science 1996, 272 (5258), 87-9. (15) Bettinger, H. F. J. Phys. Chem. B 2005, 109, 6922-6924. (16) Hashimoto, A.; Suenaga, K.; Gloter, A.; Urita, K.; Iijima, S. Nature 2004, 430, 870. (17) Urita, K.; Suenaga, K.; Sugai, T.; Shinohara, H.; Iijima, S. Phys. ReV. Lett. 2005, 94, 155502/1-15550.2/4. (18) Xu, C. H.; Fu, C. L.; Predraza, D. F. Phys. ReV. B 1993, 48 (18), 13273-13279.

J. Phys. Chem. C, Vol. 111, No. 3, 2007 1365 (19) El-Barbary, A. A.; Telling, R. H.; Ewels, C. P.; Heggie, M. I.; Briddon, P. R. Phys. ReV. B 2003, 68, 144107/1-14410.7/7. (20) Li, L.; Reich, S.; Robertson, J. Phys. ReV. B 2005, 72, 184109/118410.9/10. (21) Rossato, J.; Balerie, R. J.; Fazzlo, A.; Mota, R. Nano Lett. 2005, 5 (1), 197-200. (22) Lu, A. J.; Pan, B. C. Phys. ReV. B 2005, 71, 165416/1-16541.6/6. (23) Tien, L. G.; Tsai, C. H.; Li, F. Y.; Lee, M. H., Phys. ReV. B 2005, 72, 245417/1-24541.7/6. (24) Ding, F. Phys. ReV. B 2005, 72, 245409/1-24540.9/7. (25) Berber, S.; Oshiyama, A. Physica B 2006, 376-377, (272-275.). (26) Kim, G.; Jeong, B. W.; Ihm, J. Appl. Phys. Lett. 2006, 88, 193107/ 1-19310.7/3. (27) Carlsson, J. M.; Scheffler, M. Phys. ReV. Lett. 2006, 96, 046806/ 1-04680.6/4. (28) Thrower, P. A.; Mayer, R. M. Phys. Status Solidi A 1978, 47 (1), 11-37. (29) Lee, G.-D.; Wang, C. Z.; Yoon, E.; Hwang, N.-M.; Kim, D. Y.; Ho, K. M. Phys. ReV. Lett. 2005, 95, 205501/1-20550.1/4. (30) Mulcahy, M. F. R.; Young, B. C. Carbon 1975, 13, 115-24. (31) Lamoen, D.; Persson, B. N. J. J. Chem. Phys. 1998, 108, 33323341. (32) Sorescu, D. C.; Jordan, K. D.; Avouris, P. J. Phys. Chem. B 2001, 105, 11227-11232. (33) Incze, A.; Pasturel, A.; Chatillon, C. Surf. Sci. 2003, 537, 55s63. (34) Jelea, A.; Marinelli, F.; Ferro, Y.; Allouche, A.; Brosset, C. Carbon 2004, 42, 3189-3198. (35) Ferro, Y.; Allouche, A.; Marinelli, F.; Brosset, C. Surf. Sci. 2004, 559, 158-168. (36) Volpe, M.; Cleri, F. Surf. Sci. 2003, 544, 24-34. (37) Ferro, Y.; Marinelli, F.; Allouche, A. Chem. Phys. Lett. 2003, 368, 609-615. (38) Ferro, Y.; Marinelli, F.; Allouche, A. J. Chem. Phys. 2002, 116, 8124-8131. (39) Matsubara, T.; Maseras, F.; Koga, N.; Morokuma, K. J. Phys. Chem. 1996, 100, 2573-2580. (40) Maseras, F.; Morokuma, K. J. Comput. Chem. 1995, 16, 11701179. (41) Dapprich, S.; Komaromi, I.; Byun, K. S.; Morokuma, K.; Frisch, M. J. J. Mol. Struct. 1999, 461-462, 1-21. (42) Frauenheim, T.; Seifert, G.; Elstner, M.; Hajnal, Z.; Jungnickel, G.; Porezag, D.; Suhai, S.; Scholz, R. Phys. Status Solidi B 2000, 217, 41. (43) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Iyengar, S. S.; Millam, J. M.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Ehara, M.; Toyota, K.; Hada, M.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Kitao, O.; Nakai, H.; Honda, Y.; Nakatsuji, H.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Adamo, C.; Jaramillo, J. J.; Cammi, R.; Pomelli, C.; Gomperts, R.; Stratmann, R. E.; Ochterski, J.; Ayala, P. Y.; Morokuma, K.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, Revision C.1; Pittsburgh, PA, 2004. (44) Mokrushin, W.; Bedanov, V.; Tsang, W.; Zachariah, M.; Knyazev, V. ChemRate, Version 1.20; National Institute of Standards and Technology: Gaithersburg, MD, 2003. (45) Dias, J. R. THEOCHEM 2002, 581, 59-69. (46) Klein, E.; Lukes, V. J. Mol. Struct. 2006, 767, 43-50. (47) Nimlos, M. R.; Filley, J.; McKinnon, J. T. J. Phys. Chem. A 2005, 109 (43), 9896-9903. (48) Ellison, M. D.; Good, A. P.; Kinnaman, C. S.; Padgett, N. E. J. Phys. Chem. B 2005, 109, 10640-10646. (49) Zheng, G.; Wang, Z.; Irle, S.; Morokuma, K. J. Am. Chem. Soc., in press (50) Irle, S.; Lischka, H. J. Chem. Phys. 1995, 103 (4), 1508-1522. (51) Rettner, C. T.; Ashfold, M. N. R. Dynamics of Gas-Surface Interactions. In Kinetics of Surface Reactions; Rettner, C. T., Ashfold, M. N. R., Eds.; The Royal Society of Chemistry: Cambridge, England, 1991; Chapter 5.