New interatomic potential parameters for molecular dynamics simulations of rare-earth (RE = La, Y, Lu, Sc) aluminosilicate glass structures: exploration of RE3+ field-strength effects.

Sets of self-consistent oxygen-rare earth (RE = La, Y, Lu, Sc) interatomic potential parameters are derived using a force-matching procedure and utilized in molecular dynamics (MD) simulations for exploring the structures of RE2O3-Al2O3-SiO2 glasses that feature a fixed molar ratio n(Al)/n(Si) = 1 but variable RE contents. The structures of RE aluminosilicate (AS) glasses depend markedly on the RE(3+) cation field strength (CFS) over both short and intermediate length-scales. We explore these dependencies for glasses incorporating the cations La(3+), Y(3+), Lu(3+) and Sc(3+), whose CFSs increase due to the concomitant shrinkage of the ionic radii: R(La) > R(Y) > R(Lu) > R(Sc). This trend is mirrored in decreasing average RE(3+) coordination numbers (Z(RE)) from Z(La) = 6.4 to Z(Sc) = 5.4 in the MD-derived data. However, overall the effects from RE(3+) CFS elevations on the local glass structures are most pronounced in the O and {Al([4]), Al([5]), Al([6])} speciations. The former display minor but growing populations of O([0]) ("free oxygen ion") and O([3]) ("oxygen tricluster") moieties. The abundance of AlO5 polyhedra increases significantly from ≈10% in La-based glasses to ≈30% in their Sc counterparts at the expense of the overall dominating AlO4 tetrahedra, whereas the amounts of AlO6 groups remain <5% throughout. We also discuss the Si([4])/Al([p]) (p = 4, 5, 6) intermixing and the nature of their oxygen bridges, where the degree of edge-sharing increases together with the RE(3+) CFS.

Owing to their chemical inertness and favorable physical properties, e.g., high glass transition temperatures (T g ), high microhardness (H V ) and low thermal expansion coefficients, lanthanide-bearing AS glasses are explored as model systems for storage of long-lived actinides, 19 whereas Y-Al-Si-O glasses are exploited in radiotherapy for in situ cancer treatment. 9,20 There is a strong correlation between many physical/mechanical features of RE 2 O 3 -Al 2 O 3 -SiO 2 glasses and the RE 3+ cation fieldstrength (CFS); the CFS is given by z/R 2 , where z and R denote the charge and radius of the ion, respectively. However, the detailed structure-property relationships of RE AS glasses are poorly understood, such as why the T g -value of the glass strongly depends on the RE 3+ CFS 10,11,[15][16][17][18] but is almost independent of its RE-Al-Si composition. [9][10][11]16,17 While such peculiarities inarguably stem from the presence of RE 3+ ions and most likely originate from their bearings on structural features over both short and intermediate length-scales t1 nm, the currently incomplete structural insight hampers unambiguous rationalizations of many structure-property correlations of RE 2 O 3 -Al 2 O 3 -SiO 2 glasses. We have recently demonstrated that the RE 3+ effects on such relationships may be subtle, for example that the well-known enhanced microhardness for increasing RE 3+ CFS originates from the accompanying elevated amounts of high-coordination AlO 5 / AlO 6 polyhedra, 16,21 whose populations are in turn dictated by the RE 3+ CFS. [12][13][14][15][16][17][18] Hence, it is desirable to reach a more detailed structural picture that encompasses (i) the local RE 3+ environments (coordination numbers and RE-O distances) and their distribution across the structure, as well as (ii) the speciations and intermixing of the various SiO 4 and AlO p (p = 4, 5, 6) groups associated with the glass network. Data obtained from X-ray and neutron diffraction techniques 18,22 constitute the primary information source for feature (i), whereas the hitherto most detailed insight into the network structures stem from 17 O, 27 Al, and 29 Si solid-state nuclear magnetic resonance (NMR) spectroscopy applied to AS glasses incorporating the diamagnetic (and thereby accessible to NMR) La 3+ , Y 3+ , Lu 3+ and Sc 3+ ions [12][13][14][15][16][17][18][23][24][25] that are listed in the order of decreasing ionic radii and thereby increasing CFS. Besides including both end-members of the lanthanide series, this ion-set features a large CFS-span from 2.8 Å À2 for La 3+ to 5.4 Å À2 for Sc 3+ , thereby allowing for the probing of structural variations across a wide range of field-strengths.
As the experimentally derived structural understanding is insufficient, molecular dynamics (MD) simulation is an attractive method for gaining complementary insight. However, such reports of RE 2 O 3 -Al 2 O 3 -SiO 2 glasses are hitherto limited to the Y system, 16,[26][27][28] except in our recent modeling work on Lu 16 and La 25 AS glasses. The present contribution serves the dual purpose of (i) outlining the approach for deriving the RE-O (RE = La, Y, Lu) interionic potential parameters utilized by us in ref. 16 and 25, as well as introducing new parameters for the Sc-O pair; (ii) exploring the relationships between the local structural features and the nature of the RE 3+ ion (RE = La, Y, Lu, Sc) in two series of RE-Al-Si-O glasses; RE(2.21) and RE (2.45). Their members feature equal Al and Si molar contents but distinct glass network polymerization degrees 29 of r = 2.21 and r = 2.45, respectively, with r = n O /(n Al + n Si ), where n E and x E denote the stoichiometric amount and molar fraction of species E, respectively (the latter out of the cations): x E = n E /(n Si + n Al + n RE ), with E = Si, Al, or RE.
The glass compositions correspond to those reported by Iftekhar 11,[15][16][17]24,25 This paper is organized as follows: Section 2 describes our general MD simulation procedures, whereas Section 3 outlines our implementation of the three-step force-matching procedure of Ercolessi and Adams, 30

Computational details
The classical MD simulations emulating a melt-quench procedure were performed for an NVT ensemble in a cubic box with periodic boundary conditions. 31 The computations were carried out with the DLPOLY3 package 32,33 for E3300 atoms and a cubic box-length in the range 3.33-3.54 nm. These parameters were adjusted slightly to match the experimental density and composition of each RE AS glass, using the previously reported values for the La, 11,24 Y, 16 Lu 16 and Sc 17 bearing samples. The experimental densities ranged between 2.93 g cm À3 for Sc(2.21) and 5.05 g cm À3 for Lu (2.45). 15 The melt-quench simulation started from a structure equilibrated at 3500 K for 100 ps, followed by a 10 ps step-wise decrease of 10 K ps À1 down to 300 K, where equilibration was performed during 200 ps, of which the trajectory during the last 150 ps was used for the structural analysis. For each glass composition, this protocol was completed twelve times with distinct initial ion configurations, each generated by a random distribution subject to the constraint of a 50 pm minimum distance in any ion-pair. The average value and uncertainty of each reported modeled structural feature was derived from these distinct trajectories.
For two ion species a and b with indices j and k separated at a distance R jk a-b , the corresponding interaction energy included both long-range Coulombic where z a represents the ionic charge), and short-range Buckingham terms. The latter is given by except for the RE-O (RE = La, Y, Lu, Sc) interionic potentials derived herein that solely employed the repulsive Born-Mayer contribution: i.e., for which C 0. The sets {A a-b ,r a-b ,C a-b } are listed in  34 where e is the elementary charge.
A modified Buckingham potential was implemented to circumvent strong attractions at small interionic distances. 35 Eqn (3) and (4) were evaluated for all ion-pairs up to 0.8 nm, whereas the Coulombic interactions were calculated by a smoothed particle mesh Ewald summation 33 with a 1.2 nm real-space cut-off and an accuracy of 10 À6 . These parameters provide sufficient accuracy and are standard for MD implementations for similar glass systems (e.g., see ref. 2, 27 and 28). The equations of motion were integrated in time-steps of 2 fs by the velocity Verlet integrator approach, and a Berendsen thermostat with a 1.0 ps relaxation time constant. 31 These parameters and methods were employed throughout all MD simulations, except for the optimization procedure of the new RE-O potential parameters in Section 3, which employed a smaller set of 113 atoms and box-lengths spanning 1.12-1.15 nm, whereas 0.53 nm and 0.55 nm cut-off radii were used for the Buckingham and Coulombic interactions, respectively. The restricted set of atoms originates from computer memory constraints in the ab initio calculations outlined below.
To verify that the sets of interionic potential parameters may also reproduce the experimental glass densities, we additionally performed NPT simulations for four distinct initial configurations of each composition. The modeled densities were consistently slightly lower than their experimental counterparts, but deviated only by t3% and t5% within each RE(2.21) and RE(2.45) series, respectively. The overall largest discrepancy (5.3%) was observed for the Y(2.45) composition, whereas essentially perfect matches resulted for the two La glass compositions (o1% relative error). The sets of simulations revealed a negligible spread within AE0.01 g cm À3 for each glass; incidentally, this value is equal to the experimental uncertainty. [15][16][17] 3 Optimization of RE-O interionic potential parameters Current computer resources limit ab initio calculations to systems comprising less than a few hundred atoms over t10 ps, and simulations of larger ensembles over long timescales can only be achieved by classical MD approaches 31 that require a set of parameters for approximating the interaction energy in each ion-pair of a targeted structure. An ''ideal'' interaction potential should reproduce the experimentally assessed properties of the system, as well as providing reliable predictions over widely spanning sample compositions. Several strategies exist for constructing interatomic potentials, such as the reverse Monte Carlo, 36,37 electronic density, 38 and effective medium 39 options.
We employed the iterative force-matching procedure of Ercolessi and Adams 30 to obtain the parameters A a-b and r a-b of the Born-Mayer potential [eqn (4)]; their values were selected so as to provide the closest possible match to ab initio-derived forces. Such a protocol is convenient for applications to amorphous phases where elastic constants or lattice parameters cannot be defined. It was carried out for the RE(2.21) composition of each RE-Al-Si-O (RE = La, Y, Lu, and Sc) glass. Every iteration step (n) involves three stages: ( (2) The ion coordinates of each MD-derived glass structure were utilized in a density functional theory (DFT) full-potential linearized augmented plane wave (FP-LAPW) 43,44 calculation to provide the interionic forces [ % F(ab initio)]. The generalized gradient approximation of Perdew, Burke and Ernzerhof 45 was used with the Wien2k package. 46 One k-point of the irreducible Brillouin zone was computed by employing basis functions that obeyed the criterion K max R min = 7.0, where K max represents the magnitude of the largest basis vector in the reciprocal space and R min is the smallest muffin-tin radius in the simulated cell. The latter values were fixed throughout all computations at 1.83 a.u., 1.66 a.u., 1.51 a.u., and 1.31 a.u. for the RE, Al, Si, and O species, respectively. 46 We verified that the product K max R min = 7.0 provides converged forces. The electronic band-states were separated by À8.1 Ry to include all contributions down to the 2p and 2s orbitals for Si and Al, respectively. The force convergence criteria for the self-consistent field (SCF) procedure were 0.5 mRy per a.u. During the last SCF-cycle, the total forces were calculated using Pulay's correction. 47 (3) The DLPOLY3 package 32,33 was employed to compute a set of forces { % F j (A RE-O ,r RE-O )}. A refined set of interionic potential parameters resulted from a least-square force-fitting procedure that minimized the function where the index j runs over all RE 3+ sites in the structure and ||Á Á Á|| denotes the magnitude (norm) of the vector-difference. The Born-Mayer parameters were varied until the best match was found in eqn (5). This provided a new set of values, {A (n+1) RE-O ,r (n+1) RE-O }, extracted at the minimum of w 2 (A RE-O ,r RE-O ) and subsequently provided as input to the MD simulation stage (1) of the next, i.e., (n + 1)th iteration step.
Steps (1)-(3) of the protocol were repeated three times, resulting in the sets {A (n) RE-O ,r (n) RE-O } with n = 1, 2, and 3. Table 1 lists the {A (3) RE-O ,r (3) RE-O } values that were employed in all remaining calculations, as well as in our recent work. 16,25 The parameters obtained at each iteration stage were used to compute the potential energy [U RE-O (R RE-O )] from eqn (4), as well as the respective radial distribution function [RDF; denoted g RE-O (R RE-O )]. The latter was evaluated for an ensemble of E3300 atoms, as described in Section 2. Fig. S1 of the ESI † shows the results for the RE(2.21) glass composition and RE = {La, Y, Lu, Sc}.

Short-range structures of RE 2 O 3 -Al 2 O 3 -SiO 2 glasses
The general validity of the new potential parameters for reproducing short-range structural features of RE 2 O 3 -Al 2 O 3 -SiO 2 glasses was previously confirmed over large RE-Al-Si compositional ranges within each La, 25 Y 16 and Lu 16 system. Focussing on the structural dependency on the identity of the RE 3+ ion, we evaluate here the predictions of the MD simulations for the series of RE(2.21) and RE(2.45) specimens against experimental results from either solid-state NMR 16,17,25 (Si and Al coordinations) or diffraction 18,22 techniques; the latter data on cationoxygen distances were reported for Si-richer and RE-poorer compositions and may therefore only serve as general guides.
4.1.1 Si and Al coordinations. The MD simulations of the present glass compositions, as well as others previously reported, 16,17 reveal that Si 4+ is exclusively (>98.5%) present in tetrahedral coordination (Si [4] ), which is in excellent agreement with experimental solid-state 29 Si NMR data. [15][16][17]23,24 For this reason, we henceforth do not explicitly indicate the coordination number of Si. All Al 3+ ions (Z99.9%) coordinate either four (Al [4] ), five (Al [5] ) or six (Al [6] ) oxygen atoms, also in full accordance with 27 Al NMR results obtained from the present RE AS compositions as well as related ones. [12][13][14][15][16][17][18]23,25 Fourfold Al 3+ coordinations dominate the three coexisting AlO 4 , AlO 5 and AlO 6 polyhedra. The amounts of higher-coordination species (mainly AlO 5 ) grow primarily for increasing RE 3+ CFS and secondly for decreasing x Si , in good agreement with experimental findings. 13,14,16,17,25 These trends are witnessed in the average Al coordination number ( % Z Al ) and fractional populations {x [p] Al } of the {AlO p } polyhedra listed in Table 2 for the two RE(2.21) and RE(2.45) series of specimens. As the CFS elevates between the La 3+ and Sc 3+ ions, the modeled % Z Al values range over 4.14-4.43 for the RE(2.45) series, whereas its RE(2.21) counterpart displays nearly equal average Al coordination numbers E4.24 in the Y, Lu and Sc bearing structures. Minor variations are observed for each RE AS structure when x Si diminishes, i.e., for increasing the r-value of the glass. The elevation of % Z Al stems primarily from a growth of the AlO 5 population, whereas that for the AlO 6 groups remains low throughout (x [6] Al o 0.05). Al values are plotted against the RE 3+ CFS in Fig. 1. While excellent agreement is observed for all La and Y AS glasses, for the cases of Lu and Sc the MD simulations underestimate the AlO 5 populations (and thereby % Z Al ), notably so for the Sc(2.21) glass. Despite the largest relative discrepancy between the MD and NMR derived % Z Al -values being only E4% [for Sc(2.21)], the pronounced experimental trend of a monotonic elevation of % Z Al for increasing RE 3+ CFS is not reproduced well by the simulations of the Lu and Sc members, particularly for the r = 2.21 branch. Given the significantly faster quench-rate in the calculations relative to the physical samples, one expects the modeled x [5] Al and x [6] Al populations to be comparable to, or even higher, than their NMR-derived counterparts, as explained in detail in ref. 25. Hence, the frequently observed underestimation of % Z Al from the MD simulations is surprising (see Fig. 1). We have no explanation for these effects but refer to Jaworski et al. 25 for further discussions.  Al ) obtained by MD simulations and NMR. The latter data are given within parentheses and are reproduced from our recent work. 16 Table 3 lists a collection of average interionic Here we focus on the various cation-oxygen distances.

Crystalline model structures
The new RE-O potentials were further evaluated by energy minimizations at zero absolute temperature by applying the GULP program 48 to the set of crystalline structures shown in . Each representative XRD-derived structural parameter-set (obtained from the inorganic crystal structure database 49 ) was used as the initial input to the minimization procedure with the interionic potentials evaluated up to 2.0 nm. The space-group symmetry of each structure was preserved throughout.
The resulting lattice parameters (a, b, c), cell volume (V) and average RE-O distance ( % R RE-O ) within the unit cell are listed in Table 4 and compared with their experimental counterparts obtained using single crystal XRD. An acceptable agreement is observed: the relative deviations among the (a, b, c) parameters and the volume of the unit cell generally stay below 2% and 4%, respectively. For both the RE 2 O 3 and aluminate structures, the data derived from the new Y-O and Lu-O potential parameters accord very well with their experimental counterparts, where the largest discrepancies in the cell-lengths, volumes and average RE-O distances amount to 0.8%, 2.2% and 1.6%, respectively, while the agreement is overall significantly lower for La and Sc.  (Table 4), the globally largest discrepancy among all structures is observed for La 2 O 3 , for which the calculated values of V and % R La-O are underestimated by 7.3% and 2.4%, respectively. However, the three Sc-bearing structures manifest the overall weakest agreement between calculated and experimental lattice parameters, where discrepancies over the ranges of 1.2-2.1%, 3.8-5.8% and 2.0-3.1% are observed for the values of (a, b, c), V, and % R Sc-O , respectively. Nevertheless, our new RE-O potential parameters that were optimized directly on amorphous phases, display a good transferability to crystalline structures, except for La 2 O 3 and the Sc-based phases. For instance, the observed deviations between calculated and experimental data are not substantially larger than the analogous results of Pedone et al. 50 obtained by a Morse-based potential that was optimized directly on crystalline oxides and subsequently evaluated on several multicomponent structures. Table 4 Table 1); hence, the validity of our RE-O Born-Mayer parameters is confirmed by the comparable quality in the structural predictability observed both in their presence and absence.  Table 4). Fig. 3(a) plots the average number of oxygen species coordinated by each RE 3+ ion ( % Z RE ) against its corresponding Shannon-Prewitt radius, 42 the latter value relevant for REO 6 groups. Indeed, all MD-derived mean coordination numbers scatter around % Z RE B 6; they elevate either for increasing RE 3+ radius (i.e., decreasing CFS) or for growing r-value of the glass. The latter observation may be rationalized from the dependence of % Z RE on the glass composition, where both % Z RE and % Z Al increases for decreasing silica content (i.e., increasing x RE ), as demonstrated in ref. 16 52 as well as with several MD reports on RE-bearing (alumino)silicate glasses. 6,7,[26][27][28] Interestingly, the calculations predict significant fractions 0.09 r x [4] Sc r 0.17, suggesting that Sc 3+ might partially assume a glass-network forming role in the guise of ScO 4 groups, as reported previously for the similarly-sized Mg 2+ ion. [53][54][55] Fig. 4(a) displays a typical Sc [4] environment in the Sc(2.21) glass. The potential network-forming capability of Sc requires further experimental studies, but may constitute the origin of previously reported anomalies in 29 Si NMR shifts and unexpectedly low glass transition temperatures observed from Sc AS glasses. 15,17 5.2 Oxygen environments 5.2.1 Oxygen coordinations. We identify the coordination number p of an oxygen species O [p] in the RE-Al-Si-O glass according to its total number of bonds to Si and/or Al. Hence, O [1] and O [2] represent non-bridging oxygen (NBO) and BO species, respectively. 56 [3] (''oxygen tricluster'') moieties in the present (as well as previously reported 16,[25][26][27] ) RE AS glasses that involve high-CFS trivalent modifier ions is the primary distinction between their corresponding O speciations in mono-or di-valent based glasses of otherwise comparable compositions; compare for instance our results with previous modeled data from Ca AS glasses. 57,58 Whereas both the O [0] and O [3] populations grow slightly as the RE 3+ CFS increases, a Results of using the interionic potential parameters of Table 1 Table 1). they also depend on the RE 3+ content, with x [0] O increasing and x [3] O decreasing as the r-value of the glass elevates from 2.21 to 2.45, as revealed in Fig. 5. The minor growth of the O [0] populations for an increase in either the CFS or amount of the RE 3+ ions suggests an enhanced RE-O-RE association.

Rare-earth coordinations
While the abundance of free O 2À ions depends solely on the nature and amount of the glass modifiers, the O [3] population is additionally dictated by the n Al /n Si ratio of the glass composition. 16,25 The O [3] sites predominantly connect Al-centered polyhedra, notably so high-coordination AlO 5 and AlO 6 groups (see Section 5.2.2). Because the latter populations grow together with both x Al and the RE 3+ CFS, 16,25 the highest x [3] O -values (\0.15) are observed for high-CFS Al-rich glasses associated with relatively low values of r t 2.21, whereas x [3] O stay consistently below 0.10 in La-Al-Si-O structures. 25 These trends may be verified from our previous reports on glasses featuring variable RE : Al : Si contents within each La, 25 Y, 16 and Lu 16 system. Fig. 5(b)-(d) reveal that the O [3] moieties grow largely at the expense of their O [2] counterparts. Further, for a fixed RE : Al : Si stoichiometric ratio, Table 2 verifies a strong correlation between (i) % Z Al , which elevates concurrently with x [5] Al , as the AlO 6 population remains low throughout all glass structures; (ii) % Z O and (iii) x [3] O , all of which increase with the RE 3+ CFS. However, whereas the % Z Al / % Z O /x [3] O correlations are

Trends in cation-O [p] contacts.
We now examine the various cation-oxygen contacts, initially from the viewpoint of each Si 4+ , Al 3+ , and RE 3+ cation. For simplicity, we treat all RE [p] coordinations collectively and only summarize the general trends observed throughout all MD-generated structures.
Si exhibits the least preference for any particular O [p] coordination, meaning that the relative number of Si-O [p] bonds roughly reflect each respective x [p] O -value in Table 2. Conversely, RE-O [1] contacts are strongly preferred over their RE-O [2] counterparts, approximately by a factor of 3 compared to the predictions from a statistical intermixing based on the [2] O ratio. Fig. 6 plots the various RE-O [p] populations. The relative preferences among the various cations to coordinate the less abundant O [1] and O [3] species constitute their main distinctions: the affinity for X-O [1] bond formation decreases along the series RE(10) c Si(4) > Al [4] (2) > Al [5] \ Al [6] (1), (6) with the numbers within parentheses roughly conveying the relative preference of species X to coordinate O [1] compared to that of Al [6] that is least prone to accommodate NBO ions. Typically, \60% of all NBO ions are located at the SiO 4 groups, whereas t30% and t15% reside on AlO 4 and AlO 5 /AlO 6 polyhedra, respectively. In contrast, the affinity of a cation X to coordinate O [3] moieties manifests the reversed trend: besides relatively sparse RE-O [3] contacts, the preference for forming X-O [3] bonds scale roughly as 1 : 3 : 5: 7-10 for Si : Al [4] : Al [5] : Al [6] . As expected from bond-valence sums, O [3] moieties coordinate at most one Si atom and up to three Al [p] species, where generally at least one is a high-coordination AlO 5 /AlO 6 polyhedron, as illustrated in Fig. 4(c). The strong preference for Al [p] -O [3] linkages implies a slight Al-Al self-association across the structure, as also noted previously by Christie and Tilocca. 27 Despite the ability of Si and Al [p] to accommodate each O [p] species being independent of the nature of the RE element in the glass, the corresponding abundance of each cation-oxygen contact scales with the relative {x [1] O , x [2] O , x [3] O } values, where the x [2] O /x [3] O ratio is dictated both by the RE 3+ CFS and content, as highlighted above.
We next focus on the various Si, Al, and RE cation constellations around each O [p] site, denoted here as O [p] [q] , where q represents the total coordination number when also accounting for the q-p RE-O contacts. For both the O [1] [q] and O [2] [q] environments, the net coordination number q ranges between 2 and 4, implying that each O atom coordinates 0, 1, or 2 RE 3+ cations. In contrast, the free O 2À anions predominantly coordinate 3 (major contribution; 50-80%) or 4 (minor; 10-50%) RE species, where the abundance of the latter grows for increasing RE 3+ content (x RE ), but diminishes markedly across each RE(r) series as the CFS increases. For instance, O [0] [4] groups constitute E50% and E20% of the O [0] speciation in the La(2.45) and Sc(2.45) structures, respectively. The O [3] moieties predominantly constitute O [3] [3] groups devoid of linkages to any RE 3+ species, with the O [3] [4] population amounting to 5-12% and 10-18% for the RE(2.21) and RE(2.45) series, respectively.
The various SiO 4 /AlO p polyhedra manifest a non-preferential intermixing, leading to an MD-derived fractional population x(X-Y) of an XO p -YO q pair that is well-approximated by its corresponding product Each of y X and y Y represents the fractional population of the given species X and Y out of the total Si and Al speciation, which for the present glass compositions (that exhibit n Si /n Al = 1) corresponds to y Si = 0.5 and y [p] Al listed in Table 2. Consequently, owing to the dominance of SiO 4 and AlO 4 groups followed by AlO 5 , the most frequently encountered pairs of structural building blocks represent SiO 4 -AlO 4 , SiO 4 -SiO 4 , AlO 4 -AlO 4 motifs, and to a lesser extent SiO 4 -AlO 5 and AlO 4 -AlO 5 . Fig. 7 plots the relative abundances of the various pair-connectivities. Note that owing to the small amounts of AlO 6 polyhedra, all AlO 6 -XO p contacts remain low throughout all structures, with the highest fractions x(Si-Al [6] ) E x(Al [4] -Al [6] ) E 0.03 observed in the Sc(2.45) glass.
We verified that the modeled x(X-Y) fractions closely obey a statistical intermixing, where the deviations between the MD-derived values and the predictions from eqn (7) are readily explained by the relative preference of each {SiO 4 , AlO p } polyhedral type to accommodate NBO ions [see eqn (6)]: the comparatively most abundant Si-NBO bond formation leads consistently to somewhat lower x(Si-Si) populations relative to those predicted by eqn (7), whereas the reluctance of the AlO p groups to accommodate NBO ions implies the opposite trend for the x(Al [p] -Al [q] ) values. In the context of La 2 O 3 -Al 2 O 3 -SiO 2 glasses, a pronounced Si-Al disorder was verified experimentally for SiO 4 -AlO 4 [ref. 24] and AlO p -AlO q [ref. 25] contacts, in the latter case by a direct connectivity-probing of the various AlO p -AlO q pairs by double-quantum 27 Al NMR spectroscopy. 59 5.3.2 Extents of corner and edge sharing. As illustrated by Fig. 8(a), the Si/Al [p] -centered polyhedra merge by sharing corners or edges that constitute >85% and 5-15% out of the entire BO ensemble, respectively. The RE AS structures are devoid of face-shared SiO 4 or AlO p groups. For increasing RE 3+ CFS, a slight growth in the relative degree of edge-sharing is observed across each RE(r) glass series, which concerns both the O [2] and O [3] bridges [ Fig. 8(b) and (c)], particularly for the RE(2.45) members. Yet, the relative preferences of the BO species to participate in edge-sharing differ substantially: while \95% of all O [2] atoms connect polyhedral corners [ Fig. 8(b)], 45-60% of the O [3] counterparts are involved in edge-sharing [ Fig. 8(c)]. Note that while both the amounts of O [3] moieties (x [3] O )  [4] , Al [5] , Al [6] } speciation. Grey and dotted traces reveal the respective y(X-O [3] -Y) and y(X-O [2] -Y) contributions (where {X, Y} = {Si, Al}) to the total ADF (black trace).
and their affinity for edge-sharing increases concurrently with the RE 3+ CFS, the x [3] O values decreases with the r-value, as discussed in Section 5.2.1.
Striking differences further emerge when comparing the preferences among the various {SiO 4 , AlO 4 , AlO 5 , AlO 6 } groups to connect through corners or edges: the strong tendency of O [3] moieties to interlink the cations by edge-sharing, coupled with the inherently distinct affinities of the cations to form X-O [3] bonds, implies that the participation in edge-shared polyhedra grows along the series Si t Al [4] { Al [5] o Al [6] . Altogether, the net abundance of edge-shared XO p -YO q polyhedra scales roughly as Si-Si(0%) E Si-Al [4] (o1%) o Al [4] -Al [4] (o3%) o Si-Al [5] (E5%) o Si-Al [6] (E10%) o Al [4] -Al [5] (10-20%) o Al [4] -Al [6] (20-25%) { Al [5] -Al [5] (E50%) o Al [5] -Al [6] (60-70%) o Al [6] -Al [6] (\75%), where the numbers within parentheses reflect the edge-sharing contribution, with the remaining constituting corner-shared polyhedra of the respective type. However, whereas a majority of the AlO 5 -AlO 6 and AlO 6 -AlO 6 linkages occur through edgesharing, the total amounts of such structural motifs remain very low, as discussed above. Analogously to the relative preferences of each Al [p] -O [q] and Si-O [q] contact, no statistically ascertained variation is observed among the relative tendencies for corner/ edge-sharing of the SiO 4 /AlO p groups when the RE 3+ CFS alters. Fig. 4 illustrates a collection of structural fragments typically encountered in the RE AS glass structures. Edge-shared XO p -YO q polyhedra translate into bond-angles y X-O-Y E 901 that may be compared with their wider cornershared analogs peaking in the ranges of 110-1301 for Al-O-Al/Si and 140-1601 for Si-O-Si linkages (see Table 3). Hence, the abundance of edge-shared polyhedra may also be inspected from the plots of the MD-derived angle-distribution functions (ADFs) shown in Fig. 9 for each type of Si-O [p] -Si, Si-O [p] -Al and Al-O [p] -Al contact. The ADF amplitudes of the latter indeed increases markedly around y Al-O-Al = 901, which is most transparent for the structures involving Sc (and Lu; not shown) that exhibit the largest AlO 5 /AlO 6 populations and thereby the highest extent of edge-sharing.

Concluding remarks
New RE-O interatomic potential parameters were optimized for RE = {La, Y, Lu, Sc} by a force-matching procedure 30 involving ab initio derived forces in MD-generated RE-Al-Si-O glass structures. Overall they performed well both for reproducing diffraction-derived lattice parameters of crystalline RE-bearing oxide/aluminate structures, as well as short-range features in RE AS glasses, such as previously reported RE-O distances 18,22 and {Al [4] , Al [5] , Al [6] } speciations. 16,17,25 These self-consistent RE-O potential parameters are expected to be useful for exploring the structures of other amorphous as well as well-ordered RE-bearing oxide-based structures.
The MD-derived average Al coordination numbers ( % Z Al ) typically reproduce their experimental counterparts within o2% deviation, with the largest discrepancy (3.5%) observed for the Sc(2.21) glass. Whereas Sc 3+ mainly assumes Sc [5] coordinations and non-negligible amounts of Sc [4] species (up to E20%), all other (larger) RE 3+ ions reveal {RE [p] } speciations peaked at p = 6, with La [7] and Lu [5] constituting the second most abundant coordinations for La 3+ and Lu 3+ , respectively. The simulated structures exhibit essentially randomized connectivities among the various SiO 4 and AlO p (p = 4, 5, 6) network (associated) polyhedra, a feature attributed to the high charge of the trivalent RE 3+ ions, whose comparatively strong RE 3+ -O bonds (relative to their M + -O and M 2+ -O counterparts) perturb any Si/Al ordering and also promote otherwise energetically disfavored structural motifs, 56 such as free O 2À ions, Al [4] -O-Al [4] bridges, and Al-O [1] contacts. These trends are primarily reflected by the following structural alterations when the RE 3+  O values also increase with the RE 3+ content. (ii) Each of % Z Al and % Z O increases by elevating Al [5] and O [3] populations, respectively, at the expense of the dominating Al [4] and O [2] counterparts.
(iii) Whereas >85% of all SiO 4 and AlO p polyhedra connect through corners, the degree of edge-sharing increases steadily from E5% in the La-based glasses to E15% in the Sc analogs. Owing to the strong affinity for Al [5] -O [3] and Al [6] -O [3] bond formation combined with the pronounced tendency of O [3] moieties to create edge-sharing, large fractions of all AlO 5 and AlO 6 linkages to other groups occur via edge-sharing.
Altogether, the structures of RE 2 O 3 -Al 2 O 3 -SiO 2 glasses exhibit a pronounced disorder both in their distributions of O [p] , Al [p] and RE [p] coordinations and how the various structural motifs combine. Effectively, (i) and (ii) imply partial O [1] -O [0] and O [2] -O [3] net conversions as the RE 3+ CFS grows, thereby amounting in more complex O [p] speciations compared with the sole presence of O [1] and O [2] moieties expected from traditional theories of glass structures. 56 The present study corroborates and reinforces previous experimental and computational findings that the structures of RE 2 O 3 -Al 2 O 3 -SiO 2 glasses are strongly influenced by the RE 3+ CFS over both short and medium ranges. [12][13][14][15][16][17][18][24][25][26][27][28] The large CFS-span across the set {La 3+ , Y 3+ , Lu 3+ , Sc 3+ }, coupled with the overall monotonic (and occasionally approximately linear; see Fig. 5 and 6) trends of several structural features against the field-strength, suggest the possibility of roughly assessing the corresponding feature in another RE-Al-Si-O glass that exhibits a similar cation composition but a distinct RE 3+ ion, based on its CFS-value relative to those considered herein.