Molecular dynamics of vitreous silica — Variations in potentials and simulation regimes

Journal of Non-Crystalline Solids 426 (2015) 103–109

Contents lists available at ScienceDirect

Journal of Non-Crystalline Solids journal homepage: www.elsevier.com/ locate/ jnoncrysol

Molecular dynamics of vitreous silica — Variations in potentials and simulation regimes Ondrej Gedeon Department of Glass and Ceramics, University of Chemical Technology, Technicka 5, CZ-166 28 Prague, Czech Republic

a r t i c l e

i n f o

Article history: Received 30 April 2015 Received in revised form 23 June 2015 Accepted 2 July 2015 Available online 14 July 2015 Keywords: Molecular dynamics; Silica glass; Structure

a b s t r a c t Two potential sets, and three thermodynamic regimes were utilised in MD simulations to prepare the corresponding structures of vitreous silica. Standard structure descriptors as RDF's, angle distributions, and primitive ring distributions showed subtle differences among glasses. Structural defects as overcoordinations and dangling bonds revealed higher sensitivity to both used potential sets and TD conditions. Overcoordination is largely driven by the density while the amount of dangling bonds is given by the type of the potential. Connectivity that balances under- and over-coordination is determined rather by TD conditions than type of potential used. Ring distribution as well as radial distribution function of rings shows robustness to alternating conditions of simulation. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Experimental knowledge of glass structure is still insufficient despite huge progress in the last decades. The key problem is i) available techniques are mostly constrained to description of nearest or second nearest (with a decreasing precision) neighbours, and ii) there is(are) no universal descriptor(s) enabling uniquely and unambiguously describe the structure. Direct information about glass structure is offered by NMR (Nuclear Magnetic Resonance), IRS (Infra Red Spectroscopy), RS (Raman Spectroscopy), EXAFS (Extended X-Ray Absorption Fine Structure), ND (Neutron Diffraction), XRD (X-Ray Diffraction), and some other techniques while indirect knowledge about structure is supplied by various atomic-based models. The second crux is rather the theoretical one. Translational periodicity of crystal solid has opened rich realm of consequences of material properties coming out from this symmetry. The enormous number of data needed to describe the position of every atom shrunk to knowledge of the positions of atoms in the elemental cell and three translational vectors. Similar structural homomorphism is lacking for non-crystalline solids. Therefore a few “standard structural descriptors” are used; among them radial distribution function, RDF, or more detailed partial pair distribution function, PP RDF, are mostly used. Nevertheless, RDFs are one-dimensional projection of the spatial distribution of atoms and therefore cannot adequately catch structural anisotropy of non-crystalline solid. On the other hand RDF is the most common experimental result; it yields coordination numbers and bond lengths — information of uppermost significance for chemistry. For silica glasses Qn-units, tetrahedron formed by central

E-mail address: [email protected].

http://dx.doi.org/10.1016/j.jnoncrysol.2015.07.006 0022-3093/© 2015 Elsevier B.V. All rights reserved.

silicon with four vertices determined by oxygen atoms, with n giving the number of bridging oxygen, BO, are natural structural units. Its distribution (e.g. RS, NMR) or connectivity (NMR) is experimentally attainable and provides complemental information. Nevertheless, in glass community it is believed that glass essence is coded in medium range order (MRO) [1]. Rings or loops (various terminologies and definitions are used) are candidates of proper descriptors of MRO. They provide not geometrical but topological picture of glass structure. Primitive rings [2] or minimal ring [3], introduced as a loop that cannot decomposed into two shorter loops or as a loop with no “short cut” are mostly used. Other definitions can be found, e.g. strong rings, very strong rings [4], and basic rings [5]. Molecular dynamics (MD), atomistic simulation method, mimics the behaviour of atomic motion under given thermodynamic restrictions. Preparation of glass in MD starts at very high temperature to enable perfect mixing of atoms (ions) and to reach thermodynamic equilibrium. Gradual decrease of temperature (step by step or continuous) causes solidification of the studied system with a structure close to the experimental glass. Using classical physics approach, interatomic (mostly two- and sometimes more-particle) potentials enter into Newton's equations of motions as given (e.g. from ab-initio calculations and/or from fitting of experimental data). Treatment of thousands of particles is expected to be able to eliminate periodic boundary conditions conventionally imposed upon MD calculations [6]. Time window, on the other hand, accessible for performing the computational cooling is of a few orders shorter than in reality resulting in a high transformation temperature obtained for the simulated glass. Even though it is questionable if it is a real glass transition what is observed during simulation [7] the final computer product is a solid with its structure strongly reminding real glass. The direct comparison of simulated and real

104

O. Gedeon / Journal of Non-Crystalline Solids 426 (2015) 103–109

glass is even impaired by the fact that glass properties depend on cooling rate for both computer [8,9] and real glass [10]. Hence, glass prepared by MD should be tagged by the cooling regime used. The aim of the paper is to explore two different empirical potentials utilised for MD simulations of vitreous silica and to study the influence of various thermodynamical conditions on both the standard structural quantities as RDF, angle distribution, and ring distribution and not commonly presented structural defects. 2. MD simulation Vitreous silica was simulated by classical MD using two sets of pairwise potentials, BKS [11] and CHIK [12] (designation of both potentials comes from the first letters of author's family names), both having the same functional form
 C zi z j e2 1 Ei j ¼ Ai j exp −Bi j r − 6 þ 4πε0 r r

ð1Þ

atom within the sphere of radius D, it is considered as non-bridging (NBO), if two silicon atoms are within the sphere it is taken as bridging (BO), if more than two Si are inside the sphere, it is designated as overcoordinated (OO). At high temperatures free oxygen (FO) can be also identified as those with no Si within the sphere. Many of the bonds seem to be artificial and do not correspond to the real chemical bonds but they are rather a result of the temporary spatial packing. Therefore the following algorithm was used for the decomposition of 3-coordinated oxygen atoms. If OO is connected to at least one 5- or more- coordinated Si, the longest bond among those leading to the highest coordinated Si is removed. Remove of bond leading to the at least 5-coordinated Si (not to the 4- or less coordinated Si) is required so as the removal of the bond could not lead to the formation of the undercoordinated Si and therefore so as one type of structural defect would have not been replaced by another one. Below is the scheme, showing the result of the algorithm applied for OO with upperscripts m ≤ l ≤ k displaying the corresponding coordination. Sil

Sil

parameters of which are summarized in Table 1. The first two terms in Eq. (1) describe short-range interaction while the last one corresponds to the electrostatic long-range Coulomb potential. The Coulombic part was evaluated by summing over atoms up to 12 Å in the real space and beyond this border it was treated by Ewald's method [13] in dual space. All calculations were performed in the simulation cube comprising 7,200 atoms (2400 Si and 4800 O) with periodic boundaries. DL_Poly simulation code was utilised [14]. For each potential set the system was constrained to three different thermodynamical boundaries: one run was performed under a constant pressure (NpE ensemble) and two runs under constant volume (NVE ensembles). NVE ensembles were adjusted to i) the experimental density of vitreous silica (2.20 g cm−3), and to ii) the density of quartz (2.65 g cm−3) [15]. The system with the latter density was prepared with idea that, first, BKS potential was adjusted to quartz data and, second, the previous experience has shown a significantly smaller amount of structural defects in potassium silicate glass at quartz density [16]. Glasses prepared under the abovementioned conditions will be hereafter designated with e (experimental density), q (quartz density), and p (ambient pressure). All simulations started at 5000 K by mixing and equilibrating of the system during 25 ps. The cooling procedure consisted of 100° consecutive temperature steps; each step comprised at first the numeric control of the adjusted temperature (2500 time steps) followed by simulation obeying thermodynamic requirements (either constant volume or pressure) for the next 7500 time steps. The leap-frog algorithm was used for the numerical integration of Newton's equations of motion with an integration time step of 1 fs. Hence, the entire cooling procedure (down to 100 K) corresponds to the effective cooling rate of 1013 Ks−1. The snapshots of the simulated structures were analysed (150 structures at each temperature, i.e. 0.15 ps was the time interval between two snapshots and 7.5 ps structure evolution was recorded and analysed at each temperature) and the average quantities obtained from these snapshots are further presented. MD did not establish chemical bonds but they had to be determined indirectly by the distance criterion as follows: the bond between oxygen and silicon is established if |rSi − rO| ≤ D (D = 2.0 Å) was used. If oxygen atom has one silicon

O3

O2

Sik

Sim

Sik-1

Sim

Of course, this algorithm eliminates not only the number of OO but proportionally decreases the number of overcoordinated Si. Once bonds among atoms are established for simulated glass, distribution of Q-units, evaluation of structural defects, and primitive ring distributions were evaluated. To study spatial correlation of the rings, ring radial distribution function, RRDF, was introduced

RDFRðr Þ ¼

 1 1X

δ r− r i −r j 4πr 2 ρ N i; j

ð2Þ

where ρ is volume density of rings, N is the number of rings, and summation runs over positions of rings defined as

r¼

M X

rk

ð3Þ

k¼1

where k runs over positions of (silicon) atoms forming the ring.

Table 1 Parameters of BKS and CHIK potential, see Eq. (1). BKS

A B C zSi zO

CHIK

Si–Si

Si–O

O–O

Si–Si

Si–O

O–O

0 0 0 2.4 −1.2

18,003.76 4.87 133.54

1388.77 2.76 175

3150.46 2.851 626.75 1.91 −0.955

27,029.42 5.1586 148.1

659.59 2.59 26.84 Fig. 1. Pair potential energy of Si–O for BKS (dark red) and CHIK (green) potential sets.

O. Gedeon / Journal of Non-Crystalline Solids 426 (2015) 103–109

3. Results 3.1. Si–O pair potential energy Fig. 1 presents pair potential energy of the Si–O given by BKS and CHIK potential prescriptions. The charges associated with the simulated atoms are smaller for CHIK potential so that CHIK potential can be comprehended as a more covalent in comparison with BKS one. Note that the CHIK parameters also include Si–Si short range repulsive contribution that may compensate the smaller repulsive part coming out from the lower charges located on Si. Both potentials belong to the “strong attractive” ones [17]. The position of the minimum for BKS Si–O potential is 1.38 Å and that for CHIK potential is placed a little higher (1.45 Å). The difference between the minimum of the pair Si–O potential and an average equilibrium Si–O position in simulations (1.60–1.62 Å) shows significant contributions from other atoms. BKS pair potential is also deeper, i.e. higher energy is needed for a bond breakage with straightforward consequences for some properties, e.g. diffusion [17].

3.2. RDFs Fig. 2 documents influence of the both potential sets and thermodynamical conditions on the short-range structure. All PP RDF(Si–O)s are pretty close to each other and only additional zooming is able to unveil the tiny differences. The most probable Si–O distances for all regimes fall into the interval 1.60–1.63 Å set by neutron and X-Ray diffraction experiments [18], see Table 3. The Si–O distances seem to be determined rather by the potential used than by TD conditions/final pressures/glass densities. Not surprisingly, higher density glasses (affix q), see Table 2, shift the first Si–O peak to the lower values; and BKS generates glass with a shorter average bond length of Si–O. An interesting feature reveals the comparison of p-glasses (simulated under constant pressure) with V-glasses (simulated under constant volume).

105

Table 2 Average pressures and densities of the simulated glass at 300 K for CHIK and BKS potential sets. Densities of q- and e-glasses were set for isochoric simulations; densities of p-glasses were obtained as a result of isobaric simulations.

p [GPa] ρ [g cm−3]

CHIK-p

CHIK-q

CHIK-e

BKS-p

BKS-q

BKS-e

– 2.33

3.5 2.65

−1.1 2.20

– 2.54

0.3 2.65

−2.1 2.20

While Si–O average bond length is the same for p- and q-glasses for BKS (the simulations finished at the fairly close densities) in case of CHIK resemblance can be found between p- and e-glasses. Similar features can be also observed for the second peak of Si–O; here not only the positions but also the widths of the peaks reflect higher densities of the q-glasses. 3.3. Angle distribution Fig. 3 shows the angle distribution of O–Si–O and therefore it characterizes the shape of the silicon tetrahedron consisted of one silicon surrounded with four oxygens. The most probable angles, see Table 4, are in concordance with experimental values (108.6 ± 4.2)° [19]. The effect of the higher density (CHIK-q, BKS-q, and BKS-p), see Table 2, is reflected in the small shift of the average angle to the lower values. While CHIK-p and CHIK-e distributions are practically identical, the differences among all BKS glasses are well discernible. The angle distribution of Si–O–Si angles, Fig. 4 is expected to be strongly correlated with the first peak of PP RDF(Si–Si). In silica network Si–O–Si angle possesses the largest freedom and therefore it should be the most sensitive to the changes in density and/or TD conditions of cooling. Comparing Si–O–Si distributions with O–Si–O ones, it is seen that Si–O–Si curves are more separated even though larger statistical noise is visible since the number of O–Si–O angles are three times higher than the amount of Si–O–Si's in ideally connected silica glass. Data in the table confirm that CHIK-q mode is shifted from the values of the pair

Fig. 2. PP RDF (Si–O) for BKS (dark red) and CHIK potentials (green) at 300 K. Solid, dotted, and broken lines correspond to e-, q-, and p-glasses, respectively. Left inset figure is a zoom of the first peak and the right inset magnifies the second peak.

106

O. Gedeon / Journal of Non-Crystalline Solids 426 (2015) 103–109

Table 3 Si–O bond lengths determined as modes from PP RDF (Si–O).

Si–O [Å]

CHIK-e

CHIK-p

CHIK-q

BKS-e

BKS-p

BKS-q

1.63

1.63

1.62

1.62

1.61

1.61

(CHIK-e, CHIK-p) to the lower value and BKS-e mode is shifted from the pair (BKS-p, BKS-q) to the higher value. A review of experimental data shows large scattering; modes are ranging from 141° to 160° and FWHM from 18° to 46° [20]. Distributions for CHIK potential follow the increasing density, i.e. the higher density yields the lower average Si–O–Si angle. BKS potential reveals closeness of p- and-q glasses again. 3.4. Structural defects The chemical bonds must be determined in MD indirectly. The distance criterion with cut-off equal to 2 Å was used. However, change of the cut-off may significantly alter the connectivity of the network related with structural defects. Fig. 5 demonstrates the measure of robustness of the determined structure at ambient temperature in relation with the variation of the cut-off distance. In ideal silica glass all oxygen atoms are bridging (two bonds with two silicon atoms) and each silicon atom is bounded with four oxygen atoms and therefore the ideal connectivity generates Q4-units only. Nevertheless, both in real glass as well as in MD glass some amounts of bond defects are always present. The defects influence the distribution of Qn-units as well as network topology (connectivity, ring distribution). Although it is generally believed that the higher amount of bond defects found in MD glass in comparison to real glass is related to the higher cooling rate, little attention has been given to TD conditions under which MD glass was prepared and/or to potential sets employed. Table 5 clearly documents that the number of defect increase with the increasing density but it is also seen that CHIK potential set generates slightly higher amounts of structural defects under the same TD conditions. It is worth noting that due to the algorithm used for removing of OO, the number of Q6 units were also reduced as the high coordinated silicon may be coupled with overcoordinated oxygen. Alternative insight into the bond defect distributions gives their distinguishing into the particular oxygen defects: free oxygen (FO), dangling oxygen (DO), and overcoordinated oxygen (OO); overcoordinated silicon (OSi) was added for sake of completeness. The amounts of FO defects were zero below at about 2200 K and negligible even at high temperatures for all glasses. Fig. 6 confirms that CHIK potential generates higher amount of defects in comparison with BKS potential. Overcoordination is clearly correlated with density; the

Fig. 3. Distribution of O–Si–O angles for glasses at 300 K as obtained from simulations using BKS (dark red) and CHIK potentials (green). Solid, dotted, and broken lines correspond to e-, q-, and p-glasses, respectively.

Table 4 Modes, M, and full widths at half maximum, FWHM, of distributions of O–Si–O and Si–O–Si angles.

M(O–Si–O) FWHM(O–Si–O) M(Si–O–Si) FWHM(Si–O–Si)

CHIK-e

CHIK-p

CHIK-q

BKS-e

BKS-p

BKS-q

108.3 13.9 145 37

108.2 14.0 146 37

107.5 19.7 142 39

108.2 16.3 152 35

107.3 15.0 147 37

107.4 14.4 148 39

lowest amount of both OO and OSi is for the lowest density (e-glasses) while the highest abundancy is for the highest density (q-glasses). OSi is the dominant structural defect at ambient temperature; CHIK generates systematically more OSi than BKS under the same TD conditions. Both OO and OSi increase while DO decreases connectivity of glass network. Interestingly, BKS structure generates less overcoordination (both for Si and O) and very small amount of dangling oxygen atoms. The above discussed relations among structural defects as OSi, OO, and OSi, Qn distribution, and network connectivity, C, can be mathematically grounded. Number of bonds, NB, leading from oxygen is the same as the number of bonds leading from silicon and can be evaluated by means of coordination NB ¼

X

iN Si ðiÞ ¼

i

X

iN O ðiÞ

ð4Þ

i

where i is the coordination, and NSi(i) or NO(i) are the numbers of icoordinated Si or O in glass. On the other hand Qn relates directly to the connectivity. The presence of one Qn means the presence of one silicon atom with n connections/bonds established via different oxygen atoms with other Si atoms (not necessarily all different if two silicon tetrahedra share an edge or a face). Number of bonds can be also expressed by means of Qn by the following equation (here we limit oxygen atoms to up to three-coordinated ones for sake of simplicity) NB ¼

X

nQ n þ NðDOÞ−3NðOOÞ:

ð5Þ

n

Introducing connectivity, C, as

C¼

1 X nQ n NSi n

ð6Þ

Fig. 4. Distribution of Si–O–Si angles for glasses at 300 K as obtained from simulations using BKS (dark red) and CHIK potentials (green). Solid, dotted, and broken lines correspond to e-, q-, and p-glasses, respectively.

O. Gedeon / Journal of Non-Crystalline Solids 426 (2015) 103–109

Fig. 5. Distribution of defects for BKS p-glass. D is the cut-off distance used in the distance criterion for the establishing of Si–O bond. Similarly weak dependence was found for other glasses. Presented data are normalised to the number of Si.

107

Fig. 6. Distribution of defects for CHIK (green) and BKS (dark red) potentials at 300 K. Data are normalised to the number of Si.

where NSi is the number of silicon atoms, one gets NB ¼ N Si C þ NðDOÞ−3NðOOÞ:

ð7Þ

Assuming 3–6 coordination of Si and 1–3 coordination of O, we have the following equation coupling coordination with defects OO, OSi, and DO with connectivity C 3NSi ð3Þ þ 4NSi ð4Þ þ 5NSi ð5Þ þ 6NSi ð6Þ ¼ ¼ NðDOÞ þ 2NO ð2Þ þ 3NðOOÞ ¼ NSi C þ NðDOÞ−3NðOOÞ:

ð8Þ

For example, the second equality leads to a simple relation for connectivity expressed via oxygen coordination C¼

2NO ð2Þ þ 6NðOOÞ : N Si

ð9Þ

Expressing the number of two-coordinated oxygen as NO ð2Þ ¼ 2N Si −N ðOOÞ−NðDOÞ

4. Discussion ð10Þ

we have finally an equation relating connectivity with oxygen defects C ¼4þ4

N ðOOÞ NðDOÞ −2 : NSi NSi

ð11Þ

Connectivity calculated by Eq. (11) is presented in Table 6. Despite relatively large differences in amounts of particular defects connectivity is grouping glasses into pairs (CHIK-e, BKS-e), (CHIK-p, BKS-p), and (CHIK-q, BKS-q). The lowest (but the closest to the ideal one) is obtained for e-glasses while the highest is revealed by q-glasses. The connectivity seems to be determined more by density/TD conditions than the type of potential used. 3.5. Rings distribution There is no unique definition of the ring [3]. Nevertheless, primitive rings, intuitively introduced as rings with no shortcut, are mostly used Table 5 Distributions of Q-units (in %) of silica glasses obtained under various conditions at 300 K. Values are averages from 150 snapshots. Negligible amount (less than one unit) of Q7-units was found for BKS-q glass.

Q3 Q4 Q5 Q6

as a descriptor for glass topology. Fig. 7 compares the distributions of primitive rings. The differences between the ring topologies are subtle without any straightforward conclusions. The spatial distribution of rings, RRDF, see Fig. 8, shows the strong spatial correlation up to about 8 Å with at least four well distinguished peaks. All RRDF's are, maybe also due to the low statistics, pretty close and reasonably indistinguishable. Hence, medium range order represented by the ring distribution seems to be very robust to changes in TD conditions or potentials. The first peak at 1.6 Å seems to be unreasonable positioned at too small value. However, it arises mostly from ring clusters, i.e. rings sharing some common part. The existence of clusters comes out from the definition of the primitive ring and is related with high connectivity of glass network. As a result primitive rings in silica glass network are not disjunctive and a portion of them are mutually overlapping yielding the first peak at the short distance.

CHIK-e

CHIK-p

CHIK-q

BKS-e

BKS-p

BKS-q

0.9 96.0 3.1 0.0

0.9 95.1 4.0 0.0

0.5 93.0 6.1 0.3

0.2 97.9 1.9 0.0

0.1 96.8 3.0 0.1

0.0 95.5 4.0 0.5

MD offers a detailed picture of atomic motion based just on atomic potentials and the simple laws of Newton mechanics. The usage of statistical physics should provide a correct bridge between full information about the system provided by MD and macroscopic view represented by thermodynamics. Numerical cooling mimics glass quenching; the obtained structure and/or properties of the simulated system can be directly compared with experiments. However, there are strong limitations for such comparison from both theoretical and experimental sides. RDFs, attainable from experiments are just one-dimensional projections of the more complicated 3D space arrangements and the experimental information about atomic arrangement in glass beyond the second neighbour lacks, even if there is any, the required precision. On the other hand, potentials used in MD are empirical, cooling rates are extremely high, and connectivity of the network (bond establishment) is set indirectly, often just on base of interatomic distances. Two potential sets are used in the paper; both reveal the same functional form but differ by the values of the particular parameters. CHIK potential set utilises lower charges located on the atoms decreasing so the ionic contribution. BKS potential of Si–O pair is significantly deeper and wider and its minimum is shifted to a slightly lower value. Nevertheless, the addition of contribution of potentials from other atoms Table 6 Network connectivity of the simulated glasses calculated by the Eq. (11). Table 3. Si–O bond lengths determined as modes from PP RDF (Si–O).

Connectivity

CHIK-e

CHIK-p

CHIK-q

BKS-e

BKS-p

BKS-q

4.011

4.016

4.026

4.009

4.016

4.025

108

O. Gedeon / Journal of Non-Crystalline Solids 426 (2015) 103–109

Fig. 7. Relative distribution of primitive rings in silica glass for CHIK (green) and BKS (dark red) potentials. The distribution was evaluated at 300 K.

results in the practically identical RDF's. Fig. 2 well documents PP RDF of Si–O is nearly independent on TD conditions under which the vitreous silica was prepared. It is given by a combination of the following reasons i) minima of Si–O pair potentials are close each other for both potentials; ii) coulombic contribution into energy coming from beyond the nearest neighbour is largely isotropic; iii) a-SiO2 structure is dominantly determined by (SiO4)4− tetrahedra that are in turn given rather by steric effects; and iv) RDF is one-dimensional projection of the structure smoothing anisotropy around atoms. Angle distribution seems to be more sensitive to the TD conditions. Glasses with higher densities generate wider O–Si–O distribution as a result of less spatial symmetry in (SiO4)4− tetrahedra and their maxima are shifted to lower angles due to the higher spatial requirements imposed by the denser atomic arrangement. Larger differences are comprehensibly seen for Si–O–Si angle distribution as this angle is the most flexible part of the vitreous silica structure and therefore it can be most easily accommodated to the requirements imposed by TD. The distribution seems to be mainly determined by glass density again but two limiting distributions CHIK-q and BKS-e glasses demonstrate the increased sensitivity of the imposed TD conditions to higher (CHIK) or lower (BKS) densities for the particular potential sets. Even though RDF's and angle distributions reveal only subtle differences all presented figures document that (BKS-p, BKS-q) and (CHIK-p, CHIK-e) are forming couples with remarkable close results if compared to the third glass (e-glass for BKS and q-glass for CHIK sets). It confirms

Fig. 8. Radial distribution functions of primitive rings, RRDF's, for CHIK (green) and BKS (dark red) glasses prepared under various TD conditions.

the previous findings that BKS potential, originally optimised for quartz, bears an imprint of the origin of its derivation [7]. The specific values of defects found in simulated glasses depend on how the chemical bonds are established. The distance criterion takes into account the distances between silicon and oxygen atoms only. The amounts of overcoordinated atoms are expected to increase and the number of dangling bonds are expected to decrease with an increasing cut-off distance. Nevertheless, if PP RDF(Si–O) reveals a region within which the function is equal to zero, the number of defects arising from Si–O arrangement should be independent on a variation of the cut-off chosen from this interval. At low temperatures (300 K) the PP RDF(Si–O) exhibits a region around 2 Å where the distribution function is nearly equal to zero. Direct evaluation, see Fig. 5, shows only a small variations in the number of defects so that the used cut-off of 2 Å may be taken as reasonable. At higher temperatures, the minimum beyond the first peak on RDF(Si–O) never reaches zero and the structure must be analysed with a more care. Overcoordination was partially removed by the algorithm suggested for OO. Bonds are determined from snapshots and some of the established bonds/defects are metastable only and quickly decay [21]. Therefore their partial elimination by the used algorithm (that did not introduce any new defects) makes the final topology of glass closer to the real glass. The above presented defects are based on coordination, i.e. on geometrical positions of the atoms in the vicinity of the chosen one. Coordination/number of defects is nevertheless closely related with Qn distribution and connectivity as shown by Eqs. (4)–(11). In ideal structure of vitreous silica all Qn-units have to be Q4's only and both Q3- and Q5-units should be taken as structural defects. Table 5 shows that Q5's significantly prevail over Q3's in all glasses resulting in network connectivity higher than four. The number of Q4's are generally higher for BKS potential, indicating this way BKS potential generates lower number of structural defects in comparison with CHIK potential under the same TD conditions. The evaluated connectivity of the glass network is ranging from 4.009 to 4.026, the smallest value is obtained for BKS-e glass and the highest one for CHIK-q glass. Despite striking differences in particular defects, the connectivity of glass networks was finally determined by rather TD than type of the potential. Lower amount of DO of BKS glass was balanced by the appropriately smaller amount of OO according to Eq. (11) so that despite larger amount of oxygen defects in CHIK glasses connectivity of BKS and CHIK glasses are practically the same under the same TD conditions. Isobaric regimes show that connectivity cannot be simply related to density; large differences between densities of CHIK-p and BKS-p glasses resulted in the same connectivity. Coordination, Qn distribution, and connectivity can be alternated on base of an additional criterion of breakage of surplus connections, e.g. connections from OO or OSi. Such remove of BO would improve (toward ideal) connectivity but would also lead to the formation of additional DO so that overcoordination (amounts of Q5's and Q6's) would be decreased at the cost of increase of DO (number of Q3's). Formally, using additional bond scissors could lead to the balanced disproportionation of Q4's into Q5's and Q3's. Even though such disproportionation is a natural precursor for NBO migration and is therefore present in glass naturally, indicated approach is questionable and cannot be justified by experimental findings. The amount of defects should influence distribution of rings in the silica network, commonly used for the description of topology of medium range order. Larger amount of OO and OSi should intuitively lead to the higher number of small rings but Fig. 7 shows that this conclusion is not valid and, moreover, no simply courses explaining the changes in ring distributions can be found. It rather follows the ring distribution, and therefore medium range order topology, is robust to changes in density and/or TD regime and the alternated conditions of glass preparation can be fully accommodated by the vitreous silica topology without any significant change. Except for the purely topological description of MRO via ring distribution, RRDF adds a geometrical

O. Gedeon / Journal of Non-Crystalline Solids 426 (2015) 103–109

feature to ring relations. RRDF compromises four well recognized peaks documenting the ordering of rings up to 8 Å. However, partially due to the higher statistical scattering of data no conclusions about the differences among glasses can be drawn. On the contrary, RRDF rather confirms the robustness of MRO. 5. Conclusions Two potential sets, BKS and CHIK, and three themodynamic regimes were utilised in MD simulations to prepare the corresponding structures of vitreous silica. Standard structure descriptors as RDF's, angle distributions, and primitive ring distributions showed only subtle differences among glasses. Structural defects as overcoordinations and dangling bonds revealed higher sensitivity to both used potential sets and TD conditions. CHIK potential generates a larger number of defects under the same simulating conditions but connectivity of glasses are rather dependent on TD conditions than on potential set. Overcoordination is largely determined by the density while amount of dangling bonds is mainly determined by the potential used. Medium range order, represented by ring distribution and RDF of rings is robust and accommodates various requirements imposed on the system. Acknowledgement This work was supported by the Grant Agency of the Czech Republic through the grant No 15-12580S. References [1] P.H. Gaskell, Structure and properties of glasses — how far do we need to go? J. NonCryst. Solids 222 (1997) 1–12. [2] X.L. Yuan, A.N. Cormack, Efficient algorithm for primitive ring statistics in topological networks, Comput. Mater. Sci. 24 (2002) 343–360. [3] L. Guttman, Ring structure of the crystalline and amorphous forms of silicon dioxide, J. Non-Cryst. Solids 116 (1990) 145–147.

109

[4] K. Goetzke, H.J. Klein, Properties and efficient algorithmic determination of different classes of rings in finite and infinite polyhedral networks, J. Non-Cryst. Solids 127 (1991) 215–220. [5] O. Gedeon, M. Liska, Rings in covalent glass and an evaluation of configurational entropy associated with rings, J. Non-Cryst. Solids 360 (2013) 41–48. [6] N.T. Huff, E. Demiralp, T. Cagin, W.A. Goddard, Factors affecting molecular dynamics simulated vitreous silica structures, J. Non-Cryst. Solids 253 (1999) 133–142. [7] O. Gedeon, M. Liska, J. Machacek, Molecular dynamics and silicate glass — mutual challenges, Phys. Chem. Glasses Eur. J. Glass Sci. Technol. B 48 (2007) 388–393. [8] K. Vollmayr, W. Kob, K. Binder, Cooling-rate effects in amorphous silica: a computersimulation study, Phys. Rev. B 54 (1996) 15808–15827. [9] B.M. Lee, S. Munetoh, T. Motooka, Molecular dynamics study of velocity distribution and local temperature change during rapid cooling processes in excimer-laserannealed silicon, Comput. Mater. Sci. 37 (2006) 198–202. [10] R. Bruning, K. Samwer, Glass-transition on long-time scales, Phys. Rev. B 46 (1992) 11318–11322. [11] B.W.H. Vanbeest, G.J. Kramer, R.A. Vansanten, Force–fields for silicas and aluminophosphates based on ab initio calculations, Phys. Rev. Lett. 64 (1990) 1955–1958. [12] A. Carre, J. Horbach, S. Ispas, W. Kob, New fitting scheme to obtain effective potential from Car–Parrinello molecular-dynamics simulations: application to silica, EplEurophys. Lett. 82 (2008). [13] P.P. Ewald, Die Berechnung optischer und elektrostatischer Gitterpotentiale, Ann. Phys. 64 (1921) 25. [14] I.T. Todorov, W. Smith, K. Trachenko, M.T. Dove, DL_POLY_3: new dimensions in molecular dynamics simulations via massive parallelism, J. Mater. Chem. 16 (2006) 1911–1918. [15] O.V. Mazurin, M.V. Streltsina, T.P. Shvaiko-Shvaikovskaya, Handbook of Glass Data, Elsevier Science Publisher, Amsterdam, 1983. [16] O. Gedeon, M. Liska, Partial volumes and volume fluctuations in potassium-silicate glass: molecular dynamic simulations, J. Non-Cryst. Solids 351 (2005) 1682–1687. [17] T.F. Soules, G.H. Gilmer, M.J. Matthews, J.S. Stolken, M.D. Felt, Silica molecular dynamic force fields — a practical assessment, J. Non-Cryst. Solids 357 (2011) 1564–1573. [18] A.C. Wright, Neutron-scattering from vitreous silica. 5. The structure of vitreous silica — what have we learned from 60 years of diffraction studies, J. Non-Cryst. Solids 179 (1994) 84–115. [19] H.F. Poulsen, J. Neuefeind, H.B. Neumann, J.R. Schneider, M.D. Zeidler, Amorphous silica studied by high-energy X-ray-diffraction, J. Non-Cryst. Solids 188 (1995) 63–74. [20] W.J. Malfait, W.E. Halter, R. Verel, Si-29 NMR spectroscopy of silica glass: T-1 relaxation and constraints on the Si–O–Si bond angle distribution, Chem. Geol. 256 (2008) 269–277. [21] J. Machacek, O. Gedeon, M. Liska, Relaxation of structural units in MD simulated silicate glasses, Phys. Chem. Glasses-B 47 (2006) 266–270.