First-Principles Molecular Dynamics of Monomethylhydrazine and

Subscriber access provided by UNIV AUTONOMA DE COAHUILA UADEC

Chemical and Dynamical Processes in Solution; Polymers, Glasses, and Soft Matter

First-Principles Molecular Dynamics of Monomethylhydrazine and Nitrogen Dioxide Yulun Han, Erik K. Hobbie, and Dmitri S. Kilin J. Phys. Chem. Lett., Just Accepted Manuscript • DOI: 10.1021/acs.jpclett.9b00674 • Publication Date (Web): 12 Apr 2019 Downloaded from http://pubs.acs.org on April 14, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

First-Principles Molecular Dynamics of Monomethylhydrazine and Nitrogen Dioxide Yulun Han†, Erik K. Hobbie‡, and Dmitri S. Kilin†* †Department

of Chemistry and Biochemistry, North Dakota State University, Fargo, North Dakota 58102

‡Department

of Physics, North Dakota State University, Fargo, North Dakota 58102

ABSTRACT

The exploration of chemical reactions preceding ignition is essential for the development of ideal hypergolic propellants. Unexpected reaction pathways of a hypergolic mixture comprised of monomethylhydrazine and nitrogen dioxide are predicted through a cooperative combination of (i) spin-unrestricted ab initio molecular dynamics (AIMD) and (ii) wave packet dynamics of protons. Ensembles of AIMD trajectories reveal a sequence of reaction steps for proton transfer and rupture of C-N bond. The details of proton transfer are explored by wave packet dynamics on the basis of ab initio potential energy surfaces from AIMD trajectories. The possibility of spontaneous ignition of this hypergolic mixture at the room temperature is predicted as a quantized feature of proton transfer dynamics.

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 21

Hypergolic propellants have been widely applied in rockets and aircraft.1-2 Spontaneous ignition occurs when two components of a hypergolic propellant come into contact, without the requirement of an external ignition source.3 Traditional hypergolic propellants consist of hydrazine (N2H4) and its derivatives as fuels, e.g. monomethylhydrazine (MMH or CH3NHNH2), and dinitrogen tetroxide/nitrogen dioxide (N2O4/NO2) mixtures as oxidizers. However, questions have arisen with respect to the safety of these hypergolic propellants because of their corrosivity, toxicity, and carcinogenicity.4-5 As a result, there is considerable interest in finding alternatives to hydrazine and its derivatives as hypergolic fuels.6-17 To facilitate the development of ideal propellant fuels and to optimize combinations of fuels and oxidizers, it is important to understand the chemical reactions that precede ignition in the case of traditional hypergolic propellants. Experimental investigation of chemical reactions of hypergolic propellants is difficult because reactive transient intermediates are involved. On the other hand, theoretical approaches have been extensively adopted to study reactions between MMH and N2O4/NO2.18-23 Detailed chemical kinetic models have also been proposed.24-27 For example, Kanno et al. explored the kinetics and mechanisms of reactions of MMH and NO2 based on transition state theory (TST) and Rice–Ramsperger–Kassel–Marcus (RRKM) theory by assuming asymmetrical Eckart potential energy barriers.27-29 Nevertheless, the barrier height for the widely accepted initial reaction CH3NHNH2 + NO2 → CH3NNH2 + HONO was specified to be ~ 0.4 eV.22-23 It is far from clear why the mixture of MMH and N2O4/NO2 ignites hypergolically at the room temperature. Traditional computational chemistry rests on the classical path approximation (CPA) and the point charge approximation. The position of each nucleus is deterministic. However, a more detailed treatment is needed for chemical processes involving hydrogen atoms. One needs to


Page 3 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


quantize the degrees of freedom related to hydrogen motion because of light mass and uncertainty in position and momentum. The effect of quantum tunneling established in reaction rate theory plays a crucial role in hydrogen abstraction or proton transfer reactions.30-35 The balanced description of the quantum and classical subsystems and their interactions in molecular dynamics simulations is a delicate matter.34-38 This work seeks to bring the quantum dynamics of chemical reactions from the theory domain to practical computations. To describe the reaction of MMH and NO2, we explore the propagation of quantum wave packets on the BornOppenheimer potential energy surface (PES) by investigating the connection of CPA based ab initio molecular dynamics (AIMD) and the quantization of nuclear motion. Specifically, we use spin-unrestricted density functional theory (DFT)39-43 based AIMD to identify a sequence of reaction steps. Nuclear wave packet dynamics on the basis of ab initio PES provides insight into proton tunneling through the barrier. This work helps resolve the controversy between low temperature reaction in experiment and high barrier height in computation. A deeper understanding of the reaction mechanism of this simple system is of critical importance to shedding light on the reaction pathways of more complex systems. The reported work involves two stages: (i) modeling the sequence of elementary reactions using CPA based AIMD and (ii) refining reaction rates using nuclear wave packet dynamics. Stage (i) enables the exploration of possible scenarios in the form of an atomistic sequence of reaction events. However, those reactions involving proton motion are expected to occur at higher rates than that predicted by AIMD. Stage (ii) enables the modeling of quantum dynamics associated with proton-tunneling-induced reactant-to-product conversion. i. Ab initio molecular dynamics. AIMD calculations are performed at 2000 K for 1 ps using a time step of 1 fs. It should be noted that the AIMD trajectory generated species with excessive



Page 4 of 21

kinetic energies. To accurately assess the energy evolution, we apply a post-processing technique after the completion of AIMD calculations.44-47 The post-processing technique involves three steps. At each instant of the AIMD trajectory, we (a) extract all species from the system, (b) remove species not engaging in hydrogen transfer, and (c) perform single point energy calculations for the rest of species. Here, the removal of spectators in step (b) is a practical implementation to ensure that thermal bond elongation or contraction events in these species do not contribute to the energy of the system of interests. We also avoid image effects arising from periodic boundary conditions by adopting this post-processing technique. Note that the original simulation cell is densely packed with reactant molecules, which are polarized and serve as explicit solvent species rather than a polarized continuum. At the stage of AIMD the environment is treated in an explicit way. Prevailing interactions between reactant molecules are characterized by collision rather than continuous exposure. At a chosen instant of time, the specific reactant molecules MMH and NO2 collide, while other molecules experience free translational motion. A slice of the potential energy surface 𝑉′(𝑡) along such trajectory allows to determine the barrier height for the first reaction step as (1a)

𝐸𝑎 = max (𝑉′(𝑡)) ―𝑉′(𝑡 = 0)

ii. Wave packet dynamics. The details of wave packet dynamics described in the supporting information explain how we compute the kinetics of product concentration Pproduct(𝑡). A quantitative assessment of the reaction can be performed based on reaction yield defined as the final probability of product formation, 𝑌 = Pproduct(∞). The rate constant 𝑘 can be obtained by fitting the time dependence to an exponential function, Pproduct(𝑡 ― 𝑡′) = Pproduct(∞)[1 ― 𝑒 ―𝑘(𝑡 ― 𝑡 )] ′

𝑘=

Pproduct(∞) ∞ ∫ 𝑡′𝑑𝑡

[Pproduct(𝑡 ― 𝑡′) ― Pproduct(∞)]

(1b) (1c) 4

ACS Paragon Plus Environment

Page 5 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where 𝑡′ represents the time at which the wave packet starts climbing the barrier. The wave packet approaches the barrier when 0 < 𝑡 < 𝑡′. The numerical experiments are performed for various initial temperatures corresponding to different initial momenta 𝑃0(𝑇). The computed reaction rate constants are processed based on a modified Arrhenius equation 𝐵

ln (𝑘) = 𝐴′ ― 𝑇 + 𝐶 ln (𝑇)

(2a)

where 𝐴′, 𝐵, and 𝐶 are parameters obtained by fitting the rate constants, 𝑘(𝑇).48-51 In terms of the modified Arrhenius equation, the activation energy at a given temperature is 𝐸𝑎(𝑇) = 𝑅𝐵 + 𝑅𝐶𝑇

(2b)

where 𝑅 is the gas constant, 𝐵 and 𝑇 are in units of K, and 𝐶 is a constant. The main technical goal of the paper is to compare how much the nominal barrier height 𝐸𝑎 differs from 𝐸𝑎(𝑇), where 𝐸𝑎 is defined as the energy difference between reactant and least stable intermediate along the reaction path in Eq. (1a). We aim to prove the hypothesis that at low temperatures 𝐸𝑎(𝑇) < 𝐸𝑎 due to quantum tunneling through the barrier. iii. Computational details. Quantum mechanical calculations are carried out using spinunrestricted DFT on the basis of Kohn-Sham (KS) orbitals39 within the Vienna Ab initio Simulation Package (VASP).40-43 The generalized gradient approximation (GGA)52-53 for the exchange-correlation Perdew-Burke-Ernzerhof (PBE)53 functional is adopted. Projected augmented wave (PAW)54 potentials in a plane-wave basis are used. All atoms are treated as point charges. The atomic model consisting of one MMH molecule and eight NO2 molecules is placed into a cubic unit cell with an edge length of 10 Å under periodic boundary conditions. The minimum distance between MMH and NO2 is about 1.94 Å. Prior to AIMD calculations, the



Page 6 of 21

geometry optimized model is heated at 2000 K for 300 fs. The AIMD results are visualized using VESTA software.55 Along one AIMD trajectory, we observe three sequential steps of hydrogen transfer leading to CH3NN and 3 HONO. Figure 1 shows the energy diagram of the initial hydrogen transfer reaction. The least stable intermediate is found for the configuration at 240 fs, where one H is ejected from MMH. From 240 to 253 fs it is observed that the distance between the ejected H and donor N increases from 1.30 to 1.50 Å, while the distance between the ejected H and acceptor O decreases from 1.56 to 1.12 Å. The barrier height (𝐸𝑎) is 0.55 eV, i.e. the energy difference between the reactant and the least stable intermediate. As a reference, this amount of energy corresponds to a temperature of 6400 K. Note that AIMD calculations are performed at 2000 K. It is able to overcome the barrier via local fluctuation and redistribution of energy between degrees of freedom in noncollinear system. The energy diagrams related to other hydrogen transfer reactions are shown in supporting information. It is found that all three hydrogen transfer reactions are exothermic. Figure 2a shows the number of N2, NO, NO2, and HONO as a function of time during one AIMD trajectory. At ~ 250, ~ 580, and ~ 620 fs, one observes a stepwise decrease in the number of NO2 molecules due to three sequential hydrogen transfer reactions. A spike at ~ 630 fs is found for the curve corresponding to NO2, indicating that the amount of NO2 increases and then drops in a short period of time. This is because HONO produced after the third hydrogen transfer reaction is not stable. The transferred H is ejected and later re-captured by NO2. After ~ 630 fs, the number of NO2 molecules maintains at 5. HO-NO bond breaking and bond reforming events are abundant for HONO species, which explains the swapping behavior of curves corresponding to HONO and NO. At ~ 780 fs, one observes an increase in the number of N2 due to the rupture of C-N bond in


Page 7 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


CH3NN. The generated CH3 species then combines with OH to form CH3OH during the final stages of the trajectory. Figure 2b shows the total number of species presented in the AIMD trajectory. The dominant feature corresponds to NO2. Features corresponding to HONO, NO, and HO are relatively strong as well. The reaction of CH3 and OH takes place during the last stages of the trajectory, and thus the feature corresponding to CH3OH is weak. In this model with a small simulation cell, we express the results in terms of the absolute number of species. These results could be extrapolated to experimentally reasonable observables by recasting realistic concentrations of products and intermediates. To get more insight into hydrogen transfer reactions, we conduct Bader charge analysis5658

for configurations along the reaction pathway, as detailed in the supporting information. It is

found that hydrogen transfer reactions are indeed proton transfer reactions. Another finding is that neutral HONO is formed upon the completion of proton transfer. We have hitherto neglected the quantum nature of the proton. Based on AIMD the barrier height is 0.55 eV for the first proton transfer reaction i.e. the rate-controlling step in the initiation of the hypergolic system. It remains a puzzle how MMH and NO2 react at the room temperature to overcome such a high barrier. Here, we present a hypothesis to explain the seemingly paradox of high barrier but low activation temperature. According to this hypothesis, quantization of proton degree of freedom is expected to facilitate quantum tunneling through the barrier 𝐸𝑎. The modeled reactions proceed through the collision of reactants. We focus on the first proton transfer reaction and make a quantum correction to the translational motion of the proton during the collision. The product of the considered reaction is an intermediate for a full sequence of reactions. Figure 3a-3c show representative wave packet dynamics at 300 K computed by equations S1-S4 on the basis of first proton transfer PES. Specifically, snapshots of the wave



Page 8 of 21

packet dynamics P(𝑋,𝑡) are presented in Figure 3a. At the initial time the Gaussian wave packet is set at 𝑋0 = ― 1.0 atomic unit (au) of length with nuclear momentum 𝑃0. The position of potential curve ranges from 0.0 to 1.3 au of length. With the increase of time, we observe the evolution of the nuclear wave packet toward the potential curve. Figure 3b illustrates wave packet probability distribution dynamics as a function of time and position. The splitting of wave packet is observed at ~ 240 au of time. The splitting indicates part of wave packet passes through the potential curve, while the other reflects back, contributing to formation of products and keeping the survived reactant. Such splitting is a signature of quantum superposition effect, which is not possible in classical and semi-classical approaches. Figure 3c illustrates probability of product formation versus time. The proton transfer yield Pproduct(∞) is calculated as 0.18. The rate constants 𝑘 and yields Pproduct(∞) for other temperatures are shown in supporting information. It is found that 𝑘 and Pproduct(∞) increase with increasing temperatures. In Figure 3d, we show the plot of logarithm of rate constants against inverse of temperatures and the fitting by the modified Arrhenius equation using equation 2, with factors 𝐴′ = ― 3.65, 𝐵 = ― 34.09 and 𝐶 = 0.29. The data in Figure 3d fit equation 2 well. It is found that the computed activation energy 𝐸𝑎 (𝑇) according to the modified Arrhenius equation ranges from ~ 0.01 eV at 300 K to ~ 0.05 eV at 2000 K. In this range, 𝐸𝑎(𝑇) is much smaller than 𝐸𝑎. Note 𝐸𝑎(𝑇) is interpreted as the activation energy for proton tunneling. A conceptual explanation of the hypergolic reaction at low temperatures is provided as follows. At 300 K, for example, the kinetic energy of the system is ~ 0.03 eV. This amount is not sufficient to overcome 𝐸𝑎. However, it is greater than 𝐸𝑎(𝑇) such that the proton can be detected on the far side of the barrier through tunneling. In an ensemble of reactants, the first successful exothermic reaction will generate and dissipate enough kinetic energy into nearby molecules to trigger a chain reaction. 8 ACS Paragon Plus Environment

Page 9 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Pproduct(∞) is a technical observable from numerical experiments. The experimental equilibrium value would be obtained by multiple pass summation, away from the goal of this work. In realistic combustion chamber, the product gases experience expansion. Here the fixed equilibrium geometry of products, including their intermolecular distances, is also an approximation. Pproduct(∞) ranges from 0.18 at 300 K to 0.34 at 2000 K. Although Pproduct(∞) is qualitatively different from 0, it is still far below 1. This indicates only a small fraction has tunneled through the barrier, while the major fraction has reflected. The low yields Pproduct(∞) are due to an artifact in the modeling, in which we explore a first-order correction originating from a single pass of the wave packet through the barrier. A one-pass scattering model is chosen intentionally such that molecules separate after the collision. Quick proton transfer occurs in addition to relatively slow collision and mutual scattering of reactants. During the considered reaction, the overall distance between MMH and NO2 is approximately constant. On a longer time scale, they approach before the proton transfer reaction and depart from each other after the reaction. By performing multiple-pass modeling, one would expect the yields are close to 1. An alternative to multiple-pass modeling is to replace the flat potential by a repulsive one. The reflected fraction would experience repulsion from nitrogen and will turn around for performing multiple passes to complete the reaction. On the other hand, we should bear in mind that quasi-one-dimensional tunneling models tend to significantly overestimate tunneling effects. Tunneling correction through a quasi-one-dimensional barrier can be performed by perturbation theory taking 𝜓reactant and 𝜓product wavefunction as plane waves and express matrix element of potential as 〈𝑉〉 = ∫𝜓reactant𝑉(𝑥)𝜓productd𝑥 and rate 𝑘~|〈𝑉〉|2𝜌(𝐸) by the golden rule. However, the wave packet dynamics is more general. It would be a trivial calculation for a standardized form of the potential



Page 10 of 21

such as the Eckert potential. For a non-analytical and numerical form of the potential, however, this problem demands for numerical solution. We also calculate independent AIMD trajectories starting from different initial nuclear configurations and corresponding wave packet dynamics to assess the level of correction originating from average over ensemble of trajectories (see supporting information). Independent trajectories provide qualitatively similar potentials with some variation in barrier heights. One of the observables to represent influence of the ensemble average is the cumulative reaction yield from all elements of an ensemble. As a reference point, the yield is compared to the cumulative reaction yield from classical mechanics analysis. In classical analysis, a reaction yield takes either 0 if the initial kinetic energy of reactant is less than the barrier height or 1 if the initial kinetic energy of reactant is greater than or equal to the barrier height. Both cumulative reaction yields from classical mechanics and wave packet dynamics increase with increasing temperature. The cumulative reaction yields from classical mechanics are zero for temperatures (kinetic energies) below the lowest barrier height of the ensemble and increase in a stepwise fashion with the increase of temperature. The cumulative reaction yields from wave packet dynamics are nonzero for temperatures (kinetic energies) below the lowest barrier height of the ensemble. The nonzero yields are due to the quantum nature of proton. One expects that an average over larger ensemble will make classical curve smoother and both classical and quantum curves may converge in the limit of infinite temperature. In summary, the reported work is important from practical and methodological aspects. From a practical aspect: we have demonstrated AIMD can be used to study chemical reactions of MMH with NO2. The wave packet dynamics based on ab initio PES demonstrates that proton tunneling occurs, reducing the minimal thermal energy required to activate such reactions. Proton


Page 11 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tunneling is important for reactions between MMH and NO2 due to the low barrier height, short reaction path, and exothermic process. The finding of nonzero proton tunneling probability in the initial reaction step towards MMH combustion addresses a main challenge of contemporary rocket fuel science. It provides a conceptual explanation of the hypergolic reaction at the room temperature despite predicted high barrier height. From a methodological aspect: the analysis of quantum effect rests on parameters and technical details extracted from first-principles calculations. Most AIMD calculations treat nuclei classically. The presence of classically forbidden regions excludes several reaction pathways from modeling even if they occur experimentally due to quantum tunneling. We perform AIMD at high temperatures to identify the most probable sequence of reaction pathways. Quantum corrections are then applied to evaluate tunneling reactions at low temperatures. This work provides approximation enabling quantum corrections to CPA dynamics by the choice of preferred reaction coordinate. It demonstrates the qualitative correctness of coupled classical and quantum treatments of reaction trajectory.

AUTHOR INFORMATION Corresponding Author *E-mail: [email protected] Notes The authors declare no competing financial interests. Funding Sources This research has been supported by NSF CHE-1413614 and DOE DE-SC0001717 for methods development. ACKNOWLEDGMENT



Page 12 of 21

Authors thank DOE BES NERSC facility for computational resources, allocation award #91202, “Computational Modeling of Photo-catalysis and Photo-induced Charge Transfer Dynamics on Surfaces” supported by the Office of Science of the DOE under contract no. DE-AC0205CH11231. DSK thanks David Micha, Sergei Tretiak, Oleg Prezhdo, Svetlana Kilina, and John Hershberger for discussions. Supporting Information Available: Details about AIMD and nuclear wave packet dynamics, Bader charge analysis, preparation of PES, additional data for AIMD and wave packet dynamics.


Page 13 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


FIGURES

Figure 1. Energy diagram from one AIMD trajectory at 2000 K for the reaction of CH3NHNH2 + NO2 → CH3NNH2 + HONO. Atomistic snapshots of partially filled circles are indicated. The brown, pink, blue, and red spheres represent C, H, N, and O, respectively. Interatomic distances in Å are given for select bonds.



Page 14 of 21

Figure 2. (a) Number of species as a function of time in one AIMD trajectory at 2000 K. (b) Total number of species in the AIMD trajectory at 2000 K. The histograms are expanded to show species with low frequency.


Page 15 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. (a-c) Representative wave packet dynamics on the basis of first proton transfer PES from one AIMD trajectory at 2000 K. (a) Snapshots of wave packet dynamics. The red, green, and blue solid lines represent wave packet in the initial state P(𝑋,𝑡 = 0), the final state P(𝑋,𝑡 = 𝑡end), and potential curve 𝑉(𝑋) according to equation S1a, respectively. The offset 〈𝐸〉 represents the total energy of the wave packet. The red, yellow, and green background bands correspond to 𝑉′(𝑋ini) for 𝑋begin < 𝑋 ≤ 𝑋ini, 𝑉′(𝑋) for 𝑋ini < 𝑋 < 𝑋fin, and 𝑉′(𝑋fin) for 𝑋fin ≤ 𝑋 < 𝑋end, respectively. (b) Wave packet dynamics as a function of position and time according to equation S4a. (c) Probability of product formation versus time according to equations S4b-c. Panels (a-c) correspond to the temperature at 300 K. (d) Calculated rate constants versus temperatures and the fitting by a modified Arrhenius equation using equation 2 with 𝐴′ = ― 3.65, 𝐵 = ― 34.09 and 𝐶 = 0.29. The plot is expanded to show data points ranging from 300 to 2000 K.



Page 16 of 21

REFERENCES (1) Hurlbert, E.; Applewhite, J.; Nguyen, T.; Reed, B.; Baojiong, Z.; Yue, W. Nontoxic Orbital Maneuvering and Reaction Control Systems for Reusable Spacecraft. J. Propul. Power 1998, 14, 676-687. (2) Catoire, L.; Chaumeix, N.; Pichon, S.; Paillard, C. Visualizations of Gas-Phase NTO/MMH Reactivity. J. Propul. Power 2006, 22, 120-126. (3) Eiceman, G. A.; Salazar, M. R.; Rodriguez, M. R.; Limero, T. F.; Beck, S. W.; Cross, J. H.; Young, R.; James, J. T. Ion Mobility Spectrometry of Hydrazine, Monomethylhydrazine, and Ammonia in Air with 5Nonanone Reagent Gas. Anal. Chem. 1993, 65, 1696-1702. (4) Wang, S.; Thynell, S. T.; Chowdhury, A. Experimental Study on Hypergolic Interaction between N,N,N′,N′-Tetramethylethylenediamine and Nitric Acid. Energy Fuels 2010, 24, 5320-5330. (5) Pichon, S.; Catoire, L.; Chaumeix, N.; Paillard, C. Search for Green Hypergolic Propellants: Gas-Phase Ethanol/Nitrogen Tetroxide Reactivity. J. Propul. Power 2005, 21, 1057-1061. (6) Zhang, Y.; Gao, H.; Joo, Y.-H.; Shreeve, J. n. M. Ionic Liquids as Hypergolic Fuels. Angew. Chem., Int. Ed. 2011, 50, 9554-9562. (7) Chambreau, S. D.; Schneider, S.; Rosander, M.; Hawkins, T.; Gallegos, C. J.; Pastewait, M. F.; Vaghjiani, G. L. Fourier Transform Infrared Studies in Hypergolic Ignition of Ionic Liquids. J. Phys. Chem. A 2008, 112, 7816-7824. (8) Zhang, X.; Shen, L.; Luo, Y.; Jiang, R.; Sun, H.; Liu, J.; Fang, T.; Fan, H.; Liu, Z. Synthesis and Ignition Properties Research of 1,5-Diazabicyclo[3.1.0]Hexane Type Compounds as Potential Green Hypergolic Propellants. Ind. Eng. Chem. Res. 2017, 56, 2883-2888. (9) Chambreau, S. D.; Koh, C. J.; Popolan-Vaida, D. M.; Gallegos, C. J.; Hooper, J. B.; Bedrov, D.; Vaghjiani, G. L.; Leone, S. R. Flow-Tube Investigations of Hypergolic Reactions of a Dicyanamide Ionic Liquid Via Tunable Vacuum Ultraviolet Aerosol Mass Spectrometry. J. Phys. Chem. A 2016, 120, 80118023. (10) Carlin, C. M.; Gordon, M. S. Ab Initio Investigation of Cation Proton Affinity and Proton Transfer Energy for Energetic Ionic Liquids. J. Phys. Chem. A 2016, 120, 6059-6063. (11) McCrary, P. D.; Barber, P. S.; Kelley, S. P.; Rogers, R. D. Nonaborane and Decaborane Cluster Anions Can Enhance the Ignition Delay in Hypergolic Ionic Liquids and Induce Hypergolicity in Molecular Solvents. Inorg. Chem. 2014, 53, 4770-4776. (12) McCrary, P. D.; Chatel, G.; Alaniz, S. A.; Cojocaru, O. A.; Beasley, P. A.; Flores, L. A.; Kelley, S. P.; Barber, P. S.; Rogers, R. D. Evaluating Ionic Liquids as Hypergolic Fuels: Exploring Reactivity from Molecular Structure. Energy Fuels 2014, 28, 3460-3473.


Page 17 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(13) Russo, M. F.; Bedrov, D.; Singhai, S.; van Duin, A. C. T. Combustion of 1,5-Dinitrobiuret (DNB) in the Presence of Nitric Acid Using ReaxFF Molecular Dynamics Simulations. J. Phys. Chem. A 2013, 117, 9216-9223. (14) Sun, R.; Siebert, M. R.; Xu, L.; Chambreau, S. D.; Vaghjiani, G. L.; Lischka, H.; Liu, J.; Hase, W. L. Direct Dynamics Simulation of the Activation and Dissociation of 1,5-Dinitrobiuret (HDNB). J. Phys. Chem. A 2014, 118, 2228-2236. (15) Jiao, N.; Zhang, Y.; Liu, L.; Shreeve, J. M.; Zhang, S. Hypergolic Fuels Based on Water-Stable Borohydride Cluster Anions with Ultralow Ignition Delay Times. J. Mater. Chem. A 2017, 5, 13341-13346. (16) Li, X.; Wang, C.; Li, H.; Nie, F.; Yin, H.; Chen, F.-X. Bishydrobis(tetrazol-1-yl)Borate (BTB) Based Energetic Ionic Liquids with High Density and Energy Capacity as Hypergolic Fuels. J. Mater. Chem. A 2017, 5, 15525-15528. (17) Dane, H.; Yulun, H.; Dmitri, K. A Computational Study of the Combustion of Hydrazine with Dinitrogen Tetroxide. Journal of Nanotoxicology and Nanomedicine (JNN) 2017, 2, 12-30. (18) Frank, I.; Hammerl, A.; Klapötke, T. M.; Nonnenberg, C.; Zewen, H. Processes During the Hypergolic Ignition between Monomethylhydrazine (MMH) and Dinitrogen Tetroxide (N2O4) in Rocket Engines. Propellants, Explos., Pyrotech. 2005, 30, 44-52. (19) Nonnenberg, C.; Frank, I. Formation and Decay of Tetrazane Derivatives-a Car-Parrinello Molecular Dynamics Study. Phys. Chem. Chem. Phys. 2008, 10, 4383-4392. (20) Nonnenberg, C.; Frank, I.; Klapötke, T. M. Ultrafast Cold Reactions in the Bipropellant Monomethylhydrazine/Nitrogen Tetroxide: CPMD Simulations. Angew. Chem., Int. Ed. 2004, 43, 45864589. (21) Liu, Y.; Zybin, S. V.; Guo, J.; van Duin, A. C. T.; Goddard, W. A. Reactive Dynamics Study of Hypergolic Bipropellants: Monomethylhydrazine and Dinitrogen Tetroxide. J. Phys. Chem. B 2012, 116, 14136-14145. (22) Liu, W.-G.; Wang, S.; Dasgupta, S.; Thynell, S. T.; Goddard, W. A.; Zybin, S.; Yetter, R. A. Experimental and Quantum Mechanics Investigations of Early Reactions of Monomethylhydrazine with Mixtures of NO2 and N2O4. Combust. Flame 2013, 160, 970-981. (23) McQuaid, M. J.; Ishikawa, Y. H-Atom Abstraction from CH3NHNH2 by NO2: CCSD(T)/6311++G(3df,2p)//MPWB1K/6-31+G(D,P)

and

CCSD(T)/6-311+G(2df,P)//CCSD/6-31+G(D,P)

Calculations. J. Phys. Chem. A 2006, 110, 6129-6138. (24) Catoire, L.; Chaumeix, N.; Paillard, C. Chemical Kinetic Model for Monomethylhydrazine/Nitrogen Tetroxide Gas Phase Combustion and Hypergolic Ignition. J. Propul. Power 2004, 20, 87-92. (25) Agosta, V. D.; Seamans, T. F.; Vanpee, M. Development of a Fundamental Model of Hypergolic Ignition in Space-Ambient Engines. AIAA Journal 1967, 5, 1616-1624. 17 ACS Paragon Plus Environment


Page 18 of 21

(26) Kanno, N.; Tani, H.; Daimon, Y.; Terashima, H.; Yoshikawa, N.; Koshi, M. Computational Study of the Rate Coefficients for the Reactions of NO2 with CH3NHNH, CH3NNH2, and CH2NHNH2. J. Phys. Chem. A 2015, 119, 7659-7667. (27) Kanno, N.; Terashima, H.; Daimon, Y. U.; Yoshikawa, N.; Koshi, M. Theoretical Study of the Rate Coefficients for CH3NHNH2 + NO2 and Related Reactions. Int. J. Chem. Kinet. 2014, 46, 489-499. (28) Garrett, B. C.; Truhlar, D. G. Semiclassical Tunneling Calculations. J. Phys. Chem. 1979, 83, 29212926. (29) Eckart, C. The Penetration of a Potential Barrier by Electrons. Phys. Rev. 1930, 35, 1303-1309. (30) Caldin, E. F. Tunneling in Proton-Transfer Reactions in Solution. Chem. Rev. 1969, 69, 135-156. (31) Hu, W.-P.; Liu, Y.-P.; Truhlar, D. G. Variational Transition-State Theory and Semiclassical Tunnelling Calculations with Interpolated Corrections: A New Approach to Interfacing Electronic Structure Theory and Dynamics for Organic Reactions. J. Chem. Soc., Faraday Trans. 1994, 90, 1715-1725. (32) Truhlar, D. G.; Garrett, B. C.; Klippenstein, S. J. Current Status of Transition-State Theory. J. Phys. Chem. 1996, 100, 12771-12800. (33) Som, S.; Liu, W.; Zhou, D. D. Y.; Magnotti, G. M.; Sivaramakrishnan, R.; Longman, D. E.; Skodje, R. T.; Davis, M. J. Quantum Tunneling Affects Engine Performance. J. Phys. Chem. Lett. 2013, 4, 20212025. (34) Micha, D. A.; Stodden, C. D. An Eikonal Treatment of Electronically Diabatic Photodissociation: Branching Ratios of CH3I. J. Phys. Chem. A 2001, 105, 2890-2896. (35) Cohen, J. M.; Micha, D. A. Angular Distributions in Electronically Adiabatic Hyperthermal Collisions. An Eikonal Approach. J. Chem. Phys. 1993, 98, 2023-2031. (36) Cohen, J. M.; Micha, D. A. Electronically Diabatic Atom–Atom Collisions: A Self‐Consistent Eikonal Approximation. J. Chem. Phys. 1992, 97, 1038-1052. (37) Micha, D. A. A Self‐Consistent Eikonal Treatment of Electronic Transitions in Molecular Collisions. J. Chem. Phys. 1983, 78, 7138-7145. (38) Micha, D. A.; Thorndyke, B. The Quantum–Classical Density Operator for Electronically Excited Molecular Systems. Adv. Quantum Chem. 2004, 47, 293-314. (39) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133-A1138. (40) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B 1993, 47, 558561. (41) Kresse, G.; Hafner, J. Ab Initio Molecular-Dynamics Simulation of the Liquid-Metal–AmorphousSemiconductor Transition in Germanium. Phys. Rev. B 1994, 49, 14251-14269.


Page 19 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(42) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169-11186. (43) Kresse, G.; Furthmüller, J. Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6, 15-50. (44) Han, Y.; Anderson, K.; Hobbie, E. K.; Boudjouk, P.; Kilin, D. S. Unraveling Photodimerization of Cyclohexasilane from Molecular Dynamics Studies. J. Phys. Chem. Lett. 2018, 9, 4349-4354. (45) Han, Y.; Kilin, D. S.; May, P. S.; Berry, M. T.; Meng, Q. Photofragmentation Pathways for Gas-Phase Lanthanide Tris(isopropylcyclopentadienyl) Complexes. Organometallics 2016, 35, 3461-3473. (46) Han, Y.; Meng, Q.; Rasulev, B.; May, P. S.; Berry, M. T.; Kilin, D. S. Photoinduced Charge Transfer Versus Fragmentation Pathways in Lanthanum Cyclopentadienyl Complexes. J. Chem. Theory Comput. 2017, 13, 4281-4296. (47) Han, Y.; Rasulev, B.; Kilin, D. S. Photofragmentation of Tetranitromethane: Spin-Unrestricted TimeDependent Excited-State Molecular Dynamics. J. Phys. Chem. Lett. 2017, 8, 3185-3192. (48) Arrhenius, S. Über Die Reaktionsgeschwindigkeit Bei Der Inversion Von Rohrzucker Durch Säuren. Z. Phys. Chem. 1889, 4, 226-248. (49) Kooij, D. M. Über Die Zersetzung Des Gasförmigen Phosphorwasserstoffs. Z. Phys. Chem. 1893, 12, 155-161. (50) Trautz, M. Der Temperaturkoeffizient Chemischer Reaktionsgeschwindigkeiten. II. Die Physikalische Bedeutung Der Chemischen Reaktionsgeschwindigkeit in Gasen Und Ihre Vorausberechnung Ans Rein Thermischen Daten Der Beteiligten Stoffe. Z. Phys. Chem. 1909, 66, 496-511. (51) Laidler, K. J., The Development of the Arrhenius Equation. J. Chem. Educ. 1984, 61, 494-498. (52) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Atoms, Molecules, Solids, and Surfaces: Applications of the Generalized Gradient Approximation for Exchange and Correlation. Phys. Rev. B 1992, 46, 6671-6687. (53) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865-3868. (54) Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953-17979. (55) Momma, K.; Izumi, F. Vesta 3 for Three-Dimensional Visualization of Crystal, Volumetric and Morphology Data. J. Appl. Crystallogr. 2011, 44, 1272-1276. (56) Tang, W.; Sanville, E.; Henkelman, G. A Grid-Based Bader Analysis Algorithm without Lattice Bias. J. Phys.: Condens. Matter 2009, 21, 084204. (57) Sanville, E.; Kenny, S. D.; Smith, R.; Henkelman, G. Improved Grid-Based Algorithm for Bader Charge Allocation. J. Comput. Chem. 2007, 28, 899-908.



Page 20 of 21

(58) Henkelman, G.; Arnaldsson, A.; Jónsson, H. A Fast and Robust Algorithm for Bader Decomposition of Charge Density. Comput. Mater. Sci. 2006, 36, 354-360.


Page 21 of 21 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC


First-Principles Molecular Dynamics of Monomethylhydrazine and

Recommend Documents