Twin Induced Sensitivity Enhancement of HMX versus Shock: A

Oct 25, 2013 - ADVERTISEMENT .... The simulations show that the twin enhances the shock sensitivity remarkably, in agreement with our ... is helpful t...
0 downloads 0 Views 4MB Size
Download PDF
Article pubs.acs.org/JPCC

Twin Induced Sensitivity Enhancement of HMX versus Shock: A Molecular Reactive Force Field Simulation Yushi Wen,† Xianggui Xue,† Xiaoqing Zhou,† Feng Guo,‡ Xinping Long,† Yang Zhou,† Hongzhen Li,† and Chaoyang Zhang*,† †

Institute of Chemical Materials, China Academy of Engineering Physics (CAEP), P.O. Box 919-327, Mianyang, Sichuan 621900, China ‡ Institute of Atomic and Molecular Physics, Sichuan University, Chengdu, Sichuan 610065, China S Supporting Information *

ABSTRACT: We report the early events of a twinned HMX crystal, as well as a perfect one for a comparison purpose, shocked with various velocities in the range of 6−10 km/s for 50 ps by molecular dynamic simulations using the ReaxFF reactive force field and the multiscale shock technique. The simulations show that the twin enhances the shock sensitivity remarkably, in agreement with our recent experimental observation. That is, it exhibits a lower shock initiation stress, higher decomposition velocities, more temperature, and stress increases under the same shock conditions, relative to the perfect crystal. In addition, we find the twinning plane is hottest and the temperature decreases in terms of the distance apart from it after shock loaded earlier, suggesting possible hot spots preferred there.

1. INTRODUCTION Safety is one of the two most concerns of explosives and attracts continuous attention, due to its practical importance and many complicated sciences involved in it. Crystal quality, as well as other factors like molecular stability and surface or interfacial integrity, can influence explosive safety remarkably. The safety of an explosive is usually evaluated by sensitivity, the response of an explosive to an external stimulation: the lower the sensitivity, the higher the safety. With respect to explosive crystals, the higher quality usually leads to the lower sensitivity. One of the requirements of high quality crystals is the imperfections contained as few as possible. Twin is a common kind of imperfection of explosive crystals and has been verified to be an evident reason for increasing shock sensitivity experimentally.1 However, the details of this sensitivity enhancement remain unclear. Actually, the sensitivity mechanism involved in the evolution of explosives against stimuli is very complicated. For example, the detonation reactions, included in the sensitivity mechanism, of the simplest model explosive nitromethane have not been understood completely yet. Currently, the molecular reactive force field method like ReaxFF2 becomes a main tool to explore condensed phase reactions on large scales, and further the sensitivity mechanisms. For instance, ReaxFF has been successfully applied to study the chemical responses of energetic materials 1,3,5,7-tetranitro-1,3,5,7-azacyclo-octane (HMX), 1,3,5-trinitro-1,3,5-triazacyclohexane (RDX), 1,3,5triamino-2,4,6-trinitrobenzene (TATB), pentaerythritol tetrani© 2013 American Chemical Society

trate (PETN), and nitromethane undergoing shockwave compression and heating.2−12 This work aims to reveal the details of the twin induced enhancement of shock sensitivity of HMX, a famous explosive and the so-called king of explosives. Thousands of atoms are required to describe an HMX twin adequately with respect to its structure, much beyond usual quantum chemical calculations and preferred to the ReaxFF method.13 That is, the chemical responses of a twinned HMX (TH for short) as well as a perfect HMX crystal (PH for short) to shock with various velocities (Us) from 6 to 10 km/s were studied using the ReaxFF method combined with the multiscale shock technique (MSST).14 The simulation results show that TH is more sensitive to shock than PH, which is helpful to understand the twin-induced enhancement mechanism of shock sensitivity of HMX. It also suggests the methods employed, such as ReaxFF force field and MSST, can be extended to explore some complicated mechanisms like those involved in crystal twins.

2. METHODOLOGIES 2.1. Modeling. A super cell with a twinning plane and 320 HMX molecules (8960 atoms) was built for our simulations (see the details in S1 of the Supporting Information). The simulations on this scale regarding atomic amount should be Received: July 23, 2013 Revised: October 20, 2013 Published: October 25, 2013 24368

dx.doi.org/10.1021/jp4072795 | J. Phys. Chem. C 2013, 117, 24368−24374

The Journal of Physical Chemistry C

Article

Figure 1. Structural comparison of the perfect (a) and twinned (b) HMX crystals. The red dash shows the twinning plane, and white arrows point to the shock loading orientations. The carbon, hydrogen, oxygen, and nitrogen atoms are represented in gray, white, red, and blue, respectively.

motions and interatomic forces, which is derived from the virial theory.23 This method has been applied to nanoscale systems successfully.24

workable, since the computational cell of the multiscale shock technique14,15 follows a lagrangian point through the shock wave which enables a simulation of a system undergoing a shock wave with far fewer atoms than normally required by the direct method of inducing a shock wave in a very large computational cell.16 The structural comparison of PH and TH, as well as the loading orientations of shock, is shown in Figure 1. The principle of ReaxFF can be referred elsewhere,17 and its improved version ReaxFF_lg7 was employed (see the validation of the method in S2 of the Supporting Information). The shockwave compression was loaded along the (001) face direction using the MSST, just vertical to the twinning plane. MSST proposed by Reed et al.14 has been successfully applied to energetic materials against shock.17−19 The molecular dynamic (MD) simulations lasted for 50 ps with a time step of 0.1 fs, and the dynamic trajectories were collected per 1 ps for analyzing. Both ReaxFF and MSST are implemented in the LAMMPS software package,20,21 with which all MD simulations in this work were carried out. 2.2. Shock Sensitivity Comparison of TH and PH. Explosive sensitivity is usually assessed by the response to external stimuli. As to the shock sensitivity comparison, the results derived from the ReaxFF-MD simulations, such as decay rates of reactants, temperature, stress increases, and the shock initiation stress, were employed. To get chemical species and their distributions in the dynamic evolution, we used FindMole proposed by us22 to analyze dynamic trajectories. FindMole is a FORTRAN program whose theory for identifying a molecule or a fragment is based on Strachan et al.’s proposal.2 Then, the reaction kinetics of the HMX decay can be found based on the chemical species analyses, with an assumption that the shock induced reactions follow first order decay. This assumption was also carried out in dealing with the thermal decay of nitromethane.12 For first order reaction, the reactant decay rate k is k=

1 a ln t a−x

N

∑i = 1 ∑j < i rij ⊗ fij NkBT P= + V dV

(3)

where N is the number of atoms in the system, kB is the Boltzmann constant, T is the temperature, d is the dimensionality of the system, V is the system volume, f ij is the interatomic force on atom i by atom j, and rij is the position vector from atoms i to j. The temperatures of total systems and specified regions were calculated using the formula of the kinetic energy (KE)−temperature relationship, KE = (d/2) NkBT. The shock initiation stress P was calculated using eq 4 P = ρ0 UU s p

(4)

where ρ0, Us, and Up are the initial density, the lowest shock velocity initiating reactions, and the corresponding particle velocity, respectively.

3. RESULTS AND DISCUSSION 3.1. Comparison of the Decay Rates of PH and TH. First, we pay attention to the HMX decay versus loading time of shock with various Us. As illustrated in Figure 2, overall, to

(1)

where t, a, and x are the time and initial molecular amount and molecular amount at t of reactants, respectively. Equation 1 can be converted into eq 2, in which C is a constant. Subsequently, k can be found by fitting eq 2. ln(a − x) = −kt + C

(2) Figure 2. Evolution of the amount of HMX molecules (NHMX). The blue and red are of perfect and twinned HMX (PH and TH), respectively.

The internal stresses were calculated by LAMMPS using eq 3, equal to the sum of momentum fluxes induced by atomic 24369

dx.doi.org/10.1021/jp4072795 | J. Phys. Chem. C 2013, 117, 24368−24374

The Journal of Physical Chemistry C

Article

the Us content of 7−10 km/s, the higher Us causes the faster decomposition for both TH and PH as expected. And TH (curves in red) decays much faster than PH (curves in blue) at any given Us. TH can be always analyzed completely in 50 ps, whereas, for PH, its decomposition can be finished only at Us = 10 km/s after 16 ps, and Us below 9 km/s cannot lead to a full decomposition in the time scale of our simulations. Furthermore, as presented in Figure 2, the Us values initiating TH and PH decomposition are predicted to be between 6 and 7 and 7−8 km/s after shock for 50 ps, respectively, suggesting the easier ignition of TH. Regarding the decay rate k in Table 1, TH possesses a higher k at any shock compression than PH. All of these show that TH is more sensitive to shock than PH, in agreement with our experimental observation.1

exhibit the evolutions of cluster amounts and the maximum sizes. Some conclusions from the figures can be drawn as follows: (1) It seems that the cluster formation is always a result in shock compression induced reactions. There is no cluster formation for PH when Us = 7 km/s, which cannot initiate PH decay within the simulation time scale. Figure 3 shows no cluster formation, no reaction. As pointed out previously, this is rationally attributed to the fact that the shortened intermolecular distances caused by compression make cluster formation advantageous. (2) The condition for an easier reaction results in easier cluster formation. For example, TH is much easier to form clusters than PH at any given Us, and a higher Us leads to the easier cluster formation for both PH and TH than a lower one. Figure 3 shows a slower cluster formation for PH, after about 1 ps of shock, relative to the immediate formation for TH when Us is within 8−10 km/s. And for PH, the formation occurs at 2.5 ps when Us is equal to 8 km/s; it is shortened to about 1 ps when Us increases to 9 or 10 km/s. This time decrease for cluster formation by increasing Us can also be found for TH. (3) A higher Us leads to a bigger cluster in the simulation time scale. This is supported by the results shown in Figure 4. For PH, the maximum MWs of the clusters decrease from about 56 000 to 1400 and 600 when Us decreases from 10 to 9 and 8 km/s, while, for TH, it decreases from 60 000 to 20 000 and 4000, and (4) the cluster amounts and the maximum MWs fluctuate usually during reaction processes, implying the unceasing combinations and dissociations of these clusters. The fluctuation can be seen if most reactions can be completed in our simulation time scale, 50 ps, corresponding to the cases of Us = 8, 9, and 10 km/s for TH and Us = 10 km/s for PH. Comparing Figures 3 with 4, one can find that, on the time axes, the maximum MWs are just located on the points corresponding to local minimum cluster amounts, suggesting lots of smaller clusters are combined to bigger ones. Also, the stability of the biggest clusters in these cases can be evaluated by their life length: the biggest cluster composed of about 60% of total reactants can last for about 10 ps when Us = 10 km/s for both PH and TH, about 7 ps for the biggest cluster caused by Us = 9 km/s shock on TH, and longer for that induced by Us = 8 km/s shock on TH. As a whole, the formed clusters tend to small fragments as time proceeds. For example, Figure 5 exhibits this tendency: the cluster size decreases sharply after Us = 9 km/s shock loading on TH for about 15 ps, and the small fragments increase rapidly simultaneously. In contrast to Figure 3, for TH shocked with the same velocity, we can find that the

Table 1. Fitted k (in 1012 s−1) of Decay of HMX Shocked with Various Us (in km/s) Us

PH

TH

6 7 8 9 10

∼0 ∼0 0.003 0.021 0.278

∼0 0.038 0.115 0.754 0.939

As for the complex reaction details induced by shock compression, Ge et al. have described the initial decomposition mechanism dependent on Us using a self-consistent charge density-functional tight-binding (SCC-DFTB) method.4 Zhou et al. performed ReaxFF-MD in conjunction with the compressive shear reactive dynamics (CS-RD) method to study the shock sensitivity of β-HMX and found the initial chemical reaction is N−NO2 bond dissociation.8 In our ReaxFF-MD simulations, we will pay much attention to clusters because their formation in shocked HMX has never been mentioned, while some small fragments will be discussed due to their important contributions to the decay mechanism. At initial stage during shock loading, many molecular clusters formed by HMX decay are observed as reported formerly on shocked or heated energetic CHNO materials.3,6,19,22 Some snapshots of them are illustrated in Figure s3 in the Supporting Information. We discuss the evolution of these clusters in terms of their amounts and the maximum sizes. We defined here a fragment as a molecular cluster provided its molecular weight (MW) above 296, the MW of reactant HMX. Figures 3 and 4

Figure 3. Evolution of cluster amounts. 24370

dx.doi.org/10.1021/jp4072795 | J. Phys. Chem. C 2013, 117, 24368−24374

The Journal of Physical Chemistry C

Article

Figure 4. Evolution of the maximum molecular weight of clusters.

NH 3 , similar to other simulated or experimental results.4,8,10,25−29 Because the case of TH against Us = 8 km/s shock is much closer to test, that is, moderate pressure and temperature caused, we select the simulated products of this case when HMX is shocked for 50 ps in Figure 5f to compare with Donald L. Ornellas’ experimental observation.29 As listed in Table 2, the simulated and the observed are qualitatively Table 2. Product Comparison of Shocked HMX from and Heavily Confined Charges (a) ReaxFF-MD Simulations (b) (Fragment in mol/mol HMX)

Figure 5. Evolution of the maximum MWs of clusters and small fragment amounts when Us = 9 km/s for TH.

first decrease of the cluster amount is completed after 15 ps and the second after 35, just corresponding to the rapid increases of small fragments and reaction equilibrium except from N2 fluctuation. Next, we are concerned about some small fragments, whose evolution in the entire simulation time scale is illustrated in Figure 6. That the NO2 fragments appear sometime and disappear soon in Figure 5a suggests no evident reaction in this case, while parts b−h of Figure 6 exhibit obvious chemical reactions due to many new species caused. NO2, NO, and HONO are always involved in the initial products, and NO2 is dominant for both PH and TH undergoing various shocks, suggesting that the break of N−NO2 bonds dominates the HMX decay, while the competitive hydrogen transfers exist. These two kinds of initiation reactions for gaseous and condensed HMX have been mentioned already.4,8,10,25−28 A DFTB simulation showed that the competition between the N−NO2 bond cleavage and the hydrogen transfer is dependent on shock velocity.4 However, this cannot be discriminated in our ReaxFF-MD simulations. Subsequent reactions continue and numerous new fragments appear. In this period, the contents of N2, H2O, CO2, CO, NO2, NO, H, H2, OH, NH2, NH3, HONO, O2, HCN, and other small fragments fluctuate. As time proceeds, these fragments and clusters tend to form final products, mainly including N2, H2O, CO2, CO, H2, and

products

a

b

N2 H2O CO2 CO C(s) NH3 H2 CH4 HCN C2H6

3.68 3.18 1.92 1.06 0.97 0.395 0.30 0.039 0.008 0.001

2.46 1.33 0.75 0.22 0 (0.03, Us = 10 km/s) 0.10 0.62 0 (0.003, Us = 9 km/s) 0.003 not observed

consistent with each other. In Figure 6e, g, and h, active H and OH become dominant, due to the water dissociation at high temperatures. Presumably, they will be recombined to water as temperature decreases. 3.2. Comparison of the Stress Increase of PH and TH. Internal stress increase is also an important response of crystalline HMX against shock, and the faster increase denotes the higher shock sensitivity. Figure 7 shows the stress evolution, implying the evolutions of structures and reactions. For instance, for TH at Us = 10 km/s, first, the stress increases rapidly to about 90 GPa and drops sharply to 75 GPa, corresponding to a yielding; then, the stress waves, suggesting chemical reactions; and finally, the stress decreases abruptly, denoting a rapid volume expansion. As a whole, the stress increases with the Us increasing. And the stress of TH is higher than that of PH at any same Us, suggesting TH is more sensitive to shock than PH. After loading for about 1 ps, the largest stress of PH reaches 2, 23, 32, 42, and 70 GPa, corresponding to the Us of 6, 7, 8, 9, and 10 km/s, respectively, while those of TH reach 18, 28, 38, 55, and 24371

dx.doi.org/10.1021/jp4072795 | J. Phys. Chem. C 2013, 117, 24368−24374

The Journal of Physical Chemistry C

Article

Figure 6. Evolution of chemical fragment amounts (Nfrag) within 50 ps.

The decrease is accompanied by the system volume change, which is expanded from 55−64% to 102−106%, relative to the original one. At the same time, the HMX molecules are converted into small stable molecules such as CO2, N2, and water. Additionally, high stress cannot be generated for PH when Us is below 6 km/s, whereas about 18 GPa of stress is produced for TH at the same Us, suggesting that the twins in HMX crystals can change their mechanics obviously, and further TH is more sensitive to shock than PH. 3.3. Comparison of the Temperature Increase of PH and TH. The faster temperature increase denotes the higher sensitivity. The temperature rises gradually during the shock loading, as shown in Figure 8. At any Us within the given range, both the temperature and the heating rate of TH are higher than those of PH, respectively. The temperature of PH is much below 2000 K when Us is below 9 km/s, while it rises to 7000 K at Us = 10 km/s within 50 ps. The temperature of TH is not more than 2000 K with Us below 7 km/s, but it rises to 5500, 9000, and even more than 10 000 K as Us increases to 10 km/s. That TH exhibits higher temperature than PH suggests the faster heat release and the higher sensitivity to shock for TH.

Figure 7. Stress evolution versus the loading time of shock with different velocities.

75 GPa, higher than those of PH, respectively. The high stresses of TH under 8, 9, or 10 km/s of shock keep for several tens of picoseconds, then gradually and finally rapidly decrease. 24372

dx.doi.org/10.1021/jp4072795 | J. Phys. Chem. C 2013, 117, 24368−24374

The Journal of Physical Chemistry C

Article

velocity (Us−Up) relationship by MSST shock simulations. The relationship is very important because, if we know nothing more than the relationship, the Rankine−Hugoniot equations will specify precisely the values of pressure, volume, energy, and enthalpy during shock compression.30 As illustrated in Figure 10, the shock velocity (Us) as a function of particle velocity

Figure 8. Temperature evolution against the loading time of shock with different velocities.

We have also compared the temperatures of the twinning plane and its adjacent layers with a thickness of 0.54 nm, just the thickness of a molecular layer parallel to the plane, whose histories of early 5 ps are illustrated in Figure 9. The figure

Figure 10. Us and Up Hugoniot plots derived from experiments and calculations.

(Up) obtained from ReaxFF-MD simulations is compared with the experimental results from LASL.31 The sound velocity obtained from the ReaxFF-MD simulations is 3.6 (for PH) and 3.55 (for TH) nm/ps, only 4.4 and 5.8% larger than the experimental ones of LASL.31 It is comprehensible that compare to PH, the defect inside TH slower its sound velocity. Then, the shock stresses are calculated by assigning values of ρ0, Us, and Up. The initial density ρ0 is 1.906 g/cm3, the lowest Us initiating reactions are 7 and 8 km/s for TH and PH, and the corresponding Up are 2.2 and 2.3 km/s for TH and PH, respectively. Accordingly, the shock initiation stress of TH is predicted to be 29.4 GPa, lower than that of PH, 35.1 GPa, in qualitative agreement with the experimental observation.1

4. CONCLUSION In summary, to reveal the details of twin induced enhancement of shock sensitivity, we established twinned and perfect HMX models and investigated their responses to shock using a combined method of MSST and ReaxFF. As a result, it shows the faster decomposition and heat release, the higher stress, and the higher temperature of the twinned HMX crystal in contrast to the perfect one, showing easier reaction and further higher sensitivity to shock. Also, we find the twinning plane is hottest and the temperature decreases regarding the distance apart from it after being shock loaded, suggesting possible hot spots preferred there. The Us−Up data obtained confirms that the Up of TH is significantly higher than that of PH under the same Us, which indicates higher kinetic energy and hot-spot-forming ability of the atoms in TH. A root reason for twins enhancing the sensitivity should be due to the molecules located on the twinning plane and its neighbors possessing a higher free volume of atoms.

Figure 9. Temperature evolution of various molecular layers versus the loading time at Us = 8 km/s.

indicates a temperature gradient around the plane: the temperature decreases from the plane to the perfect layers far away. The most temperature gap of about 600 K takes place after an 8 km/s shock loaded for 1 ps. The higher temperature around the twinning plane may be attributed to the higher free volume space there. The twinning structure (see Figure s1, Supporting Information) leads to more free volume space for most atoms at the twinning interface. Those atoms may be easier to be accelerated to break the bonds and form hot spots. This may also be a root for the higher shock sensitivity of TH, because the structural difference between TH and PH is dominantly due to the plane and region around the twinning plane. We can also observe from the figure that the temperature gap is eliminated generally after 5 ps, possibly attributed to compression induced compactness and thermal transmission. 3.4. Comparison of the Shock Initiation Stress of PH and TH. Most directly, the shock initiation stress represents the shock sensitivity, and the lower shock initiation stress, the higher sensitivity. To get the shock initiation stresses of TH and PH from eq 4, we first obtained their shock velocity−particle

■

ASSOCIATED CONTENT

S Supporting Information *

Establishment of the models, validation of ReaxFF to this work, and snapshots of molecular clusters formed at the primary stage 24373

dx.doi.org/10.1021/jp4072795 | J. Phys. Chem. C 2013, 117, 24368−24374

The Journal of Physical Chemistry C

Article

(25) Manaa, M. R.; Fried, L. E.; Melius, C. F.; Elstner, M.; Frauenheim, T. J. Phys. Chem. A 2002, 106, 9024. (26) Lewis, J. P.; Glaesemann, K. R.; VanOpdorp, K.; Voth, G. A. J. Phys. Chem. A 2000, 104, 11384. (27) Chakraborty, D.; Muller, R. P.; Goddard, W. A., III. J. Phys.Chem. A 2001, 105, 1302. (28) Behrens, R. J. J. Phys. Chem. 1997, 94, 6706. (29) Ornellas, D. L. J. Phys. Chem. 1968, 72, 2390. (30) Dlott, D. D. Annu. Rev. Phys. Chem. 2011, 62, 575. (31) Marsh, S. P. LASL Shock Hugoniot Data; University of California Press: 1980; p 595.

against shock. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Author Contributions

Drs Yushi Wen and Xianggui Xue are both the first authors. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We greatly appreciate the financial support from the Science and Technology Fund of CAEP (2011A0302014 and 2013B0101010), the National Natural Science Foundations of China (11072225, 21173199, and 11302199), and Science and technology innovation Fund of Institute of Chemical Materials, CAEP(No. kjcx-201305).

■

REFERENCES

(1) Li, H. Z.; Xu, R.; Kang, B.; Li, J. S.; Zhou, X. Q.; Zhang, C. Y.; Nie, F. D. J. Appl. Phys. 2013, 113, 203519. (2) Strachan, A.; van Duin, A. C. T.; Chakraborty, D.; Dasgupta, S.; Goddard, W. A., III. Phys. Rev. Lett. 2003, 91.098301. (3) Zhang, L. Z.; Zybin, S. V.; van Duin, A. C. T; Dasgupta, S.; Goddard, W. A., III; Kober, E. M. J. Phys. Chem. A 2009, 113, 10619. (4) Ge, N. N.; Wei, Y. K.; Ji, G. F.; Chen, X. R.; Zhao, F.; Wei, D. Q. J. Phys. Chem. B 2012, 116, 13696. (5) Zybin, S. V.; Goddard, W. A., III; Xu, P.; van Duin, A. C. T; Thompson, A. P. Appl. Phys. Lett. 2010, 96, 081918. (6) Han, S. P.; van Duin, A. C. T; Goddard, W. A., III; Strachan, A. J. Phys. Chem. B 2011, 115, 6534. (7) Liu, L. C.; Liu, Y.; Zybin, S. V.; Goddard, W. A., III. J. Phys. Chem. A 2011, 115, 11016. (8) Zhou, T. T.; Zybin, S. V.; Liu, Y.; Huang, F. L.; Goddard, W. A., III. J. Appl. Phys. 2012, 111, 124904. (9) Nomura, K.; Kalia, R. K.; Nakano, A.; Vashishta, P. Appl. Phys. Lett. 2007, 91, 183109. (10) Zhou, T. T.; Huang, F. L. J. Phys. Chem. B 2011, 115, 278. (11) Shan, T.; Wixom, R.; Mattsson, A.; Thompson, A. J. Phys. Chem. B 2013, 117, 928. (12) Rom, N.; Zybin, S. V.; van Duin, A. C. T.; Goddard, W. A.; Zeiri, Y.; Katz, G.; Kosloff, R. J. Phys. Chem. A 2011, 115, 10181. (13) Nomura, K. I.; Kalia, R. K.; Nakano, A.; Vashishta, P.; van Duin, A. C. T; Goddard, W. A., III. Phys. Rev. Lett. 2007, 99, 148303. (14) Reed, E. J.; Fried, L. E.; Joannopoulos, J. D. Phys. Rev. Lett. 2003, 90, 235503. (15) Reed, E. J.; Fried, L. E.; Manaa, M. R.; Joannopoulos, J. D. In Chemistry at Extreme Conditions; Manaa, M. R., Ed.; Elsevier: New York, 2005; pp 297−326. (16) Manaa, M. R.; Reed, E. J.; Fried, L. E.; Galli, G.; Gygi, F. J. Chem. Phys. 2004, 120, 10146. (17) Fidel, C. M; Amar, M. K.; Michael, F. R., Jr.; van Duin, A. C. T.; Jonathan, P. M. Combust. Flame 2012, 159, 1272. (18) Reed, E. J.; Manaa, M. R.; Fried, L. E.; Glaesemann, K. R.; Joannopoulos, J. D. Nat. Phys. 2008, 4, 72. (19) Manaa, M. R.; Reed, E. J.; Fried, L. E.; Goldman, N. J. Am. Chem. Soc. 2009, 131, 5483. (20) Plimpton, S. J. J. Comput. Phys. 1995, 117, 1. (21) Aktulga, H. M.; Fogarty, J. C.; Pandit, S. A.; Grama, A. Y. Parallel Comput. 2012, 38, 245. (22) Guo, F.; Cheng, X. L.; Zhang, H. J. Phys. Chem. A 2012, 116, 3514. (23) Haile, J. M. Molecular Dynamics Simulation: Elementary Methods; Wiley: New York, 1992. (24) Huang, P. H.; Lai, H. Y. Phys. Rev. B 2008, 77, 125408. 24374

dx.doi.org/10.1021/jp4072795 | J. Phys. Chem. C 2013, 117, 24368−24374