ReaxFF Reactive Force Field Development and Applications for

Apr 12, 2010 - Ammonia borane (AB) has attracted significant attention due to its high hydrogen content (19.6% by mass). To investigate the reaction m...
1 downloads 17 Views 3MB Size
Download PDF
J. Phys. Chem. A 2010, 114, 5485–5492

5485

ReaxFF Reactive Force Field Development and Applications for Molecular Dynamics Simulations of Ammonia Borane Dehydrogenation and Combustion Michael R. Weismiller, Adri C. T. van Duin,* Jongguen Lee, and Richard A. Yetter Department of Mechanical and Nuclear Engineering, PennsylVania State UniVersity, UniVersity Park PennsylVania 16802 ReceiVed: January 6, 2010; ReVised Manuscript ReceiVed: March 21, 2010

Ammonia borane (AB) has attracted significant attention due to its high hydrogen content (19.6% by mass). To investigate the reaction mechanism associated with the combustion of AB, a reactive force field (ReaxFF) has been developed for use in molecular dynamics (MD) simulations. The ReaxFF parameters have been derived directly from quantum mechanical data (QM). NVT-MD simulations of single- and polymolecular AB thermolysis were conducted in order to validate the force field. The release of the first equivalent H2 is a unimolecular reaction, and MD simulations show an activation energy of 26.36 kcal mol-1, which is in good agreement with experimental results. The release of the second H2 is also a unimolecular reaction; however, the release of a third H2 requires the formation of a B-N polymer. Similar simulations were conducted with a boron and oxygen system, since the oxidation of boron will be an integral step in AB combustion, and show good agreement with the established mechanism for this system. At low temperatures, boron atoms agglomerate into a cluster, which is oxidized at higher temperatures, eventually forming condensed and gas phase boron-oxide-species. These MD results provide confidence that ReaxFF can properly model the oxidation of AB and provide mechanistic insight into the AB dehydrogation and combustion reactions. 1. Introduction Due to political and environmental reasons, the prospect of a “hydrogen economy” has been discussed and sought after for years. Hydrogen is widely accepted as a “clean” fuel, since it does not release any CO2, which has been widely attributed to climate change. Furthermore, the reaction of H2 + 1/2O2 f H2O has a combustion enthalpy of 141.9 kJ g-1, which gives it the highest heating value per unit mass of any naturally occurring molecule.2 This makes hydrogen particularly attractive as a fuel for space propulsion, since the weight of propellant is critical to an effective rocket-propelled spacecraft. Nonetheless, there exist large barriers to the widespread use of hydrogen as a fuel. Perhaps the largest barrier to a hydrogenbased economy is the ability to safely store quantities of hydrogen that are appropriate for the various applications in which it would be required. At standard temperatures and pressures, hydrogen is a very low-density gas (0.08988 g L-1), meaning that any practical application would require that gaseous H2 be stored at high pressures. The use of such highpressure gas presents a safety risk and requires the use of heavywalled steel or composite tanks, which are both heavy and expensive. Although cryogenic liquid hydrogen is currently used in some rocket propulsion applications (e.g., Space Shuttle main engine), cryogenic storage is expensive. Additionally, future manned lunar and Mars missions may require that propellants be stored for long periods in low-earth orbit, making cryogenic storage unattractive. One promising solution to these issues is the solid-state chemical storage of hydrogen, using materials such as metal hydrides, since they exhibit both a high gravimetric and volumetric density when compared with compressed or cryogenic hydrogen.4 Figure 1 (with data from Zuttel1-3) shows the * To whom correspondence should be addressed. E-mail: acv13@ psu.edu.

Figure 1. Volumetric and gravimetric hydrogen density of some selected hydrides and other hydrogen storage methods (assembled from Zuttel1-3).

volumetric hydrogen density versus gravimetric hydrogen density of various substitute forms of storing hydrogen, including various metal hydrides, liquid hydrocarbons, and hydrogen adsorbed onto carbon nanotubes. Also shown in Figure 1 is hydrogen storage in compressed and liquid forms for various pressures and storage container materials. The volumetric density is defined as the amount of hydrogen produced in mass relative to the storage volume. The gravimetric density is defined as the mass of hydrogen produced relative to the storage mass. As can be seen from Figure 1, ammonia borane (NH3BH3) exhibits high values for both of these quantities. This has prompted many recent studies into ammonia borane (AB) as a regenerative onboard hydrogen storage medium for transportation.5,6

10.1021/jp100136c  2010 American Chemical Society Published on Web 04/12/2010

5486

J. Phys. Chem. A, Vol. 114, No. 17, 2010

Previous AB studies have focused on its low temperature and pressure thermolysis,7,8 low temperature thermolysis at elevated pressures,9 and catalytic10-17 and noncatalytic hydrolysis,18 with the objective being hydrogen release at relatively low rates, appropriate for its use with proton exchange membrane fuel cells. However, little study has been made of the decomposition and oxidation of AB at high temperatures and pressures. Such study is important to the evaluation of AB as a potential fuel or fuel additive for chemical rocket propulsion systems. Because of its high hydrogen content (19.6% by weight), AB has the potential to significantly increase the specific impulse (impulse per weight of the propellant) when used as a fuel in either a composite solid propellant rocket system (solid fuel and solid oxidizer suspended in a binder), or hybrid rocket system (solid fuel combusted with liquid or gaseous oxidizer), while providing easier storage than cryogenic liquid propellant rocket engines. Since little is known about the AB oxidation process, an atomistic description of the system at relatively long time scales is highly desirable. The computational cost for ab initio quantum mechanics (QM) simulations of moderately sized systems (50-100 atoms) at long time scales are far too high. For this reason, a ReaxFF reactive force field, which retains nearly the accuracy of QM but significantly reduces the computational cost, has been developed to simulate this process. Classic force fields (FF) allow for the molecular modeling of large systems, yet require the explicit definition of all bonds, so they are unable to model chemical reactions. ReaxFF models chemical bonding with bond orders rather than explicit bonds, which allows for continuous breaking and forming of bonds. ReaxFF force field parameters are derived predominantly from QM, so it may be directly applied to novel systems that may not have been extensively studied experimentally. This makes it a very attractive tool for exploring AB oxidation kinetics. Previously published studies have used ReaxFF to investigate reactive processes in many different systems, including organic reactions,19-21 shock induced chemistry and thermal decomposition of RDX,22,23 decomposition of energetic materials used in improvised explosive devices,24 thermal decomposition of polymers,25 BiMoOx heterogeneous catalysts,26 fuel cell catalysts and membranes,27 crack propagation in silicon crystals,28 dissociation of H2 on Pt surfaces,29 storage of H2 in Mg nanoclusters,30 catalytic formation of carbon nanotubes,31 tribology of metal-metal oxide interfaces,32 and the decomposition of hydrazines.33 In this article, we first describe the development of the reactive force field. Subsequently, the results of single and multimolecule NVT-MD gaseous phase simulations of the thermolysis of AB are presented in order to demonstrate the ability of the force field to describe the chemistry. Finally, we discuss an initial application of the ReaxFF method to a boron/ oxygen system, since boron oxidation is an important step in the oxidation of AB. 2. Computational Methods 2.1. QM Description. The DFT calculations used to derive the force field parameters were performed using the B3LYP functional,34 which combines exact HF exchange with the Becke generalized exchanged functional of Becke35 and Lee, Yang, and Parr36 correlation functional. All electrons were included on the B, H, O, and N atoms using a modified variant of the Pople 6-311G** basis set.37,38 All calculations were performed using the Jaguar 6.5 program package.39 2.2. Force Field Development. The QM data were used to develop the force field parameters by employing the force field

Weismiller et al. optimization routines in the ReaxFF program. These include a single-parameter-based parabolic extrapolation method,40 using successive parameter optimization cycles to resolve the parameter correlation. 2.3. MD Settings. To study the effectiveness of the force field in modeling the thermolysis of AB, two different approaches were used. First, a single AB molecule was placed in a periodic cube with a side length of 0.75 nm. Subsequently, a system of 50 AB molecules was studied inside of a periodic cube with a side length of 2.5 nm. All simulations were performed with a constant number of atoms (N), in a constant volume (V), with a specified temperature (T), otherwise known as NVT conditions. The temperature was controlled using a Berendsen thermostat41 with a 0.1 ps damping constant. The simulation of the systems studied must accurately model flexible molecules with flexible bonds that exhibit translational, rotational, torsional, and vibrational motion. Generally, this requires a time step 1 order of magnitude smaller than the shortest motion or approximately 0.5-1.0 fs. However, with ReaxFF, a smaller time step is required because the charges and bond orders are allowed to change at every time step. For high-temperature (g2500 K) ReaxFF MD simulations, a time step of 0.1 fs allows for efficient coverage of the phase space and collisions and reactions to occur smoothly. For temperatures below 2500 K, a time step of 0.25 fs is sufficiently small. The single molecule simulations were run for multiple unique starting configurations, in order to statistically assess the decomposition. The atom coordinates and velocities for the various configurations were taken from different time steps of a simulation run at 300 K. The initiation time for the release of one H2 molecule from the AB molecule was found for each of the different configurations over a temperature range of 1000-1750 K. This process is expected to be largely unimolecular, so the single molecule simulations are the most efficient method of investigation. The 50 AB molecule simulations were run through a range of temperatures at a constant heating rate, by varying the thermostat from 300-5000 K. From this simulation, a rate of H2 release per unit change in temperature can be calculated, as well as the temperature at which the greatest amount of hydrogen is released. Furthermore, the H2NBH2 molecule was studied in the same manner in an attempt to study the second and third equivalent H2 release from the AB molecule without the system being flooded with H2 molecules. A compositional analysis of the product and intermediate species was performed using a bond order cutoff of 0.2 for the recognition of molecules. The bond order cutoff has no affect on the simulations, only the interpretation of the molecular fragments. A low bond order cutoff was used mostly because the B-N bond in the AB molecule has a rather low bond order. Table 1 shows the density, starting pressure, and final pressure for each of the molecular dynamic simulations. The use of high pressure is favored in order to facilitate enough collisions, thus keeping simulation times reasonable. 3. Results and Discussion 3.1. FF Development and Parameters. Figures 2-7 and Table 2 describe the QM-data used to develop the ReaxFF parameters for AB dehydrogenation and oxidation. The following strategy was employed in the force field development: (1) QM data were generated describing the single and (if relevant) double and triple bond dissociation for all B/N/O/H combinations. These data were used to derive initial ReaxFF bond parameters.

ReaxFF Development and MD Simulations of NH3BH3

J. Phys. Chem. A, Vol. 114, No. 17, 2010 5487

TABLE 1: Respective Configurations and Pressures for MD Simulations simulation

heating cube rate [K/fs] length [nm]

temperature [K]

single AB molecule constant temp.

N/A

0.75

50 50 50 20

0.005 22 0.005 22 N/A 0.004 44

2.5 2.5 2.5 2.5

1000; 1125; 1250, 1500; 1625; 1750 AB molecules temp. ramp 300-5000 H2NBH2 molecules temp. ramp 300-5000 AB molecules constant temp. 2500 B and 15 O2 molecules temp. ramp 300-4300

(2) After obtaining initial bond parameters the training set was extended with QM data describing angular distortions in a set of small AB-related molecules. These data were used to derive the initial ReaxFF angular parameters. (3) Subsequently, the training set was extended with reaction barriers for key reaction steps (H2 release from AB, dimerization of H2B-NH2) and reaction energies associated with H2 release from AB and with AB oxidation.

starting pressure [Mpa]

final pressure [MPa]

32.727; 36.817; 40.909, 49.090; 53.181; 57.272 13.254 13.250 110.417 9.275

65.454; 73.634; 81.818; 98.180; 106.362; 114.544 583.187 406.464 234.090 34.196

(4) Final parameter values were determined by training all parameters against the full training set. Since the ReaxFF results for N/O/H bonds and angles were already presented earlier22,23 we will focus here on the boronrelated aspects of the force field development. Figures 2-5 compare the ReaxFF and QM results for B-H bond dissociation (Figure 2), B-N bond dissociation (Figure 3), B-O bond dissociation (Figure 4), and B-B bond dissociation (Figure 5). A challenging aspect of boron chemistry involves the relative stability of 4-coordinated boron compounds (BH4, H3B-NH3, H3B-OH2), which play a key part in AB-chemistry. We observe in Figures 2-5 that ReaxFF properly describes the dissociation energies for these compounds, although there are significant differences between the ReaxFF and QM shape of the bond

Figure 2. QM and ReaxFF results for B-H bond dissociation in BH2 and BH4.

TABLE 2: QM- and ReaxFF Energies (kcal/mol) for Reactions Related to Hydrogen Release from AB and for AB Oxidation reaction

QM

ReaxFF

H3B-NH3 f H2B-NH2 + H2 H2B-NH2 f HBdNH + H2 HBdNH f BN + H2 H3B-NH3 f H3B + NH3 H3B-NH2 f H3B + NH2 H2B-NH3 f H2B + NH3 H3B-NH3 f H3B-NH3 (eclipsed) 2BH3 + O2 f 2HOBH2 2BH3 + O2 f 2O-BH3 B(OH)3 f HOB ) O + H2O HOB ) O f BO2 + H2 BO2 f BO + 1/2O2 BO3 f BO2 + 1/2O2 BO f B + O BO2 + BO2 f B2O3 + 1/2O2 B2O3 + H2O f 2HOB ) O 4H3B-NH3 + 9O2 f 2B2O3 + 12H2O + 2N2 4HB ) NH + 5O2 f 2B2O3 + 4H2O + 2N2 H3B-NH3 + O2 f HB ) NH + 2H2O

-5.1 29.5 132.8 32.8 46.9 36.1 2.3 -179.5 7.03 53.1 62.6 72.5 -26.7 190 -62.4 -7.91 -926.8

-3.6 44.7 97.5 29.8 65.0 53.0 2.21 -188.4 5.3 52.4 51.4 91.6 -23.1 201.8 -36.9 0.04 -941.3

-633.9

-585.4

-73.2

-86.6

Figure 3. QM and ReaxFF results for B-N bond dissociation in H3B-NH3, H2B-NH2, and HBdNH.

Figure 4. QM and ReaxFF results for B-O bond dissociation in H3B-OH2, H2B-OH, and HBdO.

5488

J. Phys. Chem. A, Vol. 114, No. 17, 2010

Weismiller et al.

Figure 5. QM and ReaxFF results for B-B bond dissociation in H2B-BH2.

Figure 7. QM and ReaxFF reaction energies and barriers for H2 release from H3B-NH3(a) and for H2B-NH2 dimerization (b).

Figure 6. QM and ReaxFF results for B-N-B angle distortion in H2B-NH-BH2 and for N-B-N angle distortion in H2N-BH-NH2.

dissociation curves for these overcoordinated boron species. For more traditional, 3-coordinated, boron compounds ReaxFF provides an accurate reproduction of the entire QM bond dissociation curve. Figure 7 and Table 2 demonstrate that the deviations between the ReaxFF and QM shapes of the bond dissociation curves do not affect the ability of ReaxFF to reproduce both the barriers (Figure 7) and reaction energies (Table 2) for H2 release from AB and for AB oxidation. In Figure 7b, which shows the energy curve for H2B-NH2 dimerization, the QM calculations show a point of inflection at a B-N bond distance of approximately 2.2 Å, which is not apparent in the ReaxFF curve. This feature might be associated with a significant configurational change in the dimer as it approaches the transition state for the dimerization reaction. Until a B-N distance of 2.2 Å, the approach follows a fairly straightforward trajectory; between 2.2 and 2.0 Å. However, a partial B-H(B) bond is formed between the two H2B-NH2 fragments (as described in Nguyen et al.42), which lowers the dimer energy and leads to the appearance of a shoulder as shown in Figure 7b. Some other notable discrepancies, which can be found in Table 2, are for the HBdNH f BN + H2 reaction, which shows a disagreement of 26.6% between QM and ReaxFF, and the BO2 + BO2 f B2O3 + 1/2O2 reaction, with a discrpancy of

Figure 8. Sequence of images showing the release of H2 from AB. (green ) boron, blue ) nitrogen, white ) hydrogen).

40.8%. However, most of the reaction energies are within a few kcal/mol, indicating that this force field should be suitable for describing AB chemistry. 3.2. AB Single Molecule Dissociation. Figure 8 shows the sequence of a single AB releasing one H2 molecule. These snapshots come from a simulation performed at 2000 K. As the molecule heats up, the hydrogen bonds begin to stretch significantly. This eventually leads to one of the hydrogen atoms bonded to the boron forming a bond with a hydrogen atom formally bonded to the nitrogen. The resulting H2 breaks off from the rest of the molecule, leaving H2NBH2. Figure 9 is a plot of the reaction rate constant (k) versus the inverse of temperature. The reaction rate constant was calculated by running single AB molecule simulations at several different configurations for temperatures between 1000 and 1750 K and observing the amount of time it took to release one H2. An exponential curve was fit to the results to calculate the activation

ReaxFF Development and MD Simulations of NH3BH3

Figure 9. Reaction rate constant (k) versus inverse temperature for single AB molecule, constant temperature simulations.

Figure 10. Rate of hydrogen release per change in temperature for 50 AB and 50 H2NBH2 temperature ramp simulations.

energy according to an Arrhenius model. The activation energy calculated (Ea ) 26.36 kcal mol-1 or 110.3 kJ mol-1) is very close to the value of 117 kJ mol-1 reported by ref 43 and somewhat lower than the value of 184 kJ mol-1 reported in ref 44. 3.3. AB Poly Molecular Dissociation. Figure 10 shows the rate of hydrogen release per increase in temperature from the simulations of 50 AB molecules and 50 H2NBH2 molecules, respectively. In the case of AB, the release of hydrogen peaks around 1250 K and remains substantial until a temperature of approximately 3000 K, where the H2 release appears to reach equilibrium. For the H2NBH2 simulation, the release of hydrogen peaks around 2150 K, and equilibrium is reached around 3000 K. In the case of AB, we observe the release of the first and second equivalent hydrogen, which is why the release starts at lower temperatures and remains substantial over a longer range of temperatures. In the H2NBH2 simulation, we are able to observe where the release of the second hydrogen molecule occurs. The release of a third equivalent hydrogen requires the formation of a B-N structured polymer. Although some polymerization is realized in the simulations, the time scales required for this process are much longer. It should be noted

J. Phys. Chem. A, Vol. 114, No. 17, 2010 5489 that experimental studies have shown AB to vigorously release hydrogen at 120 °C45 and, if provided ample time, release hydrogen at temperatures as low as 75 °C.7 However, to make the computing time reasonable for the present simulations, the temperature and heating rates needed to be relatively high, so that hydrogen release could be observed for the very short reaction times (∼1 ns) of the simulations. Figure 11 shows the number of molecules, the average molecular weight, the average molecular weight of species with a B-N bond, and the number of B-N bonds versus temperature for the AB temperature ramp simulation. Initially, the average molecular weight and number of B-N bonds both drop quite rapidly. This can be attributed to two phenomena: the release of the H2 from AB molecules, and the breaking of the B-N bond to form NH3 and BH3 molecules. Around 1200 K, the NH3 and BH3 molecules begin to join again, and the number of B-N bonds rises again. Also, the slope in the average molecular weight curve becomes much flatter at this point, which is also approximately where the H2 release was observed to peak. The increase in B-N bonds and the average molecular weight corresponds to the formation of polymers, an example of which is shown in Figure 11. The formation of similar polymeric aminoborane has been found experimentally by Wolf et al.7 through IR-spectroscopy, elemental analysis, and XRD analysis. The number of B-N bonds and the average B-N species molecular weight peaks around 2500 K. As the temperature is raised further, the polymers tend to break apart, as can be observed from the decrease in B-N species molecular weight and the number of B-N bonds. Panels a-d in Figure 12 show snapshots of the AB temperature ramp simulations. Figure 12a is the starting geometry. Figure 12b shows the atom coordinates at a temperature of 1083 K. At this point in the simulation, a number of AB molecules have released one equivalent hydrogen, and some have disassociated to form NH3 and BH3 molecules. Figure 12c is a snapshot at 2649 K. At this temperature, most of the AB molecules have given up their first equivalent hydrogen, as can be observed from the large quantities of H2 and H2NBH2 in the volume of interest. Some of the molecules have given up a second H2 to form HNBH, some of which have then begun to form B-N polymers. All of the NH3 and BH3 that was formed at the lower temperatures and pressures have now recombined. Figure 12d is a snapshot from the end of the simulation, where the temperature has reached 5000 K. There is some more polymerization at this point, further releasing hydrogen. Figure 13 is taken from a constant temperature simulation performed at 2500 K, using the same starting geometry as the temperature ramp simulation (Figure 10a). This temperature was chosen because it corresponds to where the greatest amount of polymerization was observed in the temperature ramp simulation. At this temperature, the AB decomposes very quickly, releasing H2 and forming H2NBH2. The H2NBH2 then decomposes (at a much slower rate) to form both HNBH and some B-N-based polymers, steadily increasing the amount of hydrogen in the system. 3.4. B + O2 MD Results. Figure 14 shows the change in the total number of molecules, number of O2 molecules, and average molecular weight (MW) for boron containing species from a B and O2 temperature ramp MD simulation. The number of molecules steadily drops throughout most of the temperature ramp, owing first to the clustering of boron atoms, and later to the oxidation of the boron cluster. At elevated temperatures, the number of molecules begins to increase again, owing to the fragmentation of the oxidized boron cluster into gas phase

5490

J. Phys. Chem. A, Vol. 114, No. 17, 2010

Weismiller et al.

Figure 11. Average molecular weight and number of molecules plotted vs temperature for 50 AB temperature ramp simulation, and an example of the polymers that are formed.

Figure 12. Snapshots from a 50 AB temperature ramp simulation (green ) boron, blue ) nitrogen, white ) hydrogen) at (a) 300 K, (b) 1083 K, (c) 2649 K, and (d) 5000 K.

molecules. The MW of the boron species initially rises rapidly until T ) 500 K, because of the boron clustering. Once the boron atoms have all formed into a cluster, the average MW rises more slowly due to oxidation on the outside of the cluster.

The MW begins to drop around 3500 K, due to the formation of some gaseous species. Panels a-f of Figure 15 show snapshots from these B and O2 temperature ramp simulations. Figure 15a is the starting

ReaxFF Development and MD Simulations of NH3BH3

J. Phys. Chem. A, Vol. 114, No. 17, 2010 5491

Figure 13. Composition vs time for a 50 AB constant temperature simulation at 2500 K.

Figure 15. Snapshots from 20B and 15O2 temperature ramp simulation (green ) boron, red ) oxygen) at (a) 300 K, (b) 700 K, (c) 1600 K, (d) 2800 K, (e) 3600 K, and (f) 4300 K.

Figure 14. Average MW for Boron species and number of molecules for 20B and 15O2 temperature ramp.

geometry (T ) 300 K and time ) 0 fs); panels b and c of Figure 15 correspond to temperatures of 700 and 1600 K and show the formation of the boron cluster and some O2 bonding at the periphery. Figure 15d corresponds to a temperature of 2800 K, where the cluster has been even further oxidized; and Figure 15e, which corresponds to 3600 K, shows the formation of a gas phase BO2 molecule. Finally, Figure 15f, which corresponds to a temperature of 4300 K, shows the cluster has consumed all but one of the O2 molecules and has significantly decomposed, forming several gaseous B2O3 molecules. 4. Conclusions Ammonia borane (AB) is comprised of 19.6% hydrogen by mass, making it an attractive fuel for several applications. A ReaxFF force field has been developed in order to study the dehydrogenation of AB and the boron oxidation. The ReaxFF parameters, have been fitted directly from QM. To demonstrate the validity of the force field for MD, NVT simulations of AB

pyrolysis and boron oxidation were performed. It was shown that AB releases the first and second equivalent H2 via unimolecular reactions; however, further dehydrogenation requires the formation of B-N polymers. Other works42 have used QM calculations to show that there exists a pathway with a lower energy barrier in which a BH3 molecule acts as a catalyst for H2 elimination. This reaction was investigated using the ReaxFF force field and shows a similar barrier (approximately 15 kcal mol-1) to that calculated by Nguyen et al.42 In the simulations presented here, a significant quantity of BH3 is formed by bond cleavage of the B-N bond at the lower temperatures. However, the BH3 molecules tend to recombine with NH3 instead of interacting with the AB molecules. At higher temperatures, the unimolecular H2 elimination is fast, and there is not sufficient time for the bimolecular BH3 catalysis reaction to take place. The catalytic BH3 reaction may become dominant if the simulations are held at lower temperatures, and the NH3 molecules are allowed to escape from the system. Single-molecule MD simulations were used to calculate an activation energy of 26.36 kcal mol-1 for the release of the first equivalent H2, which is in relatively good agreement with experimental results.44,46 A poly molecular AB system was studied using a temperature ramp (300-5000 K) simulation, and the rate of hydrogen release was shown to peak at approximately 1250 K. An identical simulation was performed with the H2NBH2 molecule, in order to investigate the second equivalent hydrogen release, and rate of hydrogen release peaked around 2150 K. Additionally, the temperature ramp simulations showed that at low temperatures, AB may dissociate to form

5492

J. Phys. Chem. A, Vol. 114, No. 17, 2010

NH3 and BH3, but these molecules recombined as the temperature eclipsed 1200 K. Note that the pressure of the system is greatly increased at this point, due to the large release of H2 molecules into the system. Furthermore, B-N polymers were observed to form from the dehydrogenated AB molecules. This polymerization takes place over much longer time scales than the unimolecular hydrogen release. This polymerization is similar to that observed by Wolf et al.7 in experiments. The polymerization peaks around 2500 K in the simulation, since at the higher temperatures the polymers tend to break apart. The temperature-ramping (300-4300 K) MD simulations of boron and oxygen show the B atoms agglomerate into a cluster at low temperatures. As the temperature is elevated, O2 molecules increasingly attack the outside of the cluster. At approximately 3500 K, some gaseous boron oxide species begin to form. The oxidation of boron cluster ions has been studied experimentally by Hanley et al.,47 and show similar processes for boron clusters with more than six atoms. At the end of the simulation, where the temperature is 4300 K, most of the O2 molecules have been consumed and the cluster has begun to break apart, forming several gaseous B2O3 molecules. The condensed phase B2O3 has a similar structure to that proposed by Brown et al. in.48 These results provide the confidence needed to employ ReaxFF to investigate the rapid oxidation of AB at an atomistic level in future work. This modeling should provide understanding of the reaction pathways and an explanation for experimental observations at the mesoscale. Acknowledgment. The authors would like to gratefully acknowledge and thank the Air Force Office of Scientific Research (AFOSR) and NASA for their sponsorship of this program under contract No. FA9550-07-1-0582. ACTvD acknowledges startup funding from KISK No. C000032472. Supporting Information Available: The force field parameters used in this work are available as Supporting Information. This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) Zuttel, A. Mater. Today 2003, 6, 24. (2) Zuttel, A. Naturwissenschaften 2004, 91, 157. (3) Zuttel, A.; Wenger, S.; Rentsch, S.; Sudan, P.; Mauron, P.; Emmenegger, C. J. Power Sources 2003, 118, 1. (4) Graetz, J. Chem. Soc. ReV. 2009, 38, 73. (5) Marder, T. B. Angew. Chem.-Int. Edit. 2007, 46, 8116. (6) Stephens, F. H.; Pons, V.; Baker, R. T. Dalton T. 2007, 2613. (7) Wolf, G.; Baumann, J.; Baitalow, F.; Hoffmann, F. P. Thermochim. Acta 2000, 343, 19. (8) Bowden, M.; Autrey, T.; Brown, I.; Ryan, M. Curr. Appl. Phys. 2008, 8, 498. (9) Baitalow, F.; Wolf, G.; Grolier, J. P. E.; Dan, F.; Randzio, S. L. Thermochim. Acta 2006, 445, 121. (10) Basu, S.; Brockman, A.; Gagare, P.; Zheng, Y.; Ramachandran, P. V.; Delgass, W. N.; Gore, J. P. J. Power Sources 2009, 188, 238. (11) Xu, Q.; Chandra, M. 2007, 446, 729.

Weismiller et al. (12) Chandra, M.; Xu, Q. J. Power Sources 2006, 159, 855. (13) He, T.; Xiong, Z. T.; Wu, G. T.; Chu, H. L.; Wu, C. Z.; Zhang, T.; Chen, P. Chem. Mater. 2009, 21, 2315. (14) Metin, O.; Ozkar, S. Energ. Fuel. 2009, 23, 3517. (15) Chandra, M.; Xu, Q. J. Power Sources 2006, 156, 190. (16) Chandra, M.; Xu, Q. J. Power Sources 2007, 168, 135. (17) Xu, Q.; Chandra, M. J. Power Sources 2006, 163, 364. (18) Diwan, M.; Diakov, V.; Shafirovich, E.; Varma, A. Int. J. Hydrog. Energy 2008, 33, 1135. (19) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. J. Phys. Chem. A 2001, 105, 9396. (20) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. J. Phys. Chem. A 2008, 112, 1040. (21) Chenoweth, K.; Van Duin, A. C. T.; Dasgupta, S.; Goddard, W. A. J. Phys. Chem. A 2009, 113, 1740. (22) Strachan, A.; Kober, E. M.; van Duin, A. C. T.; Oxgaard, J.; Goddard, W. A. J. Chem. Phys. 2005, 122, 10. (23) Strachan, A.; van Duin, A. C. T.; Chakraborty, D.; Dasgupta, S.; Goddard, W. A. Phys. ReV. Lett. 2003, 91, 4. (24) van Duin, A. C. T.; Zeiri, Y.; Dubnikova, F.; Kosloff, R.; Goddard, W. A. J. Am. Chem. Soc. 2005, 127, 11053. (25) Chenoweth, K.; Cheung, S.; van Duin, A. C. T.; Goddard, W. A.; Kober, E. M. J. Am. Chem. Soc. 2005, 127, 7192. (26) Goddard, W. A., III; van Duin, A.; Chenoweth, K.; Cheng, M. J.; Pudar, S.; Oxgaard, J.; Merinov, B.; Jang, Y. H.; Persson, P. Top. Catal. 2006, 38, 93. (27) Goddard, W. A., III; Merinov, B.; van Duin, A.; Jacob, T.; Blanco, M.; Molinero, V.; Jang, S. S.; Jang, Y. H. Mol. Simul. 2006, 32, 51. (28) Buehler, M. J.; van Duin, A. C. T.; Goddard, W. A. Phys. ReV. Lett. 2006, 96, 4. (29) Ludwig, J.; Vlachos, D. G.; van Duin, A. C. T.; Goddard, W. A. J. Phys. Chem. B 2006, 110, 4274. (30) Cheung, S.; Deng, W. Q.; van Duin, A. C. T.; Goddard, W. A. J. Phys. Chem. A 2005, 109, 851. (31) Nielson, K. D.; van Duin, A. C. T.; Oxgaard, J.; Deng, W. Q.; Goddard, W. A. J. Phys. Chem. A 2005, 109, 493. (32) Zhang, Q.; Cagin, T.; van Duin, A.; Goddard, W. A.; Qi, Y.; Hector, L. G. Phys. ReV. B 2004, 69, 11. (33) Zhang, L. Z.; van Duin, A. C. T.; Zybin, S. V.; Goddard, W. A. J. Phys. Chem. B 2009, 113, 10770. (34) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (35) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (36) Lee, C. T.; Yang, W. T.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (37) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; Defrees, D. J.; Pople, J. A. J. Chem. Phys. 1982, 77, 3654. (38) Harihara., Pc; Pople, J. A. Chem. Phys. Lett. 1972, 16, 217. (39) Jaguar 6.5; Schrodinger, LLC: Portland, OR, 2005. (40) Van Duin, A. C. T.; Baas, J. M. A.; Graaf, B. J. Chem. Soc. Faraday T. 1994, 90, 2881. (41) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (42) Nguyen, M. T.; Nguyen, V. S.; Matus, M. H.; Gopakumar, G.; Dixon, D. A. J. Phys. Chem. A 2007, 111, 679. (43) Brown, C. M.; Jacques, T. L.; Hess, N. J.; Daemen, L. L.; Mamontov, E.; Linehan, J. C.; Stowe, A. C.; Autrey, T. Physica B 2006, 385, 266. (44) Gutowska, A.; Li, L. Y.; Shin, Y. S.; Wang, C. M. M.; Li, X. H. S.; Linehan, J. C.; Smith, R. S.; Kay, B. D.; Schmid, B.; Shaw, W.; Gutowski, M.; Autrey, T. Agnew. Chem. Int. Ed. 2005, 44, 3578. (45) Hu, M. G.; Geanangel, R. A.; Wendlandt, W. W. Thermochim. Acta 1978, 23, 249. (46) Bowden, M.; Kemmitt, T.; Shaw, W.; Hess, N. J.; Linehan, J.; Gutowski, M.; Schmid, B.; Autrey, T. Mater. Res. Soc. Symp. Proc. 2006, 927, 9. (47) Hanley, L.; Anderson, S. L. J. Chem. Phys. 1988, 89, 2848. (48) Brown, R. C.; Kolb, C. E.; Rabitz, H.; Cho, S. Y.; Yetter, R. A.; Dryer, F. L. Int. J. Chem. Kinet. 1991, 23, 957.

JP100136C