Car–Parrinello Molecular Dynamics + Metadynamics Study of High

Article pubs.acs.org/JPCC

Car−Parrinello Molecular Dynamics + Metadynamics Study of HighTemperature Methanol Oxidation Reactions Using Generic Collective Variables Shaohui Zheng and Jim Pfaendtner* Department of Chemical Engineering, University of Washington, Seattle, Washington 98195-1750, United States S Supporting Information *

ABSTRACT: We used Car−Parrinello molecular dynamics (CPMD) and metadynamics in conjunction with the recently introduced social permutation invariant collective coordinates to study the mechanism of high-temperature methanol oxidation. Using a set of biased MD trajectories, we collected specific elementary reactions that arise during the simulations and assembled their connectivity in a small reaction network. A subset of the reaction network generated with metadynamics is compared to a consensus reaction network generated from many literature sources, and the many similarities indicate that this approach may be a useful way to enumerate bimolecular radical reactions in complex systems. We also demonstrated some intrinsic similarities to atomic contact maps used in our metadynamics approach and the reaction matrix/operators that are found in common mechanism generation algorithms. Extending the capabilities of these new generic collective variables for the study of complex reaction networks can help overcome limitations of enhanced sampling methods in the study of chemical reactions.

1. INTRODUCTION Theoretical and computational tools to help to predict and understand the kinetics and thermochemistry of complex reacting systems are a widely investigated research subject due to the ubiquitous importance of reactions in disparate fields of study. The use of ab initio molecular dynamics (AIMD) based approaches, frequently the Car−Parrinello (CPMD)1 scheme, to achieve this has been greatly facilitated through the application of so-called “enhanced sampling” methods, which can dramatically improve sampling of phase space and thereby facilitate the calculation of reaction rates and thermodynamics. A widely used approach is the metadynamics (MetaD) method,2,3 which enhances exploration of phase space through the construction of a time-dependent bias potential that is based on the system’s exploration of coarse descriptors known as collective variables (CVs). The MetaD approach has numerous applications to the study of chemical reactions.4,5 As evidenced by the wide range of applications, the combination of CPMD and MetaD has been successfully applied to different systems such as azulene-to-naphthalene reactions, the Beckmann rearrangement, CO2 + H2O reactions, water splitting driven by solar energy with transition metal catalysts, and base-induced elimination reactions.2,3,5−8 We note that this is just a small selection out of the hundreds of possible choices, meant only to illustrate the wide range of chemistries accessible with this versatile approach. A common feature of most investigations (of this type) using CPMD+MetaD is that the users need identify, a priori, a limited subset (typically fewer than 5) of chemical reactions using their intuition about which reactions are important. This could be fine in cases when only one or two reactions take © 2014 American Chemical Society

place, but if complex mechanisms are present, then the researcher must hypothesize the specific rate limiting or interesting steps to investigate. Furthermore, detailed knowledge in the form of hypotheses about how the reactions proceed is often required in order to successfully formulate a subset of CVs that facilitate the desired enhanced sampling. This often proves to be both nontrivial and essential in order to avoid hysteresis when constructing free-energy surfaces (FESs).2,5 Potential ways to address this challenge of CV selection include path MetaD,9 which has already been applied to studies of chemical reactions,6 or the new reconnaissance MetaD + sketch-map approach,10,11 which has also been used to study reactions.12 An additional new method was recently introduced by Pietrucci and Andreoni to use MetaD to study chemical transformations.13 They developed a new set of special CVs, the social permutation invariant (SPRINT) coordinates, to provide a generic and universal discrimination between chemical compounds. When biased with MetaD, many chemical transformations (reactions) occurred, indicating the potential to overcome the aforementioned challenges of system-dependent CV selection. Inspired by these accomplishments, we set out to extend the capabilities of these newly developed SPRINT coordinates to complex reaction networks involving free-radical reactions and applied the method to study high-temperature combustion processes. If successful, this approach would be a new point of view in the generation of complex reaction networks involving Received: January 13, 2014 Revised: April 14, 2014 Published: April 21, 2014 10764

dx.doi.org/10.1021/jp500398k | J. Phys. Chem. C 2014, 118, 10764−10770

The Journal of Physical Chemistry C

Article

energy conserving. Our simulation protocol was as follows: prior to MD simulations the wave function was converged, next a 5.1 ps thermal equilibration (no MetaD is applied during equilibration) was performed followed by another wave function optimization, and after this, production runs of 24.0 ps in length were performed. We chose the 24.0 ps time length for two reasons: (1) the lowest free-energy products, carbon dioxide/monoxide and water, are produced over this time scale; (2) multiple trials showed it is long enough to observe the time evolution of interesting elementary reactions. The number of production simulations (30) was selected on the basis of analyzing the random O2/CH3OH starting configurations to ensure the majority of different possible collision pathways were sampled (see Supporting Information (SI) Figure S1). As detailed below, we tested convergence of the simulations based on the appearance of the intermediate radical species in multiple different simulations. Namely, we have compared the generated reaction networks obtained with different numbers (25 and 30) of simulations and found no substantial differences. In order to keep reactants close together and remove diffusive limitations, we added a harmonic restraint of strength 0.5 hartree/Bohr2 (9.4 × 105 (kJ/mol)/nm2) to limit the radius of gyration (calculated on all atoms in the system) to 0.42 nm. Given that most reactions in methanol oxidation have barriers and would not occur on typical CPMD simulation time scales (658 K) oxidation of methanol.14−17 There are many motivations for selecting methanol combustion as a model system for this study: (1) the chemical reactions of methanol with oxygen require crossing high-energy barriers (e.g., the commonly used value of 44.9 kcal/mol16 for ignition of methanol), necessitating the use of enhanced sampling to accelerate ab initio molecular dynamics to observe bonds’ breaking and forming; (2) although methanol is the simplest alcohol, its oxidation process to carbon monoxide and dioxide and water is quite complex (a recently published comprehensive reaction mechanism at combustion temperatures contains around 100 reactions and 20 species);16 and (3) experimental and theoretical researchers have systematically studied this system, measured and predicted intermediates and products at different temperatures, pressures, and ratios of CH3OH to O2, and suggested many comprehensive mechanisms.14−16,18,19 All of these previous research results provide a useful benchmark for the method. In the present work, using the biased MD trajectories with the new generic SPRINT coordinates, we will show that the suggested approach facilitates collection of large sets of elementary reactions and allows for meaningful assembly of the reactions into an overall mechanism. This work is organized as follows. After introducing computational details and the algorithmic framework for mechanism generation with AIMD, the accuracy of the approach for estimating reaction energies is benchmarked with a few trial calculations. Following this, we describe a set of 30 simulations of biased methanol/O2 reactions at 1000 K starting from different initial configurations. The results show different reaction pathways and illustrate the complexity of the underlying reaction potential energy surface (PES). Finally, we summarize the current state of this new approach in the context of realizing the goal of completely automated physics-based network generation.

2. COMPUTATIONAL DETAILS Using the specific methods described below, we systematically performed successive series of CPMD+MetaD calculations, and by carefully observing the results it was possible to enumerate individual unique elementary reaction steps that occurred during the simulations. The CPMD simulations used planewave density functional theory (DFT) with the code CPMD v3.15.11,20 enhanced by the PLUMED plugin.21 The Becke exchange22 and Lee−Yang−Parr correlation (BLYP)23 density functional was selected along with the Troullier−Martins24 norm conserving pseudopotentials (PP) for oxygen, carbon, and hydrogen atoms.25 We employed the local spin density approximation (LSDA) for radicals with unpaired electrons and set a multiplicity of 3 for the model system composed of a single methanol and a single oxygen molecule. As is commonly done, we replaced all of the hydrogen atoms with deuterium (isotope = 2), whose heavier mass permits a longer time step.7 The time step in CPMD was set to 6 au (0.15 fs), and the fictitious electron mass 1000 au. The simulations were performed in the canonical ensemble (NVT) in an isolated 16.0 Å cubic box at a temperature of 1000 K. Temperature was controlled with Nose−Hoover thermostats26,27 of chain length three with frequencies for the atoms and fictitious electron kinetic energy of 2700 and 10000 cm−1, respectively. The fictitious electron kinetic energies were obtained through microcanonical ensemble (NVE) simulations for all individual systems, which were also used to verify that all parameters were

N λmax νimax,sorted

Si =

νimax =

aij = 10765

1 (λmax )M

(1)

∑ aijM νjmax j

(2)

1 − (rij/r0)n 1 − (rij/r0)m

(3)

dx.doi.org/10.1021/jp500398k | J. Phys. Chem. C 2014, 118, 10764−10770

The Journal of Physical Chemistry C

Article

where Si is the SPRINT coordinate, N is the number of total atoms, aij denotes the contact matrix (described with common switching functions3 through the variables r0 (C−C, 2.65 Å; C− O, 2.65 Å; C−H, 2.22 Å; O−H, 2.22 Å; H−H, 2.22 Å), n (6), and m (12) (the parenthetical values denote those used in this work)), aM ij is the number of walks of length M between points i and j, vmax is the largest eigenvector of the contact matrix, and λmax is the largest eigenvalue. The term “sorted” simply means sorting the eigenvector from its smallest to the largest components and putting the same type of atom together in the vector. The SPRINT coordinates (Si) of all atoms in the reactions were employed as collective variables. We note for clarification that Bucko has noted that it is important to remember that in the limit of sampling all degrees of freedom in the system, the underlying PES not the FES is actually sampled, intrinsically limiting the inclusion of entropic effects in the exploration of phase space.31 However, in the context of certain fields such as combustion the sampling of the PES is thought by some32 to be an important component of detailed kinetic modeling of reaction networks. The MetaD bias potential was constructed as a summation of Gaussian hills, represented as VG(S , t ) =

∫0

t

are performed for methanol thermal dissociation both with and without an inert molecule (N2) present. This was to reflect the fact that many elementary methanol combustion reactions include an inert collision partner, presumably to attempt to model collisional energy transfer. We confirmed that the existence of N2 molecule does not have substantial influence on the C−O bond fission of methanol. 3.1. Methanol Oxidation Simulations. To study the time evolution of various species in methanol oxidation and the capabilities of our algorithm to enumerate elementary reactions, we set up a series of CH3OH + O2 simulations with 30 random O2 configurations relative to the position and orientation of CH3OH (illustrated in SI Figure S1). This effectively allowed us to consider various stochastic effects in the system since, unlike static reaction searches using transition state theory, the results of biased CPMD simulations potentially could be dependent on starting configurations via a collision trajectory. During the equilibration stages of our calculations (5.1 ps at 1000 K) we did not observe any spontaneous reactions of oxygen and methanol in any of the 30 simulations. After the equilibration step, all 30 simulations were continued using MetaD with SPRINT coordinates to evolve the atomic topology of the system using the parameters described above. The aggregate production period for these simulations was 720 ps. Similar to the initial step observed in the benchmark methanol only simulations (SI), we universally observed fission at the C−O bond; i.e., CH3OH + O2 → •CH3 + •OH + O2, in all 30 simulations with average reaction energy of 84.3 ± 7.0 kcal/mol. Following the dissociation, we observe a huge variety in different elementary reaction channels and observed species. As a first step to analyze the properties of MetaD+SPRINT simulations for complex radical reactions, we tracked the appearance of major stable species over the course of the simulations. Shown in Figure 1 is a count (over all 30

⎛ d (S (R ) − S (R(t ′)))2 ⎞ i ⎟⎟ dt ′ ω exp⎜⎜ −∑ i 2 2 σ ⎝ i=1 ⎠ i (4)

where VG is the history-dependent potential energy, t is the simulation time, ω is an energy rate of the Gaussian hill, σi is the width for the ith CV, and d is the total number of CVs. After some initial calculations to tune the parameters, we selected a Gaussian hill height of 2.8 kcal/mol with the frequency of one hill per 500 MD steps and a Gaussian width of 1.5 (dimensionless according to eq 1). To test the accuracy and tune the performance of our CPMD calculations, we performed a number of calculations with the above approach and compared the predicted changes in Kohn−Sham energy with reaction energies obtained using static electronic structure calculations in the code Gaussian,33 which are further described in the SI. To estimate this quantity from biased CPMD+MetaD simulations of the C−O bond fission energies (reported below), we took the difference in the average Kohn−Sham (potential) energy before and after the reaction event. An upper estimate of the error for this value was calculated by summing the respective standard deviation of the Kohn−Sham energy in the plateau regions representing the reactants and products. The validity of this approach is demonstrated in the SI by comparing the values to those obtained from static DFT calculations.

3. RESULTS AND DISCUSSION It is well-known that the many parameters input into CPMD simulations, such as fictitious electron mass, cell size, and cutoff of plane-wave PP, can influence the accuracy of the results.34−36 Because of this, we performed a short series of benchmark calculations to validate our chosen set of parameters and determine if the reaction energies observed during MetaD simulations biasing the SPRINT coordinates are meaningful. In the SI we demonstrate that the reaction energies observed as identified by concomitant changes in the SPRINT coordinates and the CPMD Kohn−Sham energy do indeed reflect those obtained in the typical manner through static DFT calculations (SI Table S1 and Figures S2−S6). The benchmark simulations

Figure 1. Cumulative time evolution of count numbers of selected products (H2O, CO, and CO2) and radicals with the largest count numbers (•CHO, 3•C•H2, •CH3, H•, •OH, and HO2•) in the 30 CH3OH + O2 simulations from time 17.4 to 29.0 ps.

simulations) of selected interesting species. The full data are tabulated in the SI (Table S3). The reason for looking at the aggregate simulations is to try and get an overall picture about what the MetaD trajectories are revealing about methanol oxidation given that the biased CVs are totally generic about the types of transformations that could take place. Several salient features can be observed from analysis of the evolution 10766

dx.doi.org/10.1021/jp500398k | J. Phys. Chem. C 2014, 118, 10764−10770

The Journal of Physical Chemistry C

Article

H2O → HCO2• + H2O + H• → CO2 + H2O + 2H• → CO2 + • OH + H• + H2 → CO2 + 3•O• + 2H2. Here, the superscript number means the multiplicity of the compound. In this simulation, eight radicals are generated including •CH3, •OH, CH3OO•, 3•CH2OO•, 3•O•, 3H2C(O•)2, HCO2•, and H•. And four molecules (CH2O, CO2, H2O, and H2) are observed. We next compared the graph-theoretic bond-electron (reaction) matrices adopted in many automated network generation approaches such as the NetGen/reaction mechanism generator algorithms (RMG)38,39 and the contact matrices13 used in defining the SPRINT CVs for an elementary reaction (Figure 3). Since the two approaches are both fundamentally based on atomic connectivity, we hoped an illustrative comparison between the two representations of chemical reactions would help to connect the new approach using SPRINT CVs with the classic graph-theoretic approach that has achieved great success in the modeling of complex reaction networks. The contact matrices and reaction matrices for reaction CH3OH + O2 = •CH3 + •OH + O2 are strikingly similar. For example, the connection between C1 and O1 is marked as 1 (A21 or A12) in both contact and reaction matrices. The main differences are related to handling unpaired electrons and bond type. In the reaction matrices the A33 and A34 elements (O2−O2 and O2−O3) have values of 1 and 2 to represent one unpaired electron and a double bond, respectively; in contrast A33 and A34 in the contact matrices are 0 and 1, respectively. We point out that according to eq 3 the elements of contact matrices also depend on the corresponding bond lengths and bond characters and are not integers; for clarity we round to one significant digit for resulting in values of 0 or 1 for all elements. This analysis suggests that it may be possible to use the simulations biased with the SPRINT coordinates to discover new reaction classes and their corresponding “operators” (Figure 3, left) to augment existing reaction network generation approaches. It can also be seen from Figure 3 that the chemical reactions might be recorded by monitoring the potential energy (the Kohn−Sham energy). However, this is a degenerate measure as there are multiple states with the same distinct energy. Additionally, as shown in Figure 4, it is possible to strictly monitor the progression of species through observing changes in the SPRINT coordinates. We can see that the first (also highest) energy potential energy barrier is for the bond fission of C−O in methanol as shown in Figure 3, and •CH3 radical is crucial to activate the stable 3O2 molecule. 3H2C(O•)2 and HCO2• intermediates are important for formation of CO2 as shown in Figure 3. In Figure 4, we can see that, at time 7.5 ps, all the SPRINT coordinates are in a narrow range, which shows that the reactants CH3OH + O2 are close to each other. At time around 8.0 ps, the carbon SPRINT coordinate drops first due to the C−O bond fission and then increases by 1 because of the formation of CH3OO• radical. Similarly, at time 10.9 ps, an oxygen atom flies away and its SPRINT coordinate changes close to zero; at time 15.1 ps, 3H2C(O•)2 + H2O are formed and the weak hydrogen bond connects two fragments; therefore their SPRINT coordinates’ curves are close to each other again. In fact, we can observe that the SPRINT coordinates characteristically change along with chemical reactions in all 30 simulations. Summarizing all of the reaction events, we observed in the aggregate simulations, we compiled a reaction network that shows the connectivity between all of the carbon-containing species in the system (Figure 5; also SI Table S2). As indicated

of each of the simulations. Oxidation of carbon-containing species is steadily observed, and all methanol molecules are totally consumed in all simulations by a production simulation time of 24 ps. As shown in SI Table S3, following methyl radical, the methylperoxy (CH3OO•) and 3H2C(O•)2 radicals are the next two largest carbon-containing intermediates observed and play crucial roles for subsequent reactions. In the later stages of the simulations, once methyl, methylperoxy, and 3H2C(O•)2 radicals are consumed, we see much larger increases in the number of small molecule species (CH2O), forming a variety of radicals (•CHO, 3•CH2, and H•, etc.) as well as the thermodynamically stable species H2O and CO2, whereas CH4 is always present in the second half of the simulation, although in comparatively small amounts. Both formaldehyde and formyl radical are formed in comparatively large amounts, suggesting they could be a gateway between the substrate and final products CO/CO2. To confirm this, we continued all the simulations for a short amount and observed the expected decrease in formyl radical (just as formaldehyde begins to decrease around 26 ps). We also noticed that hydroxyl radicals (regarded as important species in combustion chemistry19,37) are always present with a stable large amount (around 15) in our production runs. Meanwhile, the numbers of hydrogen atoms and •CHO steadily increases while the number of water regarded as a quenching molecule in combustion keeps growing. Also noteworthy is that the number of hydroperoxyl (HO2•) radicals is quite stable. In addition, we compiled the distribution of all 30 carbon atoms in the aggregate simulations in the form of a histogram (Figure 2). The •CHO radical is dominant, suggesting its importance as an intermediate to generate CO and CO2 products.

Figure 2. Histogram of final counts of carbon-containing species at the end of the 30 simulations. Molecules (HCOOH, CO, CO2, CH2O, and CH4) and radicals (3H2C(O•)2, HCO2•, •CHO, 3•C•H2, 4•CH, and •CH3) from 30 biased CPMD simulations of CH3OH + O2.

Figure 3 demonstrates the time evolution of Kohn−Sham energy of a representative simulation (no. 4; similar time series are given for all simulations in SI Figures S6−S34). For illustrative purposes, we detail the biased-time evolution observed in this typical simulation that eventually yields CO2. We observed the progression of elementary reactions as CH3OH + O2 → •CH3 + •OH + O2 → CH3OO• + •OH → 3• CH2OO• + H2O → CH2O + •O• + H2O → 3H2C(O•)2 + 10767

dx.doi.org/10.1021/jp500398k | J. Phys. Chem. C 2014, 118, 10764−10770

The Journal of Physical Chemistry C

Article

Figure 3. Contact matrix and reaction matrix comparison of CH3OH + O2 → •CH3 + •OH + O2 reaction, and time evolution of elementary reactions and Kohn−Sham energy of one simulation (no. 4).

this could be potentially due to many factors. In Figure 5, we show that many carbon radicals can interconvert to each other. For example, CH3OO• ↔ •CH3, and •CHO ↔ •COOH. In addition, visual inspection of the diagram shows that the most important nodes in the reaction network are CH2O (five connections), •CHO (four connections), 3•CH2 (four connections), H2COO• (four connections), and •CH3 (four connections). Four of these are radical intermediates. In Figure 5, the reaction network compiled from our CPMD+MetaD simulations is compared to a consensus reaction network published in a recent methanol combustion review paper by Skodje et al.16 Here we use the term “consensus” because methanol combustion has been well-studied with experiment and theoretical calculations. We note that a consensus network only contains a minimum set of key reactions needed to properly reproduce experimental measurements of the species vs time concentration in batch experiments. We find many similarities between the two networks; for example, both networks show •CH3, CH2OH, •CHO, and CH2O may play important roles. However, there are some significant differences: first the C2H6 molecule and the CH3O• radical are not observed in our simulations, one of which is a key node in the consensus network; second the network from our simulations contains more radical intermediates than the review network. The absence of C2H6 is easy to understand as the simulations performed in this work are with a single methanol molecule. On the other hand, the absence of CH3O• is less clear. Our current hypothesis is that by expanding the simulations to include other dominant radicals (e.g., performing CH3OH +

Figure 4. Time evolution of SPRINT coordinates of one simulation (no. 4), and six snapshots at key time steps. M−O is a designation for the two atoms initially present as molecular oxygen.

in the figure, some reactions in our simulations are observed in both directions (indicated with a bidirectional arrow), whereas others are only observed in the direction drawn. While the lack of reversibility in MetaD sampling can be problematic for reconstruction of the projection of the FES onto the CVs, in this case our primary goal is meaningful exploration of phase space (i.e., tabulating chemical reactions), and we are not concerned about the lack of reversibility in some reactions as 10768

dx.doi.org/10.1021/jp500398k | J. Phys. Chem. C 2014, 118, 10764−10770

The Journal of Physical Chemistry C

Article

Figure 5. The networks of carbon species summarized from 30 CH3OH + O2 simulations (left) and from a 2010 methanol review paper16 (right). The area and color of circles represent the connectivity of different carbon species: the larger/darker areas have more connections. A bidirectional arrow denotes that a two-way reaction is observed. •

OH simulations) and performing simulations with two carboncontaining species, we would indeed observe the missing C2H6 and CH3O• and may make substantial progress toward obtaining a fuller picture of methanol combustion. Regarding the additional radical intermediates we observed, we note that some intermediates such as CH3OO• and HCOOH are skipped in the consensus network,16 presumably due to the overall minor influences they make on methanol combustion. We stress that the specific goal of this work is not to reproduce the consensus network but to explore the properties and potential uses of large scale simulations of ab initio dynamics biased with MetaD + SPRINT CVs. Therefore, it is not surprising that the two reaction networks in Figure 5 are different, since the consensus network is not an exhaustive list of all chemical reactions taking place but a list of the dominant reactions (as defined by the characteristic reaction rates using mass action elementary rate laws). We also point out that at this time it is impossible to determine whether the extra reactions we observed are either important reactions or superfluous ones that do not contribute to the main methanol oxidation pathways (in the case that we have identified extra reactions are not included in the consensus mechanism) or whether through the implementation of the SPRINT coordinate sampling we have introduced heuristics-based parameters (i.e., Gaussian deposition rate) that have caused us to miss some reactions. Further work is needed to systematically determine the importance of new reaction channels discovered with this approach.

interest in this approach is the fact that there are no userspecified inputs about the types of reactions that can or cannot take place; instead the reactions are scanned by a nonspecific sampling of the underlying potential energy surface. Second, while at first glance the computational cost of this approach may appear daunting, we note that these calculations were completed with a modest allocation on a national HPC resource and supplanted by local campus computing. Without question, further work is merited to understand ways to improve the computational efficiency, but with current capabilities many new studies are now possible. Third, the next key step in developing the utility of biased MetaD +SPRINT CV calculations to enumerate complex reaction networks should be further algorithmic development to automate the process of identification of new reactions, products, and intermediates. We provide the proof of concept for this in Figure 3 and propose that combining our approach with other algorithms that track averages and fluctuations on the fly could help achieve this. For example, the recently introduced variant of metadynamics with adaptive Gaussians40 is an ideal candidate. As demonstrated in Figure 3, once the reactants and products have been identified, classification of the reaction type can be facilitated through the use of chemical graph operators.

4. CONCLUSION In this work, we applied and extended the application of new generic collective variables to the study of complex radical reaction networks. Biasing AIMD simulations with the metadynamics framework (specifically the SPRINT collective variables) proved to be a useful way to enumerate large numbers of elementary reactions in high-temperature methanol oxidation. In conclusion we make several salient observations about the use of these special collective variables. First, given the many similarities between the reaction network generated from simple 1 methanol +1 O2 metadynamics simulations (cf., Figure 5 left and right), we are confident that this approach is underpinned by a physically meaningful enhanced sampling of chemical reactions that could eventually lead to new capabilities for systematic generation of reaction networks. Of particularly

Text describing computational details of validation calculations and including reaction energy comparison, methanol thermal decomposition, and accompanying references, figures showing reaction energies observed as identified by concomitant changes in the SPRINT coordinates and the CPMD Kohn− Sham energy compared to static all-electron DFT calculations, simulated results for methanol thermal dissociation with and without an inert molecule present, time evolutions of Kohn− Sham energies in CH3OH, CH3OH + N2, and various methanol oxidation simulations, and configurations used in the CPMD simulation, and tables listing reaction energy comparisons, summarized reaction channels, and frequency counts of the carbon compounds of the simulations. This material is available free of charge via the Internet at http:// pubs.acs.org.

■

ASSOCIATED CONTENT

S Supporting Information *

10769

dx.doi.org/10.1021/jp500398k | J. Phys. Chem. C 2014, 118, 10764−10770

The Journal of Physical Chemistry C

■

Article

(18) Klippenstein, S. J.; Harding, L. B.; Davis, M. J.; Tomlin, A. S.; Skodje, R. T. Uncertainty driven theoretical kinetics studies for CH3OH ignition: HO2+CH3OH and O2+CH3OH. Proc. Combust. Inst. 2011, 33, 351−357. (19) Jasper, A. W.; Klippenstein, S. J.; Harding, L. B.; Ruscic, B. Kinetics of the reaction of methyl radical with hydroxyl radical and methanol decomposition. J. Phys. Chem. A 2007, 111, 3932−3950. (20) Car, R.; Parrinello, M. The unified approach to density functional and molecular-dynamics in real space. Solid State Commun. 1987, 62, 403−405. (21) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; et al. PLUMED: A portable plugin for free-energy calculations with molecular dynamics. Comput. Phys. Commun. 2009, 180, 1961−1972. (22) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic-behavior. Phys. Rev. A 1988, 38, 3098−3100. (23) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785−789. (24) Troullier, N.; Martins, J. L. Efficient pseudopotentials for planewave calculations. Phys. Rev. B 1991, 43, 1993−2006. (25) Kleinman, L.; Bylander, D. M. Efficacious form for model pseudopotentials. Phys. Rev. Lett. 1982, 48, 1425−1428. (26) Nose, S. A unified formulation of the constant temperature molecular-dynamics methods. J. Chem. Phys. 1984, 81, 511−519. (27) Hoover, W. G. Canonical dynamicsEquilibrium phase-space distributions. Phys. Rev. A 1985, 31, 1695−1697. (28) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 12562−12566. (29) Balaban, A. T.; Ciubotariu, D.; Medeleanu, M. Topological indexes and real number vertex invariants based on graph eigenvalues or eigenvectors. J. Chem. Inf. Comput. Sci. 1991, 31, 517−523. (30) Randic, M. Unique numbering of atoms and unique codes for molecular graphs. J. Chem. Inf. Comput. Sci. 1975, 15, 105−108. (31) Bucko, T. Ab initio calculations of free-energy reaction barriers. J. Phys.: Condens. Matter 2008, 20. (32) Wang, H. On potential energy landscape and combustion chemistry modeling. Combust. Flame 2013, 160, 222−223. (33) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian09, Revision D.01; Gaussian: Wallingford, CT, USA, 2009. (34) Codorniu-Hernandez, E.; Kusalik, P. G. Mobility Mechanism of Hydroxyl Radicals in Aqueous Solution via Hydrogen Transfer. J. Am. Chem. Soc. 2012, 134, 532−538. (35) Codorniu-Hernandez, E.; Kusalik, P. G. Insights into the Solvation and Mobility of the Hydroxyl Radical in Aqueous Solution. J. Chem. Theory Comput. 2011, 7, 3725−3732. (36) Agarwal, V.; Dauenhauer, P. J.; Huber, G. W.; Auerbach, S. M. Ab Initio Dynamics of Cellulose Pyrolysis: Nascent Decomposition Pathways at 327 and 600 °C. J. Am. Chem. Soc. 2012, 134, 14958− 14972. (37) Pilling, M. J. From elementary reactions to evaluated chemical mechanisms for combustion models. Proc. Combust. Inst 2009, 32, 27− 44. (38) Broadbelt, L. J.; Stark, S. M.; Klein, M. T. Computer-generated pyrolysis modelingOn-the-fly generation of species, reactions, and rates. Ind. Eng. Chem. Res. 1994, 33, 790−799. (39) Green, W. H.; Allen, J. W.; Buesser, B. A.; Ashcraft, R. W.; Beran, G. J.; Class, C. A.; Gao, C.; Goldsmith, C. F.; Harper, M. R.; Jalan, A.; et al. RMGReaction Mechanism Generator, v4.0.1, http:// rmg.sourceforge.net/, 2013. (40) Branduardi, D.; Bussi, G.; Parrinello, M. Metadynamics with Adaptive Gaussians. J. Chem. Theory Comput. 2012, 8, 2247−2254.

AUTHOR INFORMATION

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

■

REFERENCES

We are thankful for financial support from AFOSR Grant FA9550-12-1-0133 monitored by Dr. Chiping Li. S.Z. greatly appreciates useful conversations with Dr. Fabio Pietrucci. S.Z. and J.P. are grateful to Ms. Stephanie Hoffmann for a careful reading of this manuscript. This research was supported in part by the National Science Foundation through XSEDE resources and the use of computational, storage, and networking infrastructure provided by the Hyak supercomputer system, supported in part by the University of Washington eScience Institute.

(1) Car, R.; Parrinello, M. Unified approach for molecular-dynamics and density-functional theory. Phys. Rev. Lett. 1985, 55, 2471−2474. (2) Laio, A.; Gervasio, F. L. Metadynamics: A method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Rep. Prog. Phys. 2008, 71, 126601. (3) Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 826−843. (4) Iannuzzi, M.; Laio, A.; Parrinello, M. Efficient exploration of reactive potential energy surfaces using Car-Parrinello molecular dynamics. Phys. Rev. Lett. 2003, 90, 238302. (5) Ensing, B.; De Vivo, M.; Liu, Z. W.; Moore, P.; Klein, M. L. Metadynamics as a tool for exploring free energy landscapes of chemical reactions. Acc. Chem. Res. 2006, 39, 73−81. (6) Gallet, G. A.; Pietrucci, F.; Andreoni, W. Bridging Static and Dynamical Descriptions of Chemical Reactions: An ab Initio Study of CO2 Interacting with Water Molecules. J. Chem. Theory Comput. 2012, 8, 4029−4039. (7) Stirling, A.; Iannuzzi, M.; Laio, A.; Parrinello, M. Azulene-tonaphthalene rearrangement: The Car-Parrinello metadynamics method explores various mechanisms. ChemPhysChem 2004, 5, 1558−1568. (8) Boero, M.; Ikeshoji, T.; Liew, C. C.; Terakura, K.; Parrinello, M. Hydrogen Bond Driven Chemical Reactions: Beckmann Rearrangement of Cyclohexanone Oxime into ε-Caprolactam in Supercritical Water. J. Am. Chem. Soc. 2004, 126, 6280−6286. (9) Branduardi, D.; Gervasio, F. L.; Parrinello, M. From A to B in free energy space. J. Chem. Phys. 2007, 126, 054103. (10) Tribello, G. A.; Ceriotti, M.; Parrinello, M. Using sketch-map coordinates to analyze and bias molecular dynamics simulations. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 5196−5201. (11) Tribello, G. A.; Ceriotti, M.; Parrinello, M. A self-learning algorithm for biased molecular dynamics. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 17509−17514. (12) Ceriotti, M.; Tribello, G. A.; Parrinello, M. Simplifying the representation of complex free-energy landscapes using sketch-map. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 13023−13028. (13) Pietrucci, F.; Andreoni, W. Graph Theory Meets ab Initio Molecular Dynamics: Atomic Structures and Transformations at the Nanoscale. Phys. Rev. Lett. 2011, 107, 085504. (14) Held, T. J.; Dryer, F. L. A comprehensive mechanism for methanol oxidation. Int. J. Chem. Kinet. 1998, 30, 805−830. (15) Li, J.; Zhao, Z. W.; Kazakov, A.; Chaos, M.; Dryer, F. L.; Scire, J. J. A comprehensive kinetic mechanism for CO, CH2O, and CH3OH combustion. Int. J. Chem. Kinet. 2007, 39, 109−136. (16) Skodje, R. T.; Tomlin, A. S.; Klippenstein, S. J.; Harding, L. B.; Davis, M. J. Theoretical Validation of Chemical Kinetic Mechanisms: Combustion of Methanol. J. Phys. Chem. A 2010, 114, 8286−8301. (17) Tsang, W. Chemical kinetic database for combustion chemistry. Part 2. Methanol. J. Phys. Chem. Ref. Data 1987, 16, 471−508. 10770

dx.doi.org/10.1021/jp500398k | J. Phys. Chem. C 2014, 118, 10764−10770