Letter pubs.acs.org/JPCL
ReaxFF Reactive Force-Field Modeling of the Triple-Phase Boundary in a Solid Oxide Fuel Cell Boris V. Merinov,*,† Jonathan E. Mueller,†,§ Adri C. T. van Duin,‡ Qi An,† and William A. Goddard, III† †
Materials and Process Simulation Center, California Institute of Technology, 1200 East California Boulevard, m/c 139-74, Pasadena, California 91125, United States ‡ Department of Mechanical and Nuclear Engineering, Pennsylvania State University, 136 Research Building East, University Park, Pennsylvania 16802, United States S Supporting Information *
ABSTRACT: In our study, the Ni/YSZ ReaxFF reactive force field was developed by combining the YSZ and Ni/C/H descriptions. ReaxFF reactive molecular dynamics (RMD) were applied to model chemical reactions, diffusion, and other physicochemical processes at the fuel/Ni/YSZ interface. The ReaxFF RMD simulations were performed on the H2/Ni/YSZ and C4H10/Ni/YSZ triple-phase boundary (TPB) systems at 1250 and 2000 K, respectively. The simulations indicate amorphization of the Ni surface, partial decohesion (delamination) at the interface, and coking, which have indeed all been observed experimentally. They also allowed us to derive the mechanism of the butane conversion at the Ni/YSZ interface. Many steps of this mechanism are similar to the pyrolysis of butane. The products obtained in our simulations are the same as those in experiment, which indicates that the developed ReaxFF potential properly describes complex physicochemical processes, such as the oxide-ion diffusion, fuel conversion, water formation reaction, coking, and delamination, occurring at the TPB and can be recommended for further computational studies of the fuel/electrode/electrolyte interfaces in a SOFC. SECTION: Surfaces, Interfaces, Porous Materials, and Catalysis
S
tional modeling is able to contribute significantly toward overcoming these obstacles by elucidating the fundamental mechanisms at work in these undesirable processes, so that rational approaches can be used in developing strategies for overcoming them. Furthermore, the computational methods prove to be very useful for determining the mechanisms of the desired chemical reactions that occur at the interface, for instance, hydrocarbon fuel conversion. In our study, ReaxFF reactive force-field molecular dynamics (MD) were applied to model chemical reactions, diffusion, and other physicochemical processes at the fuel/Ni/YSZ interface. We used a ReaxFF potential (see Supporting Information), which combines two earlier potentials developed and validated to describe the YSZ electrolyte5 and hydrocarbon chemistry catalyzed by Ni,6,7 respectively. In ref 5, we reported the development of the ReaxFF reactive force field for the YSZ oxide-ion conducting electrolyte. It is based on substantial data derived from QM calculations on clusters and periodic systems, such as relevant pure metallic and metal alloy phases, and various bulk and surface oxide systems, including different atomic configurations for the YSZ electrolyte. To validate the use of the ReaxFF developed, we applied it to predict the oxygen ion diffusion coefficient in YSZ as a function of temperature. These values are in excellent agreement with
olid oxide fuel cells (SOFCs) are among the most efficient fuel cells with theoretical efficiencies as high as 85%. Oxideion-conducting yttria-stabilized zirconia (YSZ) ceramics are widely applied as electrolytes in SOFCs, while perovskite-type oxides with good electronic and ionic conductivities and Ni/ YSZ cermet usually serve as the cathode and anode, respectively. SOFCs operate at very high temperature (∼1000 °C) to achieve sufficiently high oxygen conductivity, which creates some advantages and disadvantages. One of the important advantages is that SOFCs do not require precious metal catalysts, making them cheaper than other types of fuel cells. Another advantage is that SOFCs can be powered with a variety fuels, including hydrogen, carbon monoxide, hydrocarbons, and combinations thereof.1,2 In principle, the high operating temperature enables internal reforming reactions of hydrocarbon fuels, which make SOFCs even more attractive. In important respects, Ni is a suitable catalytic material for hydrocarbon reforming in SOFCs because the desired reactions are readily catalyzed on its surface. Unfortunately, if hydrocarbon fuels or other carbon-containing fuels are used, a layer of carbon is usually built up on the Ni catalyst surface due to coke formation, which rapidly breaks down the efficiency of the SOFC catalyst.3 One of the most promising ways to alleviate this problem is alloying Ni with other metals such as Fe.4 Regardless of whether pure Ni or a Ni alloy is used, the different thermal expansion coefficients of the catalytic and electrolyte materials cause another problem, delamination, which occurs at the electrode/electrolyte interface. Computa© XXXX American Chemical Society
Received: September 5, 2014 Accepted: November 4, 2014
4039
dx.doi.org/10.1021/jz501891y | J. Phys. Chem. Lett. 2014, 5, 4039−4043
The Journal of Physical Chemistry Letters
Letter
experimental results,8,9 setting the stage for the use of ReaxFF to model the transport of oxygen ions through the YSZ electrolyte for SOFCs. Because ReaxFF descriptions are already available for the Ni-metal catalyst,6 we can now consider fully first-principles-based ReaxFF MD simulations of the critical functions in SOFCs, enabling the possibility of in silico optimization of these materials. The combined Ni/YSZ ReaxFF was developed by combining the YSZ and Ni/C/H descriptions and defining the missing angular terms (e.g., Ni−O−Zr) from combination rules. Ni−Zr and Ni−Y bonded interactions were ignored; Ni/Y and Ni/Zr nonbonded parameters were obtained from combination rules. Some examples of the ReaxFF versus QM energetics, which show the reliability of the force field, are available in the Supporting Information. The ReaxFF MD simulations were performed on the H2/Ni/ YSZ and C4H10/Ni/YSZ triple-phase boundary (TPB) systems at 1250 and 2000 K, respectively. The former temperature is close to the usual operating temperature of a SOFC, while the later temperature is approximately twice as high. The higher temperature was used to increase the rate at which the potential energy surface is sampled and the sampling of high energy configurations. Thus, while an elevated simulation temperature is necessary for the observation of reactive events within a computationally feasible simulation time, care must be taken in extending the conclusions drawn to systems at lower temperatures. Figure 1 shows the periodic system applied for
Figure 2 shows a snapshot of the H2/Ni/YSZ system following 150 ps of NVT dynamics at 1250 K, which now
Figure 2. Final point of the ReaxFF MD simulation on the H2/Ni/ YSZ triple-phase boundary system.
includes 5 H2O, 51 H atoms on the Ni surface, 22 Ni−O bonds (Ni oxidation), and only 1 Ni−OH group. In addition, substantial reduction of the YSZ electrolyte, which results in the formation of 11 hydroxyl groups on the YSZ surface and 3 at the Ni/YSZ interface, was observed in our model system. Partially, the reduction of the YSZ surface might be due the limited amount of diffusive oxygen available in the model that we used for the TPB. The anodic reaction in a SOFC is water formation: H2 + O2− → H2O + 2e−. We identified two possible locations where water formation takes place in our simulation. The first is the Ni−YSZ interface, where (1) H2 dissociates on Ni and forms Ni−H bonds, then (2) O2− comes from the YSZ electrolyte first to form OH, and finally (3) another H approaches OH to form H2O. The second location where we observe the water formation reaction is the Ni(111) surface. The first three steps at this location are similar to those previously described; however, once H2O has been formed it immediately leaves the Ni surface and moves to the YSZ surface, whereas the H2O molecules formed at the Ni/YSZ interface remain there for at least the duration of our simulation. The water formation mechanism observed in our simulation is in agreement with the mechanisms following from DFT studies10−12 and from solving reaction−diffusion equations on both electrode and electrolyte surfaces.13,14 Our ReaxFF simulation not only unveils the individual mechanistic steps involved in water formation but also brings to light the structural changes that occur in each of the system’s components in operando. As can be seen in Figure 2, the simulation indicates amorphization of the Ni surface and partial decohesion at the interface (so-called delamination), which have indeed all been observed experimentally.15 In fact, the delamination of the Ni/YSZ interface in operando is well known, and one of the possible reasons is the difference in thermal expansion coefficients between Ni (13.3 × 10−6/C) and YSZ (10 × 10−6/C).15 From the crystallographic viewpoint, the distinct atomic structures of Ni and YSZ may be responsible for decohesion at the interface as well. Regardless of the exact contributions of the thermal expansion coefficients and atomic crystal structures, the atomistic description given by our simulation is consistent with experimental the observations. Computational modeling may
Figure 1. Periodic system used for modeling fuel conversion at the butane/Ni/YSZ interface.
modeling fuel conversion at the Ni/YSZ interface (240 Ni, 16 Y, 48 Zr, and 148 O atoms). We used the YSZ electrolyte with dopant concentration of ∼14 mol % Y2O3 [(Y2O3)(ZrO2)6]. The system was first equilibrated at 1250 and 2000 K using ReaxFF NPT dynamics and then we carried out ReaxFF NVT MD simulations at these temperatures on the corresponding TPB systems. The number of fuel molecules and simulation temperature were varied depending on the problem studied to observe fuel conversion reactions within reasonable time period. The reaction barriers in the H2/Ni/YSZ system with 125 H2 molecules enable us to observe numerous reactions within 150 ps of simulation time at the typical operating temperature, 1250 K, while the simulation temperature and time period were 2000 K and 2 ns for the C4H10/Ni/YSZ system with 20 butane molecules. 4040
dx.doi.org/10.1021/jz501891y | J. Phys. Chem. Lett. 2014, 5, 4039−4043
The Journal of Physical Chemistry Letters
Letter
According to Purnell et al.,16,19 the mechanism of the nbutane pyrolysis includes the following elementary reactions
help to better understand the mechanisms of electrode delamination in SOFCs, but this is beyond of the scope of our paper. Having shown that our ReaxFF potential properly describes the water formation reaction at the Ni anode of a SOFC and thus allows first-principles-based predictions of the chemical processes at the H2/Ni/YSZ TPB, we now proceed to study a complicated hydrocarbon fuel, namely, butane. The hydrocarbon-fuel−oxidation reactions, which deplete butane to produce water and carbon dioxide at the anode/ electrolyte interface of a SOFC, can be put into the following general form
C4 H10 → 2C2H5•
C4 H10 →
where the last term represents carbon deposits that cause coking of the anode. To simulate these chemical reactions at the butane/Ni/YSZ TPB, we replaced the H2 fuel used in the previous simulations with 20 butane molecules. Because the reaction rate of butane is significantly slower than that of H2, we were forced to perform these simulations at an elevated temperature (2000 K) to observe significant reaction levels within a practical simulation time. At the end of our 2 ns simulation, we observed all major productsmethane, ethane, ethylene, and propylenefound experimentally for the pyrolysis of n-butane.16−18 As can be seen in Figure 3, the number of butane molecules decreases dramatically during the first 500 ps before
•
(2)
CH3• + C4 H10 → CH4 + C4 H 9•
(3)
C2H5• + C4 H10 → C2H6 + C4 H 9•
(4)
H• + C4 H10 → H 2 + C4 H 9•
(5)
C4 H 9 → C3H6 +
∑ Cn
(1)
+ C3H 7
•
mC4 H10 + 7mO2 − → 5mH 2O + mCO2 + 14me− +
CH3•
CH3•
(6)
C4 H 9• → C2H4 + C2H5•
(7)
C2H5• → C2H4 + H•
(8)
2C2H5• → C4 H10
(9)
2C2H5• → C2H4 + C2H6
(10)
2CH3•
→ C2 H6
(11)
C2H5• + CH3• → C3H8
(12)
In our simulation, we observe all of these reactions except for reactions 3 and 4. Instead, we observe reactions: CH3• + C3H5• → CH4 + C3H4 and C2H5• + C3H5• → C2H6 + C3H4, in which similar products, CH4 and C2H6, are formed. In addition, other species, such as C4H8, C4H7, C4H6, C3H5, C2H3, C2H2, C2H, CH4, CH3, CH2, H, C, C2, C3, C4, C5, C6, C7, OH, H2O, CH2O, CH3O, CH4O, C4H10O, and C4H9O, and reactions were found as well C4 H10 → C4 H 9• + H•
C2H4 → C2H3• + H•
CH3• → CH 2 + H•
C4 H 9• → C4 H8 + H•
C2H3• → C2H 2 + H•
CH3• + H• → CH4
C4 H8 → C4 H 7• + H•
C 2 H 2 → C 2 H• + H•
C4 H 9• + H• → C4 H8 + H 2 C4 H 7• → C4 H6 + H• C 2 H• → C 2 + H•
C4 H 7• + H• → C4 H6 + H 2
C4 H6 → C4 H5• + H•
C3H 7• → CH3• + C2H4
C4 H 7• + H• → C4 H6 + H 2 C4 H5• → C2H4 + H• C3H 7• → C3H5• + H 2
C3H5• + C2H5• → C2H6 + C3H4
and so on. Many of these reactions are due to hydrogen abstraction, and some of the species are intermediates for the products shown in Scheme 1. In fact, in line with a previous study,20 a tabulation of the frequency of each elementary reaction (see Supporting Information) suggests that hydrogen abstraction occurs significantly more readily and frequently than C−C bond cleavage, so that a sizable fraction of the adsorbed fuel is converted into carbon chains. Furthermore, we observe these carbon chains reacting with each other to form longer chains, which eventually clog the surface. This suggests that the preferential cleavage of C−H before C−C bonds and the facile C−C bond formation and cleavage in C chains on Ni are important factors in coking. Thus, modifying the metal catalyst by means of alloying, to increase the ease with which C−C bonds between partially hydrogenated carbon atoms take place
Figure 3. Populations of reactant, intermediate, and product species during 2 ns simulation of butane conversion at the Ni/YSZ interface.
temporarily leveling off, primarily due to the disappearance of empty surface sites and perhaps due to decreasing the number of butane molecules, which results in a lower fuel gas pressure. Indeed, fuel-decomposition reactions that take place over the entire course of the simulation show that primarily the adsorption, and not decomposition, of C4H10 is inhibited. After 2 ns of NVT dynamics, 12 of the 20 C4H10 have been consumed to form 10 H2O (4 in the gas phase), 54 OH (most of them on the YSZ surface), 18 H, 5 H2 (all in the gas phase), 2 CH4 (both gas phase), 1 C2H6 (gas phase), 1 C2H4 (gas phase), and 15 Cn (from C to C8) species. 4041
dx.doi.org/10.1021/jz501891y | J. Phys. Chem. Lett. 2014, 5, 4039−4043
The Journal of Physical Chemistry Letters
Letter
⎡ ⎛ rij ⎞ pbo2 ⎤ BOij = BOijσ + BOijπ + BOijππ = exp⎢pbo1 ⎜ σ ⎟ ⎥ ⎢⎣ ⎝ r0 ⎠ ⎥⎦
Scheme 1. Mechanism of the Butane Conversion at the Ni/ YSZ Interface
⎡ ⎡ ⎛ rij ⎞ pbo4 ⎤ ⎛ rij ⎞ pbo6 ⎤ + ⎢pbo3 ⎜ π ⎟ ⎥ + exp⎢pbo5 ⎜ ππ ⎟ ⎥ ⎢⎣ ⎢⎣ ⎝ r0 ⎠ ⎥⎦ ⎝ r0 ⎠ ⎥⎦
Overcoordination and undercoordination energy penalties are then used to enforce the correct bond order. The total system energy is a sum of a several partial energy terms; these include energies related to lone pairs, undercoordination, overcoordination, valence and torsion angles, conjugation, hydrogen bonding, and van der Waals and Coulomb interactions. Thus, the total energy can be expressed as Esystem = E bond + E lp + Eover + Eunder + Eval + Epen + Ecoa + EC2 + Etors + E H‐bond + EvdW + Ecouloumb
Coulomb and van der Waals interactions are calculated between every pair of atoms. This allows ReaxFF to describe not only covalent bonds but also ionic bonds and the whole range of intermediate interactions. Charge distributions are calculated based on geometry and connectivity using the electronegativity equalization method (EEM).25 Coulomb interactions are treated using a seventh-order spline (Taper function).26 To keep ReaxFF from erroneously predicting a strong triple bond in C2, an additional partial energy contribution is utilized as previously reported.7
so that CHx species that can migrate to the oxide surface are produced instead of Cx chains, which tend to remain on and eventually block off the metal surface, is a potential strategy for avoiding coking. We consider the CH2O, CH3O, CH4O, and similar species, which we do observe in our simulations, as precursors for CO and CO2 formation that we do not observe. The failure to observe either CO or CO2 is most likely due to the limited amount of oxygen readily available from the YSZ electrolyte, which is also needed in the competing water formation reaction (ten H2O molecules are produced). A closer look at the presence/absence of oxygen in the various reactions observed confirms that butane decomposition takes place primarily on Ni, while the oxidation and water formation would preferentially take place on the oxide’s surface and interface, where oxygen is more readily available. The spatial separation between the primary locations for decomposition and oxidation provides a feasible explanation for why additional simulation time and oxygen may be necessary to see further progress toward complete oxidation. Scheme 1 shows the mechanism of the butane conversion at the Ni/YSZ interface derived from our simulations. Many steps of this mechanism are similar to the pyrolysis of butane.19 The products (shown in bold in Scheme 1) obtained in our simulations are the same as in experiment.17,21 This indicates that the developed ReaxFF potential reproduces well such experimental data as the oxide-ion diffusion coefficient and products of fuel (hydrogen and butane) conversion, water formation reaction, coking, and delamination occurring at the TPB and can be recommended for further computational studies of the fuel/electrode/electrolyte interfaces in a SOFC.
■
■
ASSOCIATED CONTENT
S Supporting Information *
Tables of reaction counts for hydrogen abstraction and C−C bond cleavage/formation reactions observed in the butane/Ni/ YSZ simulation; ReaxFF reactive force field used in our simulation; and Figures with examples of ReaxFFvsQM energetics. This material is available free of charge via the Internet at http://pubs.acs.org.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Present Address §
J.E.M.: Institute of Electrochemistry, University of Ulm, Albert-Einstein-Allee 47, D-89081 Ulm, Germany.
Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was partially supported by the Department of Defense Multidisciplinary University Research Initiative (MURI) program administered by the Office of Naval Research under grant N00014-02-1-0665. The facilities of the Materials and Process Simulation Center used in this study were established with grants from DURIP-ONR, DURIP-ARO, and NSF-CSEM.
COMPUTATIONAL METHOD
The unique characteristic of reactive force fields is their ability to model bond breaking and formation. This allows for descriptions of such processes as chemical reactions, diffusion, phase transitions, and so on. The ReaxFF reactive force field22 uses the bond order/bond distance relationship introduced by Tersoff23 and first applied to carbon chemistry by Brenner to describe chemical reactivity.24 Bond orders, summed from σ, π, and ππ terms, are calculated instantaneously from interatomic distances as follows
■
REFERENCES
(1) Zhu, H.; Kee, R.; Janardhanan, V.; Deutschmann, O.; Goodwin, D. J. Modeling Elementary Heterogeneous Chemistry and Electrochemistry in Solid-Oxide Fuel Cells. Electrochem. Soc. 2005, 152, A2427−A2440. (2) Ni, M.; Leung, M. K. H.; Leung, D. Y. C. Ammonia-Fed Solid Oxide Fuel Cells for Power Generation-A Review. Int. J. Energy Res. 2009, 33, 943−959.
4042
dx.doi.org/10.1021/jz501891y | J. Phys. Chem. Lett. 2014, 5, 4039−4043
The Journal of Physical Chemistry Letters
Letter
(3) Clarke, S.; Dicks, A.; Pointon, K.; Smith, T.; Swann, A. Catalytic Aspects of the Steam Reforming Of Hydrocarbons in Internal Reforming Fuel Cells. Catal. Today 1997, 38, 411−423. (4) Hecht, E. S.; Gupta, G. K.; Zhu, H.; Dean, A. M.; Kee, R. J.; Maier, L.; Deutschmann, O. Methane Reforming Kinetics within a NiYSZ SOFC Anode Support. Appl. Catal., A 2005, 295, 40−51. (5) van Duin, A. C. T.; Merinov, B. V.; Jang, S. S.; Goddard, W. A., III. ReaxFF Reactive Force Field for Solid Oxide Fuel Cell Systems with Application to Oxygen Ion Transport in Yttria-Stabilized Zirconia. J. Phys. Chem. A 2008, 112, 3133−3140. (6) Mueller, J. E.; van Duin, A. C. T.; Goddard, W. A., III. Development And Validation of ReaxFF Reactive Force Field for Hydrocarbon Chemistry Catalyzed by Nickel. J. Phys. Chem. C 2010, 114, 4939−4949. (7) Nielson, K. D.; van Duin, A. C. T.; Oxgaard, J.; Deng, W. Q.; Goddard, W. A. Development of the ReaxFF Reactive Force Field for Describing Transition Metal Catalyzed Reactions, with Application to the Initial Stages of the Catalytic Formation of Carbon Nanotubes. J. Phys. Chem. A 2005, 109, 493−499. (8) Kilo, M.; Argirusis, C.; Borchardt, G.; Jackson, R. Oxygen Diffusion in Yttria Stabilised Zirconia - Experimental Results and Molecular Dynamics Calculations. Phys. Chem. Chem. Phys. 2003, 5, 2219−2224. (9) Minh, N. Q. Ceramic Fuel Cells. J. Am. Ceram. Soc. 1993, 76, 563−588. (10) Cucinotta, C. S.; Bernasconi, M.; Parrinello, M. Hydrogen Oxidation Reaction at the Ni/YSZ Anode of Solid Oxide Fuel Cells from First Principles. Phys. Rev. Lett. 2011, 107, 206103. (11) Shishkin, M.; Ziegler, T. Direct Modeling of the Electrochemistry in the Three-Phase Boundary of Solid Oxide Fuel Cell Anodes by Density Functional Theory: A Critical Overview. Phys. Chem. Chem. Phys. 2014, 16, 1798−1808. (12) Ammal, S. C.; Heyden, A. Combined DFT and Microkinetic Modeling Study of Hydrogen Oxidation at the Ni/YSZ Anode of Solid Oxide Fuel Cells. J. Phys. Chem. Lett. 2012, 3, 2767−2772. (13) Goodwin, D. G.; Zhu, H.; Colclasure, A. M.; Kee, R. G. Modeling Electrochemical Oxidation of Hydrogen on Ni−YSZ Pattern Anodes. J. Electrochem. Soc. 2009, 156, B1004−B1021. (14) Yao, W.; Croiset, E. Modelling and Ni/Yttria-Stabilized-Zirconia Pattern Anode Experimental Validation of a New Charge Transfer Reactions Mechanism for Hydrogen Electrochemical Oxidation on Solid Oxide Fuel Cell Anodes. J. Power Sources 2014, 248, 777−788. (15) Minh, N. Q. Science and Technology of Ceramic Fuel Cells; Elsevier: Amsterdam, 1995. (16) Purnell, J. H.; Quinn, C. P. Pyrolysis of n-Butane. Proc. R. Soc. London, Ser. A 1962, 270, 267−284. (17) Sheng, C. Y.; Dean, A. M. Importance of Gas-Phase Kinetics within the Anode Channel of a Solid-Oxide Fuel Cell. J. Phys. Chem. A 2004, 108, 3772−3783. (18) Torok, J.; Sandler, S. Kinetics of Pyrolysis of n-Butane. Can. J. Chem. 1969, 47, 3863−3869. (19) Hughes, D. G.; Marshall, R. M.; Purnell, J. H. Rate Constants for Initiation of n-Butane Pyrolysis and for Recombination of Ethyl Radicals. J. Chem. Soc., Faraday Trans. 1974, 70, 594−599. (20) Mueller, J. E.; van Duin, A. C. T.; Goddard, W. A., III. Application of the ReaxFF Reactive Force Field to Reactive Dynamics of Hydrocarbon Chemisorption and Decomposition. J. Phys. Chem. C 2010, 114, 5675−5685. (21) Mallinson, R. G.; Braun, R. L.; Westbrook, C. K.; Burnham, A. K. Detailed Chemical Kinetics Study of the Role of Pressure in Butane Pyrolysis. Ind. Eng. Chem. Res. 1992, 31, 37−45. (22) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396−9409. (23) Tersoff, J. Empirical Interatomic Potential for Carbon, with Applications to Amorphous-Carbon. Phys. Rev. Lett. 1988, 61, 2879− 2882.
(24) Brenner, D. W. Empirical Potential for Hydrocarbons for Use in Simulating the Chemical Vapor-Deposition of Diamond Films. Phys. Rev. B 1990, 42, 9458−9471. (25) Mortier, W. J.; Ghosh, S. K.; Shankar, S. Electronegativity Equalization Method for the Calculation of Atomic Charges in Molecules. J. Am. Chem. Soc. 1986, 108, 4315−4320. (26) van Duin, A. C. T.; Strachan, A.; Stewman, S.; Zhang, Q.; Xu, X.; Goddard, W. A. ReaxFF(SiO) Reactive Force Field for Silicon and Silicon Oxide Systems. J. Phys. Chem. A 2003, 107, 3803−3811.
4043
dx.doi.org/10.1021/jz501891y | J. Phys. Chem. Lett. 2014, 5, 4039−4043
