Development of a ReaxFF Reactive Force Field for ... - ACS Publications

Article Cite This: J. Phys. Chem. C 2018, 122, 19613−19624

pubs.acs.org/JPCC

Development of a ReaxFF Reactive Force Field for NaSiOx/Water Systems and Its Application to Sodium and Proton Self-Diffusion Seung Ho Hahn,† Jessica Rimsza,§ Louise Criscenti,§ Wei Sun,∥ Lu Deng,∥ Jincheng Du,∥ Tao Liang,‡ Susan B. Sinnott,‡ and Adri C. T. van Duin*,†,‡ Department of Mechanical Engineering and ‡Department of Materials Science and Engineering, The Pennsylvania State University, University Park, Pennsylvania 16802, United States § Geochemistry Department, Sandia National Laboratories, Albuquerque, New Mexico 87185, United States ∥ Department of Materials Science and Engineering, University of North Texas, Denton, Texas 76203, United States Downloaded via KAOHSIUNG MEDICAL UNIV on October 8, 2018 at 20:46:12 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

S Supporting Information *

ABSTRACT: We present the development of a ReaxFF reactive force field for Na/Si/O/H interactions, which enables reactive molecular dynamics simulation of the sodium silicate−water interfaces. The force field parameters were fitted against various quantum mechanical calculations, including equations of state of different NaSiOx crystalline phases, energy barriers of a sodium cation’s transport within the sodium silicate crystal structure, interactions between the hydroxylated silica surface and sodium cation−water systems, and dissociation energies of [NaOH−n(H2O)] (n = 1−6) clusters. After the optimization process, we validated the force field capability through calculating the structures of sodium silicate crystals and glasses and transport properties of sodium ions and protons within the amorphous structures. The force field was also applied to validate the dissociation behavior of sodium hydroxides within the bulk water. Our results with the developed force field are relevant to detailed chemical dissolution mechanisms, which involve (a) the interdiffusion process of sodium ions from glasses and protons from water, (b) subsequent ionic self-diffusion of sodium ions from the subsurface region to vacancy sites at the glass−water interface, and (c) sodium ions interaction with water after leaching from the amorphous sodium silicate system. of improved and abundant computational resources.11−14 In particular, molecular dynamics (MD) simulation techniques have led to a better understanding of the structural features of silicate glasses.15,16 Although investigations of bulk structures within the conventional classical MD framework are now relatively common and reliable, there exist fairly few MD studies that also considered surface phenomena related to water,13,17,18 where complicated reactive processes contribute. Within the spatial and temporal scalability of MD simulations, the limitation of conventional empirical potentials in modeling chemical reactions, specifically the water dissociation and hydrolysis reaction at the surface, can often be resolved by using reactive potentials. ReaxFF,19,20 a reactive force field that is based on a bond-order concept, has been extensively applied in MD simulations of silica materials.21−24 The connectivity-dependent interaction terms in this reactive potential are based on a calculation of bond order, ensuring that their energy contributions smoothly decay and disappear

1. INTRODUCTION The chemical reactivity of multicomponent silicate glass surfaces under humid or aqueous environment has been an important research area for years. Experimental findings showed that exposure to humidity results in a decrease in surface energy, thereby making the glass surface vulnerable to crack propagation.1−3 Also, the formation of complex water− glass interfacial structures and subsequent gel layers on the surface of silicate materials during corrosion are critical to the long-term dissolution behaviors of these materials,4,5 especially when they are used as a host to safely dispose nuclear and other wastes.6−8 As such, the presence of water affects the behaviors of glass surfaces, making their mechanical and chemical durability environment-dependent and difficult to predict. Despite many outstanding experimental endeavors to connect chemical reactivity to the surface properties,9,10 the intrinsic difficulty in obtaining accurate atomistic structural and energy data on multicomponent noncrystalline materials hampers the formulation of detailed atomistic level mechanisms. In recent years, computer simulation methods that can provide such details have been developed, utilizing the advent © 2018 American Chemical Society

Received: June 19, 2018 Revised: August 7, 2018 Published: August 8, 2018 19613

DOI: 10.1021/acs.jpcc.8b05852 J. Phys. Chem. C 2018, 122, 19613−19624

Article

The Journal of Physical Chemistry C

silicate−water description. We show that ReaxFF can reproduce the bulk structure properties of both crystalline and amorphous sodium silicate systems and provide predictions on the sodium ion pathways within and at the sodium silicate−water interface. We believe that the applicability of this force field is diverse and it will be a useful basis for the study of surface characterization of sodium-containing silicate glasses.

upon elongation of interatomic distances. The implementation of bond-order formalism distinguishes ReaxFF from other nonreactive force fields or conventional MD methods and enables the description of chemical reactions, such as proton transfer from reactive water to the silica surface or ion exchange of the alkali component with protons at the multicomponent silicate glass−water interface, within a MD framework. One of the earliest ReaxFF reactive force field that described Si/O systems and their hydrolysis reaction is the one developed by van Duin et al.25 This force field has been successful in its application to silicon/silicon oxide crystal systems and demonstrated ReaxFF’s capability for reproducing the condensation reaction of adjacent Si−OH groups, which are prevalent at the surface of the silica glasses. Fogarty et al.26 expanded this set of Si/O/H parameters for the amorphous silica−water interface description and provided excellent agreement in structural, chemical, and electrostatic properties of the hydroxylated silica surface, while retaining the bulk properties for both amorphous silica and water. This specific force field gained credibility and has been widely employed in related fields of study.27−29 In an attempt to provide a more precise description of the surface hydroxylation, Yeon et al.30 put effort into reparameterizing the Fogarty et al. force field by including the hydrolysis reaction energy barrier for strained and unstrained Si−O structures. Although focusing on these surface reactions caused this alternative force field to develop a high number of 5-fold Si defect (Q5) in bulk silica, it offered a great surface energy prediction of hydroxylated amorphous silica.21 Meanwhile, a separate force field that was originally developed for smectite clay−zeolite composites and their dynamics with confined water by Pitman et al.31 has also proved its validity in reproducing the structural, mechanical, and dynamic properties of amorphous silica.22−24 All these parameter developments and MD simulation studies show that the ReaxFF method has been successfully applied to many challenges related to both bulk and surface science of glass materials.21,32 However, there is a significant need for more force field development beyond silica−water interactions, since the use of ReaxFF with available parameters is limited to the amorphous silica−water system or to bulk silica at this stage.33,34 Consequently, further extension of the ReaxFF force field to multicomponent glasses and their surface reactivity is highly desirable. The addition of an alkali cation component, for example sodium, is the first but critical step to address this issue and enhances the versatility of ReaxFF. The presence of alkali ions in silica network-forming multicomponent glasses reduces the connectivity of the glass network and generates terminal oxygen groups (nonbridging oxygen, NBO) that are bound to modifier cations. Therefore, hydration at the surface of alkali-containing silicate glasses is preceded by the ion-exchange process, an ionic interdiffusion between alkali ions in the glass and protons from water. Further surface hydroxylation of such materials would undergo the same process of that of an amorphous silica, which involves the breaking of tetrahedrally bridged Si−O−Si network. These series of reactions that happens at the glass surface requires a reactive potential to specifically focus on the reaction mechanisms and the energetics involved for water dissociation as well as Na−NBO bond breaking and Na−OH formation. Here, we present the development and validation of the Na/ Si/O/H ReaxFF reactive force field, an extension of the Pitman Si/O/H force field parameter set31 to the sodium

2. METHODS 2.1. ReaxFF Reactive Force Field. The ReaxFF reactive force field is a bond-order-dependent potential. The connectivity between all atom pairs is therefore not predefined, rather it is updated and modified at every MD time step. This feature allows ReaxFF to describe bond formation and dissociation, which otherwise would be best studied by quantum mechanics (QM), at a comparable computational efficiency of conventional classical MD simulations. The general description of total interaction energy in ReaxFF is as follows Esys = E bond + Eval + Etor + Eover + Eunder + E lp + EvdWaals + ECoulomb

(1)

In eq 1, Ebond (bond energy), Eval (valence angle energy), Etor (torsion angle energy), Eover (over-coordination penalty energy), Eunder (under-coordination penalty energy), and Elp (lone-pair energy) are bond-order-dependent terms, meaning that the contribution of these energies vanishes upon bond dissociation, leaving only the nonbonded interactions to be taken into account afterward. Nonbonded interaction terms include ECoulomb (Coulomb energy) and EvdWaals (van der Waals energy), and they are calculated between all atom pairs in the system. For both types of interactions, ReaxFF employs shielded potentials to avoid excessive repulsions and unphysical charges between atoms that are at close distances. The electron equilibration method (EEM)35 is used for calculations of the atomic charges in ReaxFF, and this geometry-dependent charge calculation method makes ReaxFF a reactive force field that can handle the polarization of a system. 2.2. Optimization Procedure of the ReaxFF Reactive Force Field. Application of empirical force field methods for MD simulation requires the incorporation of a reliable parameter set because the quality of the force field parameter set is directly related to the accuracy of simulation results. The parameter fitting, often referred to as the force field optimization process, is performed against a relevant training set that consists of quantum mechanical (QM) calculations and/or experimental data. In the ReaxFF reactive force field method, parameter optimizations are carried out to minimize the total error represented as eq 2, by means of the oneparameter search technique. Further details of this optimization procedure can be found in van Duin et al.36 i Xi ,QM − Xi ,ReaxFF yz zz zz σ i k {

∑ jjjjj n

total error =

i=1

2

(2)

where n represents the total number of data points included in the training set, Xi,QM and Xi,ReaxFF are the values obtained from QM and ReaxFF calculations, and σi is the assigned weight value for each data point. 19614


Article


the training set. Optimizing the parameter set to these data enables ReaxFF to properly calculate the energies associated with a hydroxylated silicate surface and a nearby sodium ion. In this study, we used a system that consists of a total of 14 water molecules to fully hydrate the sodium ion and hydrogen bond with all Si−OH (silanol) groups of the Si(OH)4 monomer. Two different Na+−Si(OH)4 complexes, a monodentate and a bidentate complex, were considered for this Si(OH)4, Na+, and 14H2O systems. For the monodentate structure, a configuration where the sodium ion is aligned with one of the silanol groups, the sodium to proton-terminated oxygen distance (dNa−O) was decreased from 4.26 to 2.08 Å in increments of 0.02 Å. Energy calculations at each point were provided as energy scan data with respect to Na−O distance. For the bidentate structure, the sodium ion binds to two oxygen atoms of the Si(OH)4 monomer. As the sodium ion approaches the monomer, an energy scan was performed by decreasing the silicon to sodium distance (dSi−Na) from 4.57 to 3.05 Å by 0.02 Å. All of the density functional theory (DFT) calculations used Gaussian 0342 program with the B3LYP functional and 6-311-G(d,p) basis set. Throughout the energy scans, all four O−Si−O angles in Si(OH)4 were constrained to mimic the hydroxylated silica surface, which exhibits less flexibility in bond length and angles than a hydrated monomer. Also, the position of the Si atom was fixed throughout the calculations to localize the hydrogen bond rearrangement around the sodium ion. 2.3.4. Sodium−Water Binding Energies and Hydration of Sodium Hydroxide by Water Molecules. For the correct sodium−water description, we investigated various structures of Na(OH)−n(H2O) (where n = 1−6) clusters at the B3LYP/ 6-311-G(d,p) level of theory. The sodium to oxygen of the hydroxyl distance was constrained, to provide the energetics related to the sodium and hydroxyl pair dissociation with respect to the amount of water cluster present in the system. 2.4. ReaxFF MD Simulations. To validate our trained Na/Si/O/H force field, we carried out various MD simulations to calculate sodium silicates structures and sodium ion transport in different environments. The MD simulations were carried out using the ReaxFF implementation within the Amsterdam density functional package.43 Unless otherwise noted, the Newton’s equation of motion was solved with the Verlet algorithm at a time step of 0.25 fs. Also, the Berendsen thermostat44 with a damping constant of 100 fs is used in the NVT (constant volume and temperature) ensemble, and the NPT (constant pressure and temperature) ensemble employs the Berendsen barostat with a damping constant of 5 ps. The ReaxFF reactive force field parameters for Na/Si/O/H, which were developed and used throughout this work, are presented in Supporting Information. 2.4.1. Crystalline and Amorphous Structures of Na2Si2O5. We used a supercell of 1 × 2 × 2 Na2Si2O5 unit cells (144 atoms) to test the structural properties of crystals and amorphous sodium silicates. An NPT ensemble was applied, and the temperature was kept at 1 K to obtain the crystalline structure. For the construction of an amorphous structure, we initially melted the same supercell at 4000 K using the NVT ensemble. Subsequently, a cooling rate of 10 K/ps was applied while ensuring an internal pressure of 0 GPa, so that the simulation box can adjust during this melt-quench process. After the structure was cooled to room temperature at 300 K, we continued the simulation at this temperature for 0.5 ns with the NPT ensemble until we observed a steady density profile.

In this study, we performed and included QM calculations in the force field training set to direct the Na/Si/O/H force field development toward the description of the NaSiOx/water system. These calculations include volume−energy relations of different NaSiOx crystals, sodium ion migration energy barriers, and sodium ion interactions with various hydrated systems. The specific computation details of each data set are discussed in the following subsections, whereas the objective of such implementations is discussed and covered in Section 3. Using the QM-based training set, parameters that are related to the shielding terms for Na charge, Na−O bond and off diagonal, Si−Na off diagonal, and valence angle terms (O− Na−O, Na−O−Na, Na−O−H, and Si−O−Na) were optimized. For force field development that involves adding a new atom component (Na) into the existing force field parameter set (Si/O/H), it is important to optimize the new parameters while retaining the validity and description capability of the original set. In this regard, all Si/O/H bond and angle combinations, e.g., Si−O bond parameters, Si−O−H angle parameters, were maintained from the Pitman et al.31 values. Retaining these parameters ensures high compatibility between the new Na/Si/O/H ReaxFF reactive force field and the previously published work and also allows consistent description of silica bulk and surface properties. 2.3. Quantum Mechanical Calculations for the Na/Si/ O/H Force Field Training Set. 2.3.1. NaSiOx Crystalline Phases. The QM calculations of energy versus volume (equation of state) data of six different SiO2−Na2O binary oxide crystalline phases: Na2SiO3, Na2Si2O5-monoclinic/ orthorhombic, Na4SiO4, Na6Si2O7, and Na6Si8O19, were performed using VASP 4.6 package37 with the projector augmented wave pseudopotentials.38 Energies for 30−40 volumes around the experimental density of each crystal phase were calculated. At each volume, the cell shape and atom positions were allowed to relax with fixed volume. Plane wave kinetic energy cutoff value was 400 eV, and the general gradient approximation functional with Perdew−Burke− Ernzerhof (PBE) parameters39 were used for the exchange and correlation interactions. The Monkhorst−Pack method was used for K-points generation with a reciprocal space resolution of 0.3 Å−1. 2.3.2. Sodium Migration Energy Barrier within Crystalline Na2Si2O5. The energy barrier of a sodium ion’s diffusion to a neighboring vacancy site was calculated from the climbing image nudged elastic band (NEB) method40,41 in Vienna ab initio simulation package.37 Here, we used an orthorhombic Na2Si2O5 sodium silicate crystalline structure (space group Pna21). A 1 × 1 × 2 supercell (72 atoms) was generated and deliberately removed one sodium ion from the lattice to create a vacancy site. The migration energy for sodium vacancy diffusion was derived by the climbing image NEB, in which five structure images are used along each sodium migration pathway. Due to the symmetry of this crystal structure, we found that the sodium vacancy site can exchange positions with three neighboring sodium ions, meaning that there exist three different pathways for sodium transport. All of the energy barrier data were included in the Na/Si/O/H training set. 2.3.3. Si(OH)4 Monomer Interaction with Sodium Cation. The NaSiOx crystal QM data introduced above were used to optimize the Na/Si/O related parameters. To ensure optimization of other parameters, such as Na−O−H angle term, additional sets of data that describe the sodium ion’s interaction with a single Si(OH)4 monomer were included to 19615


Article


Figure 1. Equations of state for six different SiO2−Na2O binary oxide crystalline phases (NaSiO3, Na2Si2O5: monoclinic, Na2Si2O5: orthorhombic, Na4SiO4, Na6Si2O7, Na6Si8O19), using QM and ReaxFF methods.

Figure 2. Sodium migration energy barrier to nearby adjacent vacancy site. The green trajectory in (a)−(c) indicates that three different pathways are available for the sodium ion. Each point in the trajectory (from 0 to 6) denotes the migration steps for the energy profiles shown in (d)−(f). QM and ReaxFF results show that ReaxFF can properly describe which pathway would be favorable (Na (blue), Si (yellow), O (red)).

constant energy) equilibrations were performed afterward to collect the trajectories. For the interdiffusion study of the sodium cation and proton, we removed three sodium ions within the amorphous sodium silicate and replaced them with protons. After the system energy was minimized with the conjugate gradient scheme within the ReaxFF code, we equilibrated the system with the NVT ensemble at 800 K. Subsequently, a 100 ps trajectory was obtained for data analysis. 2.4.3. Dissociation Rate of Sodium Hydroxide (NaOH) in Aqueous Environment. An NVT simulation at 300 K with 10 sodium hydroxide (NaOH) and 782 water molecules was

Additional 0.5 ns long NVT simulations were performed at 300 K to obtain our final amorphous sodium silicate structure. 2.4.2. Sodium Cation Self-Diffusion and Interdiffusion of Na+/H+-Exchanged Amorphous Silicate. The crystal and amorphous Na2Si2O5 structures that we created (Section 2.4.1) were used to investigate the Na+ diffusivity at different temperatures. For the ionic self-diffusion of sodium cation, structures were first equilibrated at various temperatures, 300, 600, 700, 750, 800, 850, and 900 K, with the NPT ensemble for 200 ps. We found that 200 ps simulations were sufficient for the density to converge. NVE (100 ps) (constant volume/ 19616


Article


Figure 3. Sodium cation/hydrated silica monomer interaction. Sodium ion approaching (a) to the nearest Si−OH site (monodentate) and (b) to a site in between two Si−OH groups (bidentate site) of the silica monomer. The energy jumps are accompanied by the hydration shell changing as the sodium ion approached toward −OH groups of the silica monomer (Na (blue), Si (yellow), O (red), oxygen that is within the hydration shell of a sodium cation (green)).

tion of different sodium environments within the crystals, Na− O−Na and O−Na−O valence angle terms are parameterized to provide reasonable energy descriptions during the compression and expansion of crystals. Figure 1 provides the comparison of QM and ReaxFF results of the equations of state for six NaSiOx crystal structures. The calculations show that different binary oxide states can be distinguished with the ReaxFF method. Also, the energy curve shapes generated by the ReaxFF method are in reasonable agreement with the QM data. 3.1.2. Sodium Migration Energy Barrier. It is well known that alkali modifier ions, like Na+, in glass can easily be removed from the glass surface when glass is in contact with water.45,46 These highly mobile ions can also be leached out from the bulk glass, which involves a sequential hopping of a sodium ion via adjacent vacancy sites. The same mechanism is responsible for ionic diffusion within the bulk material, for example, the sodium disilicate glass is proposed as a promising Na+ conductor because of its high mobility in the solid.47 In this work, we chose the crystalline structure to investigate the activation energy of sodium migration since the structure of glass cannot be defined as one single state. Figure 2a−c illustrates possible sodium migration pathways to the nearest vacancy site in the crystalline Na2Si2O5. Figure 2d−f shows that the energy barriers calculated by QM are well reproduced by ReaxFF although some deviations are observed.

performed to obtain the NaOH dissociation rate. A cubic cell of 28.86 Å on each side was used, and all of the Na+ and OH− were randomly placed within the simulation box to ensure full dissociation state of NaOH before any contact with water molecules. For the 100 ps of MD simulation in aqueous environment, molecular analysis using a bond order cutoff of 0.3 was performed every 50 MD steps to give molecular composition of the system.

3. RESULTS AND DISCUSSION 3.1. Optimization of the Na/Si/O/H ReaxFF Force Field. 3.1.1. Equations of State for NaSiOx Crystals. As mentioned in the Section 2.3.1, we performed QM calculations on various SiO2−Na2O binary oxide crystalline phases. The volume−energy relationship (equations of state) of six crystal structures, that of NaSiO3, Na2Si2O5: monoclinic, Na2Si2O5: orthorhombic, Na4SiO4, Na6Si2O7, and Na6Si8O19, was included as references to define the ReaxFF parameters. Since these crystals have different coordination numbers, the equations of state data can lead to the parameterization of Na− O bond energies and their relationship to the bond order. In addition, the environment around sodium can be quite different for various sodium silicate systems, for example, sodium can be either near bridging oxygens (BO) or nonbridging oxygens (NBO). Although the volume−energy relationship cannot directly regulate the bond angle distribu19617


Article


Figure 4. Sodium hydroxide pair solvation by n water molecules (n = 1−6). The fitted Na/Si/O/H force field parameter set in this study gives overall good agreement of ReaxFF against the QM data.

theory level) study48 on the hydration and dissociation phenomena on sodium hydroxide (NaOH−n(H2O), n = 1− 6) reported that regardless of the state of dissociation for these small clusters, the sodium ion tended to coordinate with four H2O. Similarly, in our calculations, as the sodium ion approaches the silanol groups of Si(OH)4, the water molecules around the ion rearrange to an energetically favorable state, one or two being replaced by hydroxyl groups to maintain the optimum coordination number around the sodium ion. After the energy jump, the QM data in Figure 3 indicate that the energies gradually rise with further approach of a sodium ion to the Si(OH)4 monomer, where strong repulsion would exist between the cation and proton from the silanol group. As the sodium ion gets closer than 2.35 Å, which is close to the Na−NBO bond length in the amorphous sodium silicate structure,16 we might expect to observe hydroxyl dissociation and Na−O bond formation. However, ion exchange was not captured by either the QM or ReaxFF calculation. Since the sodium ions move into the complexes with a limited number of water molecules (mostly four to five), they may not get close enough to the Si−OH to cause proton hopping. Additional water in the system may assist our understanding in this matter, but we find that our QM data sets were sufficient to give adequate descriptions of sodium−water interactions. Altogether, the developed Na/Si/O/H ReaxFF force field can give the correct energetics related with the solvation shell change of sodium ions near the surface of silicate glasses. 3.1.4. Sodium Hydroxide Dissociation Behavior in Aqueous Environment. After the sodium ions are leached out from the amorphous sodium silicate systems, they begin to interact with water molecules. To correctly describe the sodium cations in aqueous environment, we investigated various hydrated structures of [NaOH−n(H2O)] (where n = 1−6) clusters and included their dissociation energies to the training set. Figure 4 illustrates that the fitted parameters in ReaxFF are in good agreement with the QM data. 3.2. Applications of Na/Si/O/H ReaxFF Force Field. 3.2.1. Structure Properties of the Crystalline and Amorphous Sodium Silicate Systems. To test the newly developed Na/Si/

3.1.3. Si(OH)4 Monomer Interaction with Sodium Cation. Since the Si/O/H parameters were fixed to their Pitman values31 in this work, the capability of ReaxFF to describe the formation of a hydrated silica-rich layer at the SiO2−water interface is retained in this new version of the force field. However, the existence of sodium components in the glass requires additional reactions to be defined at the surface Na +(glass) + H 2O → H+(glass) + NaOH

(3)

NaOHaq ↔ Na +aq + OH−aq

(4)

As noted in eq 3, sodium ions in multicomponent glasses can easily leach out through the ion-exchange process with protons from water molecules. In addition, these leached alkali ions can undergo subsequent interaction with water molecules, as shown in the eq 4. To account for sodium’s interaction with Si−OH groups that are present at the surfaces of multicomponent glasses, we optimized the ReaxFF reactive force field parameters using the QM data described in Section 2.3.3. Figure 3 shows the configurations of each complex used in this study and the energy comparison of QM and ReaxFF fitting results. In Figure 3a,b, there are notable energy jumps at dNa−O = 3.40 Å and dSi−Na = 3.28 Å. Our connectivity analysis based on the bond order cut off of 0.3 in ReaxFF implies that these energy changes are related to changes in the hydration shell around the sodium ion. For the monodentate case, the sodium ion interacts with six water molecules just before the energy jump as it approaches the Si(OH)4 monomer. Then, the energy drops dramatically as one H2O molecule is removed from the hydration shell and replaces with one on the monomer’s hydroxyl groups. For the bidentate complex, we find that the sodium ion coordination number with water changes from five at dSi−Na = 3.28 Å to four as the sodium ion approaches the chevron site, between two silanol groups. Although the nature of the hydroxyl in the silanol groups is different from that of the OH− anion, we may relate the drop of energy in Figure 3 to the coordination number for the sodium ion. A QM (density functional and Moller−Plesset 19618


Article


Figure 5. (a) Crystalline and (b) amorphous Na2Si2O5 structure used in our ReaxFF simulation. Amorphous structure was obtained by melting the crystal at 4000 K followed by the 10 K/ps cooling to 300 K (Na (blue), Si (yellow), O (red)).

Table 1. Lattice Constants and Bond Lengths (Si-Bridging Oxygen/Si-Nonbridging Oxygen) of Na2Si2O5 Crystalline Structurea property

ReaxFF this work

DFT PBE this work

DFT PBE47

experiment49

a/Å b/Å c/Å α = β = γ [°] Si−Obr distance/Å Si−Onbr distance/Å

9.404 5.560 8.480 90 1.639−1.649 1.591−1.602

9.513 5.622 8.489 90 1.658−1.675 1.591−1.595

9.539 5.656 8.468 90 1.658−1.674 1.591−1.595

9.441 5.580 8.356 90 1.633−1.656 1.571−1.580

a

ReaxFF results (at 1 K) are compared with the ab initio and experimental data from literature.

O/H force field’s applicability to the sodium silicate, we carried out structural studies on both crystalline and amorphous sodium silicate models. Among the six crystal configurations from our study, the Na2Si2O5 orthorhombic structure was chosen for comparison with the ab initio and experimental values that were reported in the literature.47,49 The investigated structural properties include lattice parameters and bond lengths for the crystal system and density and partial radial pair distribution functions (RDF) for the amorphous sodium silicate structure. Figure 5 compares the Na2Si2O5 crystal structure and the amorphous structure created using ReaxFF MD simulations. Table 1 lists the calculated cell axes and Si−O bond length from the ReaxFF method and shows that ReaxFF provides good agreement with both ab initio study (in this work and from reference) and experiment.47,49 For the amorphous sodium silicate, we first calculated the density of our system and list it in Table 2 along with a reported experimental

nonreactive MD simulations and (2) the charge calculations through the EEM method required at every iteration makes it computationally challenging for ReaxFF to adopt a cooling rate of 1 K/ps or below, which is typically employed by classical MD simulations. We find that current usage of the ReaxFF method is typically limited to NPT equilibration at the room temperature, after the glass structure is obtained from the melt and quench processes through classical MD simulations.33,34 Here, despite the differences that could be generated by different thermal histories, we have demonstrated the whole melt-quench process for the sodium silicate system using the ReaxFF method. We propose that such a process is feasible with ReaxFF and reproduces a final density that is in good agreement with that of the real glass. Next, we gathered the structural information on this amorphous sodium silicate glass by calculating the radial pair distribution function (RDF). Figure 6 shows the partial RDF of the amorphous sodium silicate structure from the ReaxFF MD simulation. We leave the discussion regarding the contribution of soda content to the peak position of partial RDF in sodium silicate glasses to other refs 33, 52, 53, as there exist various experimental data and theories in the literature. For simple comparison, we can observe the elongation of the Si−O bond length to 1.65 Å from the peak position of the partial Si−O RDF, which is in good agreement with experimental determined Si−O bond lengths.50 For the partial Na−O RDF, we observe two sharp peaks at 2.05 and 2.55 Å rather than a broad peak at around 2.30−2.35 Å. From the trajectory of our simulation, we found that the Na−NBO pairs are responsible for the first sharp peak, indicating that ReaxFF recognizes two different sodium− oxygen environments, NBO and BO sites, but slightly overestimates the Coulomb interactions between Na and NBO in the amorphous sodium silicates. Although ReaxFF has the capability to transfer charges at every time step by its inherent charge redistribution scheme, we believe that this transferability between one sodium and surrounding oxygens,

Table 2. Density of the Amorphous Sodium Silicate (xNa2O·(1 − x)SiO2) property 3

density [g/cm ] a

ReaxFFa

classicalb

experimenta,49

2.463

2.406

2.493 b

For Na2Si2O5; soda content x = 0.33. For x = 0.30 glass (corresponding experimental density; 2.466 g/cm3) refs 33, 51.

density and a density calculated with classical MD simulations with the nonreactive Teter potential.33,49,50 We find that our ReaxFF simulations compare well to the reported experimental density. Note that classical MD simulations of amorphous materials are usually studied using larger simulation cells that have an order of magnitude more atoms than in our study (144 atoms). The cooling procedure used to create the amorphous structure for ReaxFF can be quite challenging because (1) ReaxFF requires up to 10 times higher time step resolution than 19619


Article

The Journal of Physical Chemistry C MSD = ⟨δr 2⟩ =

D = lim

t →∞

⟨δr 2⟩ 6t

1 N

∑ [ri(t + t0) − ri(t0)]2 i

(5)

(6)

where N is the total number of atoms, ri(t) is the position of the atom. We plotted the MSD of each atom component of the sodium silicate in Figure 7. From the 100 ps long simulation, we considered only the 0 to 50 ps results for the diffusivity calculation to ensure statistical relevance. For comparison, MSD plots of crystalline Na2Si2O5 at two different temperatures 300 and 800 K and the same for amorphous sodium silicate are presented in Figure 7a−d. For the crystalline structure, the MSDs are time-independent plateaus. On the other hand, the case of amorphous sodium silicate in Figure 7c,d shows linear increase of the Na+ displacements even at the room temperature. From our QM calculation, we identified that three different pathways exist for the sodium ion’s transport. We can expect that the sodium ion’s migration toward the corner-shared SiO4 tetrahedra, which are the cases shown in Figure 2a,b, would be less likely to happen because of the high-energy barrier. The most favorable pathway for sodium migration would be the pathway in Figure 2c, where the trajectory of the sodium ion seems to be making channels along the tetrahedral SiO4 layer. Similarly, a disordered structure of the amorphous sodium silicate glasses would open the channels and offset the high-energy barriers, thereby allowing sodium ionic diffusion to be more prominent within the amorphous sodium silicate glasses. In Figure 8, we plotted the ln[DNa] with respect to a reciprocal of the temperature to compare our sodium migration behavior to the ab initio MD (AIMD) simulation.47 From the Arrhenius equation, the activation energy can directly be derived from the slope of linearly extrapolated data. Our ReaxFF MD simulation results predicted the activation energy for the sodium migration in the amorphous Na2Si2O5 structure to be 6.67 kcal/mol. Given that the lowest energy barrier for the sodium ion to migrate within the crystalline structure was 7.87 kcal/mol (Figure 2f), structural change to amorphous state seems to provide energetically favorable pathway for sodium. Overall, the ReaxFF reproduces the activation energy of sodium diffusion at a high accuracy when compared to the 6.92 kcal/mol from AIMD. We further investigated the transport properties by replacing sodium ions with protons within the amorphous sodium silicate structure (MD simulation details in 2.4.2). Figure 9 shows the MD simulation results with the ReaxFF method for the interdiffusion of a sodium ion and a proton. Figure 9a is the trajectories of a single sodium cation and single proton within the Na2Si2O5 structure at 800 K. For the sodium ions, they are separated into modes associated with the vibration at two fixed sites. This distinct observation provides an evidence that sodium migration occurs by a hopping mechanism between two adjacent sites. In case of proton, we see that the trajectory shows a prominent increase in the number of the sites where the proton resides during the simulation. As previously discussed, sodium can move to adjacent site when there is energetically available pathway, an open channel, within the structure. Protons, however, would be less restricted to hop around nearby NBO sites because they exhibit smaller size than the sodium ions and require a closer connectivity to

Figure 6. Partial radial pair distribution function for the amorphous sodium silicate obtained from the ReaxFF simulation.

e.g., NBO to Na then Na to BO, is rather confined to the Na− NBO interaction at the temperature range that we have studied here. At high temperatures above 900 K, we find that the valley between two sharp peaks rises and forms a broad peak. In addition, the double peak could be an artifact of the relatively small system size (144 atoms in this work). Our ongoing study on the effect of different sizes (an order of magnitude larger number of atoms) and simulation procedure suggests that removal of the peak near 2.0 Å is feasible to a certain extent. We would like to ascribe this issue to the training set that we used for the development of the Na/Si/O/H force field. In the training set, we gave no information about the sodium silicates other than the crystal structures’ equations of state data. We gave higher priority, by means of the weight terms included in the error function eq 1, to the sodium cation−water interaction during the force field optimization process. During force field training, we had a stage where we could observe a resolution of this disrupted partial Na−O RDF. However, we excluded that particular set as our final fitted force field, since we observed a deterioration with the sodium−water description, which is critical for the study of the glass−water interface. Sodium− oxygen interactions then may be more adapted to the Na+− OH−[(H2O)] description, resulting in stronger interactions with the Na−NBO sites. 3.2.2. Transport Properties: Ionic Diffusion of Na+ within the Crystalline and Amorphous Sodium Silicates. We carried out MD simulations with the ReaxFF method to describe the transport behavior of sodium ions for the validation of Na/Si/ O/H force field. Such transport phenomena can be measured by means of diffusion coefficient. The diffusion coefficient D of the ionic species in the material can be determined from the time evolution of mean square displacement (MSD), where the diffusion coefficient D can then be derived from the slope of the MSD plot. 19620


Article


Figure 7. Mean square displacement (MSD) plots of crystalline Na2Si2O5 at two different temperatures (a) 300 K and (b) 800 K and of amorphous sodium silicate at (c) 300 K and (d) 800 K. We observe the linear increase of the sodium ions in the amorphous sodium silicate, which indicates that the long-range disordering of Si−O network opens up the pathway for the sodium ions to be mobile within the amorphous structure.

through the MD simulations, by placing 10 sodium hydroxides (NaOH) into the bulk water configurations. Figure 10 illustrates the time evolution of sodium ion species within the simulation box. Number of species in Figure 10 indicates how many of sodium hydroxide pairs are dissociated during the MD simulation. From the fluctuation that is observed in this figure, we may conclude that there are constant proton transfers via water dissociation and NaOH formation as well as dissociation. From time-averaged value, our force field predicts ∼67% dissociation rate within this configuration. Literature reports a valency and ionic radius correlation to pKa value for the metal hydroxide.54 Direct comparison with this correlation to our results would be challenging, but 71−88% dissociation rate derived from this correlation for the NaOH solution agrees well with our predictions. Overall, we find that the developed Na/Si/O/H force field shows good description capability of Na−OH related bond formation and dissociation in the aqueous sodium hydroxide solution.

Figure 8. Arrhenius plot of sodium ion diffusivity in the amorphous Na2Si2O5 structure. The ReaxFF MD simulation (red triangle up solid) predicts 6.67 kcal/mol for the energy barrier for the sodium migration, whereas the ab initio MD study47 (black circle) reported value is 6.92 kcal/mol for the temperature of 773−973 K and 3.92 kcal/mol within 573−723 K range.

4. CONCLUSIONS Reactive force field (ReaxFF) parameters were developed for Na/Si/O/H interactions, which enables the reactive molecular dynamics simulations on sodium silicate−water systems. We first created a training set that consists of energetic and structure data from various quantum mechanical (QM) calculations that would give adequate descriptions of our target systems and then used these data sets as our reference points for the force field fitting process. Included QM data are

oxygen atoms. Faster diffusion of proton can therefore be attributed to the size difference between the two types of ions. 3.3. Sodium Hydroxide Dissociation Behavior in Aqueous Environment. The force field description regarding sodium ions in water environment was also validated 19621


Article


Figure 9. (a) Sodium ion and proton transport trajectories and (b) MSD plot of sodium and proton in the amorphous Na2Si2O5 at 800 K.

These results are relevant for modeling the surface reactivity of alkali-containing silicate glasses and therefore indicate that the ReaxFF method is a computationally efficient tool for evaluating detailed chemical dissolution mechanisms. Such mechanisms include the interdiffusion process of sodium ions from glasses and protons from water, subsequent ionic selfdiffusion of sodium ions from subsurface region to vacancy sites at the glass−water interface, and sodium ions interaction with water after they are leached out from the amorphous sodium silicate system.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.8b05852. ReaxFF reactive force field parameters for Na/Si/O/H (PDF)

Figure 10. Number of dissociated sodium ions during 100 ps of NVT simulation with 10 sodium hydroxide (NaOH) and 782 water molecules (H2O) in a cubic periodic box of dimension 28.86 Å. The ReaxFF predicts ∼67% sodium dissociation rate in this pH condition.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: +1-814-863-6277.

(a) energy versus volume descriptions of different NaSiOx crystalline phases, (b) energy barriers of a sodium cation’s transport within the sodium silicate crystal structure, (c) interactions between hydroxylated silica surface and Na+/water systems, and (d) [NaOH−n(H2O)] (n = 1−6) cluster structures and their dissociation energies. After we have obtained the optimized force field parameter set, we validated our force field and applied it to form the amorphous sodium silicate structure and probe the diffusion behavior of sodium ions and protons within the ReaxFF MD simulations. Both structural properties and transport properties indicated that our force field gives good agreement with experimental values and ab initio MD work from the literature. The main feature of the transport phenomena within the amorphous sodium silicate is that hopping mechanisms take place for both sodium ions and protons. Although the decrease in diffusion activation energy induces sodium ions to move easily to their adjacent vacancy sites, protons show a faster diffusion due to size difference. This eventually allows them to penetrate easily to the channels formed by long-range disordering of the Si−O network in the amorphous structures. The force field was also applied to validate the dissociation behavior of sodium hydroxides within the bulk water, and we found that the reactive MD simulation could appropriately address the proton-transfer process in aqueous sodium hydroxide solutions.

ORCID

Seung Ho Hahn: 0000-0002-1722-5631 Jessica Rimsza: 0000-0003-0492-852X Louise Criscenti: 0000-0002-5212-7201 Jincheng Du: 0000-0003-4805-7498 Tao Liang: 0000-0002-5867-9939 Adri C. T. van Duin: 0000-0002-3478-4945 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by NSF DMR grant #1609107. Computations for this research (Sections 2.3.2,, and 2.4) were performed on the Pennsylvania State University’s Institute for Cyber Science Advanced Cyber Infrastructure (ICS-ACI). This work was also supported by the Laboratory Directed Research and Development (LDRD) program of Sandia National Laboratories. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. J.D. acknowledges support of NSF DMR grant #1508001. 19622


Article


■

(22) Yu, Y.; Wang, B.; Wang, M.; Sant, G.; Bauchy, M. Revisiting Silica with ReaxFF: Towards Improved Predictions of Glass Structure and Properties via Reactive Molecular Dynamics. J. Non-Cryst. Solids 2016, 443, 148−154. (23) Krishnan, N. M. A.; Wang, B.; Le Pape, Y.; Sant, G.; Bauchy, M. Irradiation-Driven Amorphous-to-Glassy Transition in Quartz: The Crucial Role of the Medium-Range Order in Crystallization. Phys. Rev. Mater. 2017, 1, No. 053405. (24) Yu, Y.; Krishnan, N. M. A.; Smedskjaer, M. M.; Sant, G.; Bauchy, M. The Hydrophilic-to-Hydrophobic Transition in Glassy Silica Is Driven by the Atomic Topology of Its Surface. J. Chem. Phys. 2018, 148, No. 074503. (25) 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. (26) 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, No. 174704. (27) 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. (28) Macià Escatllar, A.; Ugliengo, P.; Bromley, S. T. Modeling Hydroxylated Nanosilica: Testing the Performance of ReaxFF and FFSiOH Force Fields. J. Chem. Phys. 2017, 146, No. 224704. (29) Mahadevan, T. S.; Du, J. Evaluating Water Reactivity at Silica Surfaces Using Reactive Potentials. J. Phys. Chem. C 2018, 122, 9875− 9885. (30) 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. (31) Pitman, M. C.; van Duin, A. C. T. Dynamics of Confined Reactive Water in Smectite Clay−Zeolite Composites. J. Am. Chem. Soc. 2012, 134, 3042−3053. (32) Rimsza, J. M.; Deng, L.; Du, J. Molecular Dynamics Simulations of Nanoporous Organosilicate Glasses Using Reactive Force Field (ReaxFF). J. Non-Cryst. Solids 2016, 431, 103−111. (33) 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. (34) 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. (35) 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. (36) van Duin, A. C. T.; Baas, J. M.; van de Graaf, B. Delft Molecular Mechanics: A New Approach to Hydrocarbon Force Fields. Inclusion of a Geometry-Dependent Charge Calculation. J. Chem. Soc., Faraday Trans. 1994, 90, 2881−2895. (37) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B 1993, 47, 558−561. (38) Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953−17979. (39) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (40) Henkelman, G.; Uberuaga, B. P.; Jónsson, H. A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths. J. Chem. Phys. 2000, 113, 9901−9904. (41) Henkelman, G.; Jónsson, H. Improved Tangent Estimate in the Nudged Elastic Band Method for Finding Minimum Energy Paths and Saddle Points. J. Chem. Phys. 2000, 113, 9978−9985. (42) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C. et al. Gaussian 03, revision B.05; Gaussian, Inc.: Pittsburgh, PA, 2003.

REFERENCES

(1) Orowan, E. Fracture and Strength of Solids. Rep. Prog. Phys. 1949, 12, 185−232. (2) Wiederhorn, S. M. Influence of Water Vapor on Crack Propagation in Soda-Lime Glass. J. Am. Ceram. Soc. 1967, 50, 407− 414. (3) Surdyka, N. D.; Pantano, C. G.; Kim, S. H. Environmental Effects on Initiation and Propagation of Surface Defects on Silicate Glasses: Scratch and Fracture Toughness Study. Appl. Phys. A 2014, 116, 519−528. (4) 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 2017, 121, 11534−11543. (5) Rimsza, J. M.; Du, J. Nanoporous Silica Gel Structures and Evolution from Reactive Force Field-Based Molecular Dynamics Simulations. npj Mater. Degrad. 2018, 2, No. 18. (6) 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. (7) Vienna, J. D.; Ryan, J. V.; Gin, S.; Inagaki, Y. Current Understanding and Remaining Challenges in Modeling Long-Term Degradation of Borosilicate Nuclear Waste Glasses. Int. J. Appl. Glass Sci. 2013, 4, 283−294. (8) 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, No. 6360. (9) Davis, K. M.; Tomozawa, M. An Infrared Spectroscopic Study of Water-Related Species in Silica Glasses. J. Non-Cryst. Solids 1996, 201, 177−198. (10) Amma, S.; 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. (11) Tsuneyuki, S.; Tsukada, M.; Aoki, H.; Matsui, Y. FirstPrinciples Interatomic Potential of Silica Applied to Molecular Dynamics. Phys. Rev. Lett. 1988, 61, 869−872. (12) Xiao, Y.; Lasaga, A. C. Ab Initio Quantum Mechanical Studies of the Kinetics and Mechanisms of Silicate Dissolution: H+(H3O+) Catalysis. Geochim. Cosmochim. Acta 1994, 58, 5379−5400. (13) 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. 2003, 119, 6418−6422. (14) Tielens, F.; Gervais, C.; Lambert, J. F.; Mauri, F.; Costa, D. Ab Initio Study of the Hydroxylated Surface of Amorphous Silica: A Representative Model. Chem. Mater. 2008, 20, 3336−3344. (15) Huang, C.; Cormack, A. N. The Structure of Sodium Silicate Glass. J. Chem. Phys. 1990, 93, 8180−8186. (16) Yuan, X.; Cormack, A. N. Local Structures of MD-Modeled Vitreous Silica and Sodium Silicate Glasses. J. Non-Cryst. Solids 2001, 283, 69−87. (17) 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. (18) 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. (19) 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. (20) 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, No. 15011. (21) Rimsza, J. M.; Jones, R. E.; Criscenti, L. J. Surface Structure and Stability of Partially Hydroxylated Silica Surfaces. Langmuir 2017, 33, 3882−3891. 19623


Article

The Journal of Physical Chemistry C (43) ADF, SCM Theoretical Chemistry; Vrije Universiteit: Amsterdam, The Netherlands, 2018. http://www.scm.com. (44) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (45) Douglas, R. W.; El-Shamy, T. M. M. Reactions of Glasses with Aqueous Solutions. J. Am. Ceram. Soc. 1967, 50, 1−8. (46) Lopes, M.; Shelby, J. E. Introduction to Glass Science and Technology; The Royal Society of Chemistry: Cambridge, U.K., 2005. (47) 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. (48) Kumar, A.; Park, M.; Huh, J. Y.; Lee, H. M.; Kim, K. S. Hydration Phenomena of Sodium and Potassium Hydroxides by Water Molecules. J. Phys. Chem. A 2006, 110, 12484−12493. (49) Fleet, M. E.; Henderson, G. S. Epsilon Sodium Disilicate: A High-Pressure Layer Structure [Na2Si2O5]. J. Solid State Chem. 1995, 119, 400−404. (50) Du, J.; Cormack, A. N. The Medium Range Structure of Sodium Silicate Glasses: A Molecular Dynamics Simulation. J. NonCryst. Solids 2004, 349, 66−79. (51) Tilocca, A.; de Leeuw, N. H.; Cormack, A. N. Shell-Model Molecular Dynamics Calculations of Modified Silicate Glasses. Phys. Rev. B 2006, 73, No. 104209. (52) Wright, A. C.; Clare, A. G.; Bachra, B.; Sinclair, R. N.; Hannon, A. C.; Vessal, B. Neutron Diffraction Studies of Silicate Glasses. Trans. Am. Crystallogr. Assoc. 1991, 27, 239−254. (53) Zotov, N.; Keppler, H. The Structure of Sodium Tetrasilicate Glass from Neutron Diffraction, Reverse Monte Carlo Simulations and Raman Spectroscopy. Phys. Chem. Miner. 1998, 25, 259−267. (54) Davies, C. W. 280. The Electrolytic Dissociation of Metal Hydroxides. J. Chem. Soc. 1951, 1256−1258.

19624


Development of a ReaxFF Reactive Force Field for ... - ACS Publications

Recommend Documents