Development of Water Reactive Potentials for Sodium Silicate Glasses

Subscriber access provided by UNIV OF LOUISIANA

B: Fluid Interfaces, Colloids, Polymers, Soft Matter, Surfactants, and Glassy Materials

Development of Water Reactive Potentials for Sodium Silicate Glasses Thiruvilla S. Mahadevan, Wei Sun, and Jincheng Du J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.9b02216 • Publication Date (Web): 29 Apr 2019 Downloaded from http://pubs.acs.org on May 5, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Development of Water Reactive Potentials for Sodium Silicate Glasses Thiruvilla S Mahadevan, Wei Sun and Jincheng Du* Department of Materials Science and Engineering, University of North Texas, Denton, Texas (*Corresponding author email: [email protected]) Abstract: Molecular Dynamics (MD) simulations provide important insights into atomistic phenomena and are complement to experimental methods of studying glass-water interaction and glass corrosion. For simulations of glass-water systems using MD, there is a need to for a reactive potential that is capable not only to describe the bulk and surface glass structures but also reactions between glass and water. An important aspect of the glass water interaction is the dissociation of water and its interaction with glass components that can result in the dissolution and alteration in the structure of glass. These phenomena can be efficiently simulated using “Reactive” potentials that allow for dissociation of water while properly describe the bulk physical properties of water. We demonstrate a method to develop parameters for simulations of sodium silicate glasses and their interactions with bulk water. The developed parameter set were used to simulate sodium silicate glasses of different compositions and the local structure of the simulated glass is in good compliance with experimentally obtained structural information. We also demonstrate that the parameter set predicts an accurate value for the hydration number and dissociation reactions of NaOH in water. Based on these results, we posit that these simple and computationally efficient reactive potentials can be used to for further studies of water induced structural modifications in sodium silicate glasses.

1. Introduction Silicate glasses are amongst the most widely used industrial materials1,2 and research into the effect of their atomic structure and composition on their physical properties is of considerable importance. Degradation of properties of silicate glasses in the presence of water from the environment is also an important aspect with applications in geological and nuclear materials3,4. Molecular Dynamics (MD) simulations are a useful tool5 that is widely applied for the study of atomistic phenomena in glasses and can be complementary to experimental studies in the understanding of glass corrosion at the nanoscale6. Sodium silicates are an architype of wide range of practical glass compositions and serve as an excellent model system for understanding the atomic structure and property relations, as well as mechanistic understanding of the interactions of glass with water from the enviornment7–9. MD simulations studies of glass-water interactions has seen important progress in the recent decade especially with the development of “Reactive” potentials that have accounted for the dissociation of water and the formation of silanols10–12. A general rule in the simulation of reactive potentials is that a parameterization of a given type of material characterized by the local environment and bond order cannot be directly applied as the parameter set for the same element in a different environment. This has been demonstrated with 1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the ReaxFF, a bond order based reactive potential that takes into account of charge transfer hence capable of modeling of chemical reactions, where some of the parameters for Hydrogen and Oxygen as applicable to hydro-carbons are not directly suitable to be used in the simulation of bulk water13,14, as a result of which additional fitting is needed to improve the parameterizations15. Thus, while dissociative water interacting with silicates has been studied from nearly three decades back16, it is only in the recent past that sodium has been parameterized to be studied with the ReaxFF simulation of silicate glass17,18. While the ReaxFF is framework that provides capability of simulation of reactions and heterogeneous interfaces of materials including their interactions with the environment, the bond order based potential has a complex potential from with relatively large number of parameterization required for the addition of every new element. Another analytically simpler potential was recently developed for the simulation of silica glass-water interactions19, with the results that compares well with the ReaxFF potentials. The MahadevanGarofalini Force Field or MGFF was based on the concept of diffuse change, i.e. varying charge of an ion as a function of the distance, as had been observed through quantum mechanical calculations20. The framework was initially developed for description of liquid water and was later expanded to study water-silica interaction by introducing three-body terms to describe covalent nature of the silicate network. Due to the simpler analytical form of the two body interactions, this force field can be easily applied in a tabulated for in the LAMMPS5 code. These features together with Wolf summation21 makes MGFF about ten times faster than ReaxFF simulation of the same system size hence computationally advantageous to study larger systems and longer times. In an earlier study of salt induced damage to cement structures22, Na+ and Cl- ions in aqueous environment had been parameterized with the MGFF. However, in that study the charges on Na+ and Cl- were set at +1e and -1e. A charge of +1e, if used in sodium silicate glass would require a charge of -1e on the oxygen to maintain charge neutrality. MGFF had already parameterized the charges on Si, O and H to be 1.808e, -.904e and 0.452e respectively and this necessitates the charge on the Na+ to be 0.452e. Thus, in this paper, we demonstrate the method for parameterization of additional glass components starting with Na+, while maintaining a charge of 0.452e on the Na+ ion and allowing the appropriate dissociation and reaction of NaOH in water. The goal is to develop a set of self-consistent potential within the MGFF framework to enable simulations of water interactions with sodium silicate glasses, which can be further expanded to multicomponent glasses. One of the primary aims of the work was to demonstrate parameterization of sodium to be used with the MGFF framework for silica-water. While there are several established potentials based on the Buckingham form23,24 that can be used for simulation of sodium silicates, the incorporation of a polarizable water model that can react appropriately with NaOH was emphasized in this work. The reactivity of Na+ ions in bulk water was explicitly established by considering NaOH dissociation in water. 2. Computational procedure and parameter development: Molecular dynamics simulations were carried out using the Mahadevan-Garofalini force field (MGFF) 25 which has been previously used for studying water-silica interactions with a

2 ACS Paragon Plus Environment

Page 2 of 25

Page 3 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


dissociative water model that allowed for reactions12,19. In brief, the analytical form of the potential function is as follows: 𝐸2𝐵𝑜𝑑𝑦 = 𝐸𝑞𝑞 + 𝐸𝑞𝑑𝑞𝑑 + 𝐸𝑞𝑞𝑑 + 𝐸𝑞𝑑𝑞 + 𝐸𝑟𝑒𝑝 + 𝐸𝑑𝑖𝑠𝑝 (1) where 𝐸𝑞𝑞 =

𝑞𝑖𝑞𝑗

𝑟𝑖𝑗

()

𝑟𝑖𝑗 erfc

(2)

𝛽

( )erfc( )

𝑞𝑖𝑞𝑗

𝐸𝑞𝑞𝑑 = 𝐸𝑞𝑑𝑞 = 4𝑟𝑖𝑗erf

𝑟𝑖𝑗

𝑟𝑖𝑗

2𝜉𝑖𝑗

𝛽

( )erfc( )

𝑞𝑖𝑞𝑗

𝐸𝑞𝑑𝑞𝑑 = 16𝑟𝑖𝑗erfc

𝑟𝑖𝑗

𝑟𝑖𝑗

2𝜉𝑖𝑗

𝛽

[with 𝑧

erfc(𝑧𝑖𝑗)

𝐸𝑟𝑒𝑝 = 𝐴𝑖𝑗𝑟𝑒𝑝 𝐸𝑑𝑖𝑠𝑝 =

(3)

𝑧𝑖𝑗

(4) 𝑟𝑖𝑗

𝑖𝑗

= 2𝜉𝑖𝑗 𝑟

]

(5)

― 𝐶𝑖𝑗6

(6)

𝑟6𝑖𝑗

The energy terms based on charges Eq’s and Eqd’s were designed so as to make the charge on atomic species vary as a function of the distance from the central atom. In contrast to the dynamic charge equalization method (QEq) used in ReaxFF17, the E(qd) terms make a “diffuse”, albeit static representation of charges, which allows for faster calculations. Such diffuse charge terms were originally designed for molecular, non-dissociable water that simulates the properties of bulk water over a wide range of temperatures and pressures20. According to the original work, having the diffuse terms for the charge allowed for the charge on an atom to vary smoothly between q and and q-qd as a function of the distance from it. Erep and Edisp are the repulsion and dispersion terms that determine the short range interactions. The long range interactions due to the Coulombic terms are handled by the Wolf summation method and hence the attenuation term of erfc(rij/) is factored into the charge based terms21. The functional form of the potential can be parameterized to adapt it to simulations of ionic and some covalent atomistic simulations. In the original application, the dissociative form of this potential was parameterized only for water and silica12. However, the framework of the potential allows additional elements for the potential to be used for multicomponent glasses. Besides the pair potential shown above, the MGFF also contains a 3-body potential to describe partial covalent bonding in the system which is given by the following analytical form: 𝐸3𝐵𝑜𝑑𝑦(𝑟𝑖𝑗,𝑟𝑖𝑘,𝜃𝑗𝑖𝑘) = 𝜆𝑗𝑖𝑘exp

(

𝛾𝑖𝑗 𝑟𝑖𝑗 ― 𝑟𝑖𝑗

0

+𝑟

𝛾𝑖𝑘

𝑖𝑘

―

)(cos 𝜃

𝑟𝑖𝑘0

2

𝑗𝑖𝑘

― cos 𝜃0𝑗𝑖𝑘)

(7)

where  and  are the parameters to be fitted, r0’s are the cut-off distances beyond which the 3body potential goes to zero and 0 is the target angle for the configuration of atoms j-i-k. The three body term is used mainly to regulate the angular configurations in non-linear molecules and acts to impose an energy penalty when the configuration strays beyond the target angle. The values of the threebody parameters that were ultimately derived from this work are given in table 2. 3 ACS Paragon Plus Environment


In the present work, we attempt to parameterize Na+, such that it can be used in simulation of sodium silicate glasses, sodium ion in solution, as well as sodium silicate glass – water interactions. The focus of this paper is to evaluate the sodium silicate glass structure of the potential and sodium dissolution in aqueous solutions. For the parameterization of Na+, we categorized the list of parameters into two categories: (a) Parameters that determine the charge based, long range, interactions - reduced charges q and qd, the charge diffusivity term ij that determines the rate at which the charge changes as a function of the separation distance, the Wolf damping term  (b) Parameters that determine the short range repulsion-dispersion terms Arep, r and C6. As mentioned earlier, in order to maintain charge neutrality in sodium silicate glass systems, the charge on Na+ was fixed as the same as that of H+. The ionic nature of Na+ was captured by reducing the value of Na-Na to a smaller value of 4Å compared to the value of 24Å for the other elements, which allowed for the charge to reduce faster as a function of distance from the Na+ ion. Values of Na-X, where X is one of the other elements is obtained as the geometric mean of Na-Na and Na-X. The wolf damping term is common to all pairs of elements and does not change from the original value of 4.46Å. The short range parameters are critical in determining the bond distances and energies. Force matching 26,27 is a useful technique in determining or refining these parameters, this is only possible when there is significant amount of data from ab initio molecular dynamics (AIMD) simulations. In the present work, we used empirical fitting to develop the sodium related parameters. Here, the structural data from previous simulations by Du and Cormack using the partial charge Teter potential with slight modification, which has been shown to give accurate description of the sodium silicate glass structure and properties in wide composition ranges28, and equation of state data from DFT calculations of several sodium silicate crystal phases, were used to fit the parameters to obtain the correct structures as a function of glass composition. Since our emphasis on studying sodium silicate glasses is to understand the structure, dynamics and reactivity with water, reliable parameters sets for all the interactions were obtained in a reasonable time by fitting the Na-O interactions first, followed by the Na-Na, Na-Si and Na-H interactions. Figure 1 shows a plot of the Na-O and Na-Na potential as a function of separation distance with the obtained parameter set - the O-Si and O-H interaction are also shown for comparison. rNa-O and Arep Na-O were parameterized so as to obtain the minimum for Na-O at around 2.3Å, which in turn would ensure Na-O distances in simulated glass to be in the same range. In the previous work, Si-O distance had been fixed at around 1.5Å such that the combined effect of multiple atoms surrounding the central pair of Si-O resulted an eventual Si-O distance of 1.615Å. The C6 parameter of Na-O interaction was adjusted to obtaine the correct density of silicate glass as compared to values obtained with Teter/Du potentials. Fitting the two body parameters for the 28% Na2O-72%SiO2 resulted in a good compliance of the structure of the melt quenched glass to that generated by melt quench with the Teter-Du potential. Once a match for bulk sodium silicate was obtained as compared with the Teter-Du 4 ACS Paragon Plus Environment

Page 4 of 25

Page 5 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


potential, the Na-H pair potentials was parameterized. This was done to fit to the ab initio hydration number of ~5 water molecules around Na in bulk water29. This match was obtained by making slight changes to the Na-O parameters, as well as the 3-body term to maintain the Na-O-H angle at around 120˚. All the parameters, including the newly developed sodium related parameters, for Si, O, Na and H interactions are given in Table 1. All MD simulations were performed using the LAMMPS code5 and DFT calculations were carried out using the VASP package 30,31. The system sizes for MD simulations were always chosen to be larger than 20 Å, which is twice the cut off for the interactions. Further details of particular simulations are given along with the results in the next section 3. Results and discussion 3.1 Structures of sodium silicate glasses: sodium environments The validation of the structure is shown in the comparison of the Na-O and Na-Na pair distribution functions simulated using MGFF and modified Teter potential for the 28% Na2O system shown in figure 2. A cubic box containing 2424 atoms with ~28%Na2O in SiO2 was melted at 6298K and quenched to 298K at ~4K/ps and then equilibrated under NPT conditions using the both potentials to obtain the final structure. In general, all simulations were done with a timestep of 1 femtosecond and a box size which was larger than twice the interaction cut-off. The final NPT structure allowed for the system to evolve to the correct equilibrated density. As mentioned in section 2 on parameterization, trials with several parameter sets were used to run these simulations. The Na-O related parameters were adjusted by comparing the Na-O pair distribution function and the location of Na-O minima in the potential curve similar to Figure 1 to judiciously choose the next trial value of parameters. The resulting final set of parameters yielded a near perfect match for the Na-O pair distribution function and a sufficiently accurate Na-Na pdf as shown in Figure 2. Once the above reasonable structure was established for that particular composition of sodium silicate, trial runs with the fitted potentials were carried out on 6 sodium silicate compositions xNa2O.(100-x) SiO2 with x varying from 0-50%. These systems were melt-quenched starting with a random Na2O – SiO2 network and finishing with an NPT at 298 K to allow the system to evolve to the correct density, followed by an NVT run at 298K to collect trajectories to calculate the glass structures and properties. A plot of the densities of the system as a function of the Na2O content is shown in Figure 3 as compared to the experimental density32 from literature. As can be seen, the densities of the simulated glasses are fairly accurate at low Na2O compositions but are slightly higher than the experimental value for compositions with higher than 40 mol% sodium content. For example, at composition with 50 mol% Na2O the simulated density is higher than experiment by about 6%, which is well within the range obtained through simulations with other potentials18,33. The structures of the sodium silicate glass at the different compositions were analyzed by studying the Na-O and Na-Na pair distribution functions g(r)’s as has been done with other potentials28,34 – the results are shown in figure 4a and b. The g(r)’s were averaged over the trajectories through the final 200ps of during the NVT runs after the melt-quench cycle and 5 ACS Paragon Plus Environment


subsequent equilibration under NPT at 298K. The gNa-O(r) shows a narrowing and increase of the first peak as well as a slight shift to lower distances from 2.37A to 2.35A with increasing Na2O content. This has been observed experimentally35 as well as through simulations34. The oxygen coordination shell around Na also increases with increasing Na2O content from 4.6 to 5. The gNa-Na(r) shown in Figure 4b also shows behavior consistent with other MD simulation The structural features that appear as minor peaks at around 4.5Å are only present for higher Na2O content glass – at 10% Na2O these appear as abroad shoulder of the first peak. In general, a higher first peak at low Na2O content corresponds to higher cluster of Na-Na pairs which, in with the present potential seems lower than observed with the Teter-Du potential. At higher concentrations, clustering and medium range ordering of sodium ions at longer distances results in the presence of smaller peaks up to about 11Å. results28,34.

The structure was further analyzed based on the bond angle distribution of O-Na-O angles. These angles were measure based on an O-Na-O pairs with Na-O distances lower than 3Å. This distance was chosen based on the Na-O g(r) which indicates that all first neighbor Oxygens of Sodium are covered by 3Å. The bond angle distribution is similar to what has been observed through simulations with other MD force fields28,33. Oxygen atoms in the silicate network can be classified as bridging if they connect the vertices of two SiO4 tetrahedra or non-bridging, if they are bound to only one silica. Na is a network modifier in the glass structure36 and hence can occupy positions near bridging or non-bridging oxygens. Depending on the type of oxygen close to the sodium atom, the O-Na-O bond angle can be either a sharp peak around 60˚ followed by a broad peak around 90˚ for BO-Na-BO, or a peak at 90˚ and a small shoulder at around 60˚ in the case of NBO-Na-NBO. Sodium ions close to NBO’s and BO’s have features that are a combination of both the 60˚ angle and the 90˚ peaks. This distinction in the peaks becomes apparent when the overall bond angle is decomposed into its contributing components as shown in figure 5a for the silicate with 30% Na2O. The cut-off distance to recognize bridging and non-bridging oxygens was set at 2.1Å based on gSi-O(r). Figure 5b shows the bond angle distribution for the three different compositions. All the 3 distributions shown in figure 5b have been normalized to the sodium composition to allow a direct comparison between the three. At low sodium compositions, there is a higher proportion of bridging oxygen. Addition of sodium disrupts the silica network and results in the formation of more non-briding oxygen. Thus, the 60˚ peak in the BAD is lowered relative to the 90˚ peak with increasing Na2O content in the silicate. 3.2 Structures of sodium silicate glasses: Si-O network structures In order to understand the Silicate network structures, we calculated gSi-O(r) and angles in the Si-O network. Figure 6(a) shows a plot of gSi-O(r) and 6(b) shows the distribution of Si-O-Si and O-Si-O bond angles for the 3 different compositions. The O-Si-O bond angle has a mean of around 109 degrees and the Si-O-Si angle is distributed around 155 degrees. As can be seen, gSiO(r) or the bond angles do not show a major change indicating no significant change in the overall structure of the Si-O network with composition. The inset of figure 6(a) shows a slight decrease in the Si-O bond distance from 1.62Å to 1.59Å and a corresponding narrowing of the O-Si-O 6 ACS Paragon Plus Environment

Page 6 of 25

Page 7 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


distribution and lowering of Si-O-Si angles. These changes can be attributed to the formation of NBO’s which increase with the Na2O concentration. Figure 6(c) shows that the fraction of NBO’s increases non linearly with composition as is expected37. The gSi-O(r) values for the 50% Na2O composition was examined in further detail to de-convolute the BO and NBO peaks. This calculation showed that the Si-NBO distance is slightly lower at 1.59Å compared to the Si-BO peak at 1.625Å. This has also been observed in other simulations38 and provides an explantion for the shift in Si-O bond distance and bond angle distributions with increasing Na2O content. The breakup of the amorphous silica network structure in glasses by the addition of modifiers has been studied experimentally by EXAFS and the nature of modified structure is well known36. Accordingly the continuous random network formed by Si-O4 tetrahedra is disrupted by channels of modifier ions (in this case Na+) and these channels are potential pathways for diffusion of these ions. As a result of this network modification, the number of bridging oxygens reduce in a direct correlation with the quantity of the modifier ion as shown in figure 6(a). To further elucidate these alkali channels, the bulk silicate glasses were run under NVT at 298K and 398K for 400ps and 200ps. A qualitative, visual representation of the 3D network for different sodium silicate glasses is shown in the Figure 7. The position of the sodium atoms at the end of the 600ps were used to create a channel structure using a “quicksurf” rendering with Visual Molecular Dynamics (VMD) simulation software39. The path traced by sodium atoms over the 600ps simulation shows that Na+ diffusion occurs either inside or close to the sodium channels and the Si-O network is rarely broken by diffusing sodium atoms. This is further demonstrated by the NBO fractions remaining relatively unchanged from figure 6(c) over the 600ps of simulation. The silicon atoms can also be classified into the Qn species according to the number of bridging oxygens attached, where n(0-4) is the number of attached bridging oxygens. The Qn species has in glass can be experimentally characterized by NMR studies40 and is a useful tool in quantifying the network structure. MD simulations have also been used to study the correlation between evolution of Qn species and glass structure41. Figure 8 is a plot of the fraction of Q-species as a function of Na2O content in the silicate. As before, the qualification for Si-O bonds was based on the pair distribution functions consistent with the calculation of bridging oxygens. The evolution of Q-species as a function of modifier composition can be explained by two different models42,43. The first is the binary model, where a higher Q species can only be lost through the formation of a Q species of immediately lower order. In effect, this means that once Q3’s is formed, through addition of a single Na+, the next additional Na+ will exist farther away from the already formed Q3. As a result only 2 types of Q species existing at a given alkali composition. In the second, random network model, increasing alkali content can result in any of the Q-species forming and thus, different Q-species can exist at a given alkali composition. This implies that additional Na+ ions do not avoid existing Q3 species. Experimental studies have shown that the existence of all Q-species with increasing alkali content indicating a tendency towards the random network model, with the caveat that some of the Q-species show a higher composition than is predicted by the random network model. It has been noted43 that that faster quench rates in alkali silicate glasses can result in the glass being more of a continuous random network. MD simulations with Teter potential has been shown to have a lower concentration of Q3 species 7 ACS Paragon Plus Environment


consistent with the random distribution model especially at higher cooling rates44. In the present study, the variation of Q-species fraction as a function of composition shows results consistent with those predicted for the random network model, which also agrees well with the higher quenching rates used in MD simulations. As expected, the NBO concentration increased steadily with increased alkali content resulting in disruption of the silicate network. The binary model, by definition stipulates that additional alkali ions necessarily are distributed farther away from existing ones resulting in the absence of segregation of alkali rich areas. The implication of the fractions of Q-species following the random network model is that siloxane bonds are broken at random depending on the presence of Na+ ions and do not necessarily wait for the formation of a consistent number of Q-species throughout the system as required by the binary model.

3.3 Vibrational density of states of sodium silicate glasses While the positions atoms are determined to a large extent by the location of the minima of the potential curve, the vibration modes and frequencies are determined by the curvature of the potential at these minima. Thus, as another test of validating the potential parameters, the vibrational density of states (VDOS) was calculated based on the Fourier transformation of the velocity auto correlation function (VACF). In the present work, the VDOS of melt-quenched sodium silicates for different compositions was calculated over 20000 steps of 1fs each. Figure 6 shows the VDOS curves for different compositions of sodium silicates. The typical feature to observe in VDOS obtained through MD simulations45,46 is the Si-O stretch mode at around 1100cm-1, Si-O-Si symmetric stretching at around 800cm-1 and the Si-O-Si bond rocking at 400cm1. The addition of sodium also results in the shift of the 1100cm-1 peak to lower frequencies. This red shift has been observed in other alkali silicate glasses as well47,48 and can be attributed to the stretch modes vibration being dominated by the increasing non bridging oxygens. The VDOS of silicates have been studied by this method45,46,49 which show a reasonable compliance with experimental VDOS50. The presence of sodium introduces an additional peak at around 100cm-1 which49 is clearly absent in pure silica and is shown as increasing peaks with increasing sodium concentration in figure 9. 3.4 Equation of state calculation of sodium silicate crystals and comparison with QM data An additional indicator of the validity of the potential could be established by studying the structure and energy minimized volumes of the crystalline sodium silicates. Indeed, E-V curve can also be associated with the bulk modulus and hence mechanical properties through the BirchMurnaghan equation of state51–53(EOS). In view of this, the Energy-Volume equation of state was determined for various crystal systems of sodium silicates by both Density Functional Theory (DFT) calculations as well as with MD simulations. Six different crystal systems of sodium silicate were obtained from the Inorganic Crystal Structure Database 54 for this purpose. Vienna ab initio package (VASP)30,31,55 was used to carry out the DFT calculations. The exchange correlation function is treated within the generalized gradient approximation (GGA) with the PBE form56, and the electronic structure such as the interaction between ions and electrons is described by projected augmented wave pseudopotentials. The electronic wave functions are 8 ACS Paragon Plus Environment

Page 8 of 25

Page 9 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


expanded in a plane wave up to a kinetic energy cutoff of 400 eV. A Gaussian smearing of 0.05 eV was applied. The total energies as a function of the volume around the experimental volume were calculated. At each volume the cell shape and atom positions were allowed to be fully relaxed while the cell volume was maintained constant. The obtained total energy versus volume curve was fitted to the Birch53 equation of state, based on which the equilibrium volume was obtained. For the MD calculations, LAMMPS minimization of energy was done for the same crystals after replicating them to obtain systems that are at least 20Å in each dimension. Energy minimization to 10-4 Kcal/mol accuracy and box relaxation after a small change in volume and the relaxed energies were plotted as a function of the volume. Figure 10 shows the comparison of the E-V curve obtained by DFT and MD. Overall, except for Na6Si2O7, the data from MD was close the DFT data. For Na6Si2O7 also, the MD data was shifted to lower volumes – although the shift was only by about 5% and the first derivative of the curves were reasonably accurate for both the sets of data. Comparison of MD-EOS with these exact DFT calculations was also included in the most recent version of ReaxFF Sodium parameterization17. Overall, while the ReaxFF version also shows good compliance of with DFT-EOS, MGFF shows a clearer distinction between the monoclinic and orthorhombic phases of sodium disilicate. 3.5 Simulations of NaOH dissociation in water One of the significant utility of MGFF is its ability to reproduce reactions between silicates and water19 with lower computational cost as compared to ReaxFF or QM calculations. Thus, to enable description of solvation and reaction of Na+ in water, parameters for Na-H two-body and Na-O-H three-body interactions were also parameterized based on the scheme explained earlier. To test the accuracy of the parameters, the solvation and reactions of Na-OH in water was used as a benchmark. Similar tests have been conducted for the recently published reactive NaOH/Silicates parameters of the ReaxFF potential also17. A cubic box of ~29Å side was filled with 800 water molecules and 10 NaOH molecules resulting in a NaOH concentration of 0.69M. This equilibrated at 298K to study the dissociation behavior of NaOH. At low concentrations, NaOH, being an ionic solute, is expected be completely dissociated without a measurable Kb value. However, reported values of Kb for NaOH are calculated based on thermodynamic data and vary between 0.63 and 5 57,58. Based on these values the expected fraction of dissociated NaOH can vary between 60-89%. As shown in figure 11, the tested system of 800 water molecules and 10 NaOH results in an average dissociated fraction of ~81%. The large fluctuation in the values indicates on going association and dissociation of NaOH via proton transfer as has been observed with the ReaxFF potentials17. The number of water molecules surrounding the Na+ ion was also observed to be steady around 5 which is consistent with the hydration number of Sodium in water59,60. Thus, the MD model has been shown to result in a structurally consistent model for crystalline as well as amorphous sodium silicate glasses. Vibrational spectra, analyzed by Fourier transform of VACF also yielded sufficiently consistent results compared to other simulations and experiments. 4. Conclusions 9 ACS Paragon Plus Environment


We report development of parameters for the simulation of sodium silicate glasses compatible with the silica-water MGFF potentials developed a few years earlier with the aim to enable simulations of sodium silicate glass – water reactions and interactions. We established that the parameter set results reproduce a very good representation of the structure and density of sodium silicate glasses over the whole glass formation range. Addition of Na+ ions a modifier to the continuous random network of silica glass results in the formation of NBO’s and breaking up of the Si-O network, as represented in Qn distribution changes, consistent with experimental observations. Once the reactivity of Na2O in glass was analyzed, we also established the reactivity in water by studying dissociation of NaOH solutions in water. These results compare well with the recently published ReaxFF parameterization of sodium silicate and can thus be used in the study of water interaction with multicomponent glasses to understand corrosion and environmental effect on structure and properties of these glasses. Additional advantages of this framework is its higher computational efficiency and less parameters involved that would expansion of the potential set to other multicomponent silicate glass composition. Acknowledgements: This research has been made possible by funding from the Center for Performance and Design of Nuclear Waste Forms and Containers, an Energy Frontier Research Center funded by the U.S. DOE, Office of Science, Basic Energy Sciences under Award No. DE-SC0016584. The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin and the HPC Talon3 at University of North Texas for providing computational resources that have contributed to the research results reported within this paper.

References: (1)

Varshneya, A. K. Fundamentals of the Glassy State. In Fundamentals of Inorganic Glasses; Elsevier, 1994; pp 13–25.

(2)

Iler, R. K. The Chemistry of Silica (Iler, Ralph K.); John Wiley and Sons: Wilmington, Delaware, 1979.

(3)

Gin, S. Open Scientific Questions about Nuclear Glass Corrosion. Procedia Mater. Sci. 2014, 7, 163–171.

(4)

Frankel, G. S.; Vienna, J. D.; Lian, J.; Scully, J. R.; Gin, S.; Ryan, J. V; Wang, J.; Kim, S. H.; Windl, W.; Du, J. A Comparative Review of the Aqueous Corrosion of Glasses , Crystalline Ceramics , and Metals. npj Mater. Degrad. 2018, 2 (April), 15.

(5)

Plimpton, S. J. LAMMPS Molecular Dynamics Simulator http://lammps.sandia.gov/.

(6)

Gin, S.; Collin, M.; Jollivet, P.; Fournier, M.; Minet, Y.; Dupuy, L.; Mahadevan, T.; Kerisit, S.; Du, J. Dynamics of Self-Reorganization Explains Passivation of Silicate Glasses. Nat. Commun. 2018, 9 (1), 2169.

(7)

Soules, T. F. A Molecular Dynamic Calculation of the Structure of Sodium Silicate Glasses. J. Chem. Phys. 1979, 71 (11), 4570–4578. 10 ACS Paragon Plus Environment

Page 10 of 25

Page 11 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(8)

Newell, R. G.; Feuston, B. P.; Garofalini, S. H. The Structure of Sodium Trisilicate Glass via Molecular Dynamics Employing Three-Body Potentials. J. Mater. Res. 1989, 4 (2), 434–439.

(9)

Meyer, A.; Horbach, J.; Kob, W.; Kargl, F.; Schober, H. Channel Formation and Intermediate Range Order in Sodium Silicate Melts and Glasses. Phys. Rev. Lett. 2004, 93 (2), 27801.

(10)

Macià Escatllar, A.; Ugliengo, P.; Bromley, S. T. Modeling Hydroxylated Nanosilica: Testing the Performance of ReaxFF and FFSiOH Force Fields. J. Chem. Phys. 2017, 146 (22).

(11)

Fogarty, J. C.; Aktulga, H. M.; Grama, A. Y.; van Duin, A. C. T.; Pandit, S. A. A Reactive Molecular Dynamics Simulation of the Silica-Water Interface. J. Chem. Phys. 2010, 132 (17), 174704.

(12)

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

(13)

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 (5), 1040–1053.

(14)

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

(15)

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 (43), 24803–24816.

(16)

Feuston, B. P.; Garofalini, S. H. Water-Induced Relaxation of the Vitreous Silica Surface. J. Appl. Phys. 1990, 68 (9), 4830–4836.

(17)

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 NaSiO X /Water Systems and Its Application to Sodium and Proton Self-Diffusion. J. Phys. Chem. C 2018, 122 (34), 19613–19624.

(18)

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. Glas. Sci. 2017, 8 (3), 276–284.

(19)

Mahadevan, T. S.; Du, J. Evaluating Water Reactivity at Silica Surfaces Using Reactive Potentials. J. Phys. Chem. C 2018, acs.jpcc.7b12653.

(20)

Guillot, B.; Guissani, Y. How to Build a Better Pair Potential for Water. J. Chem. Phys. 2001, 114 (15), 6720–6733.

(21)

Ma, Y.; Garofalini, S. H. Modified Wolf Electrostatic Summation: Incorporating an Empirical Charge Overlap. Mol. Simul. 2005, 31 (11), 739–748. 11 ACS Paragon Plus Environment


(22)

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

(23)

Du, J.; Cormack, A. N. Molecular Dynamics Simulation of the Structure and Hydroxylation of Silica Glass Surfaces. J. Am. Ceram. Soc. 2005, 88 (9), 2532–2539.

(24)

Huang, C.; Cormack, A. N. The Structure of Sodium Silicate Glass. J. Chem. Phys. 1990, 93 (11), 8180–8186.

(25)

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

(26)

Ercolesi, F.; Adams, J. B. Interatomic Potentials from First-Principles Calculations: The Force-Matching Method. Epl 1994, 26 (8), 583–588.

(27)

Masia, M.; Guàrdia, E.; Nicolini, P. The Force Matching Approach to Multiscale Simulations: Merits, Shortcomings, and Future Perspectives. Int. J. Quantum Chem. 2014, 114 (16), 1036–1040.

(28)

Du, J.; Cormack, A. N. The Medium Range Structure of Sodium Silicate Glasses: A Molecular Dynamics Simulation. J. Non. Cryst. Solids 2004, 349, 66–79.

(29)

Rempe, S. B.; Pratt, L. R. The Hydration Number of Na+ in Liquid Water. Fluid Phase Equilib. 2001, 183–184, 121–132.

(30)

Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54 (16), 11169–11186.

(31)

Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B 1993, 47 (1), 558–561.

(32)

BANSAL, N. P.; DOREMUS, R. H. Density. In Handbook of Glass Properties; Elsevier, 1986; pp 49–100.

(33)

Pedone, A.; Malavasi, G.; Cormack, A. N.; Segre, U.; Menziani, M. C. Insight into Elastic Properties of Binary Alkali Silicate Glasses; Prediction and Interpretation through Atomistic Simulation Techniques. Chem. Mater. 2007, 19 (13), 3144–3154.

(34)

Jabraoui, H.; Vaills, Y.; Hasnaoui, A.; Badawi, M.; Ouaskit, S. Effect of Sodium Oxide Modifier on Structural and Elastic Properties of Silicate Glass. J. Phys. Chem. B 2016, 120 (51), 13193–13205.

(35)

Merzbacher, C.; White, W. Structure of Na in Aluminosilicate Glasses : A Far-Infrared Reflectance Spectroscopic Study. Am. Mineral. 1988, 73, 1089–1094.

(36)

Greaves, G. N. EXAFS and the Structure of Glass. J. Non. Cryst. Solids 1985, 71 (1–3), 203–217.

(37)

Jen, J. S.; Kalinowski, M. R. An ESCA Study of the Bridging to Non-Bridging Oxygen Ratio in Sodium Silicate Glass and the Correlations to Glass Density and Refractive Index. J. Non. Cryst. Solids 1980, 38–39, 21–26. 12 ACS Paragon Plus Environment

Page 12 of 25

Page 13 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(38)

Melman, H.; Garofalini, S. H. Microstructural Evaluation of Simulated Sodium Silicate Glasses. J. Non. Cryst. Solids 1991, 134 (1–2), 107–115.

(39)

Cross, S.; Kuttel, M. M.; Stone, J. E.; Gain, J. E. Visualisation of Cyclic and MultiBranched Molecules with VMD. J. Mol. Graph. Model. 2009, 28 (2), 131–139.

(40)

Greaves, G. N.; Sen, S. Inorganic Glasses, Glass-Forming Liquids and Amorphizing Solids. Adv. Phys. 2007, 56 (1), 1–166.

(41)

Garofalini, S. H.; Martin, G. Molecular Simulations of the Polymerization of Silicic Acid Molecules and Network Formation. J. Phys. Chem. 1994, 98 (4), 1311–1316.

(42)

Eckert, H. Structural Characterization of Noncrystalline Solids and Glasses Using Solid State Nmr. 1992.

(43)

Greaves, G. N.; Sen, S. Inorganic Glasses, Glass-Forming Liquids and Amorphizing Solids. Adv. Phys. 2007, 56 (1), 1–166.

(44)

Li, X.; Song, W.; Yang, K.; Krishnan, N. M. A.; Wang, B.; Smedskjaer, M. M.; Mauro, J. C.; Sant, G.; Balonis, M.; Bauchy, M. Cooling Rate Effects in Sodium Silicate Glasses: Bridging the Gap between Molecular Dynamics Simulations and Experiments. J. Chem. Phys. 2017, 147 (7), 74501.

(45)

Garofalini, S. H. Molecular Dynamics Simulation of the Frequency Spectrum of Amorphous Silica. J. Chem. Phys. 1982, 76 (6), 3189–3192.

(46)

Machacek, J.; Gedeon, O.; Liska, M. Normal Mode Analysis of Na2O·3SiO2glass. Phys. Procedia 2013, 48, 85–88.

(47)

Du, J.; Corrales, L. R. Structure, Dynamics, and Electronic Properties of Lithium Disilicate Melt and Glass. J. Chem. Phys. 2006, 125 (11), 1–12.

(48)

Tian, Y.; Du, J.; Han, W.; Zu, X.; Yuan, X.; Zheng, W. Thermal Conductivity of Vitreous Silica from Molecular Dynamics Simulations: The Effects of Force Field, Heat Flux and System Size. J. Chem. Phys. 2017, 146 (5).

(49)

Laurent, O.; Mantisi, B.; Micoulaut, M. Structure and Topology of Soda-Lime Silicate Glasses: Implications for Window Glass. J. Phys. Chem. B 2014, 118 (44), 12750–12762.

(50)

Haworth, R.; Mountjoy, G.; Corno, M.; Ugliengo, P.; Newport, R. J. Probing Vibrational Modes in Silica Glass Using Inelastic Neutron Scattering with Mass Contrast. Phys. Rev. B - Condens. Matter Mater. Phys. 2010, 81 (6), 98–101.

(51)

Wang, Y.; Liu, Z. T. Y.; Khare, S. V.; Collins, S. A.; Zhang, J.; Wang, L.; Zhao, Y. Thermal Equation of State of Silicon Carbide. Appl. Phys. Lett. 2016, 108 (6).

(52)

Ghiorso, M. S. An Equation of State for Silicate Melts. I. Formulation of a General Model. Am. J. Sci. 2004, 304 (8–9), 637–678.

(53)

Birch, F. Finite Elastic Strain of Cubic Crystals. Phys. Rev. 1947, 71 (11), 809–824.

(54)

Hellenbrandt, M. The Inorganic Crystal Structure Database (ICSD)—Present and Future. Crystallogr. Rev. 2004, 10 (1), 17–22. 13 ACS Paragon Plus Environment


(55)

Kresse, G.; Furthmüller, J. Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6 (1), 15– 50.

(56)

Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77 (18), 3865–3868.

(57)

Davies, C. W. 280. The Electrolytic Dissociation of Metal Hydroxides. J. Chem. Soc. 1951, No. 1256, 1256–1258.

(58)

van Eldik, R. S. Kotrlý and L. Sucha: Handbook of Chemical Equilibria in Analytical Chemistry, Ellis Horwood Ltd., Chichester 1985. 414 Seiten, Preis: £45,-. Berichte der Bunsengesellschaft für Phys. Chemie 1985, 89 (11), 1249–1249.

(59)

Hellström, M.; Behler, J. Structure of Aqueous NaOH Solutions: Insights from NeuralNetwork-Based Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2017, 19 (1), 82–96.

(60)

Imamura, T.; Ishiyama, T.; Morita, A. Molecular Dynamics Analysis of NaOH Aqueous Solution Surface and the Sum Frequency Generation Spectra: Is Surface OH – Detected by SFG Spectroscopy? J. Phys. Chem. C 2014, 118 (50), 29017–29027.


Page 14 of 25

Page 15 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Graphical abstract

Graphical Abstract shows final snapshot of the silica network represented by ball and stick model with the sodium channels represented by the yellow blobs. The red spheres are the positions of Na+ over last 100ps at 298K. This demonstrates that the diffusion of sodium ions occurs close to the Na-rich channels and the Si – rich network does not see as many Na+ ions.



Page 16 of 25

Tables

Table 1: Pair terms used in the simulations – the  value (Wolf damping factor) was 4.46Å. If the Arep value is set to 0 for any interaction, the corresponding r was set to 1.00 to avoid division by zero. Pair O-O O-H O-Si O-Na H-H H-Si H-Na Si-Si Si-Na Na-Na

ξr(Å) 0.6100 0.1998 0.3730 0.5090 1.0 0.4300 1.0 0.6400 1.0 1.0

Arep(J) 4.25e-17 2.283e-16 2.67e-16 4e-17 0 5e-16 0 7e-17 0 0

C6(J/Å-6) 4.226e-18 0 7e-18 3e-18 0 3.8e-18 0 0 0 0

Table 2: Three body terms used in the simulations Triplet(j-i-k) Si-O-Si O-Si-O H-O-H Si-O-H Na-O-H

jik J 0.099e-17 1.49e-17 3.00e-17 2.12e-17 3.47e-17

rc(Å) (ij,ik when j≠k) 2.6 3.0 1.6 2.6,1.6 3.0,1.6

(Å) (ij,ik when j≠k) 2.0 2.8 1.3 2.0,1.3 2.6,1.3


o(degrees) 109.47 109.47 100 109.47 120

Page 17 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figures:

Figure 1: Pair potential energies plotted as a function of separation distance for Na-O, H-O, Si-O and Na-Na interactions.

Figure 2: Comparing the pair distribution functions of Na-O and Na-Na of melt-quenched sodium silicate (28Na2O-72SiO2) glass with the Teter-Du and MGFF potentials.



Figure 3: Densities of simulated sodium silicate glasses as compared to experimental values32. The density was averaged over an NPT run of 100ps which resulted in a standard deviation of less than 0.2%


Page 18 of 25

Page 19 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

(b) Figure 4: (a) Na-O pair distribution functions show an increasing peak and (b) Na-Na PDF showing a decrease in peak intensity with increasing Na2O content. The higher peak at low Na2O means that the probability of finding Na-Na pairs is higher at these compositions, which indicates a tendency to segregate into alkali rich channels.



(a) (b) Figure 5:(a) Variation in BAD with composition when compared with the origin of the BAD peaks 5(b) indicates that at low Na2O content briding oxygens are more in abundance as has been indicated with other data also.


Page 20 of 25

Page 21 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(a)

(b)



(c) Figure 6: (a) Si-O pair distribution function with the inset showing details. (b) Si-O-Si and O-Si-O angles for the three different compositions (c) the fraction of NBO and 4- coordinated Si, showing a non-linear increase in NBO concentration. The grey dotted line is a visual guide showing a linear increase in NBO concentration.

Figure 7: Visual depiction of the sodium channels, The final position of Na atoms after 600ps of simulation is overlayed on the silica network shown with grey(Si) and blue (O) spheres connected by bonds The yellow blobs are a result of rendering and connecting the Na atoms depicting the channels. The red spheres are the position of the Na atoms over the final 600ps of simulation. Note that while the network gets more interconnected with increasing Na content the diffusing sodium tend to remain close to these channels.


Page 22 of 25

Page 23 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8: The evolution of Q-species component indicates that addition of Na+ ions to silicate glasses can result in these ions migrating to random locations in the network depending on the local energies. EXAFS studies indicate a similar random distribution of additional Na+ to the network although the fraction of Q3 is higher in the slower cooled glasses.



Figure 9: The VDOS of simulated sodium silicate glasses indicate the increase in the intensity at the low frequency due to Na-O vibrations and also the effect of an ~1150cm-1 shoulder for the lowered Si-O stretch.

Figure 10: The Murnaghan equation of state is a useful indicator for calculating the bulk modulus of crystals – the above graphs show that the EOS is very similar with DFT energies (dots) and MD 24 ACS Paragon Plus Environment

Page 24 of 25

Page 25 of 25 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


simulations (solid lines). There is only a 5% difference in the volume of the low energy structure of Na6Si2O7

Figure 11: Plot of the number of dissociated NaOH (Na-OH distance > 3Å) in a solution of 10 NaOH in 800 water molecules. The average fraction is about 81% indicated by the solid line.


Development of Water Reactive Potentials for Sodium Silicate Glasses

Recommend Documents