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

Abstract

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 RE₂O₃–Al₂O₃–SiO₂ 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³⁺ cation field strength (CFS) over both short and intermediate length-scales. We explore these dependencies for glasses incorporating the cations La³⁺, Y³⁺, Lu³⁺ and Sc³⁺, 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³⁺ 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³⁺ 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 AlO₅ polyhedra increases significantly from ≈10% in La-based glasses to ≈30% in their Sc counterparts at the expense of the overall dominating AlO₄ tetrahedra, whereas the amounts of AlO₆ 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³⁺ CFS.

---

1 Introduction

Thanks to their beneficial (magneto)optical properties, rare-earth (RE) bearing silicate-based glasses have found applications as laser hosts and optical amplifiers in telecommunication. For optimal performance of such devices, it is desirable to introduce high amounts of uniformly dispersed RE³⁺ ions into silica melts. Unfortunately, the RE₂O₃–SiO₂ system manifests liquid immiscibility, which leads to RE-clustering at low RE₂O₃ contents of a few mol% and to macroscopic phase-separations for RE³⁺-richer compositions.^1,2 The solubility may be increased by introducing RE³⁺ ions into a M₂O–SiO₂ based glass,^3–7 or by exploiting the ternary RE₂O₃–Al₂O₃–SiO₂ aluminosilicate (AS) system^8–18 that generally admits the formation of homogeneous amorphous phases by standard melt-quench procedures/conditions,^{8–11,15–17} unless their RE₂O₃ and/or Al₂O₃ contents are too low (typically ≲10 mol%, depending on the identity of the RE³⁺ ion^11,15).

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,¹⁹ 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₂O₃–Al₂O₃–SiO₂ glasses and the RE³⁺ cation field-strength (CFS); the CFS is given by z/R², 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³⁺ CFS^{10,11,15–18} but is almost independent of its RE–Al–Si composition.^9–11,16,17 While such peculiarities inarguably stem from the presence of RE³⁺ ions and most likely originate from their bearings on structural features over both short and intermediate length-scales ≲1 nm, the currently incomplete structural insight hampers unambiguous rationalizations of many structure–property correlations of RE₂O₃–Al₂O₃–SiO₂ glasses. We have recently demonstrated that the RE³⁺ effects on such relationships may be subtle, for example that the well-known enhanced microhardness for increasing RE³⁺ CFS originates from the accompanying elevated amounts of high-coordination AlO₅/AlO₆ polyhedra,^16,21 whose populations are in turn dictated by the RE³⁺ CFS.^12–18

Hence, it is desirable to reach a more detailed structural picture that encompasses (i) the local RE³⁺ 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₄ 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 ¹⁷O, ²⁷Al, and ²⁹Si solid-state nuclear magnetic resonance (NMR) spectroscopy applied to AS glasses incorporating the diamagnetic (and thereby accessible to NMR) La³⁺, Y³⁺, Lu³⁺ and Sc³⁺ ions^{12–18,23–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 Å⁻² for La³⁺ to 5.4 Å⁻² for Sc³⁺, 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₂O₃–Al₂O₃–SiO₂ glasses are hitherto limited to the Y system,^16,26–28 except in our recent modeling work on Lu¹⁶ and La²⁵ 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³⁺ 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²⁹ of r = 2.21 and r = 2.45, respectively, with

r = n_O/(n_Al + n_Si), (1)

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. (2)

The glass compositions correspond to those reported by Iftekhar et al.,¹⁵i.e., 0.171RE₂O₃–0.276Al₂O₃–0.553SiO₂ [RE(2.21)] or 0.238RE₂O₃–0.254Al₂O₃–0.508SiO₂ [RE(2.45)] that are associated with {x_Si = x_Al = 0.381; x_RE = 0.238} and {x_Si = x_Al = 0.340; x_RE = 0.320}, respectively. Except for Sc(2.45), which is outside the glass-forming region attainable by a regular melt-quench procedure at 1650 °C^15,17 and that has only been prepared with 20 at.% of Sc substituted by either Y or La,¹⁵ several physical properties and ²⁷Al/²⁹Si NMR parameters are reported for all RE(2.21) and RE(2.45) specimens, as well as for related RE AS compositions.^{11,15–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,³⁰ which was employed to derive new Born–Mayer RE–O interionic potential parameters for RE = La, Y, Lu, Sc. These are subsequently evaluated in Section 4, where the modeled short-range features of RE₂O₃–Al₂O₃–SiO₂ glasses—as well as the energy-optimized lattice parameters of several RE-bearing oxide and aluminate structures—are compared with experimental literature data. Section 5 explores the dependence of local and intermediate-range structures of the RE(2.21) and RE(2.45) glass series on variations of each of the four RE³⁺ ions and their respective contents, where we examine the RE and O speciations, the cation constellations at the various O environments, as well as the Si–Al intermixing and the nature of their shared bridging oxygen (BO) sites.

2 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.³¹ The computations were carried out with the DLPOLY3 package^32,33 for ≈3300 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,¹⁶ Lu¹⁶ and Sc¹⁷ bearing samples. The experimental densities ranged between 2.93 g cm⁻³ for Sc(2.21) and 5.05 g cm⁻³ for Lu(2.45).¹⁵ 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⁻¹ 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 α and β with indices j and k separated at a distance R^jk_α–β, the corresponding interaction energy included both long-range Coulombic (∼z_αz_β/[4πε₀R^jk_α–β], where z_α represents the ionic charge), and short-range Buckingham terms. The latter is given by

U^jk_α–β(R^jk_α–β) = A_α–β exp{−R^jk_α–β/ρ_α–β} − C_α–β(R^jk_α–β)⁻⁶, (3)

except for the RE–O (RE = La, Y, Lu, Sc) interionic potentials derived herein that solely employed the repulsive Born–Mayer contribution:

U^jk_α–β(R^jk_α–β) = A_α–β exp{−R^jk_α–β/ρ_α–β}, (4)

i.e., for which C ≡ 0. The sets {A_α–β,ρ_α–β,C_α–β} are listed in Table 1. The Buckingham/Born–Mayer parameters employed partial charges of −1.2e for O, 2.4e for Si, and 1.8e for Al together with all RE species,³⁴ where e is the elementary charge. A modified Buckingham potential was implemented to circumvent strong attractions at small interionic distances.³⁵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³³ with a 1.2 nm real-space cut-off and an accuracy of 10⁻⁶. 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.³¹

Table 1 Buckingham potential parameters

Pair^a	A/eV	ρ/Å	C/(eV Å^6)
a Superscripts indicate values of the partial ionic charges according to van Beest et al.^34 b Parameters from Du and Cormack.^26,40,41 c The new parameters derived herein, corresponding to the set {A^(3)RE–O,ρ^(3)RE–O} of Section 3.
O^1.2−–O^1.2− ^b	1844.7458	0.343645	192.58
Al^1.8+–O^1.2−	12 201.417	0.195628	31.997
Si^2.4+–O^1.2−	13 702.905	0.193817	54.681
Y^1.8+–O^1.2−	29 526.977	0.211377	50.477
La^1.8+–O^1.2− ^c	9537.95	0.235917	0
Y^1.8+–O^1.2−	11 019.6	0.224173	0
Lu^1.8+–O^1.2−	6657.09	0.230360	0
Sc^1.8+–O^1.2−	3095.04	0.244357	0

---

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 ≲3% and ≲5% 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 (<1% relative error). The sets of simulations revealed a negligible spread within ±0.01 g cm⁻³ for each glass; incidentally, this value is equal to the experimental uncertainty.^15–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 ≲10 ps, and simulations of larger ensembles over long timescales can only be achieved by classical MD approaches³¹ 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,³⁸ and effective medium³⁹ options.

We employed the iterative force-matching procedure of Ercolessi and Adams³⁰ to obtain the parameters A_α–β and ρ_α–β 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:

(1) For a given set of Buckingham parameters {A_α–β,ρ_α–β,C_α–β} of Si–O, Al–O and O–O interionic pairs (using the values of Du et al.^26,40,41) and Born–Mayer counterparts {A⁽ⁿ⁾_RE–O,ρ⁽ⁿ⁾_RE–O} for the RE–O pairs, a glass structure model was obtained using classical MD simulations (see Section 2) at the nth iteration stage. For the case of RE = Y, the initial set {A⁽⁰⁾_Y–O,ρ⁽⁰⁾_Y–O} was derived by mapping the Buckingham parameters of Du²⁶ onto the Born–Mayer form [{A_Y–O,ρ_Y–O,C_Y–O} → {A⁽⁰⁾_Y–O,ρ⁽⁰⁾_Y–O}] using a numerical fitting, whereas the A⁽⁰⁾_RE–O starting values were scaled for the other RE–O contacts according to the ratios of their Shannon–Prewitt⁴² ionic radii (R_RE/R_Y).

(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⁴⁵ was used with the Wien2k package.⁴⁶ One k-point of the irreducible Brillouin zone was computed by employing basis functions that obeyed the criterion K_maxR_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.⁴⁶ We verified that the product K_maxR_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.⁴⁷

(3) The DLPOLY3 package^32,33 was employed to compute a set of forces {F̄_j(A_RE–O,ρ_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³⁺ 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⁽ⁿ⁺¹⁾_RE–O,ρ⁽ⁿ⁺¹⁾_RE–O}, extracted at the minimum of χ²(A_RE–O,ρ_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⁽ⁿ⁾_RE–O,ρ⁽ⁿ⁾_RE–O} with n = 1, 2, and 3. Table 1 lists the {A⁽³⁾_RE–O,ρ⁽³⁾_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 ≈3300 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}.

4 Validation of the interionic potential parameters

4.1 Short-range structures of RE₂O₃–Al₂O₃–SiO₂ glasses

The general validity of the new potential parameters for reproducing short-range structural features of RE₂O₃–Al₂O₃–SiO₂ glasses was previously confirmed over large RE–Al–Si compositional ranges within each La,²⁵ Y¹⁶ and Lu¹⁶ system. Focussing on the structural dependency on the identity of the RE³⁺ 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 cation–oxygen 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⁴⁺ is exclusively (>98.5%) present in tetrahedral coordination (Si^[4]), which is in excellent agreement with experimental solid-state ²⁹Si NMR data.^{15–17,23,24} For this reason, we henceforth do not explicitly indicate the coordination number of Si. All Al³⁺ ions (≥99.9%) coordinate either four (Al^[4]), five (Al^[5]) or six (Al^[6]) oxygen atoms, also in full accordance with ²⁷Al NMR results obtained from the present RE AS compositions as well as related ones.^{12–18,23,25}

Fourfold Al³⁺ coordinations dominate the three coexisting AlO₄, AlO₅ and AlO₆ polyhedra. The amounts of higher-coordination species (mainly AlO₅) grow primarily for increasing RE³⁺ 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³⁺ and Sc³⁺ 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 ≈4.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₅ population, whereas that for the AlO₆ groups remains low throughout (x^[6]_Al < 0.05).

Table 2 Results of MD Simulations and ²⁷Al MAS NMR^a

Glass	Z̄ RE	Z̄ Al	x ^[4]Al	x ^[5]Al	x ^[6]Al	Z̄ O	x ^[0]O	x ^[1]O	x ^[2]O	x ^[3]O
a MD-derived RE^3+ and Al^3+ average coordination numbers (Z̄RE, Z̄Al) and fractional populations of AlOp groups (x^[p]Al) obtained by MD simulations and NMR. The latter data are given within parentheses and are reproduced from our recent work.^16,17,25 The populations are associated with uncertainties of ±0.010 (MD) and ±0.015 (NMR). The rightmost columns list the mean O coordination numbers (Z̄O) and the fractional population x^[p]O of each O^[p] (p = 0, 1, 2, 3) coordination, where only bonds to Si and Al are considered. b Results obtained by using the Y–O potential parameters of Du.^26
La(2.21)	6.39	4.13(4.14)	0.872(0.877)	0.123(0.091)	0.005(0.022)	1.84	0.009	0.229	0.677	0.085
Y(2.21)	6.18	4.24(4.25)	0.773(0.776)	0.209(0.195)	0.018(0.029)	1.86	0.010	0.226	0.655	0.109
Lu(2.21)	5.66	4.23(4.31)	0.782(0.718)	0.201(0.250)	0.017(0.032)	1.86	0.012	0.232	0.641	0.115
Sc(2.21)	5.13	4.25(4.40)	0.767(0.649)	0.217(0.307)	0.016(0.044)	1.86	0.017	0.231	0.629	0.123
La(2.45)	6.42	4.14(4.14)	0.871(0.884)	0.122(0.089)	0.007(0.027)	1.66	0.022	0.350	0.576	0.052
Y(2.45)	6.29	4.30(4.30)	0.729(0.746)	0.246(0.210)	0.025(0.044)	1.70	0.025	0.333	0.564	0.078
Lu(2.45)	5.77	4.29(4.34)	0.739(0.698)	0.237(0.260)	0.024(0.042)	1.69	0.031	0.329	0.559	0.081
Sc(2.45)	5.38	4.43	0.621	0.331	0.048	1.72	0.034	0.317	0.544	0.105
Y(2.21)^b	6.51	4.22(4.25)	0.795(0.776)	0.191(0.195)	0.014(0.029)	1.86	0.006	0.227	0.669	0.098
Y(2.45)^b	6.60	4.28(4.30)	0.738(0.746)	0.239(0.210)	0.023(0.044)	1.69	0.019	0.337	0.578	0.066

---

Table 2 compares the modeled Z̄_Al and {x^[p]_Al} results with their experimental counterparts obtained from solid-state ²⁷Al NMR in ref. 16, 17 and 25. The x^[p]_Al values are plotted against the RE³⁺ 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₅ 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 ≈4% [for Sc(2.21)], the pronounced experimental trend of a monotonic elevation of Z̄_Al for increasing RE³⁺ 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.²⁵ for further discussions.

Fig. 1 Fractional populations as obtained from MD simulations (black symbols) and ²⁷Al NMR experiments^16,17,25 (grey symbols), plotted against the RE³⁺ CFS for the as-indicated RE(2.21) and RE(2.45) AS glasses for (a) p = 4; (b) p = 5 and p = 6.

4.1.2 Cation–oxygen distances. Fig. 2 displays a selection of RDFs [i.e., g_α–β(R_α–β)] for the four RE(2.21) glasses (RE = La, Y, Lu, Sc), whereas Table 3 lists a collection of average interionic distances (R̄_α–β), as well as the R-value where g_α–β(R_α–β) is maximized (R^max_α–β). Here we focus on the various cation–oxygen distances.

Fig. 2 RDFs of the (a) Si–O, (b) Al–O, (c) RE–O, and (d) O–O pairs for each of the RE(2.21) glasses with RE = {La, Y, Lu, Sc}, the identity of which is identified by the legend in (a).

Table 3 Average interatomic distances and bond angles^a

Pair α–β (R)	R cut /pm^b	RE(2.21)					RE(2.45)
		La	Y	Y^c	Lu	Sc	La	Y	Y^c	Lu	Sc
a The distances R^maxα–β (in pm) and bond angles θ^maxα–β–γ (in degrees) were estimated from the maximum of the radial distribution function (RDF) and the angular distribution function (ADF), respectively. The corresponding average distance (R̄α–β) or bond angle ([small theta, Greek, macron]α–β–γ) over each respective distribution is given within parentheses. b The cutoff distance Rcut corresponds to the first minimum of the RDF, which varied depending on the RE^3+ cation identity {La, Y, Lu, Sc} as follows: {315, 303, 289, 275} pm for the RE–O pair; {434, 415, 405, 395} pm for Si–RE; {465, 434, 425, 400} pm for Al–RE and {495, 480, 465, 430} pm for RE–RE. c Results obtained by the potential parameters of Du.^26
O–O	395	260.3(310.0)	260.7(309.1)	261.1(308.9)	260.5(309.1)	260.9(309.0)	260.3(312.0)	261.1(310.1)	262.1(309.7)	261.1(309.9)	261.3(308.5)
Si–O	250	159.3(160.1)	159.3(160.2)	159.1(160.3)	159.3(160.2)	159.1(160.2)	159.3(160.0)	159.1(160.1)	159.1(160.2)	159.1(160.1)	158.9(160.1)
Al–O	250	176.7(180.3)	176.9(181.9)	176.9(181.8)	176.9(181.9)	176.9(182.1)	176.5(180.5)	177.1(182.3)	176.9(182.6)	176.9(182.4)	177.5(183.6)
Al^[4]–O	250	175.9(178.3)	175.9(178.3)	175.6(178.6)	175.9(178.4)	175.6(178.5)	176.1(178.0)	176.1(178.0)	175.9(178.2)	175.9(178.1)	175.6(178.1)
Al^[5]–O	250	179.9(190.9)	179.9(190.3)	179.6(190.8)	179.9(190.5)	179.6(190.4)	180.9(190.6)	180.9(190.0)	180.9(190.2)	179.6(190.0)	180.1(189.6)
Al^[6]–O	250	185.6(199.5)	181.9(198.8)	186.1(199.4)	182.4(198.4)	182.9(198.5)	181.6(200.2)	184.9(198.3)	185.9(198.7)	185.4(198.3)	185.4(197.7)
RE–O		234.3(250.9)	224.3(239.6)	228.9(240.9)	216.3(229.3)	204.9(217.6)	233.7(249.7)	224.7(238.9)	228.9(240.3)	216.5(228.8)	205.5(217.9)
Si–Si	338	312.9(309.3)	312.5(308.9)	312.9(308.8)	312.7(308.7)	313.1(308.7)	313.5(309.5)	311.3(308.3)	312.9(308.5)	312.7(308.4)	312.1(308.8)
Si–Al	352	320.3(316.2)	318.9(315.6)	318.1(315.8)	320.3(315.5)	321.1(315.6)	315.9(316.0)	315.9(315.3)	318.1(315.6)	314.9(315.3)	316.7(315.1)
Si–RE		363.5(366.2)	352.9(355.1)	354.1(353.1)	349.5(349.5)	342.1(343.0)	360.9(365.2)	351.7(354.1)	351.9(353.0)	346.9(348.6)	338.7(342.1)
Al–Al	370	314.7(317.6)	315.1(316.9)	314.3(317.5)	312.9(316.4)	315.5(315.8)	315.1(318.6)	314.7(317.1)	313.9(318.4)	313.1(316.7)	314.9(316.0)
Al–RE		359.5(370.3)	352.5(356.1)	313.3(354.5)	343.8(349.3)	338.7(339.1)	360.5(370.3)	350.5(356.1)	356.7(354.9)	346.7(349.7)	338.7(338.7)
RE–RE		370.7(405.9)	352.3(390.3)	352.3(391.5)	337.9(378.0)	360.9(360.0)	369.5(405.2)	353.5(389.4)	352.9(390.6)	340.3(377.8)	355.1(359.3)
Triplet α–β–γ (θ)
O–Si–O		108(109)	108(109)	108(109)	109(109)	108(109)	108(109)	108(109)	108(109)	108(109)	108(109)
O–Al–O		105(109)	104(108)	103(108)	104(108)	105(108)	105(109)	103(108)	102(108)	103(108)	101(108)
Si–O–Si		149(150)	148(150)	149(150)	150(150)	151(150)	151(150)	148(149)	147(149)	148(149)	146(149)
Si–O–Al		131(138)	131(136)	129(137)	129(136)	128(136)	132(138)	129(135)	129(135)	128(135)	127(133)
Al–O–Al		120(121)	120(119)	120(120)	120(118)	120(118)	121(123)	119(119)	120(120)	120(118)	121(116)

---

Si⁴⁺ exhibits the highest CFS of all cations: it arranges well-ordered SiO₄ environments, reflected in a sharp peak of the g_Si–O(R_Si–O) function [see Fig. 2(a)] and consequently essentially equal values of R^max_Si–O and R̄_Si–O, amounting to 159 pm and 160 pm, respectively (Table 3). These Si–O distances are independent of the RE³⁺ identity and display good agreement with the value R_Si–O = (162 ± 3) pm reported in ref. 18 and 22 by deconvolution of X-ray and neutron diffraction data into Gaussian peaks for Si-rich and RE-poor (RE = La, Y, Sc) AS glasses.

The Al³⁺ ion that features the second largest CFS also manifests essentially identical RDF peak maxima ≈177 pm among the various RE-bearing glasses; see Fig. 2(b) and Table 3. However, for a given RE AS glass, minor differences are now observed between R^max_Al–O and R̄_Al–O: the latter values (180–184 pm) are slightly larger than their R^max_Al–O counterparts, with the R̄_Al–O–R^max_Al–O separation increasing together with the RE³⁺ CFS along the series La³⁺ < Y³⁺ < Lu³⁺ < Sc³⁺. This trend is expected from the concomitant population-elevations of the high-coordination (AlO₅, AlO₆) polyhedra and the accompanying lengthened Al^[p]–O distances for increasing p; compare the Al^[4]–O, Al^[5]–O and Al^[6]–O bond-lengths with their aggregate “Al–O” value in Table 3. The values R^max_Al–O ≈ 177 pm and R̄_Al–O ≈ 182 pm accord well with the experimental diffraction data of ref. 18 and 22 that yielded Al–O distances in the range 180–184 pm (±3 pm) for La, Y and Sc glasses.

We next consider the oxygen distances to the RE³⁺ species, i.e., the lowest field-strength cation in each glass structure. The first peak-maximum of the respective g_RE–O(R_RE–O) function [Fig. 2(c)] becomes more diffuse as the CFS decreases along the series Sc³⁺ > Lu³⁺ ≳ Y³⁺ > La³⁺. This translates both into progressively lengthened R̄_RE–O distances [from 218 pm for the Sc(2.21) structure to 251 pm for La(2.21)] and R^max_RE–O [205 pm for Sc(2.21) to 234 pm for La(2.21)] in Table 3, as well as into steadily increased R̄_RE–O–R^max_RE–Odifferences [from 13 pm for Sc(2.21) to 17 pm for La(2.21)] as the RE³⁺ CFS diminishes. The R^max_RE–O-values of Table 3 compare very well with their corresponding diffraction-derived data,^18,22 where the latter values are stated within parentheses (or ranges of values in the case of several data-points): 234 pm (234–238 pm)^18,22 for La; 224 pm (220–226 pm)^18,22 for Y, and 205 pm (212 pm)¹⁸ for Sc. The largest discrepancy (≈3%) is observed for Sc, although it may partially originate from widely differing cation compositions between our Sc(2.21) glass and that analyzed by Sadiki et al.;¹⁸ the latter corresponds to r = 1.87, x_Al/x_Si = 0.58 and x_Sc ≈ 0.03 in our sample nomenclature. While both the Y–O interionic parameters derived herein and those of Du²⁶ provide R^max_Y–O values in good accordance with experimental diffraction data, slightly better predictions are observed from the present {A_RE–O,ρ_RE–O} set. Experimental Lu–O bond-lengths are not reported from Lu–Al–Si–O glasses. Yet, as shown below in Section 4.2, our interionic Lu–O potential parameters provide excellent predictions of R̄_Lu–O-values and lattice-parameters for some crystalline Lu-bearing phases.

4.2 Crystalline model structures

The new RE–O potentials were further evaluated by energy minimizations at zero absolute temperature by applying the GULP program⁴⁸ to the set of crystalline structures shown in Table 4; they correspond to the La₂O₃, Y₂O₃, Lu₂O₃, and Sc₂O₃ oxides, some RE-bearing aluminates (LaAlO₃, Y₃Al₅O₁₂, Lu₃Al₅O₁₂, ScAlO₃) and scandium disilicate (Sc₂Si₂O₇). Each representative XRD-derived structural parameter-set (obtained from the inorganic crystal structure database⁴⁹) 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.

Table 4 Parameters for energy-optimized structures^a

Oxide	a/pm	b/pm	c/pm	V/nm^−3	R̄ RE–O /pm^b
a Results of using the interionic potential parameters of Table 1 in lattice-energy calculations that assumed fixed (α, β, γ) unit-cell angles, and leading to the calculated (Calc.) cell dimensions (a, b, c) and volume (V). Values within parentheses represent relative deviations from the XRD-derived experimental (Exp) data. b Average RE–O distance, except for the data on α-quartz and α-Al2O3 that reports R̄Si–O and R̄Al–O, respectively. c Obtained from the Y–O potential parameters of Du.^26 d The β-value was optimized to β = 104.08°, which may be compared with the XRD-derived result of β = 102.77°. e Obtained by using the Si–O or Al–O potential parameters of ref. 26, 40 and 41 (see Table 1).
La2O3 (γ = 120°)
Exp^60	393.8	—	613.6	0.0824	253.5
Calc.	382.5 (−2.9%)	—	603.2 (−1.7%)	0.0764 (−7.3%)	247.5 (−2.4%)
LaAlO3 (γ = 120°)
Exp^61	536.5	—	1311.0	0.3268	268.2
Calc.	537.0 (0.09%)	—	1309.2 (−0.14%)	0.3270 (0.05%)	269.1 (0.34%)
Y2O3
Exp^62	1060.7	—	—	1.1935	228.6
Calc.	1060.3 (−0.04%)	—	—	1.1921 (−0.11%)	228.4 (−0.09%)
Calc.^c	1059.3 (−0.13%)	—	—	1.1887 (−0.40%)	228.1 (−0.22%)
Y3Al5O12
Exp^63	1200.6	—	—	1.7307	237.7
Calc.	1208.4 (0.65%)	—	—	1.7645 (2.0%)	239.7 (0.84%)
Calc.^c	1201.5 (0.07%)	—	—	1.7344 (0.22%)	237.8 (0.04%)
Lu2O3
Exp^64	1039.1	—	—	1.1220	223.6
Calc.	1032.3 (−0.65%)	—	—	1.1001 (−2.0%)	222.5 (−0.49%)
Lu3Al5O12
Exp^65	1190.6	—	—	1.6877	233.0
Calc.	1199.4 (0.74%)	—	—	1.7256 (2.2%)	236.8 (1.6%)
Sc2O3
Exp^66	984.9	—	—	0.9554	212.1
Calc.	1002.7 (1.8%)	—	—	1.0082 (5.5%)	216.4 (2.0%)
ScAlO3
Exp^67	493.7	523.2	720.5	0.1861	226.8
Calc.	504.0 (2.1%)	526.0 (0.54%)	742.7 (3.1%)	0.1969 (5.8%)	232.8 (2.7%)
Sc2Si2O7
Exp^68	650.3	849.8	468.2	0.2523	212.4
Calc.	658.2 (1.2%)	851.1 (0.15%)	482.1 (3.0%)	0.2619 (3.8%)	219.0 (3.1%)
α-SiO2 (γ = 120°)
Exp^69	490.2	—	540.0	0.1124	161.3
Calc.^e	493.1 (0.59%)	—	544.2 (0.78%)	0.1146 (2.0%)	158.9 (−1.5%)
BKS^34	494.0 (0.78%)	—	544.8 (0.89%)	0.1151 (2.5%)	160.1 (−0.74%)
CHIK^51	504.5 (2.9%)	—	552.0 (2.2%)	0.1217 (8.3%)	161.8 (0.31%)
α-Al2O3 (γ = 120°)
Exp^70	476.0	—	1299.6	0.2551	191.3
Calc.^e	476.7 (0.15%)	—	1316.2 (1.3%)	0.2591 (1.6%)	192.5 (0.63%)

---

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₂O₃ 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 also compares the results of using the Y–O interionic potential parameters of Du et al.²⁶ with those derived herein. Overall, the two options display comparable qualities in their structural predictions; the present {A_Y–O,ρ_Y–O} parameters perform favorably for the case of Y₂O₃, whereas those of ref. 26 give better agreement with the experimental lattice parameters of Y₃Al₅O₁₂.

The La–O and Sc–O interionic potentials offer the weakest predictability. While the LaAlO₃ lattice parameters calculated by the new La–O potential give excellent agreement with the experimental data (Table 4), the globally largest discrepancy among all structures is observed for La₂O₃, 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₂O₃ 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.⁵⁰ obtained by a Morse-based potential that was optimized directly on crystalline oxides and subsequently evaluated on several multicomponent structures. Table 4 also verifies acceptable values of the energy-minimized lattice parameters of α-quartz and α-Al₂O₃, which originated from solely employing the O–O, Si–O, and Al–O interionic parameters of Du and Cormack^26,40,41 (see 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. Moreover, the herein utilized Si–O and O–O parameters from ref. 40 and 41 compare favorably with the Beest–Kramer–Santen³⁴ (BKS) and Carré–Horbach–Ispas–Kob⁵¹ (CHIK) potentials when evaluated for the case of α-SiO₂ (see Table 4).

5 Cation field-strength effects on RE₂O₃–Al₂O₃–SiO₂ glass structures

5.1 Rare-earth coordinations

Fig. 3(a) plots the average number of oxygen species coordinated by each RE³⁺ ion (Z̄_RE) against its corresponding Shannon–Prewitt radius,⁴² the latter value relevant for REO₆ groups. Indeed, all MD-derived mean coordination numbers scatter around Z̄_RE ∼ 6; they elevate either for increasing RE³⁺ 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 and 25. Table 2 verifies our claimed correlation between Z̄_RE and Z̄_Al for the present glasses.

Fig. 3 (a) MD-derived average RE³⁺ coordination numbers (Z̄_RE) plotted against the respective Shannon–Prewitt ionic radius (assuming sixfold coordinations) for each RE(2.21) and RE(2.45) series. (b and c) Fractional populations x^[p]_RE plotted versus the coordination number p for the (b) RE(2.21) and (c) RE(2.45) series of glasses.

The underlying distribution of RE^[p] populations (x^[p]_RE) are conveyed by Fig. 3(b) and (c) for the RE(2.21) and RE(2.45) series. Coordination numbers are observed in the range 4 ≤ p ≤ 8, although the fourfold and eightfold populations are only significant for the smallest (Sc³⁺) and largest (La³⁺) cations, respectively. Except for Sc³⁺ ions that primarily form ScO₅ polyhedra, p = 6 coordinations dominate the {RE^[p]} distributions of the RE³⁺ cations, where REO₇ polyhedra constitute the second most abundant species for La³⁺ and Y³⁺, whereas >30% of the (small) Lu³⁺ cations assume fivefold coordinations. The MD-modeled Z̄_RE-values around six accord semi-quantitatively with the X-ray absorption spectroscopy results of Z̄_RE = 6 ± 0.5 observed throughout the lanthanide-ion series when RE³⁺ is present in a low amount (5 wt% RE₂O₃) in a Na₂Si₂O₅ based glass,⁵² as well as with several MD reports on RE-bearing (alumino)silicate glasses.^6,7,26–28

Interestingly, the calculations predict significant fractions 0.09 ≤ x^[4]_Sc ≤ 0.17, suggesting that Sc³⁺ might partially assume a glass-network forming role in the guise of ScO₄ groups, as reported previously for the similarly-sized Mg²⁺ ion.^53–55Fig. 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 ²⁹Si NMR shifts and unexpectedly low glass transition temperatures observed from Sc AS glasses.^15,17

Fig. 4 Some typical fragments in the MD-derived glass structures. (a) A ScO₄ tetrahedron (O atoms of which are set in black) in relation to the aluminosilicate network. (b) A SiO₄ tetrahedron accommodating one NBO ion and sharing corners with each of a SiO₄, AlO₄ and AlO₅ group. (c) A “tricluster moiety” featuring an AlO₆ polyhedron connected to SiO₄ and AlO₄ groups via corner and edge-sharing, respectively.

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.⁵⁶ They dominate the total O speciation, as follows from Table 2 that lists the fractional population {x^[p]_O} of each O^[p] moiety (p = 0, 1, 2, 3) observed in the modeled data. The relative amounts of the BO and NBO environments depend predominantly on the r-value of the AS glass and hence its RE³⁺ content; nearly constant x^[1]_O fractions of ≈0.23 and 0.32–0.35 are indeed observed for the series of RE(2.21) and RE(2.45) specimens, respectively.

The MD-predicted non-negligible amounts of O^[0] (“free O²⁻ ion”) and significant O^[3] (“oxygen tricluster”) moieties in the present (as well as previously reported^16,25–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³⁺ CFS increases, they also depend on the RE³⁺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³⁺ ions suggests an enhanced RE–O–RE association.

Fig. 5 Fractional populations of O^[p] species with (a) p = 0, (b) p = 3, and (c and d) p = 2, plotted against the RE³⁺ CFS for the as-indicated RE(2.21) and RE(2.45) series of glasses. Note that whereas the lower and upper values vary among the plots in (a–d), their spans of vertical ranges are equal throughout. Grey lines represent best fits.

While the abundance of free O²⁻ 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₅ and AlO₆ groups (see Section 5.2.2). Because the latter populations grow together with both x_Al and the RE³⁺ 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 ≲ 2.21, whereas x^[3]_O stay consistently below 0.10 in La–Al–Si–O structures.²⁵ These trends may be verified from our previous reports on glasses featuring variable RE : Al : Si contents within each La,²⁵ Y,¹⁶ and Lu¹⁶ 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₆ population remains low throughout all glass structures; (ii) Z̄_O and (iii) x^[3]_O, all of which increase with the RE³⁺ CFS. However, whereas the Z̄_Al/Z̄_O/x^[3]_O correlations are evident within both RE AS r-branches, their concomitant elevations with the CFS is only transparent along the RE(2.45) series.

5.2.2 Trends in cation–O^[p] contacts. We now examine the various cation–oxygen contacts, initially from the viewpoint of each Si⁴⁺, Al³⁺, and RE³⁺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 x^[1]_O/x^[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) ≫ 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₄ groups, whereas ≲30% and ≲15% reside on AlO₄ and AlO₅/AlO₆ 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₅/AlO₆ 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.²⁷

Fig. 6 Variation of the RE–O^[p] fractional populations associated with (a) p = 0 or (b) p = 1 and p = 2 for each as-indicated series of RE(2.21) and RE(2.45) glasses. Solid and open symbols in (b) represent data for RE–O^[1] and RE–O^[2] contacts, respectively. Grey lines represent best-fit data.

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³⁺ 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³⁺ cations. In contrast, the free O²⁻ 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³⁺ content (x_RE), but diminishes markedly across each RE(r) series as the CFS increases. For instance, O^[0]_[4] groups constitute ≈50% and ≈20% 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³⁺ 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 NBO environments are dominated by O^[1]_[2] constellations involving Si and RE, as well as O^[1]_[3] sites where O coordinates two RE cations together with either of Si or Al^[p]. The O^[1]_[3] species represent altogether 60–70% of the entire O^[1] population. The relative abundance of each O^[1]_[3]–Si/Al^[p] moiety is reflected by the relevant fractional population of the cation species (see Table 2), weighted by the respective affinity for X–NBO bond formation given by eqn (6). Furthermore, there are minor contributions from numerous distinct O^[1]_[4] structural fragments, in total representing 3–18% of the total O^[1] speciation. For increasing RE³⁺ CFS, the q = 2 populations grow and those of q = 4 diminish, while the reverse trends are observed when the RE³⁺ content elevates from x_RE = 0.24 [RE(2.21)] to x_RE = 0.32 [RE(2.45)].

The same general CFS/x_RE trends are observed in the number of O^[2]–RE contacts; the O^[2]_[4] moieties remain ≲10% throughout all structures, whereas the lower formal charge of the O^[2] atoms relative to their O^[1] counterparts implies overall higher O^[2]_[2] contributions, dominated by Si–O^[2]–Si and Si–O^[2]–Al^[4] fragments. Most Si–O^[2]–Si bridges constitute O^[2]_[2] groups devoid of contacts with RE³⁺. Conversely, roughly equal amounts of O^[2]_[2] and O^[2]_[3] motifs contribute to the overall dominating Si–O^[2]–Al^[4] linkages. The O^[2]_[3] species involve O^[2] coordinated by one RE³⁺ cation together with any two members of the set {Si, Al^[4], Al^[5], Al^[6]}. As expected from the accumulated negative charges associated with Al^[4]–O^[2]–Al^[4] bridges,⁵⁶ all such motifs require one RE³⁺ cation for stabilization, reflected in their sole presence as O^[2]_[3] (as opposed to O^[2]_[2]) groups. O^[2]_[3] moieties are also encountered for the majority of all interpolyhedral bonds involving either of AlO₅ or AlO₆, such as Si–O^[2]–Al^[5]/Al^[6] and Al^[4]–O^[2]–Al^[5]/Al^[6] bridges. Generally, the amounts of the various O^[2]_[3] groups involving Si–O–Al^[p] and Al^[p]–O–Al^[q] linkages with (p, q) = (4, 5, 6) roughly obey a statistical partitioning based on the products between the relative fractions {x_Si, x_Alx^[p]_Al} associated with each pair of Si/Al^[p] species, as discussed in detail below. Throughout all glass structures, ≈75–80% of all O^[2] sites constitute O^[2]_[2](Si, Si), O^[2]_[2](Si, Al^[4]), O^[2]_[3](Si, Al^[4]), O^[2]_[3](Al^[4], Al^[4]) and O^[2]_[3](Si, Al^[5]) moieties, where the two interlinked cations are indicated within parentheses.

5.3 SiO₄–AlO_p and AlO_p–AlO_q interpolyhedral connectivities

5.3.1 Al–Si intermixing. The MD-generated glass structures exhibit networks built by interconnected SiO₄ and AlO₄ tetrahedra, as well as network-associated AlO₅ and AlO₆ polyhedra. These groups link through O^[2] and O^[3] bridging oxygen species, both of which we henceforth collectively refer to as “BO” atoms.

The various SiO₄/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 = y_Alx^[p]_Al = 0.5x^[p]_Al, with x^[p]_Al listed in Table 2. Consequently, owing to the dominance of SiO₄ and AlO₄ groups followed by AlO₅, the most frequently encountered pairs of structural building blocks represent SiO₄–AlO₄, SiO₄–SiO₄, AlO₄–AlO₄ motifs, and to a lesser extent SiO₄–AlO₅ and AlO₄–AlO₅. Fig. 7 plots the relative abundances of the various pair-connectivities. Note that owing to the small amounts of AlO₆ polyhedra, all AlO₆–XO_p contacts remain low throughout all structures, with the highest fractions x(Si–Al^[6]) ≈ x(Al^[4]–Al^[6]) ≈ 0.03 observed in the Sc(2.45) glass.

Fig. 7 Fractional populations of the SiO₄–SiO₄, SiO₄–AlO₄, SiO₄–AlO₅, AlO₄–AlO₄, AlO₄–AlO₅, and AlO₅–AlO₅ interpolyhedral pairs encountered in each (a) RE(2.21) and (b) RE(2.45) glass structure. Each Al^[p] coordination is for brevity labeled by its respective integer “p”. Data are not shown for pairs involving AlO₆ groups as their populations remain <5% throughout.

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₄, 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₂O₃–Al₂O₃–SiO₂ glasses, a pronounced Si–Al disorder was verified experimentally for SiO₄–AlO₄ [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 ²⁷Al NMR spectroscopy.⁵⁹

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₄ or AlO_p groups. For increasing RE³⁺ 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) and their affinity for edge-sharing increases concurrently with the RE³⁺ CFS, the x^[3]_O values decreases with the r-value, as discussed in Section 5.2.1.

Fig. 8 Fraction of bridging oxygen species participating in edge-shared SiO₄/AlO_p polyhedra, plotted against the RE³⁺ CFS and considering (a) the total BO speciation (i.e., sum over O^[2] and O^[3]), or solely each individual (b) O^[2] and (c) O^[3] population.

Striking differences further emerge when comparing the preferences among the various {SiO₄, AlO₄, AlO₅, AlO₆} 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 ≲ Al^[4] ≪ Al^[5] < Al^[6]. Altogether, the net abundance of edge-shared XO_p–YO_q polyhedra scales roughly as

Si–Si(0%) ≈ Si–Al^[4](<1%) < Al^[4]–Al^[4](<3%) < Si–Al^[5](≈5%) < Si–Al^[6](≈10%) < Al^[4]–Al^[5](10–20%) < Al^[4]–Al^[6](20–25%) ≪ Al^[5]–Al^[5](≈50%) < Al^[5]–Al^[6](60–70%) < Al^[6]–Al^[6](≳75%), (8)

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₅–AlO₆ and AlO₆–AlO₆ linkages occur through edge-sharing, 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₄/AlO_p groups when the RE³⁺ 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 θ_X–O–Y ≈ 90° that may be compared with their wider corner-shared analogs peaking in the ranges of 110–130° for Al–O–Al/Si and 140–160° 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 θ_Al–O–Al = 90°, which is most transparent for the structures involving Sc (and Lu; not shown) that exhibit the largest AlO₅/AlO₆ populations and thereby the highest extent of edge-sharing.

Fig. 9 (a) Si–O–Si, (b) Si–O–Al and (c) Al–O–Al bond-angle (θ_X–O–Y) distributions for the La(2.45) (left panel) and Sc(2.45) (right) glass structures. “Al” represents the entire {Al^[4], Al^[5], Al^[6]} speciation. Grey and dotted traces reveal the respective θ(X–O^[3]–Y) and θ(X–O^[2]–Y) contributions (where {X, Y} = {Si, Al}) to the total ADF (black trace).

6 Concluding remarks

New RE–O interatomic potential parameters were optimized for RE = {La, Y, Lu, Sc} by a force-matching procedure³⁰ 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 <2% deviation, with the largest discrepancy (3.5%) observed for the Sc(2.21) glass. Whereas Sc³⁺ mainly assumes Sc^[5] coordinations and non-negligible amounts of Sc^[4] species (up to ≈20%), all other (larger) RE³⁺ ions reveal {RE^[p]} speciations peaked at p = 6, with La^[7] and Lu^[5] constituting the second most abundant coordinations for La³⁺ and Lu³⁺, respectively. The simulated structures exhibit essentially randomized connectivities among the various SiO₄ and AlO_p (p = 4, 5, 6) network (associated) polyhedra, a feature attributed to the high charge of the trivalent RE³⁺ ions, whose comparatively strong RE³⁺–O bonds (relative to their M⁺–O and M²⁺–O counterparts) perturb any Si/Al ordering and also promote otherwise energetically disfavored structural motifs,⁵⁶ such as free O²⁻ 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³⁺ CFS increases along the series La³⁺ < Y³⁺ < Lu³⁺ < Sc³⁺:

(i) An increase in the minor but non-negligible populations of free O²⁻ ions (<3.5%) is observed, suggesting a stronger RE–O–RE self-association, where each O^[0] site coordinates 3 or 4 RE³⁺ cations. The x^[0]_O values also increase with the RE³⁺ 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₄ and AlO_p polyhedra connect through corners, the degree of edge-sharing increases steadily from ≈5% in the La-based glasses to ≈15% 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₅ and AlO₆ linkages to other groups occur via edge-sharing.

Altogether, the structures of RE₂O₃–Al₂O₃–SiO₂ 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³⁺ 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.⁵⁶

The present study corroborates and reinforces previous experimental and computational findings that the structures of RE₂O₃–Al₂O₃–SiO₂ glasses are strongly influenced by the RE³⁺ CFS over both short and medium ranges.^{12–18,24–28} The large CFS-span across the set {La³⁺, Y³⁺, Lu³⁺, Sc³⁺}, 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³⁺ ion, based on its CFS-value relative to those considered herein.

Acknowledgements

This work was supported by the Swedish Research Council (contract VR-NT 2010-4943) and the Faculty of Sciences at Stockholm University. We thank Jekabs Grins and Alexander Lyubartsev for discussions.

References

1. V. McGahay and M. Tomozawa, J. Non-Cryst. Solids, 1993, 159, 246–252.
2. N. D. Afify and G. Mountjoy, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 024202.
3. J. Coon and J. E. Shelby, Phys. Chem. Glasses, 1994, 35, 47–51.
4. J. E. Shelby, Key Eng. Mater., 1994, 94–95, 43–80.
5. T. Schaller, J. F. Stebbins and M. C. Wilding, J. Non-Cryst. Solids, 1999, 243, 146–157.
6. B. Park, H. Li and L. R. Corrales, J. Non-Cryst. Solids, 2002, 297, 220–238.
7. J. A. Johnson, C. J. Benmore, D. Holland, J. Du, B. Beuneu and A. Mekki, J. Phys.: Condens. Matter, 2011, 23, 065404.
8. A. Makishima, M. Kobayashi, T. Shimohria and T. Nagata, J. Am. Ceram. Soc., 1982, 65, C210–C211.
9. M. J. Hyatt and D. E. Day, J. Am. Ceram. Soc., 1987, 70, C283–C287.
10. J. T. Kohli and J. E. Shelby, Phys. Chem. Glasses, 1991, 32, 67–71.
11. S. Iftekhar, J. Grins and M. Edén, J. Non-Cryst. Solids, 2010, 356, 1043–1048.
12. J. T. Kohli, J. E. Shelby and J. S. Frye, Phys. Chem. Glasses, 1992, 33, 73–78.
13. T. Schaller and J. F. Stebbins, J. Phys. Chem. B, 1998, 102, 10690–10697.
14. P. Florian, N. Sadiki, D. Massiot and J. P. Coutures, J. Phys. Chem. B, 2007, 111, 9747–9757.
15. S. Iftekhar, J. Grins, P. N. Gunawidjaja and M. Edén, J. Am. Ceram. Soc., 2011, 94, 2429–2435.
16. S. Iftekhar, B. Pahari, K. Okhotnikov, A. Jaworski, B. Stevensson, J. Grins and M. Edén, J. Phys. Chem. C, 2012, 116, 18394–18406.
17. B. Pahari, S. Iftekhar, A. Jaworski, K. Okhotnikov, K. Jansson, B. Stevensson, J. Grins and M. Edén, J. Am. Ceram. Soc., 2012, 95, 2545–2553.
18. N. Sadiki, L. Hennet, P. Florian, Y. Vaills, D. Massiot and J.-P. Coutures, Mater. Tech., 2010, 98, 409–421.
19. P. J. Hayward, Radioactive Waste Forms for the Future, North-Holland, Amsterdam, 1988.
20. S. D. Chen, J. F. Hsieh, S. C. Tsai, W. Y. Lin, K. Y. Cheng and S. J. Wang, Nucl. Med. Commun., 2001, 22, 121–125.
21. B. Stevensson and M. Edén, J. Non-Cryst. Solids, 2013, 378, 163–167,  DOI:10.1016/j.jnoncrystsol.2013.06.013.
22. I. Pozdnyakova, N. Sadiki, L. Hennet, V. Cristiglio, A. Bytchkov, G. J. Cuello, J. P. Coutures and D. L. Price, J. Non-Cryst. Solids, 2008, 354, 2038–2044.
23. E. Leonova, A. S. Hakeem, K. Jansson, B. Stevensson, Z. Shen, J. Grins, S. Esmaeilzadeh and M. Edén, J. Non-Cryst. Solids, 2008, 354, 49–60.
24. S. Iftekhar, E. Leonova and M. Edén, J. Non-Cryst. Solids, 2009, 355, 2165–2174.
25. A. Jaworski, B. Stevensson, B. Pahari, K. Okhotnikov and M. Edén, Phys. Chem. Chem. Phys., 2012, 14, 15866–15878.
26. J. Du, J. Am. Ceram. Soc., 2009, 92, 87–95.
27. J. K. Christie and A. Tilocca, Chem. Mater., 2010, 22, 3725–3734.
28. J. K. Christie and A. Tilocca, J. Phys. Chem. B, 2012, 116, 12614–12620.
29. M. Edén, J. Non-Cryst. Solids, 2011, 357, 1595–1602.
30. F. Ercolessi and J. B. Adams, Europhys. Lett., 1994, 26, 583–588.
31. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford, Clarendon, 1989.
32. W. Smith and T. R. Forester, J. Mol. Graphics, 1996, 14, 136–141.
33. I. T. Todorov, W. Smith, K. Trachenko and M. T. Dove, J. Mater. Chem., 2006, 16, 1911–1918.
34. B. W. H. van Beest, H. J. Kramer and R. A. van Santen, Phys. Rev. Lett., 1990, 64, 1955–1958.
35. V. A. Bakaev and W. A. Steele, J. Chem. Phys., 1999, 111, 9803–9812.
36. R. H. Swendsen, Phys. Rev. Lett., 1979, 42, 859–861.
37. A. Lyubartsev and A. Laaksonen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 52, 3730–3737.
38. E. Lomba and J. L. Lòpez-Martín, J. Chem. Phys., 1996, 104, 5244–5250.
39. K. W. Jacobsen, J. K. Nørskov and M. J. Puska, Phys. Rev. B: Condens. Matter Mater. Phys., 1987, 35, 7423–7442.
40. J. Du and A. N. Cormack, J. Am. Ceram. Soc., 2005, 88, 2532–2539.
41. J. Du and L. R. Corrales, J. Non-Cryst. Solids, 2007, 353, 210–214.
42. R. D. Shannon, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Cryst., 1976, 32, 751–767.
43. M. Petersen, Comput. Phys. Commun., 2000, 126, 294–309.
44. D. J. Singh and L. Nordstrom, Planewaves, Pseudopotentials and the LAPW Method, Springer-Verlag, Berlin, 2006.
45. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
46. P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka and J. Luitz, WIEN2k, An Augmented Plane Wave + Local Orbitals Program for Calculating Crystal Properties, Karlheinz Schwarz, Techn. Universität Wien, Wien, Austria, 2001.
47. P. Pulay, Mol. Phys., 1969, 17, 197–204.
48. J. D. Gale and A. L. Rohl, Mol. Simul., 2003, 29, 291–341.
49. G. Bergerhoff and I. D. Brown, Crystallographic Databases, International Union of Crystallography, Chester, 1987.
50. A. Pedone, G. Malavasi, M. C. Menziani, A. N. Cormack and U. Segre, J. Phys. Chem. B, 2006, 110, 11780–11795.
51. A. Carré, J. Horbach, S. Ispas and W. Kob, Europhys. Lett., 2008, 82, 17001.
52. M. R. Cicconi, G. Giuli, E. Paris, P. Courtial and D. B. Dingwell, J. Non-Cryst. Solids, 2013, 362, 162–168.
53. C. A. Angell, P. A. Cheeseman and R. R. Kadiyala, Chem. Geol., 1987, 62, 83–92.
54. M. Guignard and L. Cormier, Chem. Geol., 2008, 256, 111–118.
55. D. R. Neuville, L. Cormier, V. Montouillout, P. Florian, F. Millot, J.-C. Rifflet and D. Massiot, Am. Mineral., 2008, 93, 1721–1731.
56. M. Edén, Annu. Rep. Prog. Chem., Sect. C, 2012, 108, 177–221.
57. L. Cormier, D. Ghaleb, D. R. Neuville, J.-M. Delaye and G. Calas, J. Non-Cryst. Solids, 2003, 332, 255–270.
58. A. Pedone, E. Gambuzzi, G. Malavasi and M. C. Menziani, Theor. Chem. Acc., 2012, 131, 1147.
59. A. Y. H. Lo and M. Edén, Phys. Chem. Chem. Phys., 2008, 10, 6635–6644.
60. P. Aldebert and J. P. Traverse, Mater. Res. Bull., 1979, 14, 303–323.
61. B. C. Chakoumakos, D. G. Schlom, M. Urbanik and J. Luine, J. Appl. Phys., 1998, 83, 1979–1982.
62. M. M Faucher and J. Pannetier, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 1980, 36, 3209–3211.
63. A. Nakatsuka, A. Yoshiasa and T. Yamanaka, Acta Crystallogr., Sect. B: Struct. Sci., 1999, 55, 266–272.
64. F. Hanic, M. Hartmanová, G. G. Knab, A. A. Urusovskaya and K. S. Bagdasarov, Acta Crystallogr., Sect. B: Struct. Sci., 1984, 40, 76–82.
65. F. Euler and J. A. Bruce, Acta Crystallogr., 1965, 19, 971–978.
66. A. Bartos, K. P. Lieb, M. Uhrmacher and D. Wiarda, Acta Crystallogr., Sect. B: Struct. Sci., 1993, 49, 165–169.
67. R. J. Hill and I. Jackson, Phys. Chem. Miner., 1990, 17, 89–96.
68. Y. I. Smolin, Y. F. Shepelev and A. P. Titov, Kristallografiya, 1972, 17, 857–858.
69. G. A. Lager, J. D. Jorgensen and F. J. Rotella, J. Appl. Phys., 1982, 53, 6751–6756.
70. L. Lutterotti and P. Scardi, J. Appl. Crystallogr., 1990, 23, 246–252.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c3cp51726h

---