Evaluating Water Reactivity at Silica Surfaces Using Reactive Potentials

Subscriber access provided by University of Sussex Library

C: Surfaces, Interfaces, Porous Materials, and Catalysis

Evaluating Water Reactivity at Silica Surfaces Using Reactive Potentials Thiruvilla S Mahadevan, and Jincheng Du J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b12653 • Publication Date (Web): 19 Apr 2018 Downloaded from http://pubs.acs.org on April 19, 2018

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 29 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

Evaluating Water Reactivity at Silica Surfaces using Reactive Potentials Thiruvilla S Mahadevan, Jincheng Du* Department of Materials Science and Engineering, University of North Texas, Denton, Texas76207 (*Corresponding author. Email: [email protected]) Abstract: Understanding the interactions between amorphous silica surfaces and water provides insight into material degradation of glassy materials in aqueous environment such as corrosion of nuclear waste glasses. Molecular dynamics (MD) simulations of water and nanometer sized silica structures were used in this work to evaluate the reactivity of silica flat and surfaces with curvature. We compared two dissociable water/silica potentials, namely the Reactive Force Field (ReaxFF) and the Mahadevan-Garofalini water/silica force field (MGFF) that have been in development over the last decade, to study their performance in simulating bulk water as well as silica-water interactions. Significant differences in the properties of bulk water that had a bearing on the surface interactions were observed between the two types of potentials, as well in the same potential type with two parameterization for ReaxFF, suggesting need of improvement of the existing water/silica ReaxFF potentials. Our simulation results show that a majority of the silanols were formed by reaction between water and a strained siloxane bonds that mainly exist on the surface of amorphous silica, within a few nanoseconds of the simulation time scale, in agreement with previous studies. Effect of surface curvature on the reactivity with water was investigated. Our results indicate that defect concentration at the surface bears a strong correlation to the concentration of silanols (SiOH) that eventually form. We observe undercoordinated Si at the surface that are attacked by water before the hydrolysis reaction of the siloxane bonds and demonstrate possible mechanisms of water reacting with these undercoordinated Si. We also find that the method of generating surfaces in simulation determines the defect concentration and hence influences the reactivity of the amorphous silica surface.

1. Introduction Water and silica are among the most abundant compounds on the earth’s crust and surface and the study of interactions between the two have significant implications in geosciences, glass technology, nuclear sciences, biomaterials and biological systems. As a consequence, there is a scientific interest in understanding the interactions between the two at wide range of length and timescales and research in this field has been in progress for several decades now1,2. Computational modeling of these interactions has been done in a variety of methods3 including ab initio, DFT, and semi-empirical Quantum Mechanical (QM) methods in the atomistic and 1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 29

angstrom scale4–7, classical MD simulations up to multiple8–10 nanometer scale, and greater lengths through multiscale methods11–13. Several empirical potentials for classical MD simulations have been used in the study of silica14,15 that sufficiently reproduce the surface and bulk structure of amorphous silica. The potentials used for silica simulation are all atomistic and they allow for breaking of bonds between atoms and diffusion of ionic species. There is an even greater number of models for classical potentials for simulations of water16 (16 and references therein) – however most of these water potentials were designed for simulation of a narrow range of physical properties of water. An important aspect of silica-water interactions is the hydrolysis reaction that leads to the formation of silanols and dissociation of a water molecules which would control several other chemical interactions in the system1,4,17. This chemical reaction can be reproduced by either ab-initio QM calculations4,18,19 or by combining classical MD potentials and QM calculations6,12,20,21 but such simulations are limited in size due to the computational complexity involved. Classical MD simulations using reactive potentials enable study of the effect of medium range structure features and relaxation due to surface and bulk effect. They also provide trajectories up to nanometer time scale to observe diffusion associated reaction processes. As the silica-water reaction involves physisorption and chemisorption of water on the silica surfaces, hydrolysis reactions, and proton transfer that is commonly observed to mediate various water involved reactions, it is desirable for the water potential to be a flexible water molecule capable of dissociation and description of chemical reactions with silica. While early models of such dissociative water based on pairwise potentials and Coulomb interactions with partial charges8,22 were capable of reproducing some of the surface interactions, recent reactive models23,24 have been paramterized to reproduce the dynamics of hydroxylated silica surfaces25, the long range structure and properties of water as well as silica water interactions. The Reactive Force Field (ReaxFF) is capable of describing chemical reactions and finds wide applications in simulations of reactive systems and systems with heterogeneous interfaces. ReaxFF was introduced in 200123 which allowed for atomistic depiction of all the species. Oxygen and Hydrogen were parameterized from the earliest version in 2003, and were categorized as parameters for either aqueous (reactive water) or combustion (hydro carbons) reactions26. ReaxFF water was parameterized for silica water interactions by Fogarty et.al27 which was later refined and re-parameterized by Larsson28 et. al and Yeon29et. al. Several of these ReaxFF parameters have been tested earlier for Silica water interactions10,30 , water interaction in silica gels31 and have also been compared to other potentials32. Dissociative, reactive potential by Mahadevan and Garofalini24 was initially designed in 2007 studying the expansion of water in nanopores. It is based on an interaction distance dependent charge scheme similar to the water potential by Guillot and Guissani although the water molecules

2 ACS Paragon Plus Environment

Page 3 of 29 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


were rigid in their model16. This has hence been used to study several reactions in silica-water systems33 and aqueous sodium chloride with silica34. In general, while several potentials are available in literature for classical MD simulations of silica-water interactions, very few of them can reproduce the important dissociation reactions of water. This is because, most of the modeuls of water use either a rigid, or at best, a flexible, nondissociable, molecular structure. Very few of the models mentioned24,28,29,35 above are flexible and dissociable besides being able to reproduce some of the bulk properties of water and can be termed as reactive potentials. Thus, one of the purposes of this work is to evaluate the suitability of some of these potentials for simulation of large systems and complex geometries. Amorphous silica surfaces have atomic36 and nanometer scale roughness that affect its their reactivity37. Since water diffusing in silica usually occurs in constraining geometries that have surface curvatures like in silica gels38,39, there is a need to evaluate silica water interactions in these conditions. To this end, concave and convex surfaces were exposed to water and their reactivities were evaluated.

2. Overview of empirical potentials Among the available dissociable water potentials, we investigate the Reactive Force Field (ReaxFF), developed by van Duin et.al23 and the dissociable reactive potential of MahadevanGarofalini24. The details of the potential function and its implementation are given in the authors’ original articles. In brief, the ReaxFF is based on consideration of several multibody interactions and the overall system energy (for silica water systems) is described as: 𝐸𝐸𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 =

(𝐸𝐸𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 + 𝐸𝐸𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣 + 𝐸𝐸𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡 +𝐸𝐸𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 + 𝐸𝐸𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢 + 𝐸𝐸𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙 )

(𝐸𝐸𝐻𝐻−𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 + 𝐸𝐸𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣𝑣 + 𝐸𝐸𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶𝐶 )

𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵 𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼

𝑁𝑁𝑁𝑁𝑁𝑁−𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵 𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼𝐼

+

where the subscripts correspond to bond, valence-angle, torsion-angle, overcoordination, undercoordination, and lone pair for bonded interactions and Hydrogen bond, Van der Waals, Coulomb terms for non-bonded interactions. The bonded interactions are calculated only at the “intra-molecular” level or between atoms directly connected by a bond and are determined by the bond order. In ReaxFF, the uncorrected bond order is expressed as a function of the bond distance as: 𝐵𝐵𝐵𝐵′𝑖𝑖𝑖𝑖 =

�

𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡

exp(𝑝𝑝𝑏𝑏𝑏𝑏,1 ∙ �

𝑟𝑟𝑖𝑖𝑖𝑖 𝑃𝑃𝑏𝑏𝑏𝑏,2 � ) 𝑟𝑟𝑜𝑜

The above uncorrected bond order is corrected for under and over coordination and is used to calculate the bond energy as: 3 ACS Paragon Plus Environment


Page 4 of 29

𝑝𝑝

𝐸𝐸𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏 = −𝐷𝐷𝑒𝑒 ∙ 𝐵𝐵𝐵𝐵𝑖𝑖𝑖𝑖 ∙ exp �𝑝𝑝𝑏𝑏𝑏𝑏,1 (1 − 𝐵𝐵𝐵𝐵𝑖𝑖𝑖𝑖𝑏𝑏𝑏𝑏,1 )�

All the bonded interactions are calculated as a function of the bond order. The unbonded interactions also involve a charge on the atoms, which is updated using charge equilibration (QEq) methods. The details of the functions and parameters are given in references23,35,40 and a complete description is beyond the scope of this work. Overall, the potential function involves over 20 parameters for each atom type and more for bond types, angle types and dihedral angles. The parameters are calibrated based on results of QM simulations and several sets of such ReaxFF parameters exist for water. In this work, we investigate 2 of the recent parameter sets, namely by Yeon29and Larsson28. Both these parameter sets were obtained by different optimization of the set by Fogarty et. al.10,27 developed for water-silica systems. The Mahadevan-Garofalini potential (MGFF) has a much simpler analytical form describing pair and 3 body forces. While the analytical form is based on charge varying as a function of the distance16, the potential was reparametrized to simulate a dissociable atomistic water that accurately reproduced physical properties of water such as structure, density self-diffusion coefficients over a wide range of temperatures and pressures24. The total energy as a function of the distance between atoms (rij) is given as:

where

𝐸𝐸2𝐵𝐵𝐵𝐵𝐵𝐵𝐵𝐵 = 𝐸𝐸𝑞𝑞𝑞𝑞 + 𝐸𝐸𝑞𝑞𝑑𝑑 𝑞𝑞𝑑𝑑 + 𝐸𝐸𝑞𝑞𝑞𝑞𝑑𝑑 + 𝐸𝐸𝑞𝑞𝑑𝑑 𝑞𝑞 + 𝐸𝐸𝑟𝑟𝑟𝑟𝑟𝑟 + 𝐸𝐸𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 𝐸𝐸𝑞𝑞𝑞𝑞 = 𝐸𝐸𝑞𝑞𝑞𝑞𝑑𝑑 = 𝐸𝐸𝑞𝑞𝑑𝑑 𝑞𝑞 = 𝐸𝐸𝑞𝑞𝑑𝑑 𝑞𝑞𝑑𝑑 =

𝑖𝑖𝑖𝑖

𝑞𝑞𝑖𝑖 𝑞𝑞𝑗𝑗 𝑟𝑟𝑖𝑖𝑖𝑖 erfc � � 𝑟𝑟𝑖𝑖𝑖𝑖 𝛽𝛽

𝑞𝑞𝑖𝑖 𝑞𝑞𝑗𝑗 𝑟𝑟𝑖𝑖𝑖𝑖 𝑟𝑟𝑖𝑖𝑖𝑖 erf � � erfc � � 4𝑟𝑟𝑖𝑖𝑖𝑖 𝛽𝛽 √2𝜉𝜉𝑖𝑖𝑖𝑖

𝑞𝑞𝑖𝑖 𝑞𝑞𝑗𝑗 𝑟𝑟𝑖𝑖𝑖𝑖 𝑟𝑟𝑖𝑖𝑖𝑖 erfc � � erfc � � 16𝑟𝑟𝑖𝑖𝑖𝑖 2𝜉𝜉𝑖𝑖𝑖𝑖 𝛽𝛽

𝐸𝐸𝑟𝑟𝑟𝑟𝑟𝑟 = 𝐴𝐴𝑟𝑟𝑟𝑟𝑟𝑟

erfc(𝑧𝑧𝑖𝑖𝑖𝑖 ) 𝑟𝑟𝑖𝑖𝑖𝑖 �with 𝑧𝑧𝑖𝑖𝑖𝑖 = � 𝑧𝑧𝑖𝑖𝑖𝑖 2𝜉𝜉𝑖𝑖𝑖𝑖

𝐸𝐸𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 =

𝑖𝑖𝑖𝑖

−𝐶𝐶6 𝑟𝑟𝑖𝑖𝑖𝑖6

Complete details and values of the parameters are given in the original paper24. The above formulation includes the Wolf attenuation factor erfc(rij/β) applied to all the Coulumbic terms to handle the long range interactions hence provides high computational efficiency. The 3 body potential is also a function of the bond distances and angles and details of all the parameters are given in the original paper. In comparison with ReaxFF, this MGFF has 12 parameters requiring 4 ACS Paragon Plus Environment

Page 5 of 29 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


calibration for water, including the 3 body terms. Parameter sets for MGFF have also been developed for the silica-water systems and have been demonstrated to accurately predict diffusion and reactivity between water and silica surfaces41.

3. Results 3.1 Structure and properties of bulk water One of the important criteria for considering the suitability of an MD potential for water is the extent to which it reproduces the structure of the system. A simple way of evaluating the structure is through comparing the pair/radial distribution function and coordination of species with experimental data. In the case of water, the molecular structure can be obtained from the OxygenOxygen pair distribution function and to a lesser extent the other OH and HH distribution functions. The O-O pair distribution function of water simulated from the three potentials (two versions of ReaxFF and MGFF) that were studied in this work are in good agreement with experimental distribution functions obtained through scattering experiments42,43 (Figure 1). The radial distribution functions were obtained over 100 ps of an ensemble of 1000 water molecules in an NVE ensemble following an NPT to equilibrate at the right temperature (298K) and densities. While both the ReaxFF potentials showed higher peaks and valleys, these deviations smoothen out at longer distances indicating relatively similar bulk structures. If the peaks can be considered to be the a direct indicator of the location and distribution of the shell of neighboring atoms, then all the models are consistent with the experimentally verified structure with respect to molecular coordination up to the first shell of neighbors. The differences in the structure arises in the short range density of the neighboring molecules and this also results in differences in the overall density as will be shown later in this article. More significant differences between the potentials were found in the densities and diffusion coefficients of bulk water measured from simulations. Figure 2 shows the measured densities as a function of temperature compared to the experimental equation of state of water obtained from steam tables44. While the 298K densities were also measured for a 1000 molecule system under NPT and NVT conditions, the densities at higher temperature was measured only through the Gibbs ensemble MD simulations detailed in45, where bulk water is equilibrated at a given temperature with a box of vacuum added on top of it. In the present work, a simulation box of 30 nm side was first equilibrated at NPT for 100 ps at the given temperature and then further NVE equilibrated for another 100 ps to check for thermal stability. NVT equilibration for a further 100 ps with the expanded box was carried out to calculate the density profile. The experimental equation of state curve shown is from the NIST website46. At the time of this work, a new set of ReaxFF parameters for water was published35 and we have done comparison with the two existing ReaxFF parameters chosen in this work. The equation of state for this new set parameter set shows 5 ACS Paragon Plus Environment


Page 6 of 29

much better agreement to the experimental values of water while maintaining similar structural features. However, no compatible silica potential was included in the potential development which makes direct comparison of water/silica reaction impossible. Hence comparison below sticks to the two set of aforementioned water/silica ReaxFF potentials. The Mean Square displacement of Oxygen molecules was calculated for a 1000 molecule system in a cubic box of side ~3.1nm over 300 ps and averaged overs sections of 100ps with NVE conditions following NPT. The Mean Square Displacement (MSD) vs time plots are shown in the supplemental information. The box size was kept sufficiently larger than the interaction curoff distance to ensure independence of diffusion from system size47. The time scale was chosen to be sufficiently large to overcome discrepancies induced by atomic vibrations and to have the MSD measured in the linear regime. The self-diffusion coefficient of oxygen in water was found to be 2.45 x10-9m2/s for MGFF, 2.86 x10-9m2/s for Yeon et. al. and 3.78 x10-9m2/s for Larsoon et. al. ReaxFF parameters. Experimentally, the diffusion coefficient is most widely reported to be around 2.3x10-9m2/s48 at 298K, although values of 2.13x10-9m2/s49 -2.57x10-9m2/s50 were also reported. Potential type and parameters have a significant influence51 in the exact value of diffusion coefficient measured in MD simulations and hence there is a large variation in the reported values. Considering this, the diffusion coefficient obtained with MGFF (2.45x10-9m2/s) is reasonably close to the experimental value while both the ReaxFF parameter sets yielded higher coefficients than experiments. 3.2 Simulations of molecular clusters To understand the origin of the differences in the bulk properties, simulations of smaller clusters of water were done. Stabilized water clusters containing 2-9 molecules, the coordinates of which were obtained from Shields et. al.52 and a single molecule of water were placed a 31Åx31Åx31Å box and equilibrated at 298K (NVT) to obtain energies. The difference between the average energy of the clusters and that of a single molecule of water over 10000 steps was used as the average cohesive energy of the cluster. Similarly, the difference between average energy of stabilized bulk water and that of the single molecule was also calculated as the cohesive energy of water. The calculated cohesive energy for water with MGFF is -11.12 kCal/mol which is close to the theoretical value -11.23 kCal/mol expected for flexible potentials that include intramolecular interactions16,24. The cohesive energies calculated with ReaxFF were both -8.46 kCal/mol. The energies of the molecular clusters at 298 K are plotted in figure 3. In general, for systems with negative cohesive energies, larger clusters have to be more stable and this is demonstrated by the energies decreasing as a function of cluster size in all the three models. Since MD simulations only approximate the quantum effects through the potential function, the energies obtained will not match those of DFT methods although the configurations seen with MGFF have energies closest to DFT studies52.


Page 7 of 29 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


While a complete analysis of set of cluster configurations is not in scope of this work, figure 3b shows an example of the stable structure of a 5 molecule cluster with MGFF and Reax-Yeon. The starting 5 molecule cluster was taken from one form of ice with 4 tetrahedrally arranged water molecules hydrogen bonded to a central water molecule. With all MD simulations, this structure was found to be stable only at low temperatures – at higher temperatures the stable structure is a cyclically bonded ring of 5 water molecules which is consistent with DFT calculations52. This low temperature 5 water cluster helps demonstrate that the first neighbor distance between water molecules is nearly the same with all the potential sets. This is also consistent with the pair distribution function shown in figure 1 - all the first shell of neighbors are approximately at 2.76 Å from the central water molecule. However, there is a consistent difference in the way the hydrogen bond donating neighbors themselves around the central molecule. This difference results in allowing space for a higher number of first neighbors at around 2.76 Å and a lower number of neighbors at around 3.2A with both the ReaxFF parameter sets compared to MGFF. This could possibly manifest further at larger distances creating low density structures in the bulk. In line with examining small molecule clusters, the interaction between a single silicic acid molecule and a single water molecule was also studied to understand the difference in the potentials. For this simulation, a single molecule of Si(OH)4 was added to a box of side 31 Å and one molecule H2O were added to the box 11A away from the Si(OH)4. The system was equilibrated at 298 K and the water molecule was brought incrementally closer to the silicic acid molecule over 1ns. The silicon atom was fixed in its original location but all the OH’s were allowed to move around it thus allowing for all possible configurations of Si(OH)4. Figure 4 shows a plot of the normalized energy of the system as a function of the separation distance between the silicon atom and the oxygen of the water. Clearly, a minimum energy configuration that is stable over a wider range forms with MGFF while both the ReaxFF parameter sets show a wide fluctuation in their energies. When the water molecule is forced even closer to Si(OH)4 (rSi-Ow < 2 Å), the multiple hydrogen bonded configuration that forms with MGFF results in the expulsion of the OH group located furthest across from the incoming H2O and a formation of new SiOH bond with the water molecule. While this reaction occurs with all the three potentials studied, the energy barrier that the water molecule needs to overcome to get as close to Si(OH)4 is higher in both the ReaxFF parameter sets – thus indicates a lower likelihood of silanol formation in bulk silica-water systems. 3.3 Water interactions with flat silica surfaces To evaluate the interactions between silica glass and water, bulk silica glass structures of 3000 atoms in a 35.6Åx35.6Åx35.6Å box were made by a standardized melt quench procedure starting from a crystalline quartz structure as described in an earlier work41. One side (top surface) of the bulk glass was then exposed to a column of vacuum at least twice as big as the original glass box and heated to 1298K with 4 angstroms of the other side (bottom) of the glass rendered immobile. A relaxed glass surface forms when this system was cooled back to 298K. All the melt-quench and 7 ACS Paragon Plus Environment


Page 8 of 29

surface relaxation of glasses were done with MGFF since it was around 12x faster than ReaxFF as explained earlier. Layers of equally spaced water molecules (288 molecules) were added on top of the cooled surface and NVT simulations were done for at least 1ns (10,000,000 steps) to allow interaction between the water molecules and the surface. The simulations with MGFF took about 27 hours on 16 processors for 1.1ns, while the ReaxFF simulations had to be split into 3 parts of ~0.4ns each of which took 72 hours on 32 processors for 3864 atoms. This saving in time is one of the main advantages of using MGFF in simulation of the water/silica interactions. Figure 8a shows the concentration of silanols formed in the system as a function of time (Fig. 8a left) and a snapshot of the system (Fig. 8a right). The equilibrium silanol concentration is 6.1 /nm2 for MGFF, 4.1/nm2 for Yeon and 3.0/nm2 for Larsson ReaxFF. In this case we define a silanol as an oxygen that is within 2.0Å of a silicon atom and within 1.2Å of a Hydrogen atom. This definition was based on the first Si-O pair distribution peak which had reached a zero value by 2.2Å in all the glasses. Experimental observations have put the total silanol concentration on dried silica surfaces to be in the range of 4.2/nm2 – 5.7nm2 53 while silica gels and other hydrated forms of silica have shown up to 10 silanols/nm2 depending on treatment1. These numbers do not generally take into account silica that is deeper than a few nanometers underneath the surface. While MGFF produces a slightly higher concentration of silanols compared to the Zhurvalev number53, it is still in range considering that this number includes the silanols formed up to 1 nm under the surface. As is common to adsorption processes, the initial ballistic regime of rapid formation54 is common to all the three potentials, however, the ReaxFF parameter sets yield a significantly smaller silanol concentration at the surface. A small number of water molecules are directly adsorbed onto the surface forming an Si-OH2 in the initial 0.2 ns with the MGFF potential and these are shown as the dotted line in figure 5. This is in line with one of the demonstrable mechanisms of silanol formation where the water molecule is chemisorbed onto a silicon atom before dissociating to lend a hydrogen atom to a nearby oxygen. The lower density of water and the decreased stability of the Si(OH)4-H2O complex provide an explanation for the lower surface silanol density for systems simulated with ReaxFF. In the case of simulations with Larsson potentials, water molecules tended to form an agglomerate over the surface without contacting with the surface. Hence, an additional push was given to the water molecules by means of an initial momentum towards the surface. As a result the onset of the ballistic regime was earlier and of shorter duration as seen by the earlier start and steeper initial slope of the silanol concentration curve. The source of the oxygen atoms involved in the silanol can be obtained based on its index since the first 3000 atoms are from glass and the rest are from the added water. If each of the water molecule were to split into exactly one OH or H, then exactly 50% of the oxygens involved in the silanols will be from glass. However the percentage of glass oxygens involved in the silanols evolves to about 55% with MGFF and Reax-Yeon and to about 67% with Reax-Larsson. This 8 ACS Paragon Plus Environment

Page 9 of 29 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


indicates the presence of free OH- or H3O+ or silanol formation by mechanisms other than breaking of a siloxane bond and subsequent hydrolysis). Very rarely do we observe any free H+, although it is possible to find some with lifetimes shorter than 0.1 ps which happens when the hydrogen is between transfers from one O to another. The evolution of the free OH-, defined as any oxygen that is attached to exactly one hydrogen and is greater than 2.2 Å away from any silicon atom, concentration is shown in figure 5b for all the three potentials. The presence of free OH- explains the higher fraction of glass oxygens contributing to silanols. This also highlights the differences in the reactivities between the potentials. There was a negligible amount of H3O+ that also form with shorter lifetimes and only within the initial 0.2 ns. Figure 6 shows the concentration profile of the various oxygen species in the system plotted along the normal to the glass surface. The density profile was measured over the final 0.1ns of the simulation and the speciation of oxygen is defined based on the Si-O longest bond distance of 2.0 Å and O-H bond distance of 1.2 Å. By definition, the silanols are included as glass oxygens. The nominal number density of oxygen in a glass (density 2.2 g/cc) is expected around 44/nm3 and oxygen in water is 33/nm3. MGFF and Reax-Yeon both show the glass density to be the in right range while Reax-Larsson shows a slightly higher value. This is because of the initial ‘push’ given to the water molecules in the system which would otherwise only hover on top of the silica surface without approaching reacting distances. This increased density also manifests as an earlier decline in the concentration of Glass oxygens in the profile. Generally, at the hydrophilic surface, an increased density profile of water has been observed in previous studies41,55,56, which, in the present study is seen only with MGFF. The thickness of water layer is spreads for only about 5Å over the glass surface with MGFF while it spreads to about 10Å in both the Reax parameter sets. On the other hand, the penetration of whole water molecules as well as silanols is much deeper with MGFF as compared to either of the Reax potentials. Thus, the higher concentration of water over the glass and the deeper penetration of water into glass results in a higher silanol concentration with MGFF. The oxygens involved in OH- and H3O+ were too small to have an impact in the plot and have not been shown. 3.4 Reactivity of curved surfaces Amorphous silica surfaces have atomic and nanometer scale roughness that affect its reactivity37,41. The roughness at the atomic scale is captured in the “flat” surface simulations presented in the previous section. In order to understand the effect of curvature on the reactivities further, fully curved surfaces were created by melt-quenching a spherical block of h-cristobalite of radius 1.5nm with a reflecting spherical region (feature available in LAMMPS region command57) surrounding the crystal that prevents the diffusion of atoms from the molten block. The initial spherical block was obtained by separating out only the atoms within 1.5nm away from the center of a 3.5nm side cubic block. Similarly the reminder of the cubic block from the above was melt quenched with a reflecting sphere inside the block to create a spherical region of negative 9 ACS Paragon Plus Environment


Page 10 of 29

curvature inside the silica block. Drifting of the atom clusters was prevented by “freezing” of some of the atoms far away from the surface (the center of the sphere and the outer region in the block). Thus, spherical silica and silica with a spherical void of 3nm radii was obtained which could be used as surfaces of negative (concave) curvature or positive (convex) curvature. Water at the right density was filled in the region of the cubic block not occupied by silica and simulations were carried out to 1ns at 298K with the MG potential. Snapshots of the system with indicating the schematic arrangement of curved glass and water is shown in the inset of figure 7. The evolution of silanol concentration for these systems are compared with those in the flat surface system of section 3.3 in figure 7. Since the curved surfaces are in the immediate vicinity of water, there is an almost immediate onset of silanol formation in both the curved surfaces which is in contrast to the situation with the delayed onset of silanols in both the flat surfaces. However, the curved surfaces consistently show a lower reactivity to water with values of silanol concentration closer to the usually cited value1,53 of around 4.2-5.7/nm2. The reactivities of surfaces are a function of several factors including the surface energy and the defect concentration. As demonstrated by Rimsza et.al9, the calculated surface energies may vary depending on the potential being used and also connectivity of the glass. Nevertheless, the connectivity of the glass, which has a direct correlation to the quantity of defect species are a possible cause for differences in the reactivities and were examined further. Vitreous silica surfaces contain several types of defects and these defect sites are susceptible to reaction with water and formation of silanols8,21. One of the mechanisms of silanol formation in our simulations has been through the valence satisfaction of undercoordinated silicon by attachment of an OH- species or a water molecule which then donates a proton to an oxygen in proximity. Such valence satisfaction of undercoordinated silicons has been observed with MGFF41 and ReaxFF27 and also through ab-initio studies4. Figure 8 shows possible pathways for these reactions to occur that was observed with MGFF on the flat silica surface. Figure 9 shows a plot of the fraction of 3 coordinated silicons and their evolution with time for the surfaces studied. The flat silica surface shows an increased number of 3 coordinated silicons at the start, compared to either of the curved surfaces in our simulations. Since the flat surface system has 2 surfaces with only the top exposed to water, only the 3 coordinated Si found in the top half of the system are shown. Valence satisfaction of 3 coordinated Si occurs almost immediately on contact with water and this can be seen in the steep fall in their fraction. As evidenced from figure 5a water does not react immediately with the silica surface in the flat surface system because of the gap between the water cloud and the silica surface. Nevertheless, water begins to accumulate at the surface starting at around 20 ps and the fraction of 3 coordinated Si drops down to the level seen in the other curved surface systems. The distribution of undercoordinated Si (Si3) is almost exclusively close to the surface of the glass and evidence for this can be obtained by evaluating the density of various species along the 10 ACS Paragon Plus Environment

Page 11 of 29 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


axis normal to the glass surface. Figure 9 shows the distribution of Si, O and Si3 along in the top ~15 Å of the glass. Based on the way the flat surface was formed (Cleaving from the bulk), a much higher concentration of Si3 is formed at the surface but annealing resulted in the repair of some of them, the remainder of which show up at the surface. Vitreous silica surfaces are always oxygen terminated and this explains the tail of the oxygen distribution curve extending beyond that of Si. This concentration distribution curve was averaged over the first 20 ps before reaction with water after which the concentration of Si3 at this surface reduces to almost zero. While water does penetrate to around 11Å under the surface, it still does not reach the location of the deep seated Si3 (seen as the blip in the curve around 26Å) which persists all through the simulation. The calculation of concentration profile of species in the concave and convex systems is complicated because of the curvature and movement of the atoms at the surface which change the center of mass. However, these were not necessary since the fraction of Si3 does not change with time on the curved surfaces as seen in figure 10. Thus, these calculations establish that valence satisfaction of 3 coordinated Si is one of the mechanisms of silanol formation.

4. Discussions We evaluated the suitability of two sets of ReaxFF parameters and the MGFF reactive potentials for studying water reactivity at silica surfaces. These potentials were designed to reproduce the properties of bulk water and water-silica reactivities. Although all the potentials are reasonably accurate in reproducing the structure of bulk water, there were significant differences in the densities, diffusivities, cluster formation and cluster reactivities among these. Specifically, at low temperature simulations, the arrangement of the first shell of neighbors was a skewed tetrahedron with the ReaxFF parameters while it is closer to the ice tetrahedron with MGFF. This allows for a higher number of water molecules to gather around the central water molecule which is evidenced by the higher first peak for gO-O(r) with both the ReaxFF parameter sets. As a consequence, ReaxFF shows lower density and higher diffusion coefficient in bulk water. A newer version of the ReaxFF parameter sets is available at this time which has density and diffusivities closer to that of bulk water – however, this is yet to be parameterized for interactions with SiO2. Further evaluation of all the potentials was also done by studying the interaction of water with silica surfaces. The difference in density and reactivity of water molecule with Si(OH)4 manifests as a lower concentration of silanols, shallower penetration of water into the silica surface and lower density of water on top of the silica surface with both the Reax parameter sets. Several earlier studies have indicated the formation of layers of varying density55,58 at the surface which we observed as the densification of water at the surface. Silanol densities per unit surface area were observed to be slightly higher under certain conditions with the MGFF potentials. The commonly published value of 4.2-5.7/nm2 59, is based 11 ACS Paragon Plus Environment


Page 12 of 29

on the saturation of silanols on a cleaved cristobalite surface1 and does not explicitly include the penetration of water and sub-surface reactions. Including these could result in silanol concentrations upto 8/nm2 – which makes the values observed in this work well within range. While silanol formation is widely reported as being due to water attack on a strained siloxane bond8,21, valence satisfaction of undercoordinated species also contributes to silanol formation. Undercoordinated silicon being a more reactive site for water attack has been observed in earlier simulations4 as well as experimental studies17 but needs to be at the right configurations for the reactions to occur21. While these reactions were observed with all the potentials we tested, it was more prevalent with MGFF because of the higher concentration of water at the surface allowing greater access to defect sites. The absence of such point defects resulting in water molecules not dissociating has also been studied through spectroscopic measurements7,60,61 and shows the role played by surface structure in reactivity of water. In this work, we also found a correlation between the initial concentration of defective Si3s and the eventual number of silanols that are formed. Most of the defective Si3s were present at the exposed surface and are reacted upon by water very early in the simulation. Surface generation in our simulations was done by either the widely used technique of cleaving and annealing8,9,41 or by melt quenching with a reflective boundary that constrains the atoms to the required shape. In general, the defect species were higher in number with the cleave-anneal technique, although longer annealing at higher temperature relative to the system size, can remove some of these defects. Both these methods of surface preparation in simulation correspond to ways by which silica surfaces are made experimentally. Figure 11 shows the concentration of Si3s for silica surfaces prepared by both these methods in simulations. Clearly, the defect concentration at the surface is higher with the standard cleave-anneal surface. Consequently, these regions have a greater propensity to react with water and result in the higher number of silanols forming earlier as shown in the figure 11(b). This sets the path for an eventual higher silanol concentration early on in the simulation. Subsequent formation of silanols are at similar rates for both the surfaces. Our results corroborate experimental results that the method of preparation of surfaces bears a correlation to the reactivity of silica. We plan to use reactive potentials like MGFF for future studies of water under various levels of confinement in silicates.

5. Conclusions Based on our analysis of the dissociative/reactive water potentials, we observe that while ReaxFF potential describes the interactions more explicitly based on QM calculations, our results show that a simpler potential function like MGFF provides accurate description of water structure, diffusion, vapor equation of state, as well as water/silica reactivity. It is also computationally less resource consuming, hence access to longer time and larger length scale simulations, while providing similar levels of analysis capabilities. Analysis of silica surfaces reacting with water reveals further differences between the potentials and the reactivity was in general lower with the 12 ACS Paragon Plus Environment

Page 13 of 29 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


ReaxFF potentials. While the surface concentrations of silanols in glass vary in a narrow range, higher defect concentration at the surface in the form of undercoordinated silicons (Si3s) results in an initial higher concentration of silanols which these higher numbers being maintained for longer periods. It was also observed that the method of preparation of surfaces plays a significant role in the concentration of these defect species at the surface and hence has an influence on the reactivity of the surface. Comparison of the flat, concave and convex surfaces show similar reactivity and final silanol concentration if the surfaces were created in similar method. This suggests that the introduction of surface defects play an important role on surface reactivity and curvature at 1-2 nm scale does not significantly affect silica surface reactivity. The MGFF potential parameters can be used for larger scale and longer time MD simulations of water/silica systems such as diffusion in porous gel structures and in developing multicomponent glasses with development of compatible potential parameters.

Acknowledgements This research has been made possible by funding from the Center for Performance and Design of Nuclear Waste Forms and Containers, an Energy Frontier Research Center funded by the U.S. DOE, Office of Science, Basic Energy Sciences under Award # DE-SC0016584. The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin and the HPC Talon3 at University of North Texas for providing computational resources that have contributed to the research results reported within this paper.

References: 1.

Iler, R. K. The Chemistry of Silica (Iler, Ralph K.). ISBN: 978-0-471-02404-0 (John Wiley and Sons, 1979).

2.

O’Connor, T. L. & Greenberg, S. A. The Kinetics for the Solution of Silica in Aqueous Solutions. J. Phys. Chem. 62, 1195–1198 (1958).

3.

Du, J. & Rimsza, J. M. Atomistic computer simulations of water interactions and dissolution of inorganic glasses. npj Mater. Degrad. 1, 16 (2017).

4.

Walsh, T. R., Wilson, M. & Sutton, A. P. Hydrolysis of the amorphous silica surface. II. Calculation of activation barriers and mechanisms. J. Chem. Phys. 113, 9191–9201 (2000).

5.

Wilson, M. & Walsh, T. R. Hydrolysis of the amorphous silica surface. I. Structure and dynamics of the dry surface. J. Chem. Phys. 113, 9180–9190 (2000).

6.

Mischler, C., Horbach, J., Kob, W. & Binder, K. Water adsorption on amorphous silica surfaces: a Car–Parrinello simulation study. J. Phys. Condens. Matter 17, 4005–4013 (2005). 13 ACS Paragon Plus Environment


Page 14 of 29

7.

Chen, Y.-W. & Cheng, H.-P. Interaction between water and defective silica surfaces. J. Chem. Phys. 134, 114703 (2011).

8.

Garofalini, S. H. Molecular dynamics computer simulations of silica surface structure and adsorption of water molecules. J. Non. Cryst. Solids 120, 1–12 (1990).

9.

Rimsza, J. M., Jones, R. E. & Criscenti, L. J. Surface Structure and Stability of Partially Hydroxylated Silica Surfaces. Langmuir 33, 3882–3891 (2017).

10.

Rimsza, J. M., Yeon, J., van Duin, A. C. T. & Du, J. Water Interactions with Nanoporous Silica: Comparison of ReaxFF and ab Initio based Molecular Dynamics Simulations. J. Phys. Chem. C 120, 24803–24816 (2016).

11.

Trickey, S. B., Yip, S., Cheng, H.-P., Runge, K. & Deymier, P. A. A Perspective on Multi-scale Simulation: Toward Understanding Water-silica. J. Comput. Mater. Des. 13, 1–12 (2006).

12.

Cheng, H.-P. et al. Quantum, classical, and multi-scale simulation of silica–water interaction: molecules, clusters, and extended systems. J. Comput. Mater. Des. 13, 161– 183 (2006).

13.

Bromley, S. T. & Slater, B. Modelling the structures and reactivity of silica and water: from molecule to macroscale. in CECAM conference

14.

Rahman, A., Mandell, M. J. & McTague, J. P. Molecular dynamics study of an amorphous Lennard‐Jones system at low temperature. J. Chem. Phys. 64, 1564–1568 (1976).

15.

Garofalini, S. H., Halicioglu, T. & Pound, G. M. Computer simulation of amorphous surfaces: energy profile and binding energy. J. Non. Cryst. Solids 37, 411–415 (1980).

16.

Guillot, B. & Guissani, Y. How to build a better pair potential for water. J. Chem. Phys. 114, 6720–6733 (2001).

17.

Lee, J.-W., Tomozawa, M. & MacCrone, R. K. Point defect formation and annihilation in silica glass by heat-treatment: Role of water and stress. J. Non. Cryst. Solids 354, 1509– 1515 (2008).

18.

Rimola, A. & Ugliengo, P. A quantum mechanical study of the reactivity of (SiO)2defective silica surfaces. J. Chem. Phys. 128, 204702 (2008).

19.

Ma, Y., Foster, A. S. & Nieminen, R. M. Reactions and clustering of water with silica surface. J. Chem. Phys. 122, 144709 (2005).

20.

Du, M.-H., Kolchin, A. & Cheng, H.-P. Water–silica surface interactions: A combined quantum-classical molecular dynamic study of energetics and reaction pathways. J. Chem. Phys. 119, 6418–6422 (2003).

21.

Leed, E. A., Sofo, J. O. & Pantano, C. G. Electronic structure calculations of physisorption and chemisorption on oxide glass surfaces. Phys. Rev. B 72, 155427 (2005).

22.

Du, J. & Cormack, A. N. Molecular Dynamics Simulation of the Structure and Hydroxylation of Silica Glass Surfaces. J. Am. Ceram. Soc. 88, 2532–2539 (2005). 14 ACS Paragon Plus Environment

Page 15 of 29 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


23.

Chenoweth, K., van Duin, A. C. T. & Goddard, W. A. ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J. Phys. Chem. A 112, 1040–1053 (2008).

24.

Mahadevan, T. S. & Garofalini, S. H. Dissociative Water Potential for Molecular Dynamics Simulations. J. Phys. Chem. B 111, 8919–8927 (2007).

25.

Pedone, A. et al. FFSiOH: a New Force Field for Silica Polymorphs and Their Hydroxylated Surfaces Based on Periodic B3LYP Calculations. Chem. Mater. 20, 2522– 2531 (2008).

26.

SCM TEAM. ReaxFF: reactive MD with graphical interface & analysis tools. (1995). Available at: https://www.scm.com/doc/ReaxFF/General.html. (Accessed: 15th November 2017)

27.

Fogarty, J. C., Aktulga, H. M., Grama, A. Y., van Duin, A. C. T. & Pandit, S. A. A reactive molecular dynamics simulation of the silica-water interface. J. Chem. Phys. 132, 174704 (2010).

28.

Larsson, H. R., van Duin, A. C. T. & Hartke, B. Global optimization of parameters in the reactive force field ReaxFF for SiOH. J. Comput. Chem. 34, 2178–2189 (2013).

29.

Yeon, J. & van Duin, A. C. T. ReaxFF Molecular Dynamics Simulations of Hydroxylation Kinetics for Amorphous and Nano-Silica Structure, and Its Relations with Atomic Strain Energy. J. Phys. Chem. C 120, 305–317 (2016).

30.

Côté, A. S., Cormack, A. N. & Tilocca, A. Reactive molecular dynamics: an effective tool for modelling the sol–gel synthesis of bioglasses. J. Mater. Sci. 52, 9006–9013 (2017).

31.

Rimsza, J. M. & Du, J. Interfacial Structure and Evolution of the Water–Silica Gel System by Reactive Force-Field-Based Molecular Dynamics Simulations. J. Phys. Chem. C 121, 11534–11543 (2017).

32.

Macià Escatllar, A., Ugliengo, P. & Bromley, S. T. Modeling hydroxylated nanosilica: Testing the performance of ReaxFF and FFSiOH force fields. J. Chem. Phys. 146, 224704 (2017).

33.

Lockwood, G. K. & Garofalini, S. H. Proton Dynamics at the Water–Silica Interface via Dissociative Molecular Dynamics. J. Phys. Chem. C 118, 29750–29759 (2014).

34.

Webb, M. B., Garofalini, S. H. & Scherer, G. W. Molecular Dynamics Investigation of Solution Structure between NaCl and Quartz Crystals. J. Phys. Chem. C 115, 19724– 19732 (2011).

35.

Zhang, W. & van Duin, A. C. T. Second-Generation ReaxFF Water Force Field: Improvements in the Description of Water Density and OH-Anion Diffusion. J. Phys. Chem. B 121, 6021–6032 (2017).

36.

Garofalini, S. H., Mahadevan, T. S., Xu, S. & Scherer, G. W. Molecular mechanisms causing anomalously high thermal expansion of nanoconfined water. ChemPhysChem 9, 1997–2001 (2008). 15 ACS Paragon Plus Environment


Page 16 of 29

37.

Rimola, A., Costa, D., Sodupe, M., Lambert, J.-F. & Ugliengo, P. Silica Surface Features and Their Role in the Adsorption of Biomolecules: Computational Modeling and Experiments. Chem. Rev. 113, 4216–4313 (2013).

38.

Peri, J. B. & Hensley, A. L. The surface structure of silica gel. J. Phys. Chem. 72, 2926– 2933 (1968).

39.

Battisha, I. K., El Beyally, A., El Mongy, S. A. & Nahrawi, A. M. Development of the FTIR properties of nano-structure silica gel doped with different rare earth elements, prepared by sol-gel route. J. Sol-Gel Sci. Technol. 41, 129–137 (2007).

40.

van Duin, A. C. T., Dasgupta, S., Lorant, F. & Goddard, W. A. ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 105, 9396–9409 (2001).

41.

Mahadevan, T. S. & Garofalini, S. H. Dissociative Chemisorption of Water onto Silica Surfaces and Formation of Hydronium Ions. J. Phys. Chem. C 112, 1507–1515 (2008).

42.

Soper, A. K. The Radial Distribution Functions of Water as Derived from Radiation Total Scattering Experiments: Is There Anything We Can Say for Sure? ISRN Phys. Chem. 2013, 1–67 (2013).

43.

Soper, A. K. The radial distribution functions of water and ice from 220 to 673 K and at pressures up to 400 MPa. Chem. Phys. 258, 121–137 (2000).

44.

Wagner, W. & Pruß, A. The IAPWS Formulation 1995 for the Thermodynamic Properties of Ordinary Water Substance for General and Scientific Use. J. Phys. Chem. Ref. Data 31, 387–535 (2002).

45.

Panagiotopoulos, A. Z. in Observation, Prediction and Simulation of Phase Transitions in Complex Fluids 463–501 (Springer Netherlands, 1995). doi:10.1007/978-94-011-00656_11

46.

National Institute of Standards and Technology. NIST Chemistry WebBook, NIST Standard Reference Database Number 69. doi:10.18434/T4D303

47.

Yeh, I.-C. & Hummer, G. System-Size Dependence of Diffusion Coefficients and Viscosities from Molecular Dynamics Simulations with Periodic Boundary Conditions. J. Phys. Chem. B 108, 15873–15879 (2004).

48.

Holz, M., Heil, S. R. & Sacco, A. Temperature-dependent self-diffusion coefficients of water and six selected molecular liquids for calibration in accurate 1H NMR PFG measurements. Phys. Chem. Chem. Phys. 2, 4740–4742 (2000).

49.

Simpson, J. H. & Carr, H. Y. Diffusion and Nuclear Spin Relaxation in Water. Phys. Rev. 111, 1201–1202 (1958).

50.

Wang, J. H. Self-Diffusion Coefficients of Water. J. Phys. Chem. 69, 4412–4412 (1965).

51.

Wu, Y., Tepper, H. L. & Voth, G. A. Flexible simple point-charge water model with improved liquid-state properties. J. Chem. Phys. 124, 24503 (2006).

52.

Shields, R. M., Temelso, B., Archer, K. A., Morrell, T. E. & Shields, G. C. Accurate Predictions of Water Cluster Formation, (H 2 O) n =2−10. J. Phys. Chem. A 114, 11725– 16 ACS Paragon Plus Environment

Page 17 of 29 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


11737 (2010). 53.

Zhuravlev, L. T. & Potapov, V. V. Density of silanol groups on the surface of silica precipitated from a hydrothermal solution. Russ. J. Phys. Chem. 80, 1119–1128 (2006).

54.

Mahadevan, T. S. et al. Diffusion transport of nanoparticles at nanochannel boundaries. J. Nanoparticle Res. 15, 1477 (2013).

55.

Gallo, P., Ricci, M. A. & Rovere, M. Layer analysis of the structure of water confined in vycor glass. J. Chem. Phys. 116, 342 (2002).

56.

Spohr, E., Hartnig, C., Gallo, P. & Rovere, M. Water in porous glasses. A computer simulation study. J. Mol. Liq. 80, 165–178 (1999).

57.

Plimpton, S. J. LAMMPS Molecular Dynamics Simulator. J. Comp. Phys 1–19 (1995). doi:https://doi.org/10.1006/jcph.1995.1039

58.

Bourg, I. C. & Steefel, C. I. Molecular Dynamics Simulations of Water Structure and Diffusion in Silica Nanopores. J. Phys. Chem. C 116, 11556–11564 (2012).

59.

Zhuravlev, L. T. The surface chemistry of amorphous silica. Zhuravlev model. Colloids Surfaces A Physicochem. Eng. Asp. 173, 1–38 (2000).

60.

Wendt, S. et al. The interaction of water with silica thin films grown on Mo(112). Surf. Sci. 565, 107–120 (2004).

61.

Yates, J. T. Water interactions with silica surfaces: A big role for surface structure. Surf. Sci. 565, 103–106 (2004).



Page 18 of 29

TOC Graphic


Page 19 of 29 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 1: O-O pair distribution function in bulk water showing matching location of the first and second peaks. Experimental values were taken from Soper42



Page 20 of 29

Figure 2: The computed equation of state measured using Gibbs Ensemble Molecular Dynamics for MGFF, Reax-Larsson and Reax-Wei in comparison with experimental data. Results from Yeon potential had very similar values to the Larsson potential. The MGFF and recently optimized ReaxFF parameter set from Wei et.al35 show excellent agreemeng with experiments.


Page 21 of 29 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


(a)

(b)

Figure 3: (a) Cluster energies of (H2O)n , n=1-9. While all the models do predict lower energies (and hence higher stabilities) of larger clusters in general, the energies of MGFF are closer in line with previous results. (b)Stable configurations of the 5 water cluster obtained with MGFF (top) and Reax-Yeon(bot). Values for ΔH298 are from Shields52



Page 22 of 29

Figure 4: Energy of the Si(OH)4 – H2O system plotted as a function of the distance of separation between the molecules. The inset figure shows the stable, multiple hydrogen bonded configuration formed with MGFF at around 2.6A separation


Page 23 of 29 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


(a)

(b)

Figure 5: (a) Comparison of silanol concentration as a function of time on flat silica surface from simulations with the three potentials – the dotted line at the bottom is the concentration of [SiOH2] which was formed only with the MGFF potential. The inset shows the arrangement of the system with atomistic representations depicted inside the transparent blocks (b) Plot of the ratio of (OH)- and H2O over time. A small amount of H3O+ was also observed, but only in the initial 0.2ns of simulation with short lifetimes and are not shown. The high percent of OH- with Larsson parameter set correlates to the larger percent of glass oxygen involved in silanol formation



Page 24 of 29

MGFF

Reax-Yeon

Reax-Larsson

Figure 6: Plot of the density profiles of different types of oxygens along the normal to the surface (Z axis) at the top of the glass. Glass oxygens are attached to at least one Si, Water oxygens are attached to exactly 2 H and no Si and SiOH oxygens are attached to at least 1 Si and 1 H. The nominal number density of water oxygen in a glass is expected around 44/nm3 which is seen in the glass oxygens with all the three potentials. Bottom panel is for Reax-Larsson, middle is ReaxYeon and top is MGFF. All the curves have been smoothened with csplines to improve visual clarity.


Page 25 of 29 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 7: Silanol formation in “true” curved surfaces without any flat surface regions. The pictures on the side show the snapshot of the system. The silica block has been rendered to show some of the surface features



Page 26 of 29

Figure 8: Snapshot of possible reaction pathways of silanol formation from 3 coordinated Si. (a)3 (c) shows a water molecule attaching to a single Si and releasing a proton to another water molecule forming a hydronium. In (d)-(f) water molecule releases a proton to a neighboring NBO 3 and forming an OH- which then attaches to the Si . This also shows one of the possible means by which a free OH- can be generated. The grey bond-network in the background are to show the location of the atoms near the surface.


Page 27 of 29 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 9: Concentration of 3 coordinated Si (Si3) as a percent fraction of total Si in the initial 100ps of simulation for curved and flat surface systems. While these defective Si’s are usually very low in bulk silica – and are also repaired during annealing, some of them persist at the surface which are repaired almost immediately on reaction with water.



Page 28 of 29

Figure 10: Density profiles of species along the surface normal. The Si3 concentration reduces to almost zero in the top surface after the initial 200 ps. The very small number of (seen as the small bump around 26Å deep surface Si3 persists through the simulation even though water has penetrated up to 25Å.


Page 29 of 29 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 11: (a) Density of defective silicons in Melt-Quenched and Cleave-Annealed samples compared to average silicon density in the samples (b) Silanol formation kinetics for the same samples showing lower concentration for the Melt-Quenched sample


Evaluating Water Reactivity at Silica Surfaces Using Reactive Potentials

Recommend Documents