Molecular Dynamics Simulation of Amorphous SiO2, B2O3, Na2O

Jun 27, 2019 - Molecular Dynamics Simulation of Amorphous SiO2, B2O3, Na2O–SiO2, Na2O–B2O3, and Na2O–B2O3–SiO2 Glasses with Variable ...
0 downloads 0 Views 4MB Size
Article Cite This: J. Phys. Chem. B 2019, 123, 6290−6302

pubs.acs.org/JPCB

Molecular Dynamics Simulation of Amorphous SiO2, B2O3, Na2O− SiO2, Na2O−B2O3, and Na2O−B2O3−SiO2 Glasses with Variable Compositions and with Cs2O and SrO Dopants Pooja Sahu,†,‡ Antariksh Avinash Pente,† Mandar Deepak Singh,† Islam Ali Chowdhri,† Kuldeep Sharma,§ Madhumita Goswami,§ Sk. Musharaf Ali,*,†,‡ K. T. Shenoy,† and Sadhana Mohan†

Downloaded via KEAN UNIV on July 28, 2019 at 09:12:08 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.



Chemical Engineering Division and §Glass and Advanced Materials Division, Bhabha Atomic Research Center, Mumbai 91400085, Maharashtra, India ‡ Homi Bhabha National Institute, Mumbai 400094, Maharashtra, India S Supporting Information *

ABSTRACT: Selection of suitable glass composition for vitrification of high-level radioactive wastes (HLWs) is one of the major challenges in nuclear waste reprocessing. Atomic and molecular level understanding of various structural, thermodynamical, and dynamical properties of a glass matrix can help in preliminary screening and thus reduce the dependency to some extent on tedious experimental procedures. In that context, extensive molecular dynamics (MD) simulations have been performed to calculate various microscopic properties of the glass matrix. The present article demonstrates that the “Buckingham potential-included long-ranged Coulomb interaction” can be utilized to simulate the glasses of varied compositions. The proposed simulation model has been validated for a wide range of glass compositions: pure glass matrixSiO2 and B2O3; binary glass mixturesSiO2−B2O3, Na2O−SiO2, and Na2O−B2O3; ternary glassNa2O−SiO2− B2O3; and also the Cs2O- and SrO-doped matrix of sodium borosilicate. Most importantly, the MD results have been validated with those of in-house synthesized glasses. The effect of alkali addition on the density and network connectivity of the glass matrix has been explored. The results capture well the boron anomalies for varied concentrations of network formers and network modifiers. The intermediate structural ordering in glasses has been explored by calculating the partial and total structure factors. Further, the characteristic vibration density of states of constituent atoms in the glass matrix is determined. In addition, the glass structures with the addition of dopant oxides Cs2O and SrO have been examined as they are known to be prime heat-generating agents in HLWs. The results establish the structure and dynamics of the doped glass matrix to be a complex nature of the dopant’s mass, concentration, charge, and ionic radius. The present MD results might be of great academic and technological significance for further studies in the field of vitrification and prediction of effects associated with the dopant’s nature and concentration.

1. INTRODUCTION

important role in the disposal of high-level radioactive wastes (HLWs).3−5 For immobilization of HLWs, borosilicate glass is preferred by most countries, with a composition of nearly 70% glass-forming agents and about 30% waste constituents.6 Prior to melting, the glass components are added to the liquid waste

Glass is an amorphous solid with a wide spectrum of applications because of its peculiar optical, mechanical, and chemical properties.1 Being optically transparent or translucent, it is used in window panes, tables, and lenses. Optical fibers made of glass are capable of transmitting data at a very high speed in the form of light.2 On adding certain metallic salts, one can impart color to the glass and increase its aesthetic value. In today’s scientific horizon, glass plays another very © 2019 American Chemical Society

Received: April 1, 2019 Revised: June 25, 2019 Published: June 27, 2019 6290

DOI: 10.1021/acs.jpcb.9b03026 J. Phys. Chem. B 2019, 123, 6290−6302

Article

The Journal of Physical Chemistry B stream, which produces a solidified radioactive product.7 The individual components of glass play one of the three basic roles in the glass structure: network former, intermediate, or modifier.8 Constituents such as silica and boric oxide are the major network formers. In addition to the role of network formers, the physical properties of vitrified glass are also known to have clear dependence on the nature and quantity of the added alkali modifiers,9 which is demonstrated by the glass transition temperature, viscosity, optical basicity, absorption edge, and thermal-expansion coefficient.10,11 In particular, the dependence on the nature of alkali modifier originates from structural variation. Earlier studies12,13 show that the effects introduced by the nature of alkali oxides are less pronounced than those induced by the quantity of modifiers. Such a disagreement on the role of metal ions drives the parallel studies on the structure and properties of the glass matrix with various compositions.14−16 Another important feature is the non-monotonic variation of the physical properties of the borosilicate glass matrix with a change in the alkali oxide concentration, a phenomenon known as the “boron anomaly.”17 Many studies on the borate glass aimed at elucidating the nature and relative population of the borate units while explaining the boron anomaly effect. Krogh-Moe et al.18 using infrared spectroscopy, Turner and co-workers19 using NMR, and Gruenhut and MacFarlane.20 using Raman spectroscopy have made pioneering contribution in this field. The results reveal that with the lower concentration of alkali oxides, the network connection is modified because of the conversion of bridging oxygen (BO) to non-BO (NBO). However, different yields of NBO conversion have been reported from previous studies;18−20 therefore, there is a need to be revisited again. Furthermore, the primary source of heat generation in nuclear waste glass is from the ionization of gamma rays and electronic energy deposition by alpha and beta particles.21,22 In HLWs, heat generation is caused mainly by the beta decay of the short-lived fission products, 90Sr and 137Cs, which are responsible for self-heating of the vitrified glass.22 However, only a few studies13,23−25 are available for the Cs- and Srdoped glasses and thus invite further exploration. It is also necessary to study the properties, such as leaching, chemical reactivity, and mechanical disintegration, which are effective over a very long time, in general, longer than the life span of one generation.26,27 Such events can be effectively predicted within a few days or weeks using molecular dynamics (MD) simulations. Apart from this, the immobilization of waste stream into vitrified glass requires handling of the highly toxic and radioactive materials. Preparation and experimental analysis of such systems on a laboratory scale is arduous and potentially hazardous.21 Conversely, the MD simulations and computational analysis can be used to prepare and analyze such systems with the alleviation of the abovementioned drawbacks of the laboratory methods.28,29 Also certain conventional analysis techniques such as neutron and X-ray scattering are one-dimensional in approach, whereas glasslike disordered solids need three-dimensional imaging of the network.30 This is a paramount benefit of using MD where the atomic structure can be directly visualized. However, the credibility of the MD simulations depends on the transferability of the potential model and force field parameters. A variety of potential models have been applied to simulate glasses using MD simulations namely Buckingham potential, Born Mayer Huggins (BMH) potential, Morse potential, and Tersoff potential.31−33 The use of Buckingham potential

successfully reproduced the radial distribution function (RDF) profile in close concordance with the X-ray diffraction data.20 However, the iterative conversion of the BMH potential to Buckingham potential including van der Waal terms will be tedious for systems having more than two types of atoms since many pair of combinations are possible. The studies by Takada et al.34 reported MD simulations of B2O3 using Buckingham− Coulombic potential for B−B and O−O interactions and Morse−Coulombic potential for the B−O interaction. The article gives a broad description of all possible structures of B2O3, which is of very high academic use. The work is ingenious in the use of coordination-dependent potential, which shows the presence of a significant number of boroxol rings. However, in the use of transferable potential, the maximum temperature used is 1500 K, which does not allow the system to expand rapidly by retaining the desired density. Also, the considered number of atoms is very low in all simulated systems. The studies by Tilocca et al.35 establish that the minimum number of 3000 atoms is necessary in order to simulate an amorphous system.31,35 Furthermore, Kieu et al.28 performed MD simulation of sodium borosilicate (NBS) glasses using the Buckingham potential with Coulombic terms and a modification of the composition-dependent AB−O parameter. The article explains the need for a potential that depends on the composition. However, in this study, the origin of assumption about the charge ratio of B atoms in BO3 and BO4− units as well for O atoms is understated and unexplained. Subsequently, Du and Cormack11 reported MD simulations of sodium silicate glasses using the Buckingham potential with parameters developed by Tilocca et al.,36 along with an additional modification of the repulsive term B/rn for r < r0, to overcome the problem of negative infinity as r → 0 at high temperatures. Hereby, the earlier studies utilized a wide spectrum of potential and potential parameters in order to capture boron anomaly. Also, the developed potentials have been validated only for a few compositions of glass, and hence there is a demand to conduct the study over a wide range of glass compositions. In addition to these, the effect of a dopant in the base NBS glass matrix has been explored earlier partially and needs further exploration. The present article demonstrates that Buckingham potentialincluded Coulomb interaction model potential can be utilized to simulate the glasses of varied compositions: pure glass SiO2 and B2O3; binary glass mixturesSiO2−B2O3, Na2O− SiO2, and Na2O−B2O3; ternary glass mixturesNa2O−SiO2− B2O3; and Cs2O-, SrO-doped borosilicate glass matrix. The effects associated with alkali addition have been validated by studying the various compositions of sodium silicate and sodium borate. Subsequently, the boron anomalies in the glass structure have been captured with various compositions of SiO2−B2O3, Na2O−B2O3, and Na2O−SiO2−B2O3. Further, the MD results have been validated with the experimentally synthesized NBS glass matrix in our laboratories. In addition, the effect of Sr and Cs dopings has also been studied. The MD results establish the structural distribution of the glass matrix in terms of density, glass transition temperature, number of BO atoms, fraction of BO4− units, and structure factors for intermediate ordering. The vibration dynamics of the constituent atoms has been explored by studying the characteristic vibration density of states (VDOS). The MD results have been corroborated by the experimental studies conducted with the in-house-synthesized Cs2O- and SrOdoped NBS matrix. 6291

DOI: 10.1021/acs.jpcb.9b03026 J. Phys. Chem. B 2019, 123, 6290−6302

The Journal of Physical Chemistry B

2. COMPUTATIONAL DETAILS The present simulations were conducted with LAMMPS37 MD package. The interaction between atoms was simulated using the combination of Buckingham potential (for short-range interactionsexpression shown as eq 1) and Coulomb potential (for long-range interactionsexpression shown by eq 2). N

UVdW =

N

∑ ∑ Aij e

−rij / ρij

i=1 j>i

UCoul =



qB′ qO′

qiqj (2)

The parameters ρij, Aij, and Cij control the narrowness of the potential model. The separation between the ith and jth atoms, having partial charges qi and qj, respectively, is represented by rij, and ε is the permittivity. The above potential model with ρij = 1/bij (where ρij is the ionic-pair dependent length parameter and bij is the inverse of ρij; the value of ρij and bij depends on the interacting atom species i and j) is also known as the BKS potential model, named after the developers Van Beest, Kramer, and van Santen.38 A similar potential form with different coefficient values has also been used by Kieu et al.28 and others39,40 in order to simulate the glass structures of varied compositions with numerous objectives. The coefficients ρij, Aij, and Cij, which are used in the present study for the short-range interactions, are shown in Table 1.

A (kcal/mol)

ρ (Å)

C (kcal Å6 mol−1)

Si−Si28 Si−O28 Si−B28 B−B28 B−O28 O−O28 Na−O28 Na−Na45 Cs−Cs45 Cs−O45 Sr−Sr31 Sr−O31,46

19 191.2 1 041 824.56 7767.1 2785.3 101 323.7 207 621.69 2 768 285.06 218 500 92 000 738 063.951 90 000.0 335 032.651

0.29 0.161 0.29 0.35 0.17 0.265 0.17 0.23 0.384615 0.267191 0.3 0.245015

0.0 1061.2085 0.0 0.0 0.0 1955.7383 0.0 0.0 0.0 5732.52 0.0 1880.779

(3)

K = NSiO2 /NB2O3

(4)

R = NNa 2O/NB2O3

(5)

qi′ = qi − NB

(qB′ − qB) N

; i = Na, Si, O, Cs, Sr

(6)

Here, NB and N, respectively, represent the number of boron atoms and the total number of atoms in the glass matrix. All simulated glass matrices were prepared with 4500 atoms. The binary glass mixtures of borosilicate and alkali-added glasses of sodium silicate and sodium borate were prepared and simulated with various compositions. The glass matrix of NBS was prepared with a stoichiometry of 21SiO2.6B2O3.7Na2O, similar to experimentally synthesized glass in our laboratories. Further, the Cs2O- and SrO-doped matrices of NBS, that is, NBS−5Cs, NBS−10Cs, NBS−5Sr, and NBS−10Sr, were prepared and simulated. The suffix numbers 5 and 10 in NBS represent the mass composition of doping oxides Cs2O and SrO in the NBS matrix (percentage of base matrix mass). The initial structures of all these glass matrices were generated using PACKMOL program.42 The doped matrices of NBS− 10Cs and NBS−10Sr were also prepared in the laboratory (synthesis details included in Scheme S1 of the Supporting Information). All systems were initially heated at 5000 K using the canonical (NVT) ensemble for 10 ns in order to remove the memory effects of the initial glass structure. This was followed by quenching at a rate of 5 K/ps in NVT. The systems were furthermore equilibrated for 20 ns dynamics using the isothermal−isobaric ensemble (NPT) at 300 K and P = 1 atm. Within 2−5 ns of the equilibration run, all glass matrices were seen to reach the desired density by change in volume. Further, the production runs were performed for 30 ns, and the generated data was utilized for the analysis of structural and dynamical properties of simulated glass. The simulation steps involved in the glass vitrification are depicted in Figure 1. All simulation studies were conducted with a time-step of 1.0 fs. To control the temperature, the Nose’ Hoover thermostat43 was used with the velocity Verlet integration scheme.43 In order to avoid the artifacts of surface boundaries, periodic boundary conditions were applied in all three directions.43 The images from the post analysis of simulation data were produced by visual MD package.44

Table 1. Pair Potential Parameters for the Glass Matrix pair

i

and A0 = 1.49643; A1 = 0.29504; A2 = −0.2565; A3 = 0.08721; A4 = −0.01323; A5 = 0.00073; A6 = 0.00315 (if R > 0.55); and A6 = 0 (if R ≤ 0.55). qB′ and qO′ represent the modified charge on boron and oxygen atoms. Further, eqs 3 and 6 were used to find the modified charges on atoms

(1)

4πεrij

ÉÑ ÑÑ ∑ AiR + A 0ÑÑÑÑÑ ÑÑÖ i=1 5

where

Cij rij 6

ÄÅ ÅÅ Å = −ÅÅÅA 6K 2 + ÅÅ ÅÇ

Article

Particle−particle−particle mesh method41 was employed to account for the long-range interactions. The charges on the atoms were taken to be composition-dependent, as used in ref 28. The partial charge model of atoms was selected for Coulombic interaction. The initial fixed partial charges assumed on atoms O, Si, B, and Na are −0.945, 1.89, 1.4175, and 0.4725, respectively, as per studies of Kieu et al.28 It is to be noted that the alkali ions with a similar vacancy were assigned identical initial charges, as followed by the studies of Wang et al.,40 and therefore, the charge on the doping atom Cs was the same as on Na and the charge on the Sr atom was considered to be 0.945 like other divalent alkali earth metal ions. Further, the assumed partial charges were modified as per the composition of the glass matrix using the following expressions28

3. RESULTS AND DISCUSSION 3.1. Validation of Methodology: Pure Glass Matrix of SiO2 and B2O3. First, the potential model and force field parameters were validated by simulating the pure glass matrix of amorphous SiO2 and amorphous B2O3, which are major network formers in the NBS glass matrix. Pure SiO2 provides very high mechanical strength and an extremely high chemical 6292

DOI: 10.1021/acs.jpcb.9b03026 J. Phys. Chem. B 2019, 123, 6290−6302

Article

The Journal of Physical Chemistry B

the mass by volume ratio of the glass matrix; mass M of the matrix was calculated as M = ∑mi·Ni, where mi represents the mass of the ith atom and Ni is total number of ith atoms in the glass matrix, and V is the volume of the matrix) was measured to be 2.38 g/cm3 for SiO2 and 1.82 g/cm3 for B2O3. The estimated density of amorphous SiO2 is found to be 8.18% higher than the experimental value; however, it is very close to the earlier reported values.47,48 The density of amorphous B2O3 was found to be in good agreement with the experimental value.34 Furthermore, the structural analysis of these glasses was done using RDF,49 estimated using eq 7. g (r ) = ρ− 2

Figure 1. Simulation steps for glass vitrification.

∑ ∑ δ(ri)δ(rj − r) i

resistance against leaching; however, having a very high melting point, it cannot be utilized alone in the glass matrix. B2O3 is added to the glass matrix in order to lower the melting temperature, which makes the glass vitrification practically viable. Both SiO2 and B2O3 exist in the form of crystalline as well as amorphous phase in the nature. However, with the objective of glass, the amorphous phase was simulated, images for which are shown in Figure 2a. The density (estimated from

=

V N2

j≠i

∑ ∑ δ(r − rij) i

j≠i

(7)

Here, ρ and N denote the number density and the total number of atoms/molecules in the simulation system, whereas V and δ represent the volume and Dirac delta, respectively. The results show the peak positions at 1.62, 2.58, and 3.17 Å, respectively, for Si−O, O−O, and Si−Si in SiO2 and at 1.35, 2.37, and 2.57 Å corresponding to B−O, O−O, and B−B in B2O3 (RDF profiles shown in Figure S1). The coordination number of Si and B atoms was seen to be varied from 1 to 5 and 1 to 4, respectively, as shown in Figure 2b. Nevertheless, peaks were noticed at CN = 4 for Si−O and at CN = 3 for B− O, representing the tetrahedral network structure for SiO2 and the trigonal connection of B with three oxygen atoms. The same was furthermore corroborated by the angle distribution of O−Si−O and O−B−O shown in Figure 2c, where the maximum peak intensity was observed at 104.37° and 120°, respectively, reflecting the tetrahedral coordination of O to Si and the trigonal geometry of BO3 units. It should be noted that both the O−Si−O and O−B−O angle distribution profiles show a shoulder at ∼90°, which might be indicative of 5-fold coordinated Si in SiO2 and 4-fold coordinated B in B2O3 glass. The availability of 5-fold Si and 4-fold coordinated B was quantified as ∼7 and 0.28%, respectively. Furthermore, the characteristic structural factor of SiO2 and B2O3 glass was calculated by Fourier transform of pair distribution functions using eqs 8 and 9,50 which gives a better understanding of the glass structure for intermediate and long distances. Sij(k) = 1 + ρ

∫0



4πr 2⌈gij(r ) − 1⌉

sin(kr ) dr kr

(8)

where ρ is the number density and Sij(k) represents the partial structure factor for the i−j pair at frequency k, obtained by the Fourier transform of the corresponding pair correlation function gij(r). The total structure factor50 was then obtained by weightage average of the partial structure factor,50 as shown below. Figure 2. (a) Snapshot for the structure of amorphous SiO2 and amorphous B2O3 (color code: Si: white, O: red, and B: blue). Inset images show SiO4 tetrahedra, BO3 planar structure, and BO4− tetrahedral unit. (b) Coordination distribution profile and (c) angle distribution profile for SiO2 and B2O3 glass. (d) Total structure factor and (e) VDOS for SiO2 and B2O3.

S(k ) =

∑i , j cicjlil jSij(k) ∑i , j cicjlil j

(9)

where ci and li, respectively, represent the fraction and neutronscattering length of the ith atom in pair i−j. 6293

DOI: 10.1021/acs.jpcb.9b03026 J. Phys. Chem. B 2019, 123, 6290−6302

Article

The Journal of Physical Chemistry B

mode and are related to the vibration of B atoms.57 On the other hand, the peaks at higher frequency of 960 and 1180 cm−1 represent the vibration of oxygen atoms connected to trigonal coordinated boron. It has to be noted that the vibrational frequency of this mode depends on the strength of the O−B−O bond-angle interaction and therefore might deviate slightly for other potential models. Essentially, the selected potential parameters well reproduce the experimental VDOS and earlier reported data.55−57 3.2. Binary Mixtures of Borosilicate Glass. Being the base structure of the NBS matrix, it is important to study the characteristic properties of the borosilicate glass. The structural models of the borosilicate glass were prepared using different compositions of SiO2 and B2O3 in the matrix. The results in Figure 3a show that the density of glass increased from 1.75 to

As common features, one can notice that the structure factor of SiO2 in Figure 2d shows the first peak at ∼2.15−2.25 Å−1 and then a second peak at around 5.45 Å−1. Similar to this, the structure factor of B2O3 shows a first peak at 2.4−2.45 Å−1 and a negative dip at 3.0 Å−1. These peaks are supposed to be related with the local binding inside the structural units of SiOn and BOn. Another important characteristic is the presence of a negative dip at smaller wave vectors around 1.2−1.3 Å−1 for SiO2 and 1.45−1.65 Å−1 for B2O3, known as first sharp diffraction peak (FSDP), and is associated with the two connected SiO4 tetrahedrals and the two connected trigonal boron (BO3) units, respectively. The position of the pronounced peaks or negative dips is in good agreement with the results reported in various literature studies.51−53 However, in view of comparison with the experimental data, one might have difficulty in observing the FSDP from the neutron-scattering experiment, as the peaks/dips at a small k factor [1−2 Å−1] are not of much significance. Furthermore, the selected potential model for SiO2 and B2O3 was validated from the characteristic vibration spectrum VDOS [ψ(v)] determined by performing the Fourier transform of the mass weightage velocity autocorrelation function,43,54 as shown by eqs 10−12. ψ (v ) =

2 lim kT τ →∞

τ

∫−τ C(t )e−i2πvt dt

(10)

Here, T is the temperature of the system and k is the Boltzmann constant. The velocity autocorrelation function, C(t), can be expressed as N

C(t ) =

Figure 3. (a) Density and (b) BO atoms (left Y axis) and BO4− (B4) units (right Y axis) in the borosilicate glass matrix.

3

∑ ∑ ⟨mici k(t )⟩

2.14 g/cm3when there is an increase in the SiO2 concentration from 10 to 70% (by mole), which is in excellent agreement with the earlier studies reported by Doweidar et al.58 Also, it was noticed that the increase in the silica content leads to increased conversion of planar BO3 (B3) units into BO4− (B4) tetrahedrals, as the fraction of B4 was found to increase with the increase in SiO2 concentration in the glass matrix, as shown by the right Y axis of Figure 3b. The same was furthermore supported by the increasing number of BO atoms with increase in SiO2 percentage shown by left Y axis of Figure 3b. In other words, the network connection becomes stronger with higher concentration of SiO2 in the glass matrix. In the present study, the number of BO atoms were determined by counting the number of oxygen atoms which were connected to two Si/B provided that each Si/B was bonded to the four/three nearest oxygen atoms. The observed trends are in accordance with the results reported by Wang et al.40 However, the reported density and B4 units by Wang et al.40 are marginally higher than the presented data because of the presence of Na2O and CaO alkali ions in the reference borosilicate matrix,40 whereas the results in Figure 3a,b are reported with the binary mixture of SiO2−B2O3 only, without any alkali ions in the glass matrix. No significant variation in the local structure of SiO2 was observed with variation in the composition of the borosilicate matrix. The peak for Si−O was always noticed at 1.61 Å with the tetrahedral coordination of O to Si. On the other hand, the coordination environment of B was seen to be varied with the addition of SiO2 in the glass, as observed from increasing BO atoms and B4 in Figure 3b. Though the addition of SiO2 in borosilicate improves the network connectivity and hence the chemical and mechanical strength of the glass matrix; nevertheless, the very high concentration of SiO2 in the

(11)

i=1 k=1

where ci k(t ) = lim

τ →∞

1 2τ

τ

∫−τ vik(t′ + t )vik(t′) dt′

(12) k

Here, mi denotes the mass of the ith molecule, vi (t) is the kth velocity component of the ith molecule at time t, and cik(t) is represented as a velocity autocorrelation function for molecule i in the kth direction at time t. The vibrational properties play a significant role in the estimation of specific heat, heat conductance, and transparency of glass. The data shown in Figure 2e is in agreement with the studies by Kob and Ispas55 as well as Tanguy et al.56 The VDOS of SiO2 show a flat peak at 200−600 cm−1 and another peak in between ∼1050 and 1070 cm−1. The lower frequency peak is known to have originated from the bending and rocking motion of oxygen atoms. To some extent, the broadening of the lower frequency peak reflects the disordered structure of amorphous SiO2; however, most of the part comes from the inability of the potential model to accurately reproduce the VDOS spectra. The peak broadening is seen to have significantly reduced in the case of first principles studies as well as in the experimental studies of VDOS for SiO2.55,56 Nevertheless, such a flat structure of VDOS of SiO2 at the lower frequency zone is not new but has been reported by many earlier MD simulation studies.55,56 Further, the peaks at higher frequency ∼1050−1070 cm−1 in the VDOS of SiO2 capture well the modes concerning intra-tetrahedral excitations. Additionally, the VDOS spectrum of B2O3 shows four characteristic peaks at 130, 670, 960, and 1180 cm−1. The first two peaks at 130 and 670 cm−1 represent the breathing 6294

DOI: 10.1021/acs.jpcb.9b03026 J. Phys. Chem. B 2019, 123, 6290−6302

Article

The Journal of Physical Chemistry B

On the contrary, with Na2O composition higher than 40%, the added sodium ions are observed to participate as the network modifier, and then both the density and B4 fraction were reduced with the increase in the Na2O concentration. Similar observations have been reported by Shelby et al.60 Additionally, the altering behavior of sodium ions in sodiumborate glass was furthermore confirmed by the angle distribution profiles for angle O−B−O, where the peak intensity was observed to be reduced for 120°; on the contrary, it increased for ∼109.7° with the increase in the Na2O content from 0 to 40%, and the reverse was followed for a further increase in the Na2O concentration from 40 to 60%, as shown in Figure 4c. In fact, the relative change in the peak intensity of angle distribution from 120° ↔ 109.7° reflects the exchange of BO3 ↔ BO4− and vice versa. The observed anomalies in the B coordination and the altered behavior of Na2O in sodium-borate glass are in line with the previous reported data.59−61 In sequence, akin to the sodium-silicate matrix, the number of BO atoms were always found to reduce with the increase in the Na2O concentration in the sodiumborate glass, as reflected from the data shown in Figure 4d. 3.4. NBS Glass Matrix. Following the validation of methodologies for pure glasses and binary glass mixtures, the NBS glass matrix was simulated with a composition of 21SiO2.6B2O3.7Na2O, that is, 59.70% SiO2, 19.76% B2O3, and 20.53 %Na2O. The structure of the simulated NBS matrix is shown in Figure 5a. The density of the MD-simulated NBS

borosilicate matrix is avoided as it would increase the melting temperature of the glass matrix. On the contrary, too much B2O3 content will reduce the chemical strength of glass as B2O3 has a low leaching resistance, but yet it is added to reduce the melting temperature of glass for easy operation and to avoid the engineering problems of material of construction. 3.3. Effect of Alkali Addition: Sodium Silicate and Sodium Borate. Alkali addition to the glass matrix reduces the operating temperature of glass by increasing the thermal conductivity, ionic conductivity, and electrical conductivity. However, the presence of excess alkali may lead to devitrification of glass and therefore must be avoided. The same can be predicted from Figure 4a, where the number of

Figure 4. (a) Number of BO atoms in sodium-silicate glass, (b) density (left Y axis) and B4 units (right Y axis) in the sodium-borate glass matrix, (c) angle distribution profile for O−B−O in the sodiumborate glass matrix, and (d) number of BO atoms in the sodium borate glass matrix.

BO atoms was found to be reduced with the increase in the Na2O concentration in sodium-silicate glass.12 The distribution of Na 2 O in sodium-silicate glasses with Na 2 O mole contribution of 10 and 60% is shown in Figure S2. Furthermore, very interesting observations were noticed with the sodium addition to the B2O3 matrix. Initially, similar to sodium silicate, the density of sodium borate was seen to be increasing with the increase in the Na2O concentration59 till 40%; thereafter, a decrease in the density was observed with further addition of Na2O in the borate matrix, as shown by the left Y axis of Figure 4b. To be precise, the initially added sodium ions (for Na2O < 40%) were seen to reside within the cavities formed by boroxyl rings. Occupancy of sodium ions in the interstitial positions causes larger increase in mass than the volume of glass matrix, and therefore leads to an increase in density for the sodium-borate matrix with the Na2O addition for concentration smaller than 40%. At the same time, the presence of alkali ions within the boroxyl cavities provides the charge compensation to BO4− tetrahedrals and so contribute to BO3 → BO4− conversion. As a result, the fraction of B4 in the sodium-borate glass matrix was observed to be increased with increase in the Na2O moles for a concentration smaller than 40%, as shown by the right Y axis of Figure 4b.

Figure 5. (a) Snapshot for the NBS glass matrix (SiO4 tetrahedrals: blue; BO3 and BO4− units: purple; sodium ions: yellow). (b) Coordination number distribution, and (c) angle distribution in the NBS glass matrix.

glass matrix was estimated to be 2.42 g/cm3, which was in good agreement with the density (*2.4 g/cm3) of the in-house synthesized glass matrix of the same composition. Details of the synthesis are included in the Supporting Information as Scheme S1. Figure 5b shows the coordination number distribution of Si and B for oxygen atoms. Using the cutoff radius of 2.15 Å for Si−O and B−O and 3.1 Å for Na−O, the average coordination number of Si, B, and Na was estimated to be 4.0, 3.8, and 10.28, respectively. The RDF peak position for Si−O, B−O, and Na−O was noticed at 1.651, 1.48, and 2.55 Å in sequence.28 The angle distribution profiles in Figure 5c represent the tetrahedral coordination of oxygen to Si, as peak for O−Si−O was obtained at 109.4°. Two peaks were noticed for O−B−O at 109.4° as well as at 120°, indicating two 6295

DOI: 10.1021/acs.jpcb.9b03026 J. Phys. Chem. B 2019, 123, 6290−6302

Article

The Journal of Physical Chemistry B Table 2. Composition of Simulated NBS Glass M[SiO2/(SiO2 + B2O3)]

SiO2

B2O3

Na2O

R[Na2O/B2O3]

K[SiO2/B2O3]

171 187 209 226 230 237 243

0.22 0.29 0.45 0.89 1.17 1.78 3.47

0.25 0.67 1.50 4.03 5.66 9.08 18.67

266 250 237 224 215

1.01 0.67 0.50 0.40 0.33

3.01 1.67 1.00 0.59 0.33

NBS15 0.20 0.40 0.60 0.80 0.85 0.90 0.95

194 423 691 1019 1110 1208 1307

781 634 460 253 196 133 70

0.75 0.63 0.50 0.37 0.25

794 625 473 336 215

264 375 474 564 642

NBS20

Furthermore, one can notice that the density of the NBS20 matrix is marginally higher than that of the NBS15 matrix for M smaller than 0.6. This is in line with the earlier results noted for the density of sodium silicate and sodium borate, where the density was observed to be increased with the increase in the alkali concentration, as discussed in Section 3.3. Essentially, M represents the silica content of the glass matrix; therefore, the network connectivity was increased with increase in M, as reflected from the increasing BO atoms as well as increasing B4 fraction as displayed in Figure 6. In addition, the results display comparatively lower BO atoms and higher B4 fraction for NBS20 than NBS15 as expected. Particularly, the reduced network connectivity is predictable with the higher concentration of sodium in the glass matrix, also shown in Figure 4a,d for sodium silicate and sodium borate, respectively. Likewise, the data of Figure 6 supports the observation (made in Section 3.3) that the B4 fraction increases with the increase in the sodium concentration for Na2O concentration smaller than 40% (as considered in this section). Another important point to be noted from Figure 6 is the change in the trend of density, BO atoms, and B4 profiles that happen nearly at R ≥ 1. The presented results can be found in good agreement with the earlier results reported by Wang et al.40 and Barlet et al.62 3.5. Effect of Cs/Sr Doping in NBS Glass. Cs137 and Sr90 are known as the major heat-generating isotopes in the nuclear waste. If added to a glass matrix, they lead to heating of the glass matrix over a period of initial 600 years of deposition.22 The study of Cs/Sr doping is quite important from the perspective of chemical strength as glasses subjected to high temperature are easily prone to leaching. Also, too much content of Cs/Sr may lead to devitrification of the glass matrix. In addition, Cs-based glass matrix13 (Cs-pencils) can be utilized for medical applications. In India, presently Co60 is being used as a gamma source in medical applications. However, very frequent replacements are required with Co60 because of the very short half-life of 5.27 years. On the other hand, Cs137, being the major constituent of the nuclear waste and having longer half-life (∼30 years), can be a suitable candidate to replace Co60. This would not only reduce the reactivity and revenue losses but also can be utilized over a longer period of time. On the other hand, Sr90 finds potential use in nonsteroidal anti-inflammatory drugs, which are used to treat pain associated with conditions like rheumatoid arthritis, migraine, sprains of muscles and joints, gout, and in mild to

different coordination environments for B atoms in the NBS glass matrix. The peak at 109.4° corresponds to tetrahedral coordination of B with the oxygen atoms in BO4− units, and another peak at 120° represents the B atom coordinated with three oxygen atoms in the BO3 unit. The relative concentration of BO4− and BO3 was estimated to be ∼63 and 29%, respectively. The peaks for angles Si−O−Si, B−O−B, and O− Na−O were found to be rather broad in nature, with a peak position of 146°, 138°, and 151°, respectively, which is in good agreement with the available data.28 Furthermore, two set of simulations, NBS15 and NBS20, were performed with Na2O molar compositions of 15 and 20%, respectively. The corresponding values of R(Na2O/ B2O3), K(SiO2/B2O3), and M[SiO2/(SiO2 + B2O3)] are shown in Table 2. Figure 6 represents the variation in density, BO

Figure 6. Density, % BO, and % B4 units in NBS glass vs M = [SiO2/ (SiO2 + B2O3)]; results for NBS15 are shown by a connected dashed line, whereas the results of NBS20 are shown by a symbol only.

atoms, and BO4− tetrahedral units (B4) with respect to M. The results show an increase in density with the increase in M till M = 0.8; thereafter, reduced density was noticed as M approached 1, which corresponds to the pure SiO2 matrix. Similar effects were noticed with the increase in SiO2 concentration in the borosilicate matrix, as discussed in Section 3.2, and also have been reported earlier by Wang et al.40 for the NBS glass matrix. 6296

DOI: 10.1021/acs.jpcb.9b03026 J. Phys. Chem. B 2019, 123, 6290−6302

Article

The Journal of Physical Chemistry B

alkali content from 12.28% in the NBS matrix to 13.45 and 13.72% in the doped matrix with 10% Cs and 10% Sr, respectively. The results can be found in good agreement with the studies of Avramov et al.,64 where almost a linear decrease in Tg with the mole fraction of alkali modifiers has been reported. In particular, addition of alkali ions in the glass matrix leads to the opening of the network structure, which reduces the melting point of the glass matrix. In line with the alkali content, Tg of Sr-doped glass was found to be smaller than that of Cs-doped glass. However, the observed difference in Tg of Cs- and Sr-doped glass is expected from the higher number of NBO atoms (i.e., more open glass network) created during Sr doping than the Cs doping in the glass matrix. This has been discussed more clearly in the next paragraph. The network connectivity in the doped matrix was determined by counting the BO atoms and estimating the coordination environment. Figure 8a displays the percentage of

moderate fever in some cases. Therefore, the studies were extended to examine the effects of Cs and Sr doping in NBS glass, named as NBS−5Cs, NBS−5Sr, NBS−10Cs, and NBS− 10Sr. The suffix numbers 5 and 10 in NBS represent the mass composition of doping oxides Cs2O and SrO in the NBS matrix (percentage of the base matrix mass). The images of doped glasses are shown in Figure S3. The doped glasses NBS−10Cs and NBS−10Sr were also synthesized in our laboratories using the same methodology as discussed in Scheme S1. The density and glass transition temperature (Tg) of NBS and doped NBS matrix determined from the MD simulations as well as from the experiments (details are included in the Supporting Information as Schemes S2 and S3) are presented in Table 3. Table 3. Density and Glass Transition Temperature of Doped NBS Glass Matrixa glass

density (*exp) [g/cm3]

Tg (*exp) [K]

NBS NBS−5Cs NBS−10Cs NBS−5Sr NBS−10Sr

2.42 (*2.4) 2.49 2.56 (*2.5) 2.49 2.55 (*2.6)

590 (*560) 545 (*532) 540 (*526)

a

The suffix numbers 5 and 10 in NBS represent the mass composition (percentage of the base matrix mass) of doping oxides Cs2O and SrO in the NBS matrix.

In MD simulations, Tg was obtained from the plot of inverse density versus temperature, with a range of temperature corresponding from a glass state to a melt state. The data points of both the zones were fitted to a straight line using the least squares method, and the intersection point of the fitted lines to the glass zone and melt zone was used to determine the glass transition temperature,63 as shown in Figure 7. On the

Figure 8. (a) BO and (b) coordination number distribution in the NBS and doped NBS glass matrix. Coordination number profile vs cutoff radius for the (c) Cs-doped NBS glass matrix and (d) Sr-doped NBS glass matrix

BO atoms for these systems. The results show a nearly similar number of BO atoms for dopants Cs as well as Sr when added in 5% by mass. However, with the increase in the dopant mass by 10%, a sudden decrease in BO atoms was noticed with Sr, whereas for Cs, the drop in BO atoms was not that significant. In fact, the conversion of BO to NBO depends on the charge and radius of the dopant. The larger the dopant radius, the more the conversion of BO to NBO is expected. Second, the dopants with higher charge will convert more BO to NBO. For example, with the addition of one Cs atom, maximum one BO will be converted to NBO, whereas the presence of one Sr atom in the glass matrix might convert two BO atoms into NBO atoms. The point to be noted is that the ionic radius of Cs (1.69 Å) is larger than that of Sr (1.33 Å), however, the charge of Cs is smaller than that of Sr. Hereby, the ionic radius and charge on Cs and Sr interplays in the reverse direction. It seems that for a lower concentration of the dopant, the network connectivity is dominated by the ionic radius. However, for a higher concentration of the dopant (e.g., 10%

Figure 7. Glass transition temperature Tg for NBS (left Y axis) and doped NBS with 10% Cs2O (NBS−10Cs) and 10% SrO (NBS−10Sr) (right Y axis).

other hand, experimental Tg was estimated from the differential thermal analysis, as discussed in Scheme S3 of the Supporting Information. The results in Table 3 show that the glass transition temperature obtained from the MD simulations is marginally higher than that from the experimental data, which might be associated with the evaporation of alkali ions at higher temperature in the experiments. In addition, the region used for fitting of data points might affect the estimated Tg in the experiments as well as in the MD simulations. Furthermore, a decrease in Tg can be noticed while doping, which is supposed to originate from the overall increase in the 6297

DOI: 10.1021/acs.jpcb.9b03026 J. Phys. Chem. B 2019, 123, 6290−6302

Article

The Journal of Physical Chemistry B

matrix (2.15−2.25 Å−1) and preshifted compared to B2O3 (2.4−2.45 Å−1), signifying the strength of local binding inside the structural units of SiOn and BOn in NBS and doped NBS to be in between those of SiO2 and B2O3 glass matrix. Further, the decrease in the intensity of peaks in NBS and doped NBS compared to the pure glass matrix is expected from the presence of alkali ions in the NBS and doped NBS. As the alkali ions create NBO atoms in the glass matrix, the larger the number of NBO atoms, the lesser the peak intensity. A more clear understanding of the total structure can be obtained from the partial structure factors. Therefore, the partial structure factor for pairs Si−Si, B−B, Si−B, Si−Na, B− Na, Na−Na, Si−dopant, B−dopant, and dopant−dopant was estimated and shown in Figure 9b. For any considered pair, one might notice the significant difference in FSDP (1.2−1.8 Å−1) compared to other peaks in the partial structure factors of glass with different compositions. On the other hand, only marginal effects were noticed for the peaks at ∼2.15−2.5 Å−1. This shows that though connection of MOn (M = Si, B, Na Cs, and Sr) units largely depends on the composition of the glass matrix, however, the local binding inside the structural units of MOn is only marginally affected by the glass composition. The results show split of FSDP (see the partial structure of Si−Na, Si-dopant, and B-dopant) into two subpeaks or subdips with the increase in the alkali concentration or with increased NBO atoms in the system. Also, the difference in the peak positions for the partial structure factor of M−Na, M−Cs, and M−Sr (where M represents the network forming elements B or Si) shows the different connecting environment of them in the glass matrix. For the simulated glass, dopant−dopant peaks were not much pronounced, which is due to their very small concentration in the glass matrix. This is why the doping effects could not be visualized from the total structure factor. However, the marked difference in the local environment of Na, Cs, and Sr can be observed from the corresponding partial structural factors. Likewise, the characteristic VDOS of NBS and doped NBS were estimated, and the results are shown in Figure 10. The VDOS show a first peak at ∼120 cm−1, a second broad peak at ∼430 cm−1, a shoulder around ∼700 cm−1 in the lower frequency zone, and a peak at ∼1027 cm−1 in the highfrequency zone. In particular, the peaks at ∼65−940 cm−1 are due to rocking, bending, and stretching contributions of Si− O−Si. The intermediate peaks at 400−600 cm−1 are from the combined effect of symmetric stretching of bonding oxygen atoms T−OBO−T (T = B, Si) + breathing vibration of four- to six-membered rings. The shoulder at 700 cm−1 is specifically due to Si−OBO−Si stretching, whereas the peak at higher frequency (∼1027 cm−1) in the VDOS spectrum contributes from the stretching motion of non-BO atoms connected to Si. One can notice the higher peak intensity at ∼120 cm−1 for the NBS matrix compared to the pure matrix, and effects become more pronounced with the doping of NBS. This reflects the increased bending and stretching motion of Si−O−Si with alkali addition. In particular, the addition of alkali ions to the glass matrix makes the glass structure more flexible due to opening of the network. The earlier results by Fu et al.65 showed the widening of peaks at 500−700 cm−1 with the increase in the non-BO atoms in the system, and the same can be visualized from the VDOS of NBS and doped NBS compared to the pure glass matrix reported herein. In agreement to the studies by Pedone et al.,12 the alkali-added matrices of NBS and doped NBS show a very high peak

by mass), contribution of charges on the dopants becomes more significant. Successively, the probability distribution for the coordination of Na−O, Cs−O, and Sr−O in NBS and doped NBS was estimated, and the corresponding results are shown in Figure 8b. The cutoff radius was selected to be 3.75 Å for Na−O and 4.05 Å for dopant-O coordination. The results show a variation in the coordination number of alkali ions with a change in the dopant concentration. In particular, the reduced coordination of Na−O was observed with the addition of dopant in the glass matrix. The average coordination of Na−O was reduced from 10.28 to 8.5 with the doping of the glass matrix with Cs, as shown in Table 4. The average coordination of oxygen to Table 4. Coordination Number of Alkali Ions in the NBS and Doped NBS Matrixa CN [Na−O] Rcut = 3.75 Å SBN SBN-5Cs SBN-10Cs SBN-5Sr SBN-10Sr

10.28 9.48 8.47 9.51 8.77

CN [Cs−O] Rcut = 4.05 Å

CN [Sr−O] Rcut = 4.05 Å

13.08 12.82 9.04 9.06

a

The suffix numbers 5 and 10 to NBS represent the mass composition of doping oxides Cs2O and SrO (percentage of base matrix mass).

doping elements Cs and Sr was found to be ∼13 and ∼9, respectively. Only a marginal change in the coordination number of these doping agents was noticed with the change in the dopant concentration for the considered range. In contrast to network formers, the coordination of dopants seems to be quirky as no sharp peak was noticed for the coordination of Cs and Sr. In fact, the structural environment in the vicinity of these dopants was found to be highly dependent on the cutoff radius, and the same is shown in Figure 8c,d, respectively, for Cs and Sr doping in the NBS matrix. Furthermore, the characteristic total structure factor was estimated for NBS and doped NBS glass matrix and is displayed in Figure 9a. NBS glass shows a global first peak at 0.95−1.1 Å−1, whereas for doped glasses, this peak seems to split into two subpeaks. This can be found in agreement with the earlier studies of Pacaud et al.,50 where the FSDP of the NBS glass matrix was observed to split into subpeaks with the increase in the alkali concentration. As stated earlier, the FSDP is associated with two connected SiO4 tetrahedrals and/or two trigonal boron units; therefore, postshifting of FSDP of NBS and doped NBS compared to SiO2 and preshifting of NBS and doped NBS compared to B2O3 reflects that the connection between SiO4 tetrahedrals or trigonal BO3 units in alkali-added glasses is weaker than the pure SiO2 matrix but stronger than the pure B2O3 matrix. In other words, the strength of NBS and doped NBS can be expected in between those of the pure matrix of SiO2 and B2O3. Another important observation was the difference in the intensity of peaks and dips, which was found to be in the order of NBS > NBS−5Sr ∼ NBS−5Cs ∼ NBS−10Cs > NBS−10Sr. Interestingly, the order of difference can be found in accordance with the BO atoms observed (see Figure 8a). The same was true for the second peak as well. In particular, the second peak for the NBS glass matrix was noticed at 2.15 Å−1, which was postshifted to 2.35−2.45 Å−1 for doped glasses. Akin to FSDP, the second and higher peaks in NBS and doped NBS were postshifted compared to the SiO2 6298

DOI: 10.1021/acs.jpcb.9b03026 J. Phys. Chem. B 2019, 123, 6290−6302

Article

The Journal of Physical Chemistry B

Figure 9. (a) Total structure factor of the NBS matrix and doped NBS matrix and comparison with the pure glass matrix of SiO2 and B2O3. (b) Partial structure factors of pairs (i) Si−Si, (ii) B−B, (iii) Si−B, (iv) Si−Na, (v) B−Na, (vi) Na−Na, (vii) Si-dopant, (viii) B-dopant, and (ix) dopant−dopant in NBS, doped NBS matrix, and comparison with the SiO2 and B2O3 matrix.

intensity at the lower frequency of 120 cm−1. Furthermore, the peaks for pure SiO2 and pure B2O3 at the higher frequency zone seem to be merged in the case of NBS and doped NBS matrix, which signifies that the peak at 1030 cm−1 in the case of NBS and doped NBS results from the collective motion of non-BO atoms connected to Si (Si−ONBO) and BO atoms connected to trigonal-coordinated B3. The results show a slight postshift at 120 cm−1 and a preshift at 1100 cm−1 for doped NBS compared to NBS, which might be associated with the weakening of the Si−OBO−Si network and increasing Si−ONBO with doping of the glass matrix. The reported peak positions can be found in good agreement with the studies reported by Kob and Ispas55 for the NBS glass matrix. Since the glassy state of a similar NBO shares very similar peak positions, only

negligible differences were noticed for doped glasses compared to the NBS glass matrix, however, a significant change can be noticed for NBS or doped NBS matrix compared to the pure glass matrix of SiO2 and B2O3.

4. SUMMARY AND CONCLUSIONS Glass plays an important role in the disposal of HLW for which NBS glass matrix has been preferred by most of the countries. The physical properties of the glass matrix depend on the nature and concentration of the alkali modifiers. MD simulations can provide a guideline to reduce the tediousness of experimental procedures for fabrication and characterization and also reduce the associated potential hazards. The present 6299

DOI: 10.1021/acs.jpcb.9b03026 J. Phys. Chem. B 2019, 123, 6290−6302

Article

The Journal of Physical Chemistry B

the NBS glass matrix was observed to be increased with the increase in alkali oxides for a concentration smaller than 40%. Most importantly, the results show the change in the trend of network connectivity nearly at R ≥ 1. Besides, Cs and Sr being the prime heat source in HLW were considered with a mass composition of 5 and 10% in the NBS matrix. The results show a nearly similar content of BO atoms for Cs and Sr doping when added in amount of 5% by mass. However, significant differences were noticed in BO atoms and structure factors of Cs-doped and Sr-doped NBS matrices with a doping concentration of 10%, which seems to be interconnected with the charge and ionic radius of the dopant. Furthermore, the characteristic VDOS of NBS and doped NBS were estimated, and the results showed that the peak intensity at ∼120 cm−1 was increased for the NBS matrix compared to the pure matrix of SiO2 and B2O3, and effects become more pronounced while doping. This was expected from the opening of network because of the presence of alkali ions in the NBS and doped NBS glass matrix, which makes the glass network more flexible. In addition, the peaks for NBS and doped NBS at the higher frequency zone were anticipated as the collective motion of non-BO atoms connected to Si (Si− ONBO) and BO atoms connected to trigonal coordinated B3. The VDOS of the doped matrix show a slight postshift at 120 cm−1 and a preshift at 1100 cm−1 compared to NBS, which reflects the weakening of the Si−OBO−Si network and increasing Si−ONBO with doping of the glass matrix. As a whole, the VDOS spectra show only negligible differences for doped glasses compared to the NBS glass. On the other hand, a significant change could be noticed for NBS or doped NBS compared to the pure glass matrix of SiO2 and B2O3. In summary, the present article demonstrates that the Buckingham potential-incorporated Coulomb interaction model with the selected force field parameters can be utilized to simulate the glasses with varied compositions. The results have been validated with the experimentally produced glasses in our laboratories. The results well capture the boron anomalies for varied concentration of network formers and network modifiers. Most importantly, the selected potential model and force field parameters not only capture well the structural distribution of glasses but also validate the characteristic vibration spectra. In addition, the validation of simulation data with the experimental observations assures the versatility of selected models for simulation of the vitrified glass matrix with various compositions. The results might be very helpful in the prediction of effects associated with the dopant’s nature and concentration. The presented MD results provide a microscopic understanding of the glass structure and phenomena associated with the change in the glass composition. The present work will hence have great academic as well as technological significance for further guidance in the selection of suitable glass composition for vitrification.

Figure 10. VDOS of NBS and doped NBS (NBS−10Cs and NBS− 10Sr) and comparison with the pure glass matrix of SiO2 and B2O3.

MD studies conducted with a wide range of glass compositions, pure matrix to multicomponent matrix, established that the Coulomb interaction-incorporated Buckingham potential model with the recommended force field parameters can be utilized to simulate glass of varied compositions. The used potential model and force field parameters were first validated by simulating the pure glass matrix of amorphous SiO2 and amorphous B2O3, the major network formers in the glass matrix. The density, structural morphology, and dynamics of the simulated SiO2 and B2O3 glass matrices were found in agreement with the experimental observations. Furthermore, the structure models of borosilicate were prepared by varying the SiO2 concentration from 10 to 70%. Results showed that the network connection became stronger with higher concentration of SiO2, as observed from an increasing BO atoms and BO3 to BO4− conversion. Furthermore, the effect of sodium addition to the pure glass matrix of SiO2 and B2O3 was observed by varying the sodium concentration from 10 to 60% in sodium silicate and sodium borate. The results display a rise in NBO atoms with the increase in the Na2O concentration. Interesting results were observed for sodium addition in the B2O3 matrix, where it was found that the density of glass and BO3 to BO4− conversion increases with alkali addition till 40% and thereafter starts reducing with furthermore addition of Na2O. In particular, the initially (up to 40%) added sodium ions were seen to reside within the cavities formed by boroxyl rings, where they act as a charge compensator to BO4− tetrahedrals and so convert BO3 → BO4−. Once the occupancy of Na ions in the interstitial positions is fulfilled (that appears at 40% Na2O concentration), no more BO4− conversion happens, and the added Na2O now acts as a network modifier, leading to a reduced density and BO4− units in the sodiumborate glass matrix. In sequence, the NBS matrix was simulated with a composition of 21SiO2.6B2O3.7Na2O. The density and glass transition temperature of the simulated matrix was found in good agreement with the experimentally synthesized glass matrix of the same composition in our laboratories. In addition, two set of simulations, NBS15 and NBS20, were performed with varied R(Na2O/B2O3), K(SiO2/B2O3), and M[SiO2/(SiO2 + B2O3)] and fixed Na2O concentration, that is, 15 and 20% by mole, respectively. The network connectivity was observed to become stronger with the increase in the SiO2 concentration in the NBS matrix. In line with the observations noticed with sodium silicate and sodium borate, the density of



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.9b03026. Details for the experimental synthesis of NBS glass and doped NBS glasses; analysis details for the density of glass from experiments; analysis details for the glass transition temperature (Tg) from experiments; density (g/cm3) of the NBS matrix and the doped NBS matrix; 6300

DOI: 10.1021/acs.jpcb.9b03026 J. Phys. Chem. B 2019, 123, 6290−6302

Article

The Journal of Physical Chemistry B



(19) Turner, G. L.; Smith, K. A.; Kirkpatrick, R. J.; Oldfield, E. Boron-11 Nuclear Magnetic Resonance Spectroscopic Study of Borate and Borosilicate Minerals and a Borosilicate Glass. J. Magn. Reson. 1986, 67, 544−550. (20) Gruenhut, S.; MacFarlane, D. R. Molecular Dynamics Simulation of Heavy Metal Fluoride Glasses: Comparison of Buckingham and BHM Potentials. J. Non-Cryst. Solids 1995, 184, 356−362. (21) Weber, W. J.; Ewing, R. C.; Angell, C. A.; Arnold, G. W.; Cormack, A. N.; Delaye, J. M.; Griscom, D. L.; Hobbs, L. W.; Navrotsky, A.; Price, D. L.; Stoneham, A. M.; Weinberg, M. C. Radiation Effects in Glasses Used for Immobilization of High-level Waste and Plutonium Disposition. J. Mater. Res. 1997, 12, 1948− 1978. (22) Weber, W. J. Radiation and Thermal Ageing of Nuclear Waste Glass. Procedia Mater. Sci. 2014, 7, 237−246. (23) Hayakawa, S.; Ohtsuki, C.; Matsumoto, S.; Osaka, A.; Miura, Y. Molecular Dynamic Simulation of Heterogeneity and Chemical States of Fluorine in Amorphous Alkaline Earth Silicate Systems. Comput. Mater. Sci. 1998, 9, 337−342. (24) Osipov, A. A.; Eremyashev, V. E.; Mazur, A. S.; Tolstoi, P. M.; Osipova, L. M. Structure of Cesium−Borosilicate Glasses According to NMR Spectroscopy. Glass Phys. Chem. 2017, 43, 287−293. (25) Minami, T.; Tokuda, Y.; Masai, H.; Ueda, Y.; Ono, Y.; Fujimura, S.; Yoko, T. Structural Analysis of Alkali Cations in Mixed Alkali Silicate Glasses by 23Na and 133Cs MAS NMR. J. Asian Ceram. Soc. 2014, 2, 333−338. (26) Marples, J. A. C. The Preparation, Properties, and Disposal of Vitrified High Level Waste From Nuclear Fuel Reprocessing. Glass Technol.-Part A 1988, 29, 230−247. (27) Connelly, A. J.; Travis, K. P.; Hand, R. J.; Hyatt, N. C.; Maddrell, E. Composition−Structure Relationships in Simplified Nuclear Waste Glasses: 1. Mixed Alkali Borosilicate Glasses. J. Am. Ceram. Soc. 2010, 94, 151. (28) Kieu, L.-H.; Delaye, J.-M.; Cormier, L.; Stolz, C. Development of Empirical Potentials for Sodium Borosilicate Glass Systems. J. NonCryst. Solids 2011, 357, 3313−3321. (29) Guillot, B.; Sator, N. A Computer Simulation Study of Natural Silicate Melts. Part II: High Pressure Properties. Geochim. Cosmochim. Acta 2007, 71, 4538−4556. (30) Suck, J. B.; Chieux, P.; Raoux, D.; Riekel, C. Methods in The Determination of Partial Structure Factors of Disordered Matter by Neutron and Anomalous X-Ray Diffraction. Proceedings of the ILL/ ESRF Workshop: Grenoble, France, 1992. (31) Du, J. Challenges in Molecular Dynamics Simulations of Multicomponent Oxide Glasses. Molecular Dynamics Simulations of Disordered Materials; Springer Series in Materials Science: Switzerland, 2015; Vol. 215, pp 157−180. (32) Pedone, A.; Malavasi, G.; Menziani, M. C.; Cormack, A. N.; Segre, U. A New Self-Consistent Empirical Interatomic Potential Model for Oxides, Silicates, and Silica-Based Glasses. J. Phys. Chem. B 2006, 110, 11780−11795. (33) Kubicki, J. D.; Lasaga, A. C. Molecular Dynamics Simulations of SiO2 Melt and Glass: Ionic and Covalent Models. Am. Mineral. 1988, 73, 941−955. (34) Takada, A.; Catlow, C. R. A.; Price, G. D. Computer Modelling of B2O3: part II. Molecular Dynamics Simulations of Vitreous Structures. J. Phys.: Condens. Matter 1995, 7, 8693−8722. (35) Tilocca, A. Cooling Rate and Size Effects on the MediumRange Structure of Multicomponent Oxide Glasses Simulated by Molecular Dynamics. J. Chem. Phys. 2013, 139, 114501. (36) Tilocca, A.; Leeuw, N. H. d.; Cormack, A. N. Shell-Model Molecular Dynamics Calculations of Modified Silicate Glasses. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 104209. (37) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (38) van Beest, B. W. H.; Kramer, G. J.; Santen, R. A. v. Force Fields for Silicas and Aluminophosphates Based on Ab Initio Calculations. Phys. Rev. Lett. 1990, 64, 1955−1958.

glass transition temperature (Tg) of the NBS matrix and the doped NBS matrix; RDF, coordination distribution, and angle distribution for amorphous SiO2 and amorphous B2O3; structure of sodium silicate with 10% Na2O and 60% Na2O in the sodium-silicate glass matrix; and snapshot of the doped NBS glass matrix (PDF)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +91-22-25591992. ORCID

Sk. Musharaf Ali: 0000-0003-0457-0580 Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS The computer division of BARC is acknowledged for providing ANUPAM Supercomputing facility. The author(s) received no specific funding for this work.



REFERENCES

(1) Al-Azzawi, S.; Yasin, F. Optical and Chemical Properties of Glass with Special Reference to Building. Sol. Wind Technol. 1985, 2, 41− 51. (2) Gliemeroth, G. Optical Properties of Optical Glass. J. Non-Cryst. Solids 1982, 47, 57−68. (3) Kaushik, C. P. Indian Program for Vitrification of High Level Radioactive Liquid Waste. Procedia Mater. Sci. 2014, 7, 16−22. (4) Plodinec, M. J. Borosilicate Glasses for Nuclear Waste Immobilization. Glass Technol-part A 2000, 41, 186−192. (5) Harrison, M. T. Vitrification of High Level Waste in the UK. Procedia Mater. Sci. 2014, 7, 10−15. (6) Wicks, G. G. Immobilization of Hazardous and Radioactive Waste into Glass Structure Senior Advisory Scientist: USA, 1997. (7) Grover, J. R. The Solidification of High-Level Radioactive Wastes Final Report. IAEA Bulletin-20 (NUREG/CR--0895) 1979, 261. (8) Manaktala, H. K. An Assessment of Borosilicate Glass as a Highlevel Waste Form Nuclear Regulatory Commission Contract NRC-02-88005 (CNWRA 92-017): San Antonio, Texas, 1992. (9) Evaluation of Solidified High-Level Waste Forms. IAEA-TECDOC239 VIENNA, 1981. (10) Delaye, J.-M.; Ghaleb, D. Molecular Dynamics Simulation of SiO2 + B2O3 + Na2O + ZrO2 Glass. J. Non-Cryst. Solids 1996, 195, 239−248. (11) Du, J.; Cormack, A. N. The Medium Range Structure of Sodium Silicate Glasses: A Molecular Dynamics Simulation. J. NonCryst. Solids 2004, 349, 66−79. (12) 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, 3144−3154. (13) Vegiri, A.; Varsamis, C.-P. E.; Kamitsos, E. I. Composition and Temperature Dependence of Cesium-borate Glasses by Molecular Dynamics. J. Chem. Phys. 2005, 123, 014508. (14) Cormack, A. N.; Cao, Y. Molecular Dynamics Simulation of Silicate Glasses. Mol. Eng. 1996, 6, 183−227. (15) Gou, F.; Greaves, G. N.; Smith, W.; Winter, R. Molecular Dynamics Simulation of Sodium Borosilicate Glasses. J. Non-Cryst. Solids 2001, 293-295, 539−546. (16) Delaye, J. M.; Ghaleb, D. Molecular Dynamics Simulation of a Nuclear Waste Glass Matrix. Mater. Sci. Eng. B 1996, 37, 232−236. (17) Ehrt, D. Structure, Properties and Applications of Borate Glasses. Glass Technol.-Part A 2000, 41, 182−185. (18) Krogh-Moe, J. The Structure of Vitreous and Liquid Boron Oxide. J. Non-Cryst. Solids 1969, 1, 269−284. 6301

DOI: 10.1021/acs.jpcb.9b03026 J. Phys. Chem. B 2019, 123, 6290−6302

Article

The Journal of Physical Chemistry B

of Polymer-Functionalized Neutral Macrocyclic Host. ACS Appl. Mater. Interfaces 2018, 10, 20968−20982. (64) Avramov, I.; Vassilev, T.; Penkov, I. The Glass Transition Temperature of Silicate and Borate Glasses. J. Non-Cryst. Solids 2005, 351, 472−476. (65) Fu, X.; Wang, A.; Krawczynski, M. J. Characterizing Silicate Glasses with Vibrational Spectroscopy In 47th Lunar and Planetary Science Conference, 2016; p 2470.

(39) Malavasi, G.; Menziani, M. C.; Pedone, A.; Segre, U. Void Size Distribution in MD-Modelled Silica Glass Structures. J. Non-Cryst. Solids 2006, 352, 285−296. (40) Wang, M.; Anoop Krishnan, N. M.; Wang, B.; Smedskjaer, M. M.; Mauro, J. C.; Bauchy, M. A New Transferable Interatomic Potential for Molecular Dynamics Simulations of Borosilicate Glasses. J. Non-Cryst. Solids 2018, 498, 294−304. (41) Luty, B. A.; van Gunsteren, W. F. Calculating Electrostatic Interactions Using the Particle−Particle Particle−Mesh Method with Nonperiodic Long-Range Interactions. J. Phys. Chem. 1996, 100, 2581−2587. (42) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. Packmol: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30, 2157−2164. (43) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press, Oxford University, 2017. (44) Humphrey, W.; Dalke, A.; Schulten, K. VmdVisual Molecular Dynamics. J. Mol. Graph. 1996, 14, 33−38. (45) Machacek, J.; Gedeon, O. Q-species in Alkali-disilicate Glasses. Ceramics 2003, 47, 45−49. (46) Xiang, Y.; Du, J. Effect of Strontium Substitution on the Structure of 45S5 Bioglasses. Chem. Mater. 2011, 23, 2703−2717. (47) Pedone, A.; Malavasi, G.; Menziani, M. C.; Segre, U.; Cormack, A. N. Molecular Dynamics Studies of Stress-Strain Behavior of Silica Glass under a Tensile Load. Chem. Mater. 2008, 20, 4356−4366. (48) Valle, R. G. D.; Venuti, E. High-pressure Densification of Silica Glass: A Molecular-dynamics Simulation. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 3809−3816. (49) Allan, D.; McQuarrie Statistical Thermodynamics; USB, 1984. (50) Pacaud, F.; Delaye, J.-M.; Charpentier, T.; Cormier, L.; Salanne, M. Structural Study of Na2O−B2O3−SiO2 Glasses from Molecular Simulations Using a Polarizable Force Field. J. Chem. Phys. 2017, 147, 161711. (51) Ohmura, S.; Shimojo, F. Mechanism of Atomic Diffusion in Liquid B2O3: An Ab Initio Molecular Dynamics Study. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 224206. (52) Fábián, M.; Sváb, E.; Mészáros, G.; Révay, Z.; Proffen, T.; Veress, E. Network Structure of Multi-component Sodium Borosilicate Glasses by Neutron Diffraction. J. Non-Cryst. Solids 2007, 353, 2084−2089. (53) Pedesseau, L.; Ispas, S.; Kob, W. First Principles Study of a Sodium Borosilicate Glass-former I: The Liquid State. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 91, 134201. (54) Balucani, U.; Zoppi, M. Dynamics of the Liquid State; Instituto di Elettronica Quantistica CNR: Firenze: Italy, 1994. (55) Kob, W.; Ispas, S. First-principles Simulations of Glass-formers, 2016. arXiv:1604.07959. (56) Tanguy, A. Vibration Modes and Characteristic Length Scales in Amorphous Materials. JOM 2015, 67, 1832−1839. (57) Verhoef, A. H.; den Hartog, H. W. A Molecular Dynamics Study of B2O3 Glass Using Different Interaction Potentials. J. NonCryst. Solids 1992, 146, 267−278. (58) Doweidar, H. Considerations on the Structure and Physical Properties of B2O3−SiO2 and GeO2−SiO2 glasses. J. Non-Cryst. Solids 2011, 357, 1665−1670. (59) Inoue, H.; Masuno, A.; Watanabe, Y. Modeling of the Structure of Sodium Borosilicate Glasses Using Pair Potentials. J. Phys. Chem. B 2012, 116, 12325−12331. (60) Shelby, J. E. Introduction to Glass Science and Technology, 2nd ed.; Royal Society of Chemistry: Cambridge, 2005. (61) Han, Y. H.; Kreidl, N. J.; Day, D. E. Alkali Diffusion and Electrical Conductivity in Sodium Borate Glasses. J. Non-Cryst. Solids 1979, 30, 241−252. (62) Barlet, M.; Kerrache, A.; Delaye, J.-M.; Rountree, C. L. SiO2− Na2O−B2O3 Density: A Comparison of Experiments, Simulations, and Theory. J. Non-Cryst. Solids 2013, 382, 32−44. (63) Sahu, P.; Ali, S. M.; Shenoy, K. T.; Mohan, S. Structure, Dynamics, and Adsorption of Charged Guest within the Nanocavity 6302

DOI: 10.1021/acs.jpcb.9b03026 J. Phys. Chem. B 2019, 123, 6290−6302