Surface Reactivity and Leaching of a Sodium Silicate Glass under an

May 31, 2019 - Sodium ions are bound close to the surface as Na–OH ion pairs during the ..... Q3, Q4) and surface ionization corresponding to the pH...
0 downloads 0 Views 4MB Size
Download PDF
Article Cite This: J. Phys. Chem. C 2019, 123, 15606−15617

pubs.acs.org/JPCC

Surface Reactivity and Leaching of a Sodium Silicate Glass under an Aqueous Environment: A ReaxFF Molecular Dynamics Study Seung Ho Hahn and Adri C. T. van Duin* Department of Mechanical Engineering, The Pennsylvania State University, University Park, Pennsylvania 16802, United States

Downloaded via UNIV OF SOUTHERN INDIANA on July 28, 2019 at 15:47:06 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: The presence of leachable ions in multicomponent silicate glasses complicates the understanding of their surface structure and chemistry under a humid or aqueous environment, in particular when compared to pure silica−water interactions. Here, we performed large-scale ReaxFF reactive molecular dynamics (MD) simulations and investigated the structural and chemical modifications that take place at the sodium silicate glass−water interface. From our simulation results, we observed that the water interaction gradually replaces sodium ions associated with the surface nonbridging oxygen (NBO) and found that molecular water dissociates into H+ and OH− in the presence of these modifier cations. The protonation of the NBO site produces silanol (SiOH), and our results denote that protons can diffuse further into the bulk region by the discrete proton hopping between two adjacent NBO sites. Consequently, the number of silanol formation increases along the glass−water reaction, rather than reaching an equilibrium point within the nanosecond time scale of our MD simulations. In addition, sodium leaching and its mobility at the interface were studied to provide adequate information on how the ion-exchange process dictates their transport within the interface. Sodium ions are bound close to the surface as Na−OH ion pairs during the initial stages of the water interaction and exhibits a certain residence time at the surface before they are released to the aqueous media. The reactive MD simulation presented here can successfully provide physical insights needed to understand the surface chemistry with the time and length scales which would otherwise be difficult for first-principles or experimental studies to evaluate.

1. INTRODUCTION Side-by-side with a growing number of glass applications in different fields, scientific interest in understanding the water interaction with the glass surface has also increased. One of the examples would be the chemical corrosion behavior of glass by water which is of great importance for the safe usage of glass materials for the long-term disposal of radioactive waste.1−4 Another application of glass as a device for promoting the growth and repair of damaged tissue5,6 requires site-by-site assessment of the hydrophilic/hydrophobic nature of the surface because the surface modification takes place with the immersion of body fluids. The increase in wear resistance of the glass surface at high relative humidity also demands attention, as the adsorbed water layers to the sodium-leached sites are believed to be playing a crucial role for the enhancement of this mechanical property.7,8 As such, the presence of water is closely related to the changes in surface properties, and therefore, a thorough understanding of water interactions with glass surfaces is highly necessary. In experiments, the glass−water interface can be probed by different methods, and many endeavors have been made to this end. The chemical composition changes of the glass surface after water interactions can be analyzed using X-ray photoelectron spectroscopy (XPS) or nuclear magnetic resonance.9−11 The vibration spectroscopy methods that use nonlinear optical processes such as sum-frequency generation © 2019 American Chemical Society

(SFG) and second-harmonic generation allows the surface characterization of adsorbed of water layer structures and polarization of water molecules.12−14 While all of these surface probing techniques are essential to highlight the consequent changes upon water contact, critical information regarding the structural characteristics and chemical reaction mechanisms is rather difficult to obtain. In recent years, computational modeling of glass with molecular dynamics (MD) simulations has been providing fruitful information to support the findings of experiments and details to the unanswered questions, expanding its versatility from bulk property investigations to surface characterizations at an atomistic resolution. Within a MD simulation, interatomic force calculations are operated by either quantum mechanical (QM) methods or empirical force fields. QM methods are very applicable to simulate reactive processes of the glass−water interface because they can implicitly describe changes in electron behaviors. This level of theory provides first-principles level accuracy on energetics of the chemical and physical modifications.15−17 However, the high computational expense restricts ab initio MD simulation’s scalability to picoseconds and typical maximum system sizes of few hundred atoms. Received: March 28, 2019 Revised: May 28, 2019 Published: May 31, 2019 15606

DOI: 10.1021/acs.jpcc.9b02940 J. Phys. Chem. C 2019, 123, 15606−15617

Article

The Journal of Physical Chemistry C

(NaSiOx)−water interactions42 to expand the modeling capability from ab initio approaches to interatomic potential domain.43,44 Key descriptions such as sodium vacancy migration in sodium silicate crystal, Na+ interaction with hydroxylated silica surface and with water, and NaOH− n[H2O] (n = 1−6) dissociation energies were calculated from density functional theory (DFT) and introduced for parameterizing the potential to properly reproduce essential dynamical processes at the sodium silicate glass−water interface.42 In this work, we present large-scale reactive MD simulations of the sodium silicate glass−water interface. The surface reactivity was investigated by looking into the details of the hydroxylation mechanism which accompanies adsorption and dissociation of water molecules. In addition, sodium leaching and its mobility at the interface were studied from the MD trajectories to provide adequate information on how the structural modifications affect their behaviors. Although a detailed assessment of surface reactive sites and critical roles of modifier ions can be provided with ab initio models, the system size is often not large enough for statistical analysis of structural features. Furthermore, important dynamical aspects such as pH, ionic strength of the solution, and the types of cations and anions present in solution which influence the protonation−deprotonation equilibria at the phase boundary are difficult to be meaningfully evaluated within the time and length scale of the DFT calculations. Our work alleviates these limitations, providing new physical insights into atomistic resolution on the chemical dynamics associated with the glass− water interface.

Meanwhile, these temporal and spatial constraints are considerably eased for the conventional MD simulations with interatomic potentials. Over the last decades, a number of efficient interatomic potentials for silicate crystals and silicabased glasses have been developed. Such potentials include BKS potential18 and extensions of this rigid ionic model, Teter19 and Pedone20 potentials as well as shell-ion models which take account of polarization effects.21−23 A reliable description of the short- and medium-range order of bulk silicate and glasses can be obtained with these two-body interatomic potentials, and structural, mechanical, and vibrational properties are relatively well understood from MD simulations.24 Contrary to extensive modeling and simulation work devoted to the bulk, investigation on silica surfaces and aqueous interfaces had been limited primarily because of the lack of potential’s capability to reproduce the surface chemistry associated with hydroxylated groups such as silanol (Si−OH). Despite the importance of the area density and distribution of silanol groups on the surface properties, prior models disregarded the alteration of surface chemistry by assuming neutral or nonstoichiometric surfaces with silanol groups.25−28 This limited applicability was resolved with the inclusion of realistic surface models and quantitative predictions of surface properties and selective adsorption of peptides, drugs, surfactants, and synthetic polymers, as a function of the surface environment and pH are now directly comparable with experimental measurements.29−32 While full-range modeling of variable surface chemistry for all types of silica, temperature treatment, pH values, and synthesis methods have become feasible, a limitation still persists.31,33,34 Because of the inherent nature of the analytical form of the bonded harmonic force field, these force fields cannot describe any surface reactivity and chemical reaction mechanisms at the interfaces of silica-based glass. Various reactive potentials that were developed and parameterized to describe water dissociation behavior and subsequent chemisorption of dissociated species on the glass surface have been important in this extent. Consequently, a solid understanding of interfacial properties is established for the silica−water interactions with reactive potentials such as dissociative water potentials developed by Mahadevan and Garofalini35,36 and the Si/O/H parameterization in ReaxFF.37,38 However, silica−water interactions are only part of the full glass−water interface challenge because many of functional or engineering glasses are complex multicomponent glasses. In case of silicate glasses with network modifier alkali/alkaline cations, the presence of these leachable ions further complicates the glass−water interaction. For example, surface nonbridging oxygen (NBO) sites are prone to water dissociation and protonated when modifier alkali cations (M+) neutralize the remaining OH− groups. Subsequent reaction dissociates M−OH into M+ and OH−, and these ions are dissolved to the liquid media. The residual glass surface is subject to the protonation−dissociation equilibria according to the pH value and ionic strength of the aqueous solution. Until recently, no reactive potential parameterizations were available for such applicability to multicomponent glass− water systems. For this reason, the specific effects of the network modifier cations and water dissociation mechanisms at the surface of multicomponent silicate glasses had been studied mostly with first-principles methods.17,39−41 Recently, we developed a ReaxFF reactive force field for sodium silicates

2. COMPUTATIONAL DETAILS 2.1. Simulation MethodsReaxFF Reactive Force Field. ReaxFF45,46 is one of the most versatile reactive force field which has proved its applicability in various research areas in chemistry and physics of materials science domain.47−50 It is a bond-order-dependent potential, meaning that the connectivity between all atom pairs is updated depending on the atomic local environment during the MD simulation. This feature allows ReaxFF to reliably describe the chemical reaction process as well as the energetics involved. The total energy of the system in ReaxFF generally consists of contributions from bond-order-dependent terms and nonbonded interaction terms. The bond order is calculated based on the interatomic distances of all atom pairs at every iteration. Hence, the energy contributions from bond-order-dependent terms such as bond, valence angle, and torsion angle vanish upon bond dissociation, while leaving only the nonbonded interactions to be considered afterward. The van der Waals and Coulomb energies are the contributions of the nonbonded interactions, and ReaxFF employs the electronegativity equalization method51 to derive the atomic charges. Further details of the ReaxFF potential forms can be found in publications by van Duin and Chenoweth.45,52 In this work, sodium silicate glass−water reaction simulations were performed with the ReaxFF reactive force field to capture all of the dynamic and reactive situations at the corresponding solid−liquid interface. 2.2. Sodium Silicate Glass−Water Interface; System Preparation and Simulation Setup. For the system preparation, we targeted to build the sodium silicate glass− water interface by a three-step process in this work. First, an amorphous sodium silicate glass structure was generated, followed by the preparation of water system, and the last step 15607

DOI: 10.1021/acs.jpcc.9b02940 J. Phys. Chem. C 2019, 123, 15606−15617

Article

The Journal of Physical Chemistry C

constants of 0.1 and 5 ps were used for the Berendsen thermostat and barostat, respectively. In all of the reactive MD simulations presented here, we employed our recently developed Na/Si/O/H ReaxFF parameter set,42 and simulations were performed with the ReaxFF implementation within the Amsterdam Density Functional (ADF) package.56

was to model the glass surface to be in contact with an aqueous medium. To maximize the surface area, we chose to build a relatively large model of bulk glass with 6000 atoms. The amorphous sodium silicate system of composition SiO2/Na2O = 70:30 was generated through the melt-quench process. This methodology is adopted in most cases of bulk glass construction in computational chemistry and was found out to produce reliable bulk glass property, despite the relatively large quenching rate compared to that of the experiment.53,54 An initial mixture of 1400 SiO2 and 600 Na2O units were inserted in a periodic cubic box of a = 45.434 Å, having the initial density of 2.47 g/cm3. Using canonical (NVTconstant number of atoms, volume, and temperature) ensemble, the system was annealed from 4000 K. After the system temperature had reached 300 K from gradual cooling, we further equilibrated the system using NPT (constant number of atoms, pressure, and temperature) ensemble at 1 atm, 300 K condition to let the system’s density converge. The last trajectory from this step was used to perform the final equilibration with NVT ensemble for another 500 ps. Following the bulk glass construction, we modeled the water system to be combined for generating the glass−water interface. From previous studies, ReaxFF parameterization of O/H was reported to produce slight underestimation of water density.36,55 Because the O/H parameters of the ReaxFF potential used in this study also have foundation on that, we performed supplementary simulations to correctly select the dimension of the water system to be incorporated to the glass surface and to reduce any unphysical contact due to different compressibility of the systems. Therefore, the density of bulk water was investigated, and an NPT ensemble was applied to the bulk water system with 1500 water molecules for 250 ps. After we obtained the steady density for the system of this particular size, the water box was adjusted to fit the surface area of the sodium silicate glass. With the bulk glass structure obtained from the procedures above, we used this configuration to generate a free surface prior to modeling the glass−water interface. The system dimension in z-direction was expanded by 25 Å from both top and bottom sides, resulting in the vacuum region that separates the top surface from the bottom surface. While threedimensional periodicity is maintained, we can benefit from this separation by generating two free surfaces which would increase the statistics for evaluating the chemical modifications of the glass surface. Afterward, both the top and bottom surfaces were relaxed at 300 K, followed by annealing at 800 K and then cooled back to 300 K. To each of contacting surfaces, 1500 water molecules were inserted (3000 H2O in total) to fill the vacuum region. During this stage, the cell parameter c was adjusted to introduce the correct density of bulk water. Care was taken to avoid any initial contact between the atoms of the surface and those of water molecules, by allowing at least 2 Å separation. Further details of the bulk system preparation, surface construction, and relaxation process are elaborated in the Supporting Information (Figures S1−S5). Altogether, the glass−water interface system contained 15 000 atoms in total, and the glass−water reaction was carried out through MD runs with NVT equilibration at 300 K for 1 ns. To produce results with reasonable statistics, we carried out the MD run in triplicate with different random placements of water molecules in the starting configurations. Unless otherwise noted, the Newton’s equation of motion was solved with the Verlet algorithm at a time step of 0.25 fs. Damping

3. RESULTS AND DISCUSSION 3.1. Chemical Modification at the Glass−Water Interface. For water interactions with pure silica, we would observe silanol (SiOH) formation and hydrolysis of siloxane (Si−O−Si) bonds at the surface upon water contact.36,57,58 In case of sodium silicate−water interactions, the presence of leachable ions in the sodium silicates additionally requires following ion-exchange processes to be accounted for59,60 SiONa + H 2O → SiOH + NaOH

SiOH + NaOH → SiOH + Na + + OH−

SiO−···Na+ groups (surface NBO sites, sodium siloxide) that are directly exposed on the surface of sodium silicate glasses undergo these ion-exchange reactions and form NaOH. Subsequently, NaOH dissociates into Na+ and OH−. From the perspective of Na+ transport, these reaction sequences are referred as the “leaching” process that sodium ions in glass are extracted to the liquid media. Once the superficial cations are leached out, the silanol groups of the residual ion-exchanged silica surface are subject to the protonation and dissociation equilibria with respect to the pH value of the aqueous solution. It is therefore essential to evaluate the chemical modification of the surface (e.g., area density of silanol groups and of siloxide groups) and cation dissociation to the solution to understand the effects of surface chemistry and pH to the surface properties of silica and silica-based glasses. Well-established silica surface models31,32 that have fully considered different surface environments (Q2, Q3, Q4) and surface ionization corresponding to the pH can be a benchmark for the assessment of reactive MD simulation results. Our sodium silicate surface model is composed of mixed Q2/Q3 surface environment with a ratio of approximately 1:1 (see Supporting Information, Figure S6). This indicates that the surface is roughly estimated to have a silanol density of ∼7 groups per nm2 (upper range of 9.4−4.7 nm−2) at the point of zero charge. With this in mind, it is anticipated that the ion-exchange reaction would provoke chemical modification toward the estimated value. In the following, we first present a chemical modification at the sodium silicate glass−water interface during the initial stage of the glass−water reaction. The outcomes of sodium leaching and adjustments in the degree of ionization upon water contact will be discussed in Section 3.3 through qualitative comparison with the previous studies. The classification between the glass region and water region of the starting configuration (t = 0 ps, shown in Figure 1a) is clear because we deliberately placed the water molecules in the system to prevent them from making any initial close contact with the glass surface. Immediately after the start of the reactive MD simulation, the separation between two regions becomes less clear as the water molecules rapidly start to fill up the empty spaces and interact with the glass surface. Because of the ion-exchange process, sodium ions in the glass are substituted with protons in water, and Figure 1b (t = 1000 ps) clearly shows that these mobile ions in the glass have 15608

DOI: 10.1021/acs.jpcc.9b02940 J. Phys. Chem. C 2019, 123, 15606−15617

Article

The Journal of Physical Chemistry C

Figure 1. Snapshots of the glass−water interface showing (a) initial configuration of the system and (b) configuration after 1000 ps of reactive MD simulation. Color legend: Na (blue), Si (ivory), Oglass (gray), Owater (red), H (white). Note that only the top side of the contact surface is illustrated here. Figure 2. Mass density profiles of the silicon (colored in gray), sodium (blue), and hydrogen (red) atoms after 1 ns of the sodium silicate glass−water reaction. The dashed line indicates dividing boundaries between the glass region (section in the middle) and water region (left and right sections of the middle portion). The profile of sodium atoms explicitly shows that they had been leached out to the water region.

indeed leached out to the water region. However, to quantitively assess this kind of chemical modification during the earliest stages of the glass−water reaction, sharp boundaries should be defined to distinguish in which media the reaction product species are located at. In this respect, we employed a methodology to define the dividing boundary (Gibbs surface) on each side of the glass− water interface. The same concept was adopted for the silica− water interface38 as well as for the vapor−liquid water interface61 and has been reported to describe the boundary for assessing property modifications within the region. Because silicon atoms are considered to be essentially immobile, low diffusion, network formers within the sodium silicates62 and the Si-profiles from our simulations indeed suggest structural integrity of the glass surface throughout the length and time of our simulations, we determined the boundaries with the silicon atom’s number density profile such that following equation is satisfied.

Figure 3. Time evolution of different species’ number change during the glass−water reaction(a) ion-exchange between Na+ from bulk glass with H+ from water and (b) magnified portion at [0, 100 ps]. Within 1.5 ps upon water contact, water molecules rapidly start to penetrate the glass surface.

∫water region up to Gibbs surface (n(z) − n1)dz =−

∫Gibbs surface up to glass region (n(z) − n2)dz

surface area. The increase in both species is the result of the ion-exchange reaction, and unlike the silica−water interaction, there do not exist an equilibrium concentration of the main products of this reaction.38,58 While one may expect the 1:1 stoichiometric balance between the ion-exchange products (Na+/H+) at the glass−water interface, the difference between the total number of silanol formation and number of leached sodium ions during 1 ns of the glass−water reaction can be best explained by the presence of reactive sites at the surface. First of all, even though considerable simulation time has been awarded to relax the glass surface during sample preparation (noted in the Supporting Information), we find that there always are undercoordinated Si (3-coordinated Si, IIISi) atoms present at the surface. After exposure to the water media, these undercoordinated Si atoms are almost immediately consumed in the reaction where the hydroxyl from the dissociated water molecule binds to.35,36,63,64 This high affinity of water molecules to these defect sites results in higher number of silanol after 1 ns of reaction as clearly observed from Figure 3a. Another reason can be the contribution of the hydrolysis of siloxane bond at the surface. For pure silica, siloxane bonds are widely reported as one of the reactive sites, that is attacked by water molecules to produce silanol at the surface.65−68 In our simulation, the hydrolysis of the siloxane bonds in sodium silicates are not prevalent, and the contribution to overall SiOH number is negligible. However, it is worth mentioning

where n(z) is the number density of silicon atoms along the zdirection of the whole system, n1 and n2 are the randomly distributed silicon number density per slab volume in bulk water and bulk glass. With this dividing surface, any silicon atoms residing outside would be identified as being within the interface.38 For sodium ions, they would be considered as leached out from the surface when their relative position from the boundary is on the outside. Also, for protons, the boundary would be the reference point whether they are part of the glass network at the interface or they have diffused into the bulk region of the glass. Figure 2 shows the mass density profiles of the Si, Na, and H atoms in the glass−water interface model, after exposing the top and bottom glass surfaces to the aqueous environment for 1 ns. The mass density was calculated by cutting the entire system into 0.5 Å slabs along the z-direction and summing up the mass of atoms therein. The profile of the sodium ions explicitly shows that they could be leached out upon water contact and released to the solution. Likewise, the hydrogen mass density profile reflects that there had been ion-exchange process and atomic diffusion into the bulk sodium silicate. Using the dividing boundaries at both interfaces, we analyzed the changes in chemical species of the glass surface. Figure 3a presents the time evolution of averaged total number of leached sodium ions and that of silanol (SiOH) units per 15609

DOI: 10.1021/acs.jpcc.9b02940 J. Phys. Chem. C 2019, 123, 15606−15617

Article

The Journal of Physical Chemistry C

ment is established for the experimental observations69−71 and predictions from MD simulations for the silanol concentration in pure silica, a direct comparison of a silicate glass with leachable ions demands careful considerations. In the experiment, the oxygen speciation is derived from the atomic % information from XPS within the top 10 nm region of the surface, and hydrous species (mostly SiOH) concentration is determined based on the stoichiometry of the glass composition. This calculation from the compositional quantification neglects any spatial heterogeneity in composition because the dealkalization region of the glass surface is often probed up to 100 nm with XPS. On the contrary, ReaxFF MD results denote that only a handful amount of hydrous species was able to diffuse into 5−6 Å of the subsurface region. Considering the compositional gradient of the species populating at or within the subsurface region, the areal density we calculated earlier (5.79 ± 0.14 silanols per unit surface area) can be revisited with a modified region selection. From Figure 2, the mass density profile of sodium atoms denotes that there is a sodium-depleted zone located approximately 3 Å below the Gibbs surface boundary. We selected all atom components within this region and investigated the concentration of bridging oxygen (BO), NBO, SiOH, OH−, and IIISi. Figure 4b shows the plot of concentration change before water interaction and after 1 ns of MD simulation. On the basis of the updated selection, our results predict a silanol concentration of 4.8 nm−2. The areal density for the same species is measured to be 3.15−3.5 nm−2 for soda lime silica (SLS, sodium calcium silicate) glass,70,71 which is component-wise most alike to our simulated glass. The discrepancy of SiOH concentration can be explained with the difference in the surface alkali content, as our modeled glass had [Na] = 23.4 at. % at the surface. SLS glass typically possesses less amount of modifier cations,72,73 even with the sodium and calcium content combined. The high concentration of NBO (Figure 4b) before water interaction is the direct evidence of the high alkali content in our simulated glass, as the network connectivity is reduced by the large number of network modifiers. In the end, approximately 2/3 of these surface NBOs were protonated as denoted in Figure 4b. One of many advantages of the MD simulation technique is its capability to capture the atomistic origin of the exact structures that are responsible for distinct features of a material. Here, we observed that 0.73 OH− are present per unit surface area after exposure to the water environment and that they had made their way into the subsurface region. Experiment analysis with the glass composition can only probe OH species as OH groups, rather than distinguishing individual contributions of the oxygen speciation. The surface characterization by SFG and ATR−IR facilitates the understanding of the contribution to the hydrous species, but direct information, as can be obtained from MD simulations, can assist the experimental interpretation. 3.2. Reactivity of the Surface NBO Sites. In the previous section, we obtained a quantitative assessment of the main products of the glass−water reaction based on the number change of different chemical species with respect to the simulation time. Along with this chemical modification analysis, which provides useful information for examining the stoichiometry of the reaction, we investigated the MD trajectories to depict the hydroxylation mechanisms that happen at the sodium silicate glass surface, specifically focusing on the surface NBO sites.

that their contribution could be very different when sufficient amount of sodium ions is leached out. The chemical modification taking place at the interface can be further described through Figure 3b, which shows the magnified portion of Figure 3a (from t = 0 to 100 ps range). In addition to the time evolution of number of leached sodium ions and silanol units, that of undercoordinated Si atoms as well as water molecules that had penetrated beneath the surface are plotted. Within 1.5 ps of water contact, the simulation results show a rapid increase in the number of sodium ions that populate just above the surface boundaries. At the same time, water molecules quickly start to penetrate below the glass surface. This denotes that sodium−water interactions are built up and preceded prior to any protonation of the surface NBO sites. A slightly slower rate of the increase in SiOH formation during this short time period also indicates that water adsorption to the surface NBO sites develops after the sodium ions are surrounded by immersed water molecules. The details of this kind of hydroxylation mechanism which involves surface NBO sites are further elucidated in the following section (Section 3.2). Additionally, we inspected the silicate network structures before and after the glass−water interaction by examining the initial glass surface configuration and last trajectory of simulation. Figure 4a shows a part of the hydroxylated sodium

Figure 4. (a) Schematics of the sodium-leached/hydroxylated sodium silicate surface. (b) Species concentration change of the sodium silicate surface before and after (leached) water interaction. Note that the silanol concentration shown here (4.8 nm−2) is a value added up within the enclosed region, that span across 3 Å below the surface boundary (illustrated as dashed box in (a)).

silicate glass surface after removing all of the molecules near the surface region based on the bond-order cutoff value of 0.3 in ReaxFF analysis. Using these structures, we evaluated the concentration of the glass network and hydrous species. Upon 1 ns of leaching, the silanol concentration of sodium silicate glass was calculated as 5.79 ± 0.14 nm−2. This areal density is derived by including all of the silanols formed (terminal value at 1 ns, Figure 3a) during the glass−water reaction. As we mentioned earlier, pure silica exhibits an equilibrium concentration during water interaction. This equilibrium is reached when all of the surface defect sites (undercoordinated Si atoms and NBOs) on the silica surface are hydroxylated and MD simulation results reported that they are typically converged within 200−250 ps.36,38,58 While a decent agree15610

DOI: 10.1021/acs.jpcc.9b02940 J. Phys. Chem. C 2019, 123, 15606−15617

Article

The Journal of Physical Chemistry C

water and is replaced with the water solvation shell. This functions as one of water dissociation paths during the MD trajectory until the water solvation shell around sodium stabilizes and eventually allows sodium cations to diffuse further away from the surface and into the aqueous fluid. The higher population of sodium at the very proximate region of the glass surface (regions around 20 < z < 30 Å and 70 < z < 80 Å in Figure 2) therefore indicates that majority of the modifier cations remain bound at the surface, even though they may be defined as the “leached” components. This series of process is illustrated in Figure 6a,b, where the trajectory provides the same schematics of the hydroxylation of an NBO site but is associated with water dissociation and subsequent NaOH interactions.

Figure 5a−d illustrates how a surface NBO site is hydroxylated within the aqueous environment. Almost

Figure 5. Hydroxylation mechanism at the surface NBO sites. After the glass−water interaction establishes, a water molecule dissociates and protonates the NBO. The remaining dissociated water (OH−) participates as a proton recipient in consecutive proton transfer or hopping events observed with the surrounding water molecules. Color legend: Na (blue), Si (ivory), Oglass (gray), Owater (red), and H (white). Surface NBO is colored in orange, whereas the participants of the water dissociation and subsequent proton hopping are denoted in green. Note that ReaxFF only uses a single oxygen typedifferent oxygen colors are used here only to assist the visualization of the reactive events.

Figure 6. Hydroxylation of NBO sites by water dissociation and Na− OH formation. (a) Water adsorption at the NBO site, (b) NBO protonation (SiOH), and subsequent establishment of Na−OH interaction. Color legend: Na (blue), Si (ivory), Oglass (gray), Owater (red), and H (white). Surface NBO, adsorbed water, and Na+ interacting with OH− are colored in orange, green, and purple, respectively. Note that ReaxFF only uses a single oxygen and sodium atom typesdifferent oxygen and sodium colors are used here only to assist the visualization of the reactive events.

immediately after few picoseconds (typically within 3 ps, not shown in Figure 5a−d) of the simulation time, we could capture vast amount of water molecules that engage with sodium ions at the glass surface. During this development of Na+−water interaction, water molecules start to bind at the exposed NBO sites by hydrogen bonding and gradually replace SiO−Na to SiO−···H2O, as shown in Figure 5a with a dashed line. After molecular water is bound to the glass surface, the water molecule can chemisorb onto the NBO sites, resulting in the silanol formation (Figure 5b). In the meantime, the remainder of the dissociated water, a hydroxyl (OH−), requires stabilization by positive charge species from the surroundings. Our MD simulation results provide that the OH− products of the SiO−/water reaction can be neutralized from either H+ or Na+. Figure 5c shows the neutralization of OH− by proton donation from adjacent water molecules. The water molecule that initially engaged with the SiO− site is therefore converted back to molecular water, leaving a newly formed OH− nearby. This OH− transfer can further propagate from the interface to the water region, by consecutive water dissociation and proton hopping along the water chains.17 The arrows that are shown in Figure 5c,d indicate the direction of the proton jumps. While a solid reaction pathway that involves continuous water dissociation can successfully provide a neutralization of the OH− at the initial reactive site of the surface, the ReaxFF MD trajectory confirms that a second mechanism could alleviate the local negative charge of initially formed OH− groups at the corresponding water dissociated sites. In this case, the presence of Na+ is importantthese cations bind with OH− and form Na−OH at the interface. The Na−OH interaction is continued until the OH− grabs H+ from adjacent

3.3. Sodium Leaching and the Residence Time. It was noted in the Section 3.2 that the surface NBO sites associated with the sodium ions can induce water adsorption and dissociation, resulting in the formation of OH− groups interacting with nearby sodium ions. The Na−OH ion pairs typically have a fairly long lifetime until sufficient water molecules solvate the sodium ions and have them completely released to the water medium. Figure 7a presents one representative event of the sodium ion’s trajectory during the ReaxFF MD simulation. This graphical trajectory highlights three stages in the sodium leaching process. First, the sodium ions reside in the subsurface region, from which they diffuse to the glass−water interface when a channel opens up by leaching out sodium ions closer to the interface. Once the sodium ions have reached the surface region, they start to interact with water and associate with the OH− at the surface. During this ion-pair stage, sodium ions tend to reside at the surface for a certain amount of time. The existence of a plateau region at t = 215−310 ps as shown in Figure 7b explains that it is where the Si−OH···NaOH interaction is built and maintained until Na+ dissolution to the bulk water happens. It is notable that these sodium ions tend to reside at the surface for a certain time before they are completely released to the water medium. To explain this particular dynamic feature of the sodium behavior, we investigated the lifetime or “residence time” of the sodium ions. Here, the sodium residence time is calculated from the point when a sodium cation starts to enter the reference region, as shown in Figure 8a. This fixed region is set as the 5 Å range from the Gibbs 15611

DOI: 10.1021/acs.jpcc.9b02940 J. Phys. Chem. C 2019, 123, 15606−15617

Article

The Journal of Physical Chemistry C

analysis. The second type of leaching, on the other hand, was to consider only the sodium ions that have reached further away from the 5 Å reference region, after examining their final positions from the 1 ns of MD trajectories. Figure 8b shows the residence time of all of the leached sodium ions. The large number count at [900, 1000 ps] is consistent with the mass density profile in Figure 2 in that the majority of the sodium ions remain bound near the glass surface for almost the whole length of the MD simulations. Although these sodium ions are considered as being leached out from the bulk glass by the definition, it is hard to capture any relevant physical insights from this type of leached ions. On the other hand, Figure 8c presents the result of the same analysis but after ruling out the sodium ions located close to the surface. Among these leached ions, we found that one-third of sodium ions are released to solution within 200 ps of MD simulations. The shortest residence time of sodium ions was 22.50 ± 4.33 ps, whereas the longest residence time was 990.00 ± 2.87 ps. The histogram of the residence time in Figure 8c shows that the number count decays with respect to the increase in time variable. On the basis of this analysis, the release of sodium ions to the aqueous fluid can be regarded such that they are being forced out from the Na+-enriched region which is built close to the surface. This process occurs in a serial order that preceding ions are disposed to the solution and allow the rest to be followed. The time decay indicates the slow sodium leaching process, as the majority of sodium ions at the glass− water interface do not make it to the solution even after 1 ns of reactive MD simulations. The surfaces of our sodium silicate glass model consist of mixed Q2/Q3 species with sodium ions neutralizing the negative surface charge of superficial NBO sites. Ion-exchange processes result in the formation of NaOH, and surface ionization degree is expected to decrease in contact with water. At the same time, a fraction of the cation dissociates into the solution after the residence time as discussed above. In general, the amount of ionized groups depends on the surface type, pH, ionic strength, and type of cation, and inclusion of these aspects has become feasible through silica surface models with a realistic surface chemistry description. We have considered the same quantitative assessment to our simulation results by deriving the surface ionization change with respect to time (Figure 9a,b). As a result of Na+/H+ ion exchange, surface

Figure 7. (a) Representative trajectory of sodium ion and (b) its time evolution of relative position to the boundary at the top surface. Panels (c,d) show how the sodium ions exist at the surface and in the solution, respectively. Color legend: Na (blue), ONa−OH (yellow), Owater (red), and H (white). Note that Na−Owater bonds are illustrated just to emphasize their nearest coordination.

Figure 8. (a) Sodium trajectory and reference region (shaded in gray) for the calculation of residence time. The residence time histogram of (b) all of the leached sodium ions and (c) the leached sodium ions, in which their final positions at 1 ns were outside of the reference region.

dividing boundary. The end time of the residence, which would indicate that sodium ion is no longer bound at the surface, was then checked when a sodium ion had left the reference region and had migrated to the liquid medium. For the leached sodium ions, we evaluated their residence time according to these criteria and plotted it as a histogram (Figure 8b,c). In this analysis, we assorted the sodium ions into two different leaching types. The first group was chosen to include all of the sodium ions outside of the surface boundary for the

Figure 9. (a) Integral count of leached sodium ions on the sodium silicate−water interface at times 100, 250, 500, 750, and 1000 ps. z = 0 indicates one surface boundary and a second surface exists at z ≈ 54 Å. Dashed lines indicate the 5 Å-reference from surface boundaries to distinguish dissociated sodium ions to the solution. (b) Amount of dissolved Na+ per nm2 with respect to the ionization degree at the corresponding time frame. Error bars indicate the standard error from three glass−water reaction simulations. 15612

DOI: 10.1021/acs.jpcc.9b02940 J. Phys. Chem. C 2019, 123, 15606−15617

Article

The Journal of Physical Chemistry C ionization decreases to ∼52% within 100 ps of glass−water reaction. At 1 ns, the lowest degree of ionization is exhibited, and the total number of dissolved sodium cations increases from ∼0.25 to ∼0.98 nm−2. This trend is qualitatively consistent with the well-established silica surface models in that higher total charge density lead to close proximity of the majority of sodium ions at the surface and intermediately ionic surfaces (∼20%) show the maximum sodium dissociation. Within the duration of our simulation, however, an ionization degree of 20% or below could not be reached. Sodium leaching involves the dissociation of NaOH in the form of Na+ and OH−. The concentration of hydroxide (OH−) ions determines the pH value of the solution, and the surface exposed to the aqueous interface undergoes adjustments in surface chemistry. In this regard, the pH after 1 ns of glass− water interaction was computed by taking account of all of the hydroxide species in the solution. The average value of [OH−] = 0.46 ± 0.03 M which corresponds to pH ≈ 13.7 was derived based on the molecular analysis with a bond-order cutoff of 0.3. Experimental literature studies show the evidence of increase in surface ionization (negative surface charge density) toward higher pH but typically with an upper limit pH of 10.74−80 The surface ionization degree and pH of the current MD simulation results are both in the higher range and are subject to vary with longer simulations (>1 ns), which makes it difficult to directly correlate with available data from experiment or modeling. One possible way to eliminate these shortcomings would be an immersion of the reacted glass surfaces into a new solution of pH 7. A recent development of ReaxFF-reactive force field parameters for electrolytes (H+, Li+, Na+, K+, Cs+/OH−, F−, Cl−, I−)−water systems81 can also be introduced to model aqueous solution and probe surface ionization as a function of pH. Taking this into account, the current ReaxFF MD simulation arguably provides atomic-level monitoring of the ion-exchange and proton transfer reactions as well as dynamics over a nanosecond time scale. The refinement of atomistic force fields for silica in the recent years has enabled accurate computation of immersion energies of silica surfaces in water, adsorption isotherms, and contact angles to less than 5% deviation from experiment.31 ReaxFF also predicts these thermodynamic properties at a comparable agreement (see Supporting Information Section S6), which suggests that it can further be employed to understand and rationally explain the contribution of glass thermal history and surface exposal to wide range of aqueous environments on the alteration of surface properties. 3.4. Sodium and Proton Transport at the Glass− Water Interface. The MD simulations herein present another important aspect of sodium silicate glass−water interaction, which is the diffusion of protons through the glass structures. Figure 2 indicates that these protons can penetrate the surface, and we observed that the penetration depth can be as large as 6 Å during our simulation. The origin of these atoms can be tracked down from the trajectory, and we found two existing species that are mainly responsible for the diffusion within the bulk sodium silicates. First of all, the diffusion of the water molecule in its intact form was captured within our simulation duration. In this type of diffusion, water molecules could successfully find their way inward the silicate network, without having to chemically react with the network itself. As a high fraction of network modifier cations opens up the glass network, the presence of sodium ions in silicate glasses

provides possible channels which allow water molecules to penetrate beneath the surface. Figure 10a denotes the

Figure 10. Diffusion of hydrogen atoms to the bulk sodium silicate glasses (a) as in its intact form of molecular water, (b,c) by proton hopping between adjacent NBO sites(OA → OB, colored in purple).

schematic of the water diffusion that was captured in our simulation. The continuous line (colored in green) shows the trajectory path of the water molecule of interest, which had penetrated below the surface boundary (dashed line across the Figure 10a). Once the diffused water molecules have reached to a point when their mobility is limited, all of the molecules had dissociated into a proton and a hydroxyl. The oxygen speciation of OH− shown in Figure 4b confirms that the diffusion through water dissociation can be the dominant mechanism in confined spaces.38,65,82 Another type of diffusion can be found in the form of proton hopping from one NBO site to the nearest of the same kind, as portrayed in Figure 10b,c. After identifying different atomic species that are responsible for the transport within the bulk glass, we calculated the diffusion coefficients of each species and evaluated the environmental effects on diffusion. The atomic diffusion was described using the mean square displacement (MSD). The MSD was computed using the following equation 1 MSD = ⟨δr 2⟩ = ∑ [ri(t + t0) − ri(t0)]2 N i where N is the total number of ions and ri(t + t0) and ri(t0) are the displacements of the i ion at t + t0 and t0, respectively. In this work, we performed the glass−water reaction MD 15613

DOI: 10.1021/acs.jpcc.9b02940 J. Phys. Chem. C 2019, 123, 15606−15617

Article

The Journal of Physical Chemistry C

intact water to penetrate but channel size is small enough to constrain the diffusion. The water diffusion at the interface is therefore governed by structural confinement as well as the hydrogen bond network with surface silanol groups. Further diffusion of water species inside the bulk glass depends on the glass network, as observed from the trajectory path of intact water (Figure 8a). Another finding from the glass−water reaction simulation is the difference in sodium transport, from which we can observe that the mobility of sodium ions is being affected by their environment. The sodium ions interacting with OH− at the very surface diffuse slower because their diffusion occurs as whole unit (Na−OH ion pair) rather than the sodium ions alone. Once the Na−OH interaction is substituted to Na− water interaction and released to the solution (bulk water), their transport is less limited. Our simulation predicts DNa,solvation = 19.97 (±1.95) × 10−6 cm2/s for sodium ions in the water region, whereas DNa,interface = 16.26 (±1.64) × 10−6 cm2/s for sodium ions that populates within 5 Å from the surface boundary. In addition to the water and sodium diffusivity, we also evaluated proton transport. In sodium silicate glass−water interface systems, we find that protons diffuse faster than what they are predicted to do in pure silica− water systems. This feature can be explained with the number of NBO defects that exist within the bulk glass. As we mentioned previously, a high fraction of network modifier cations results in fragmentation of the glass network and produces NBOs within the bulk glass. These defects are where the proton transport occurs by discrete hopping between two adjacent NBO sites. This proton hopping, or transfer within the bulk, allows these atoms to migrate through the glass and result in an increased atomic diffusion coefficient. Pure silica glass contains almost no NBO defects in the bulk,84 and therefore, their diffusion to the bulk region may not be as fast as in sodium silicate glasses.

simulations in triplicate. The MSDs from each configuration were calculated and averaged to establish statistical significance. After we obtained the averaged MSD for each species, the slope of the MSD plot was extracted to determine the selfdiffusion coefficient (D). D = lim

t →∞

⟨δr 2⟩ 6t

To take account for different natures of atomic species, we used different protocols for distinguishing the types of atoms of interest. For the diffusivity of protons, we selected the hydrogen atoms that are in the form of SiOH, so that it could correctly derive the D value of proton hopping in the bulk glass system. In case of water diffusion, we first identified whether the water molecules had penetrated beneath the glass surface and then selected oxygen atoms of those molecules to investigate the diffusivity of water. The reason for calculating oxygen diffusivity is that the diffusion coefficient of oxygen in water (Ow) is considered to be the representative of its translational mobility. In case of sodium ion’s diffusivity, we grouped sodium ions into twosodium ions in water region and those within the interfacial layer. Figure 11 shows the MSD calculated for Na+, H+, and Ow in the sodium silicate−water interface at 300 K. Two reference

4. CONCLUSIONS In this work, reactive MD simulations with the ReaxFF potential were carried out to investigate the structural and dynamical features of the sodium silicate glass−water interface. Our nanosecond run of glass−water reactions using the ReaxFF MD framework provides sufficient information to derive various properties, which would otherwise be difficult to obtain from the short time scale of ab initio MD simulations. We primarily focused on the initial stages where the chemical modifications of the glass surface happen through the ionexchange process between sodium ions (Na+) from the glass and protons (H+) from water. From the simulation results, it is observed that the surface NBO sites are where the water reactivity is strong. With the presence of modifier cations, water molecules that are initially bound at these NBO sites dissociate into H+ and OH−. As a result, NBO sites are protonated and the remaining hydroxyls are captured by sodium ions. During this hydration of the surface, it was expected that the chemical reactions involving the Na+/H+ ion exchange would be stoichiometrically met. However, we found that small fraction of undercoordinated Si atoms are hydroxylated as well and causes a higher amount of silanol units than total number of leached sodium ions during the 1 ns of glass−water reaction. Because of the NBO protonation with the presence of sodium ions, we found that Na−OH interactions are mainly developed and populate at the glass−water interface. These

Figure 11. Average MSD of Na+, H+ and Ow (water). The dashed lines indicate the oxygen diffusivity in bulk water (colored in red, refs 38, 83) and proton diffusivity within the bulk silica (colored in green, ref 38).

lines, Ow in bulk water38,83 and H+ in bulk pure silica,38 were inserted to distinguish the transport behavior within different environments. Table 1 summarizes the diffusion coefficients Table 1. Diffusion Coefficient (D, 10−6 cm2 s−1) in the Sodium Silicate Glass−Water Interface property DNa in solvation DNa at the glass−water interface DH of SiOH DO of H2O within sodium silicate bulk

19.97 16.26 11.67 9.71

± ± ± ±

1.95 1.64 0.57 0.79

calculated from the slope of the MSDs. Based on the diffusion coefficients, we can derive several distinctive features of the transport at the interface. One of the main characteristics would be the confined water diffusivity within the glass when compared to their transport in bulk systems. The decrease in water diffusion is intimately related with the glass surface structure, where the fragmentation of the glass network allows 15614

DOI: 10.1021/acs.jpcc.9b02940 J. Phys. Chem. C 2019, 123, 15606−15617

The Journal of Physical Chemistry C

■

interactions are continued along the trajectory before sufficient water molecules solvate the sodium ions and cause them to be completely released to the water medium. To elucidate their diffusion behavior at the interface, we first investigated the residence time of the leached sodium ions to provide adequate information on how the structural modifications affect their behaviors. In addition, the transport property of different species was investigated from the MD trajectory. We calculated the diffusivity of leached Na+ in solution, leached Na+ but bound at the interface as in NaOH, hydrogen atoms of SiOH, and oxygen atoms of water that had penetrated beneath glass surface. For the diffusion of hydrogen atoms, the origin of their diffusion can be attributed to the proton hopping within two adjacent NBO sites and opening up of the glass network due to the high fraction of modifiers, which eventually allows diffusion of water molecule in its intact form. All of the detailed assessment of surface chemistry and dynamics of modifier ions from the large-scale ReaxFF reactive MD simulation suggest that it can successfully be used toward understanding complex glass−water systems which would otherwise be difficult for first-principles or experimental studies to evaluate. ReaxFF also resolves major limitations of prior empirical force fields by taking account of the reactivity, and thus, the prediction of glass surface chemical modification and accurate computation of aqueous interfacial properties have become feasible. The augmentation of this reactive potential with a surface model database31 would open up the opportunities to assess different adsorption behavior of water, biomolecules, and polymers by customizing the silicate glass surfaces with different thermal treatments and adjusting pH conditions of aqueous solutions.

■

REFERENCES

(1) Cailleteau, C.; Angeli, F.; Devreux, F.; Gin, S.; Jestin, J.; Jollivet, P.; Spalla, O. Insight into Silicate-Glass Corrosion Mechanisms. Nat. Mater. 2008, 7, 978−983. (2) Gin, S.; Guittonneau, C.; Godon, N.; Neff, D.; Rebiscoul, D.; Cabié, M.; Mostefaoui, S. Nuclear Glass Durability: New Insight into Alteration Layer Properties. J. Phys. Chem. C 2011, 115, 18696− 18706. (3) Gin, S.; Jollivet, P.; Fournier, M.; Angeli, F.; Frugier, P.; Charpentier, T. Origin and Consequences of Silicate Glass Passivation by Surface Layers. Nat. Commun. 2015, 6, 6360. (4) McGrail, B. P.; Kumar, A.; Day, D. E. Sodium Diffusion and Leaching of Simulated Nuclear Waste Glass. J. Am. Ceram. Soc. 1984, 67, 463−467. (5) Hench, L. L.; Xynos, I. D.; Polak, J. M. Bioactive Glasses for in Situ Tissue Regeneration. J. Biomater. Sci., Polym. Ed. 2004, 15, 543− 562. (6) Tilocca, A. Structural Models of Bioactive Glasses from Molecular Dynamics Simulations. Proc. R. Soc. A 2009, 465, 1003− 1027. (7) Bradley, L. C.; Dilworth, Z. R.; Barnette, A. L.; Hsiao, E.; Barthel, A. J.; Pantano, C. G.; Kim, S. H. Hydronium Ions in SodaLime Silicate Glass Surfaces. J. Am. Ceram. Soc. 2013, 96, 458−463. (8) He, H.; Qian, L.; Pantano, C. G.; Kim, S. H. Mechanochemical Wear of Soda Lime Silica Glass in Humid Environments. J. Am. Ceram. Soc. 2014, 97, 2061−2068. (9) Maekawa, H.; Saito, T.; Yokokawa, T. Water in Silicate Glass: 17O NMR of Hydrous Silica, Albite, and Na2Si4O9 Glasses. J. Phys. Chem. B 1998, 102, 7523−7529. (10) Sprenger, D.; Bach, H.; Meisel, W.; Gütlich, P. XPS Study of Leached Glass Surfaces. J. Non-Cryst. Solids 1990, 126, 111−129. (11) Shchukarev, A.; Rosenqvist, J.; Sjöberg, S. XPS Study of the Silica−Water Interface. J. Electron Spectrosc. Relat. Phenom. 2004, 137−140, 171−176. (12) Shen, Y. R. Surface Properties Probed by Second-Harmonic and Sum-Frequency Generation. Nature 1989, 337, 519. (13) Ong, S.; Zhao, X.; Eisenthal, K. B. Polarization of Water Molecules at a Charged Interface: Second Harmonic Studies of the Silica/Water Interface. Chem. Phys. Lett. 1992, 191, 327−335. (14) Barnette, A. L.; Kim, S. H. Coadsorption of N-Propanol and Water on SiO2: Study of Thickness, Composition, and Structure of Binary Adsorbate Layer Using Attenuated Total Reflection Infrared (ATR-IR) and Sum Frequency Generation (SFG) Vibration Spectroscopy. J. Phys. Chem. C 2012, 116, 9909−9916. (15) Masini, P.; Bernasconi, M. Ab Initio Simulations of Hydroxylation and Dehydroxylation Reactions at Surfaces: Amorphous Silica and Brucite. J. Phys.: Condens. Matter 2002, 14, 4133− 4144. (16) Mischler, C.; Horbach, J.; Kob, W.; Binder, K. Water Adsorption on Amorphous Silica Surfaces: A Car−Parrinello Simulation Study. J. Phys.: Condens. Matter 2005, 17, 4005−4013. (17) Tilocca, A.; Cormack, A. N. Modeling the Water−Bioglass Interface by Ab Initio Molecular Dynamics Simulations. ACS Appl. Mater. Interfaces 2009, 1, 1324−1333. (18) van Beest, B. W. H.; Kramer, G. J.; van Santen, R. A. Force Fields for Silicas and Aluminophosphates Based on Ab Initio Calculations. Phys. Rev. Lett. 1990, 64, 1955−1958. (19) Cormack, A. N.; Du, J.; Zeitler, T. R. Alkali Ion Migration Mechanisms in Silicate Glasses Probed by Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2002, 4, 3193−3197. (20) Pedone, A.; Malavasi, G.; Menziani, M. C.; Cormack, A. N.; Segre, U. A New Self-Consistent Empirical Interatomic Potential Model for Oxides, Silicates, and Silica-Based Glasses. J. Phys. Chem. B 2006, 110, 11780−11795. (21) Sanders, M. J.; Leslie, M.; Catlow, C. R. A. Interatomic Potentials for SiO2. J. Chem. Soc., Chem. Commun. 1984, 1271−1273. (22) Tilocca, A.; de Leeuw, N. H.; Cormack, A. N. Shell-Model Molecular Dynamics Calculations of Modified Silicate Glasses. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 104209.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.9b02940.

■

Article

Details of the bulk sodium silicate glass construction, glass surface construction and water system evaluation, Qn-species distribution of the bulk and surface of the glass model, surface energy of the sodium silicate glass before and after 1 ns of glass−water reaction, reliability evaluation of ReaxFF for silica−water systems, and ReaxFF reactive force field parameters for Na/Si/O/H (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Seung Ho Hahn: 0000-0002-1722-5631 Adri C. T. van Duin: 0000-0002-3478-4945 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the National Science Foundation (grant no. DMR-1609107). Computations for this research were performed on the Pennsylvania State University’s Institute for Cyber Science Advanced Cyber Infrastructure (ICS-ACI). 15615

DOI: 10.1021/acs.jpcc.9b02940 J. Phys. Chem. C 2019, 123, 15606−15617

Article

The Journal of Physical Chemistry C (23) Tilocca, A.; de Leeuw, N. H. Structural and Electronic Properties of Modified Sodium and Soda-Lime Silicate Glasses by Car−Parrinello Molecular Dynamics. J. Mater. Chem. 2006, 16, 1950−1955. (24) Pedone, A. Properties Calculations of Silica-Based Glasses by Atomistic Simulations Techniques: A Review. J. Phys. Chem. C 2009, 113, 20773−20784. (25) Feuston, B. P.; Garofalini, S. H. Oligomerization in Silica Sols. J. Phys. Chem. 1990, 94, 5351−5356. (26) de Leeuw, N. H.; Higgins, F. M.; Parker, S. C. Modeling the Surface Structure and Stability of α-Quartz. J. Phys. Chem. B 1999, 103, 1270−1277. (27) Du, J.; Cormack, A. N. Molecular Dynamics Simulation of the Structure and Hydroxylation of Silica Glass Surfaces. J. Am. Ceram. Soc. 2005, 88, 2532−2539. (28) Pedone, A.; Malavasi, G.; Menziani, M. C.; Segre, U.; Musso, F.; Corno, M.; Civalleri, B.; Ugliengo, P. FFSiOH: A New Force Field for Silica Polymorphs and Their Hydroxylated Surfaces Based on Periodic B3LYP Calculations. Chem. Mater. 2008, 20, 2522−2531. (29) Patwardhan, S. V.; Emami, F. S.; Berry, R. J.; Jones, S. E.; Naik, R. R.; Deschaume, O.; Heinz, H.; Perry, C. C. Chemistry of Aqueous Silica Nanoparticle Surfaces and the Mechanism of Selective Peptide Adsorption. J. Am. Chem. Soc. 2012, 134, 6244−6256. (30) Heinz, H.; Lin, T.-J.; Kishore Mishra, R.; Emami, F. S. Thermodynamically Consistent Force Fields for the Assembly of Inorganic, Organic, and Biological Nanostructures: The INTERFACE Force Field. Langmuir 2013, 29, 1754−1765. (31) Emami, F. S.; Puddu, V.; Berry, R. J.; Varshney, V.; Patwardhan, S. V.; Perry, C. C.; Heinz, H. Force Field and a Surface Model Database for Silica to Simulate Interfacial Properties in Atomic Resolution. Chem. Mater. 2014, 26, 2647−2658. (32) Emami, F. S.; Puddu, V.; Berry, R. J.; Varshney, V.; Patwardhan, S. V.; Perry, C. C.; Heinz, H. Prediction of Specific Biomolecule Adsorption on Silica Surfaces as a Function of PH and Particle Size. Chem. Mater. 2014, 26, 5725−5734. (33) Heinz, H. Adsorption of Biomolecules and Polymers on Silicates, Glasses, and Oxides: Mechanisms, Predictions, and Opportunities by Molecular Simulation. Curr. Opin. Chem. Eng. 2016, 11, 34−41. (34) Heinz, H.; Ramezani-Dakhel, H. Simulations of Inorganic− Bioorganic Interfaces to Discover New Materials: Insights, Comparisons to Experiment, Challenges, and Opportunities. Chem. Soc. Rev. 2016, 45, 412−448. (35) Mahadevan, T. S.; Garofalini, S. H. Dissociative Chemisorption of Water onto Silica Surfaces and Formation of Hydronium Ions. J. Phys. Chem. C 2008, 112, 1507−1515. (36) Mahadevan, T. S.; Du, J. Evaluating Water Reactivity at Silica Surfaces Using Reactive Potentials. J. Phys. Chem. C 2018, 122, 9875− 9885. (37) van Duin, A. C. T.; Strachan, A.; Stewman, S.; Zhang, Q.; Xu, X.; Goddard, W. A. ReaxFFSiO Reactive Force Field for Silicon and Silicon Oxide Systems. J. Phys. Chem. A 2003, 107, 3803−3811. (38) Fogarty, J. C.; Aktulga, H. M.; Grama, A. Y.; van Duin, A. C. T.; Pandit, S. A. A Reactive Molecular Dynamics Simulation of the SilicaWater Interface. J. Chem. Phys. 2010, 132, 174704. (39) Geneste, G.; Bouyer, F.; Gin, S. Hydrogen−Sodium Interdiffusion in Borosilicate Glasses Investigated from First Principles. J. Non-Cryst. Solids 2006, 352, 3147−3152. (40) Zapol, P.; He, H.; Kwon, K. D.; Criscenti, L. J. First-Principles Study of Hydrolysis Reaction Barriers in a Sodium Borosilicate Glass. Int. J. Appl. Glass Sci. 2013, 4, 395−407. (41) Pierce, E. M.; Frugier, P.; Criscenti, L. J.; Kwon, K. D.; Kerisit, S. N. Modeling Interfacial Glass-Water Reactions: Recent Advances and Current Limitations. Int. J. Appl. Glass Sci. 2014, 5, 421−435. (42) Hahn, S. H.; Rimsza, J.; Criscenti, L.; Sun, W.; Deng, L.; Du, J.; Liang, T.; Sinnott, S. B.; van Duin, A. C. T. Development of a ReaxFF Reactive Force Field for NaSiOx/Water Systems and Its Application to Sodium and Proton Self-Diffusion. J. Phys. Chem. C 2018, 122, 19613−19624.

(43) Yu, Y.; Wang, B.; Wang, M.; Sant, G.; Bauchy, M. Reactive Molecular Dynamics Simulations of Sodium Silicate Glasses  Toward an Improved Understanding of the Structure. Int. J. Appl. Glass Sci. 2017, 8, 276−284. (44) Dongol, R.; Wang, L.; Cormack, A. N.; Sundaram, S. K. Molecular Dynamics Simulation of Sodium Aluminosilicate Glass Structures and Glass Surface-Water Reactions Using the Reactive Force Field (ReaxFF). Appl. Surf. Sci. 2018, 439, 1103−1110. (45) 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. (46) Senftle, T. P.; Hong, S.; Islam, M. M.; Kylasa, S. B.; Zheng, Y.; Shin, Y. K.; Junkermeier, C.; Engel-Herbert, R.; Janik, M. J.; Aktulga, H. M.; et al. The ReaxFF Reactive Force-Field: Development, Applications and Future Directions. npj Comput. Mater. 2016, 2, 15011. (47) Shin, Y. K.; Kwak, H.; Zou, C.; Vasenkov, A. V.; van Duin, A. C. T. Development and Validation of a ReaxFF Reactive Force Field for Fe/Al/Ni Alloys: Molecular Dynamics Study of Elastic Constants, Diffusion, and Segregation. J. Phys. Chem. A 2012, 116, 12163− 12174. (48) Islam, M. M.; Ostadhossein, A.; Borodin, O.; Yeates, A. T.; Tipton, W. W.; Hennig, R. G.; Kumar, N.; van Duin, A. C. T. ReaxFF Molecular Dynamics Simulations on Lithiated Sulfur Cathode Materials. Phys. Chem. Chem. Phys. 2015, 17, 3383−3393. (49) Yoon, K.; Rahnamoun, A.; Swett, J. L.; Iberi, V.; Cullen, D. A.; Vlassiouk, I. V.; Belianinov, A.; Jesse, S.; Sang, X.; Ovchinnikova, O. S.; et al. Atomistic-Scale Simulations of Defect Formation in Graphene under Noble Gas Ion Irradiation. ACS Nano 2016, 10, 8376−8384. (50) Hong, S.; Sheng, C.; Krishnamoorthy, A.; Rajak, P.; Tiwari, S.; Nomura, K.-i.; Misawa, M.; Shimojo, F.; Kalia, R. K.; Nakano, A.; et al. Chemical Vapor Deposition Synthesis of MoS2 Layers from the Direct Sulfidation of MoO3 Surfaces Using Reactive Molecular Dynamics Simulations. J. Phys. Chem. C 2018, 122, 7494−7503. (51) Mortier, W. J.; Ghosh, S. K.; Shankar, S. ElectronegativityEqualization Method for the Calculation of Atomic Charges in Molecules. J. Am. Chem. Soc. 1986, 108, 4315−4320. (52) 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 2008, 112, 1040−1053. (53) Yuan, X.; Cormack, A. N. Local Structures of MD-Modeled Vitreous Silica and Sodium Silicate Glasses. J. Non-Cryst. Solids 2001, 283, 69−87. (54) Du, J.; Cormack, A. N. The Medium Range Structure of Sodium Silicate Glasses: A Molecular Dynamics Simulation. J. NonCryst. Solids 2004, 349, 66−79. (55) Zhang, W.; van Duin, A. C. T. Improvement of the ReaxFF Description for Functionalized Hydrocarbon/Water Weak Interactions in the Condensed Phase. J. Phys. Chem. B 2018, 122, 4083− 4092. (56) Baerends, E. J.; Ziegler, T.; Atkins, A. J.; Autschbach, J.; Bashford, D.; Baseggio, O.; Bérces, A.; Bickelhaupt, F. M.; Bo, C.; Boerritger, P. M.; et al. ADF2018, SCM, Theoretical Chemistry; Vrije Universiteit: Amsterdam, The Netherlands, 2018; https://www.scm. com. (57) Webb, E. B.; Garofalini, S. H. Relaxation of Silica Glass Surfaces Before and After Stress Modification in a Wet and Dry Atmosphere: Molecular Dynamics Simulations. J. Non-Cryst. Solids 1998, 226, 47− 57. (58) Yeon, J.; van Duin, A. C. T. ReaxFF Molecular Dynamics Simulations of Hydroxylation Kinetics for Amorphous and NanoSilica Structure, and Its Relations with Atomic Strain Energy. J. Phys. Chem. C 2016, 120, 305−317. (59) Luo, J.; He, H.; Podraza, N. J.; Qian, L.; Pantano, C. G.; Kim, S. H. Thermal Poling of Soda-Lime Silica Glass with Nonblocking ElectrodesPart 1: Effects of Sodium Ion Migration and Water Ingress on Glass Surface Structure. J. Am. Ceram. Soc. 2016, 99, 1221−1230. 15616

DOI: 10.1021/acs.jpcc.9b02940 J. Phys. Chem. C 2019, 123, 15606−15617

Article

The Journal of Physical Chemistry C (60) Luo, J.; Banerjee, J.; Pantano, C. G.; Kim, S. H. Vibrational Sum Frequency Generation Spectroscopy Study of Hydrous Species in Soda Lime Silica Float Glass. Langmuir 2016, 32, 6035−6045. (61) Sokhan, B. V. P.; Tildesley, D. J. The Free Surface of Water: Molecular Orientation, Surface Potential and Nonlinear Susceptibility. Mol. Phys. 1997, 92, 625−640. (62) Lei, X.; Jee, Y.; Huang, K. Amorphous Na2Si2O5 as a Fast Na+ Conductor: An Ab Initio Molecular Dynamics Simulation. J. Mater. Chem. A 2015, 3, 19920−19927. (63) Leed, E. A.; Pantano, C. G. Computer Modeling of Water Adsorption on Silica and Silicate Glass Fracture Surfaces. J. Non-Cryst. Solids 2003, 325, 48−60. (64) Walsh, T. R.; Wilson, M.; Sutton, A. P. Hydrolysis of the Amorphous Silica Surface. II. Calculation of Activation Barriers and Mechanisms. J. Chem. Phys. 2000, 113, 9191−9201. (65) Rimsza, J. M.; Du, J. Ab Initio Molecular Dynamics Simulations of the Hydroxylation of Nanoporous Silica. J. Am. Ceram. Soc. 2015, 98, 3748−3757. (66) Leed, E. A.; Sofo, J. O.; Pantano, C. G. Electronic Structure Calculations of Physisorption and Chemisorption on Oxide Glass Surfaces. Phys. Rev. B: Condens. Matter Mater. Phys. 2005, 72, 155427. (67) Rimola, A.; Ugliengo, P. A quantum mechanical study of the reactivity of (SiO)2-defective silica surfaces. J. Chem. Phys. 2008, 128, 204702. (68) 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 2016, 120, 24803−24816. (69) Zhuravlev, L. T.; Potapov, V. V. Density of Silanol Groups on the Surface of Silica Precipitated from a Hydrothermal Solution. Russ. J. Phys. Chem. 2006, 80, 1119−1128. (70) Banerjee, J.; Kim, S. H.; Pantano, C. G. Elemental Areal Density Calculation and Oxygen Speciation for Flat Glass Surfaces Using X-Ray Photoelectron Spectroscopy. J. Non-Cryst. Solids 2016, 450, 185−193. (71) Sheth, N.; Ngo, D.; Banerjee, J.; Zhou, Y.; Pantano, C. G.; Kim, S. H. Probing Hydrogen-Bonding Interactions of Water Molecules Adsorbed on Silica, Sodium Calcium Silicate, and Calcium Aluminosilicate Glasses. J. Phys. Chem. C 2018, 122, 17792−17801. (72) Amma, S.-i.; Kim, S. H.; Pantano, C. G. Analysis of Water and Hydroxyl Species in Soda Lime Glass Surfaces Using Attenuated Total Reflection (ATR)-IR Spectroscopy. J. Am. Ceram. Soc. 2016, 99, 128−134. (73) Amma, S.-i.; Luo, J.; Pantano, C. G.; Kim, S. H. Specular Reflectance (SR) and Attenuated Total Reflectance (ATR) Infrared (IR) Spectroscopy of Transparent Flat Glass Surfaces: A Case Study for Soda Lime Float Glass. J. Non-Cryst. Solids 2015, 428, 189−196. (74) Bolt, G. H. Determination of the Charge Density of Silica Sols. J. Phys. Chem. 1957, 61, 1166−1169. (75) Tadros, T. F.; Lyklema, J. Adsorption of Potential-Determining Ions at the Silica-Aqueous Electrolyte Interface and the Role of Some Cations. J. Electroanal. Chem. Interfacial Electrochem. 1968, 17, 267− 275. (76) Milonjić, S. K. Determination of Surface Ionization and Complexation Constants at Colloidal Silica/Electrolyte Interface. Colloids Surf. 1987, 23, 301−312. (77) House, W. A.; Orr, D. R. Investigation of the PH Dependence of the Kinetics of Quartz Dissolution at 25 °C. J. Chem. Soc., Faraday Trans. 1992, 88, 233−241. (78) Zhuravlev, L. T. Surface Characterization of Amorphous Silicaa Review of Work from the Former USSR. Colloids Surf., A 1993, 74, 71−90. (79) Sonnefeld, J. Determination of Surface Charge Density Constants for Spherical Silica Particles Using a Linear Transformation. J. Colloid Interface Sci. 1996, 183, 597−599. (80) Zhuravlev, L. T. The surface chemistry of amorphous silica. Zhuravlev model. Colloids Surf., A 2000, 173, 1−38. (81) Fedkin, M. V.; Shin, Y. K.; Dasgupta, N.; Yeon, J.; Zhang, W.; van Duin, D.; van Duin, A. C. T.; Mori, K.; Fujiwara, A.; Machida, M.;

et al. Development of the ReaxFF Methodology for Electrolyte− Water Systems. J. Phys. Chem. A 2019, 123, 2125−2141. (82) Bakos, T.; Rashkeev, S. N.; Pantelides, S. T. H2O and O2 Molecules in Amorphous SiO2: Defect Formation and Annihilation Mechanisms. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 195206. (83) Krynicki, K.; Green, C. D.; Sawyer, D. W. Pressure and Temperature Dependence of Self-Diffusion in Water. Faraday Discuss. Chem. Soc. 1978, 66, 199−208. (84) Feuston, B. P.; Garofalini, S. H. Empirical Three-body Potential for Vitreous Silica. J. Chem. Phys. 1988, 89, 5818−5824.

15617

DOI: 10.1021/acs.jpcc.9b02940 J. Phys. Chem. C 2019, 123, 15606−15617